My favorite prime number generator
eli.thegreenplace.net
eli.thegreenplace.net
┬─┬───────────────────────────────────┬──── ────┬
│ │ ┬─┬ ────┬─┬────────────────────── ┼───┬ ┬───┼
│ │ └─┤ ────┼─┼───────────────────┬── ┼───┼ │ ┬ │
│ │ │ ┬───┼─┼───────────────────┼── ┼─┬─┼ │ ┼ │
│ │ │ │ ─ ┼─┼─┬─────┬──── ──────┼─┬ │ ├─┘ └─┤ │
│ │ │ │ ┬ └─┤ │ ┬─┬ ┼─┬─┬ ──┬───┼─┼ ├─┘ ├─┘
│ │ │ └─┤ └─┤ └─┤ │ ├─┘ ──┼─┬─┼─┼ │ │
│ │ │ │ │ │ ├─┘ ┬─┼─┼─┼─┼ │ │
│ │ │ │ │ ├─┘ └─┤ │ ├─┘ │ │
│ │ │ │ └───┤ │ ├─┘ │ │
│ │ │ │ │ ├─┘ │ │
│ │ │ │ ├─────────┘ │ │
│ │ │ ├─────────┘ │ │
│ │ └───┤ │ │
│ │ ├───────────────────────────┘ │
│ ├───────┘ │
└─┤ │
└───────────────────────────────────────────┘
λ 1 (1 ((λ 1 1) (λ λ λ 1 (λ λ 1) ((λ 4 4 1 ((λ 1 1) (λ 2 (1 1)))) (λ λ λ λ 1 3 (2 (6 4))))) (λ λ λ 4 (1 3)))) (λ λ 1 (λ λ 2) 2)
[1] https://tromp.github.io/cl/cl.html
[2] https://github.com/tromp/AIT/blob/master/characteristic_sequ...
> On Pi day 2023 I gave an online talk about AIT and BLC based on these slides.
Perhaps watching the linked-to talk will somewhat clarify matters.
12: [3, 2]
25: [5]
49: [7]
The Sorenson paper linked in the post is also a great read, I remember reading it a couple of years ago and it is very clear. It adds many further optimizations on top of standard sieve of Eratosthenes (e.g. wheel factorization is a generalization of the idea of considering only numbers in {1, 5} mod 6, or {1, 7, 11, 13, 17, 19, 23, 29} mod 30, or …), and segmenting is what you need to do to avoid running out of memory.
The fastest such generator (using sieve-based methods) seems to be https://github.com/kimwalisch/primesieve — it implements all those optimizations and more.
However, if you just want the nth prime, it's even faster to compute π(n) by combinatorial/analytic methods, instead of a sieve: the same author has https://github.com/kimwalisch/primecount for that. (Try the ten-trillionth prime number.)
It doesn't quite produce primes directly. It gives a sequence of integers, some of which are powers of 2. Those powers of 2 are 2^2, 2^3, 2^5, 2^7, 2^11, ..., i.e., the sequence of 2 to the power of primes.
For those who have never seen FRACTRAN it is an esoteric programming language invented by John Conway. It is Turing-complete.
Here's how you execute a FRACTRAN program, with input N.
1. Step through the list of fractions in order until you find a fraction f such that fN is an integer or you run out of fractions.
2. If you run out of fractions halt.
3. Output fN.
4. Replace N with fN.
5. Goto step 1.
Wikipedia gives some sample programs and explains in detail how the heck FRACTRAN can compute: https://en.wikipedia.org/wiki/FRACTRAN
p_n = 1 + Sum_{i=1}^{2^n} Floor((n/(Sum_{j=1}^{i} Floor((Cos((j-1)!+1)/j) Pi)^2)))^(1/n))
https://en.wikipedia.org/wiki/Formula_for_primes
https://www.youtube.com/watch?v=j5s0h42GfvM
It's not fast but it really does work.
The last time this formula came up, I added some references to the Wikipedia article (current version: https://en.wikipedia.org/w/index.php?title=Formula_for_prime... about the "worthlessness" of such formulas, but it doesn't deter people from being interested in that formula… as long as they don't imagine it's useful in any way, it's fine to enjoy its perverse cleverness, I suppose.
The video does well to quote Wilf on “the disease was preferable to the cure” and to end with “Could it [C. P. Willans] be a pseudonym for someone who didn't want to tarnish their reputation by writing about useless prime-generating formulas?” :-) Dudley (who mentions a large list of other such formulas) also has a nice conclusion:
> The conclusion to be drawn from all this, I think, is that formulas for formulas' sake do not advance the mathematical enterprise. Formulas should be useful. If not, they should be astounding, elegant, enlightening, simple, or have some other redeeming value.
https://gist.github.com/scythe/5fb962722934c58c60430180beab8...
In the spirit of the blog post, I'll let you guess how it works :p
q = 101
D = {'102': [3, 2], '105': [7, 5], '121': [11], '169': [13], '289': [17], '361': [19], '529': [23], '841': [29], '961': [31], '1369': [37], '1681': [41], '1849': [43], '2209': [47], '2809': [53], '3481': [59], '3721': [61], '4489': [67], '5041': [71], '5329': [73], '6241': [79], '6889': [83], '7921': [89], '9409': [97]}
q = 101
S = 121
r = 11
D = {'102': [3, 2], '105': [7, 5]}https://github.com/llvm-mirror/libcxx/blob/78d6a7767ed57b501...
There is no good reason to use this one except in a code golf environment that includes all headers by default, which is where I learned about it.
It asymptotically speeds up the process of generating prime numbers. It has time complexity O(n/(log log n))
https://www.ams.org/journals/mcom/2004-73-246/S0025-5718-03-...
They are essentially the same in all conceivable practical cases in this Universe.
log log (atoms in universe) < 100
It is so fast, I wonder what the algorithm behind it is.
p: 100999999, the prime numbers is 2059519669
Noted in the details section of the link you provided:
"Primality testing on numbers larger than (2^31) uses the probabilistic Miller-Rabin algorithm."
I understand that not using yield would mean it is not strictly a "generator", but the spirit of this problem is generating the sequence of primes (unless I am missing something).
But one does not always end up creating a large list of prime numbers, do they?
> why not just return a large list of prime numbers?
Because that requires upfront knowledge of the filtering factor or absence thereof.
A generator means you don't care, you generate an infinite sequence and the consumer is free to do whatever they want.
> It is obvious that returning an array of prime numbers is what one would do in C/C++/Rust.
It's certainly not obvious for Rust.
Everything I replied to is the opinion you personally expressed.
And note that the author had already preempted your comment in the first place:
> When we want a list of all the primes below some known limit, gen_primes_upto is great, and performs fairly well. There are two issues with it, though:
> We have to know what the limit is ahead of time; this isn't always possible or convenient.
primes = gen_primes() for item, prime in zip(some_iterator, primes):
With the right algorithm you can solve it in Python in a max of 30 seconds. With the wrong one, you sometimes can't solve it in C++ in the lifetime of the universe.
Therefore you should program in the language that you find it easiest to express yourself in. And not in a language chosen for raw speed.
You'll see the win pretty quickly when you have to search a large range for primes.
This isn't just shifting the goal posts, you are now talking about something completely different.
With the right algorithm you can solve it in Python in a max of 30 seconds
Then you can do the same thing in C++, but when it runs the C++ program will be done 100x faster. A lot of modern C++ is as clean as python if you were being more explicit about types.
With the wrong one, you sometimes can't solve it in C++ in the lifetime of the universe.
Nothing about this makes any sense.
Therefore you should program in the language that you find it easiest to express yourself in. And not in a language chosen for raw speed.
If you continue on with python and need more speed you will run around in circles trying to get it while the C++ version is already done. Anything where speed can be a bottleneck is not fit to be written in pure python.
This pattern plays out over and over again. If something is resource intensive, a scripting language is not a long term solution.
Absolutely none of your arguments make sense in a context of a series of programming puzzles of fixed size. Back in the real world, there is a real need for exploratory programming on mathematical problems. I' prefer to use Python for that. If your problems are statistical, R is a better choice. C++ is a poor fit for exploratory programming.
C++ is a good choice for serious mathematical computing. But not always the best. https://www.hpcwire.com/2020/01/14/julia-programmings-dramat... shows the growing popularity of Julia, even for the most demanding computations.
Looking to the future, I would choose Rust over C++ for most new projects where people currently use C++. It seems bizarre to me that C++ evangelists can talk about type safety and ignore memory safety, when memory safety is a giant problem for security and correctness in any large program.
So, yeah. C++ is great for what it's great for. But if I'm going some exploratory back of the envelope calculations, I'm still reaching for Python. And yield is one of the reasons why.
It works fine for me. Are you sure this isn't just because you don't have a lot of practice with C++?
It seems bizarre to me that C++ evangelists can talk about type safety and ignore memory safety, when memory safety is a giant problem for security and correctness in any large program.
This has nothing to do with the current thread, but this isn't really a problem in modern C++. Are you you sure you are up to date with modern C++? Most programs don't look too different than python due to for iteration and auto type deduction.
I have no idea why you are trying to evangelize C++ here. I use a wide variety of languages for different purposes. For Project Euler type tasks I switched to Python about 15 years ago because of yield, and the convenience of using objects as hash keys. If you make different choices and they work for you, then have fun.
You seem to be saying there is some advantage to python because of the yield feature and I just don't think that's true. You can keep state in a class and call a method. There are a lot of disadvantages though.
triple sigh
I'd say the flexibility to not need to know how big of a prime you'll need is the main draw.
The programming flexibility is indeed the win.
def primes ():
yield 2
yield 3
f = {}
n = 5
for p in primes():
if p == 2:
continue
p2 = p*p
while n < p2:
if n in f:
double_p = f.pop(n)
m = n + double_p
while m in f:
m += double_p
f[m] = double_p
else:
yield n
n += 2
f[n] = 2*pTo do that it recursively contains a second prime generator that only produces up to the square root. Which contains one that produces up to the fourth root. And so on until you get down to 3. The work and memory for those recursive generators is a rounding error compared to the initial prime generator. The work and memory saved by not keeping unnecessary primes is a significant win.
When writing ‘fizz buzz’, would you expect to return an array of 100 strings, some of which are numbers, some ‘fizz’es and some ‘buzz’es?
Maybe I'm misunderstanding you, but Rust has native generators with yield: https://doc.rust-lang.org/beta/unstable-book/language-featur...
C++ also does using coroutines + co_yield (https://en.cppreference.com/w/cpp/language/coroutines#co_yie...) or iterators (old syntax).
I can’t say I have a favourite but I do enjoy multiplying two random prime numbers with each other and using the random result when I need a random number with added random.
To add maximum random, sometimes I like to return a string with a random Gertrude Stein quote, segfault, or make the PC speaker beep instead of returning a number. So random!
Teaching my kids coding, we also arrived at most of the other optimizations, in Janet:
(defn primes-sieve []
(fiber/new (fn []
(yield 2) # the only even prime, yield it explicitly.
(let [sieve (table/new (math/pow 2 12))]
(loop [n :range [3 nil 2]] # only odd numbers
(if-let [current-composite (sieve n)]
(do
(each prime-factor current-composite
(let [next-composite-key (+ n (* 2 prime-factor))] # add 2*prime-factor, else composite would be even
(if-let [next-composite (sieve next-composite-key)]
(array/push next-composite prime-factor)
(put sieve next-composite-key @[prime-factor]))))
(put sieve n nil))
(do
(yield n)
(put sieve (* n n) @[n])) # next composite is n^2. Anything less is composite a smaller prime, already in the sieve
)))))) import itertools as it
def erat3( ):
D = { 9: 3, 25: 5 }
yield 2
yield 3
yield 5
MASK= 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0,
MODULOS= frozenset( (1, 7, 11, 13, 17, 19, 23, 29) )
for q in it.compress(
it.islice(it.count(7), 0, None, 2),
it.cycle(MASK)):
p = D.pop(q, None)
if p is None:
D[q*q] = q
yield q
else:
x = q + 2*p
while x in D or (x%30) not in MODULOS:
x += 2*p
D[x] = pHardcoding those values is a semi-requirement of the above technique.
d = defaultdict(list)
# Magic
for i in range(2, 10000):
for y in i in d and d.pop(i) or [i]:
d[i+y] += (y,)
# Not sure how to print this in a succinct way
print(sorted(list(set().union(*d.values())))) from collections import defaultdict
d = defaultdict(list)
for i in range(2, 10000):
for y in d.pop(i, [i]):
d[i + y].append(y)
for y in sorted(y for x in d.values() for y in x):
print(y)It basically uses the priority queue to store all the composites up to a certain point and then you can assume the next one is prime.
I had some (admittedly not well written) code for it too https://gist.github.com/weskerfoot/4699275
It can also incorporates the wheel factorization optimization mentioned in this article.
Memory Usage is greatly reduced by using BitSet.
https://github.com/krnaveen14/PrimeNumbers/blob/master/src/m...
https://codegolf.stackexchange.com/a/66915/41984
OMG, it's been 8 years already!
The Genuine Sieve of Eratosthenes https://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf
even a little research into cryptography will inform you that prime numbers are used heavily in encryption algorithms.
https://stackoverflow.com/questions/439870/why-are-primes-im...
"even a little research" will tell you people shouldn't roll your own cryptography and to just use a library.
So, this still doesn't answer my question.
but to counter your point: people do have to write those libraries that you depend on. others have to validate them, and yet still others are doing their best efforts to invalidate them and attack them.
but you are correct that most people should not roll their own cryptography suites.
* https://github.com/intc/c-notes/blob/master/threads/4-thrd_e...
{(~T∊T∘.×T)/T←1↓⍳⍵} 25
{T~∘∊∘.×⍨T←1↓⍳⍵} 25
((⊢~∘∊∘.×⍨)1↓⍳) 25[1] https://en.wikipedia.org/wiki/APL_(programming_language)
First 15 integers:
⍳ 15
┌→──────────────────────────────────┐
│1 2 3 4 5 6 7 8 9 10 11 12 13 14 15│
└~──────────────────────────────────┘
1↓⍳15
┌→────────────────────────────────┐
│2 3 4 5 6 7 8 9 10 11 12 13 14 15│
└~────────────────────────────────┘
T←1↓⍳15
T∘.×T
┌→─────────────────────────────────────────────────┐
↓ 4 6 8 10 12 14 16 18 20 22 24 26 28 30│
│ 6 9 12 15 18 21 24 27 30 33 36 39 42 45│
│ 8 12 16 20 24 28 32 36 40 44 48 52 56 60│
│10 15 20 25 30 35 40 45 50 55 60 65 70 75│
│12 18 24 30 36 42 48 54 60 66 72 78 84 90│
│14 21 28 35 42 49 56 63 70 77 84 91 98 105│
│16 24 32 40 48 56 64 72 80 88 96 104 112 120│
│18 27 36 45 54 63 72 81 90 99 108 117 126 135│
│20 30 40 50 60 70 80 90 100 110 120 130 140 150│
│22 33 44 55 66 77 88 99 110 121 132 143 154 165│
│24 36 48 60 72 84 96 108 120 132 144 156 168 180│
│26 39 52 65 78 91 104 117 130 143 156 169 182 195│
│28 42 56 70 84 98 112 126 140 154 168 182 196 210│
│30 45 60 75 90 105 120 135 150 165 180 195 210 225│
└~─────────────────────────────────────────────────┘
T∊T∘.×T
┌→──────────────────────────┐
│0 0 1 0 1 0 1 1 1 0 1 0 1 1│
└~──────────────────────────┘
T,[0.5]T∊T∘.×T
┌→────────────────────────────────┐
↓2 3 4 5 6 7 8 9 10 11 12 13 14 15│
│0 0 1 0 1 0 1 1 1 0 1 0 1 1│
└~────────────────────────────────┘
~T∊T∘.×T
┌→──────────────────────────┐
│1 1 0 1 0 1 0 0 0 1 0 1 0 0│
└~──────────────────────────┘
(~T∊T∘.×T)/T
┌→────────────┐
│2 3 5 7 11 13│
└~────────────┘
&2=+/~x!/:\:x:!100