How to get a prime collection within 2^63 - 1 in Python

Question

I want to get a prime number set within 2^63 - 1 in Python,I have seen the following code on the web:

limit = 2**63 - 2
p = [True] * (limit + 1)
# p = bitarray(limit + 1)
p[0] = p[1] = False
for i in range(2, int(math.sqrt(limit) + 1)):
    if p[i]:
        for j in range(i * i, limit + 1, i):
            p[j] = False
prime = [i for i in range(limit + 1) if p[i]]
print(prime)

But when I run this program, the compiler complains that can not fit 'int' into an index-sized integer. I tried to solve the problem with bitarray, but the index of the array is still too big.

Added `import math` as first line and set a lower limit (e.g. 100) ... then the program works perfectly. With `limit = 2 ** 63 - 1` I'm getting a MemoryError ... — quant, Nov 16 '18 at 08:29
This looks like the Sieve of Eratosthenes. It's super quick, but very memory hungry. Generating primes up to 2^63 would take millions of Gigabytes of memory using the naive implementation. — DatHydroGuy, Nov 16 '18 at 08:31
what does "within 2^63-1" mean? There are too many primes < 2^63 to enumerate all of them with current technology. — President James K. Polk, Nov 16 '18 at 18:14

score 1 · Answer 1 · answered Nov 16 '18 at 08:24

1

You can use sympy:

import sympy

print(list(sympy.primerange(0,2**63-1)))

but as 2^63 is quite large this will take some time.

answered Nov 16 '18 at 08:24

T. Ewen

126
4

1

How long you think it will take to compute and print that list 216289611853439384 integers? – President James K. Polk Nov 16 '18 at 21:57
1

Longer than anyone would be able to wait, not to mention the amount of memory you would need. – T. Ewen Nov 19 '18 at 08:08

quant · Accepted Answer · 2018-11-16T08:57:13.983

1

You can use the following code. It is using the Sieve of Eratosthenes in combination with a generator function in order to reduce the memory usage of this algorithm. It is furthermore exploiting the less commonly known fact that every prim number > 4 can be written as 6*n ± 1.

import math

limit = 2 ** 63 - 1

def isPrim(n, belowPrims):
    limit = int(math.sqrt(n))
    for prim in belowPrims:
        if prim > limit: break
        if n % prim == 0: return False 
    return True

def calcPrims():
    yield 2
    yield 3
    toTest, nextN = [], 6
    while True:
        p1 = nextN - 1
        if isPrim(p1, toTest):
            yield p1
            toTest.append(p1)
        p2 = nextN + 1
        if isPrim(p2, toTest):
            yield p2
            toTest.append(p2)
        nextN += 6

for prim in calcPrims():
    if prim > limit:
        break
    print(prim)

Edit

This link here https://primes.utm.edu/notes/faq/six.html explains briefly why every prim number can be written in the form 6*n ± 1.

edited Nov 16 '18 at 08:57

answered Nov 16 '18 at 08:43

quant

2,184
2
19
29

1

@PatrickArtner Oh f**k you are right, the correct limit is 3 not 2. I updated the code accordingly. .. In case someone is interested why every prim number can be written in the form 6n+-1, this link here explains it briefly https://primes.utm.edu/notes/faq/six.html – quant Nov 16 '18 at 08:50
1

Add the link to your answer - I dont tell if you dont tell about the `3` – Patrick Artner Nov 16 '18 at 08:51
2

@quant Please use better language here. – Blue Nov 16 '18 at 08:53
If I want to generate a prime number set between 5e9 and 2 ** 63 - 1, how should your code be modified? – fanfei Nov 16 '18 at 12:29
@fanfei The Sieve of Eratosthenes algorithm requires you to know the prim numbers between 0 and 5e9 in order to calculate the prim numbers above. So there is no shortcut by not calculating these prim numbers. However you can modify the line `print(prim)` to `if prim < 5e9: print(prim)` in order to stop outputting them. – quant Nov 16 '18 at 12:41

hiro protagonist · Answer 3 · 2018-11-16T09:15:02.310

if you have a primes() generator of some kind, you could do this:

is_prime_var = 0

MAX = 1 << 5
last_p = 0
for p in primes():
    if p > MAX:
        break
    print(p, p-last_p)
    is_prime_var <<= (p - last_p)
    is_prime_var |= 1
    last_p = p
is_prime_var <<= (MAX - last_p - 1)

now the locations of the primes are stored (in reversed order) in the integer is_prime_var.

then the expression (is_prime >> (MAX-n-1)) & 1 would be 1 if n is prime; 0 otherwise:

def is_prime(n):
    return bool((is_prime_var >> (MAX-n-1)) & 1)

you could use primes() from this excellent answer as prime generator.

thers is also this answer of mine about a fast and memory-efficient sieve of eratosthenes. might also be interesting.

How to get a prime collection within 2^63 - 1 in Python

3 Answers3