[HN Gopher] My favorite prime number generator
       ___________________________________________________________________
        
       My favorite prime number generator
        
       Author : mfrw
       Score  : 137 points
       Date   : 2023-08-23 14:05 UTC (8 hours ago)
        
 (HTM) web link (eli.thegreenplace.net)
 (TXT) w3m dump (eli.thegreenplace.net)
        
       | svat wrote:
       | This is really elegant! To make the explanation more concrete,
       | here's a snapshot of the algorithm's state (the dict D) after it
       | has processed numbers up to 10 (inclusive):                   12:
       | [3, 2]         25: [5]         49: [7]
       | 
       | So each prime up to 10 is present exactly once. When we process
       | 12, we move 2 to the list for 14, and 3 to the list for 15.
       | 
       | 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 p(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.)
        
       | SeanLuke wrote:
       | Believe it or not, there exists a closed-form equation which
       | computes the nth prime number. It's Willans's formula:
       | 
       | 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.
        
       | krnaveen14 wrote:
       | Reminds me of the time where I've implemented highly optimised
       | version of Sieve of Eratosthenes in Java, just for a Code Golf
       | Challenge.
       | 
       | 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!
        
       | scythe wrote:
       | The memory usage can in fact be optimized further with a small
       | change:
       | 
       | https://gist.github.com/scythe/5fb962722934c58c60430180beab8...
       | 
       | In the spirit of the blog post, I'll let you guess how it works
       | :p
        
         | svat wrote:
         | Wow this is a great optimization! You're under-selling it :) It
         | changes the memory requirement from O(p(n)) to O(p([?]n)), e.g.
         | after processing numbers up to 100, the internal state in the
         | OP's algorithm is:                   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]}
         | 
         | while with your "small change", the state is just:
         | q = 101         S = 121         r = 11         D = {'102': [3,
         | 2], '105': [7, 5]}
        
       | tromp wrote:
       | My favorite prime number generator is
       | +-+-----------------------------------+---- ----+         | | +-+
       | ----+-+---------------------- +---+ +---+         | | +-+
       | ----+-+-------------------+-- +---+ | + |         | |   |
       | +---+-+-------------------+-- +-+-+ | + |         | |   | | -
       | +-+-+-----+---- ------+-+ | +-+ +-+ |         | |   | | + +-+ |
       | +-+ +-+-+ --+---+-+ +-+     +-+         | |   | +-+   +-+ +-+ |
       | +-+ --+-+-+-+ |       |         | |   |   |     |   | +-+
       | +-+-+-+-+ |       |         | |   |   |     |   +-+     +-+ | +-+
       | |       |         | |   |   |     +---+         | +-+   |       |
       | | |   |   |         |         +-+     |       |         | |   |
       | |         +---------+       |       |         | |   |
       | +---------+                 |       |         | |   +---+
       | |       |         | |       +---------------------------+       |
       | | +-------+                                   |         +-+
       | |           +-------------------------------------------+
       | 
       | which graphically depicts the 167-bit lambda term
       | l 1 (1 ((l 1 1) (l l l 1 (l l 1) ((l 4 4 1 ((l 1 1) (l 2 (1 1))))
       | (l l l l 1 3 (2 (6 4))))) (l l l 4 (1 3)))) (l l 1 (l l 2) 2)
       | 
       | which generates the characteristic sequence of primes [1] [2].
       | 
       | [1] https://tromp.github.io/cl/cl.html
       | 
       | [2]
       | https://github.com/tromp/AIT/blob/master/characteristic_sequ...
        
         | MrJeffUtah wrote:
         | [dead]
        
         | pixelpoet wrote:
         | Super cool to see one of the authors of the most popular Go
         | scoring methods (particularly for MCTS) posting here :) Respect
         | and thanks!
        
         | hgsgm wrote:
         | What's the "characteristic sequence"?
        
           | hidroto wrote:
           | its a list with each prime index contains a one/true value
        
         | cobaltoxide wrote:
         | Could you explain it? Those two links are Greek to me.
        
           | tromp wrote:
           | The first link starts out as:
           | 
           | > 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.
        
       | philipLutz wrote:
       | I enjoyed this, but I am a little stuck on how stuck the author
       | is on yield. If one ends up creating a large list of prime
       | numbers with this generator, why not just return a large list of
       | prime numbers? The author clearly notes how optimizing Python
       | code is "odd" when you can more easily switch languages, but
       | languages like C/C++/Rust would not have pretty code because it
       | does not have yield. It is obvious that returning an array of
       | prime numbers is what one would do in C/C++/Rust.
       | 
       | 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).
        
         | btilly wrote:
         | Try solving some Project Euler problems. You'll see the win
         | pretty quickly when you have to search a large range for
         | primes.
        
           | CyberDildonics wrote:
           | Why would that be true?
        
             | btilly wrote:
             | The list you need to process is large enough that it would
             | be prohibitive to put it all in memory at once. Furthermore
             | you do not start knowing things like how many primes to
             | look at. Therefore generating and dealing with them
             | iteratively makes for simple code and a lot fewer
             | headaches.
        
         | jameshart wrote:
         | If the problem is _outputting_ the primes, then _storing_ them
         | in a list is not a necessary part of that problem.
         | 
         | 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?
        
         | vlovich123 wrote:
         | > but languages like C/C++/Rust would not have pretty code
         | because it does not have yield
         | 
         | 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.cpprefere
         | nce.com/w/cpp/language/coroutines#co_yie...) or iterators (old
         | syntax).
        
         | masklinn wrote:
         | > If one ends up creating a large list of prime numbers with
         | this generator
         | 
         | 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.
        
           | philipLutz wrote:
           | "I can't try to rewrite it in C++ or Rust for now, due to the
           | lack of generator support; the yield statement is what makes
           | this code so nice and elegant, and alternative idioms are
           | much less convenient."
        
             | masklinn wrote:
             | That is the context of the first sentence in your comment,
             | pretty much the only part I left alone.
             | 
             | 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.
        
         | sp332 wrote:
         | This implementation runs until it crashes. So you need a way to
         | output the values as you go, or they will be lost when you run
         | out of memory. It also lets you control how many values are
         | generated if you just want to see a few without actually using
         | all the RAM.
        
         | Vexs wrote:
         | It's a bit of a convolution, but let's say we've got some
         | iterator of unknown size and we want to have a prime each
         | iteration. We can do the following in python, which saves us
         | from having to pre-compute the primes _or_ know how long
         | `some_iterator` is.
         | 
         | primes = gen_primes() for item, prime in zip(some_iterator,
         | primes):
        
       | naillo wrote:
       | This person has an incredible work ethic
        
       | rullelito wrote:
       | My favorite one:                   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()))))
        
         | paulddraper wrote:
         | Beautiful                 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)
        
       | manicennui wrote:
       | I've always enjoyed the Perl regex solution.
        
       | omoikane wrote:
       | My favorite prime number generator is the undocumented
       | __next_prime():
       | 
       | 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.
        
       | tzs wrote:
       | My favorite is this FRACTRAN program: (17/91, 78/85, 19/51,
       | 23/38, 29/33, 77/29, 1/17, 11/13, 13/11, 15/2, 1/7, 55/1) with
       | input 2.
       | 
       | 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
        
       | oxxoxoxooo wrote:
       | And another one:
       | 
       | The Genuine Sieve of Eratosthenes
       | https://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf
        
       | tzot wrote:
       | The fastest variant I know of this algorithm is the one found at
       | https://stackoverflow.com/a/3796442/6899 :
       | 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] = p
        
         | svat wrote:
         | This adds two things over the OP's initial version, both
         | mentioned in the post: wheel factorization (with 30), and the
         | "linear probing" to avoid lists. It's possible to keep
         | improving (up to a point) by using e.g. wheel with 210 instead
         | of 30, and so on, if one doesn't mind hard-coding more and
         | more.
        
         | simlevesque wrote:
         | isn't hardcoding 2 3 and 5 a bit of a cheat ? It seems like you
         | could hardcode it to yield the X next primes and it'd be
         | faster.
        
           | stouset wrote:
           | The point isn't hardcoding those numbers. The point is that
           | you don't iterate on any _multiples_ of those numbers (which
           | we know are composite) and doing allows us to eliminate 2 /3
           | of the iterations. It's an extension of enumerating over only
           | odd numbers, which eliminates 1/2 of the iterations.
           | 
           | Hardcoding those values is a semi-requirement of the above
           | technique.
        
             | simlevesque wrote:
             | Thank you, it makes more sense now.
        
       | shaftoe444 wrote:
       | My favourite prime number generator is this one using go channels
       | to create filters. Not efficient but fun.
       | 
       | https://go.dev/play/p/9U22NfrXeq
        
       | tantalor wrote:
       | My favorite prime number generator is
       | http://alpha61.com/primenumbershittingbear/
        
         | proactivesvcs wrote:
         | It doesn't need Flash Player any more? It's the year 2023 and I
         | can run the Prime Number Shitting Bear again? Hallelujah!
        
           | [deleted]
        
       | redocneknurd wrote:
       | What are you using your prime numbers for?
        
       | dvh wrote:
       | {(~T[?]T[?].xT)/T-1|[?][?]}
        
         | rgrau wrote:
         | Here's a solution in k, that uses a similar array oriented
         | approach:                   &2=+/~x!/:\:x:!100
         | 
         | I wrote a little explanation of it, in case anyone is curious
         | about how it works: https://github.com/kidd/arraylangs-
         | index/blob/master/noteboo...
        
         | jodrellblank wrote:
         | {(~T[?]T[?].xT)/T-1|[?][?]} 25
         | {T~[?][?][?].x[?]T-1|[?][?]} 25
         | (([?]~[?][?][?].x[?])1|[?]) 25
        
         | ducttapecrown wrote:
         | What language is that? AWK??
        
           | jodrellblank wrote:
           | It's an APL sieve, for an input of 15 (primes up to 15) it
           | does these steps:
           | 
           | First 15 integers:                         [?] 15
           | +-----------------------------------+         |1 2 3 4 5 6 7
           | 8 9 10 11 12 13 14 15|
           | +~----------------------------------+
           | 
           | 1 element dropped from the front of that:
           | 1|[?]15         +---------------------------------+
           | |2 3 4 5 6 7 8 9 10 11 12 13 14 15|
           | +~--------------------------------+
           | 
           | That array stored in variable T:
           | T-1|[?]15
           | 
           | A table of outer-product multiplication, T by T:
           | T[?].xT
           | +--------------------------------------------------+
           | | 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|
           | +~-------------------------------------------------+
           | 
           | The elements of T which exist (are found) in the
           | multiplication table, boolean result:
           | T[?]T[?].xT         +---------------------------+         |0
           | 0 1 0 1 0 1 1 1 0 1 0 1 1|
           | +~--------------------------+
           | 
           | Showing that together with the input, 4 6 8 are found, 2 3 5
           | are not:                         T,[0.5]T[?]T[?].xT
           | +---------------------------------+         |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|         +~--------------------------------+
           | 
           | Numbers not found in the multiplication table are prime, so
           | invert the test result so primes get 1 and composites get 0:
           | ~T[?]T[?].xT         +---------------------------+         |1
           | 1 0 1 0 1 0 0 0 1 0 1 0 0|
           | +~--------------------------+
           | 
           | Use that to 'compress' the input to keep things where the 1
           | positions are and drop things where the 0s are:
           | (~T[?]T[?].xT)/T         +-------------+         |2 3 5 7 11
           | 13|         +~------------+
           | 
           | You can play with it at https://tryapl.org/ using the bar at
           | the top to enter the symbols.
        
           | paulddraper wrote:
           | It's got funny symbols, so I assume APL. [1]
           | 
           | [1] https://en.wikipedia.org/wiki/APL_(programming_language)
        
       | Exuma wrote:
       | He said he has "used this code many times"... have you ever
       | needed to find prime numbers in your daily programming?
        
         | thrownblown wrote:
         | [flagged]
        
       | nemo8551 wrote:
       | Can I just say how surprised I am the many people have a
       | favourite prime number generator.
       | 
       | 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.
        
         | cobaltoxide wrote:
         | That doesn't seem like a good way to generate "random" numbers,
         | since it will only produce numbers with exactly two prime
         | factors.
        
           | amenghra wrote:
           | Correct. It will almost never produce an even number, so
           | there's already one bit of entropy lost. It will also never
           | produce lots of other numbers.
        
       | max_ wrote:
       | In J generating prime numbers is simply p:21, that gives you 79
       | which is the 21st prime number.
       | 
       | It is so fast, I wonder what the algorithm behind it is.
       | 
       | [0]: https://code.jsoftware.com/wiki/Vocabulary/pco
        
         | noman-land wrote:
         | How high can you go with this?
        
           | max_ wrote:
           | I have J on my android phone, and the highest I could go was
           | 
           | p: 100999999, the prime numbers is 2059519669
        
         | jonahx wrote:
         | > I wonder what the algorithm behind it is.
         | 
         | Noted in the details section of the link you provided:
         | 
         | "Primality testing on numbers larger than (2^31) uses the
         | probabilistic Miller-Rabin algorithm."
        
       | ChrisCinelli wrote:
       | The Sieve of Atkin is VERY fast at generating primes between 1
       | and n.
       | 
       | It asymptotically speeds up the process of generating prime
       | numbers. It has time complexity O(n/(log log n))
       | 
       | https://www.baeldung.com/cs/prime-number-algorithms
        
         | anonymoushn wrote:
         | Your link seems to say it has time complexity O(n)?
        
           | hgsgm wrote:
           | If you count everything carefully, it's 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
        
             | anonymoushn wrote:
             | This seems to be a matter of the above source reporting a
             | variation of the algorithm that is asymptotically slower
             | than the algorithm given in the paper *and* uses
             | asymptotically more memory. Thanks for posting the paper!
        
       | anaphor wrote:
       | There's a functional / lazy version I like that uses a priority
       | queue in order to get the next prime number
       | https://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf
       | 
       | 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.
        
       ___________________________________________________________________
       (page generated 2023-08-23 23:00 UTC)