Fastest way to list all primes below N
Warning: timeit results may vary due to differences in hardware or
version of Python.
Below is a script which compares a number of implementations:
- ambi_sieve_plain,
- rwh_primes,
- rwh_primes1,
- rwh_primes2,
- sieveOfAtkin,
- sieveOfEratosthenes,
- sundaram3,
- sieve_wheel_30,
- ambi_sieve (requires numpy)
- primesfrom3to (requires numpy)
- primesfrom2to (requires numpy)
Many thanks to stephan for bringing sieve_wheel_30 to my attention.
Credit goes to Robert William Hanks for primesfrom2to, primesfrom3to, rwh_primes, rwh_primes1, and rwh_primes2.
Of the plain Python methods tested, with psyco, for n=1000000,
rwh_primes1 was the fastest tested.
+---------------------+-------+
| Method | ms |
+---------------------+-------+
| rwh_primes1 | 43.0 |
| sieveOfAtkin | 46.4 |
| rwh_primes | 57.4 |
| sieve_wheel_30 | 63.0 |
| rwh_primes2 | 67.8 |
| sieveOfEratosthenes | 147.0 |
| ambi_sieve_plain | 152.0 |
| sundaram3 | 194.0 |
+---------------------+-------+
Of the plain Python methods tested, without psyco, for n=1000000,
rwh_primes2 was the fastest.
+---------------------+-------+
| Method | ms |
+---------------------+-------+
| rwh_primes2 | 68.1 |
| rwh_primes1 | 93.7 |
| rwh_primes | 94.6 |
| sieve_wheel_30 | 97.4 |
| sieveOfEratosthenes | 178.0 |
| ambi_sieve_plain | 286.0 |
| sieveOfAtkin | 314.0 |
| sundaram3 | 416.0 |
+---------------------+-------+
Of all the methods tested, allowing numpy, for n=1000000,
primesfrom2to was the fastest tested.
+---------------------+-------+
| Method | ms |
+---------------------+-------+
| primesfrom2to | 15.9 |
| primesfrom3to | 18.4 |
| ambi_sieve | 29.3 |
+---------------------+-------+
Timings were measured using the command:
python -mtimeit -s"import primes" "primes.{method}(1000000)"
with {method} replaced by each of the method names.
primes.py:
#!/usr/bin/env python
import psyco; psyco.full()
from math import sqrt, ceil
import numpy as np
def rwh_primes(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Returns a list of primes < n """
sieve = [True] * n
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i]:
sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
return [2] + [i for i in xrange(3,n,2) if sieve[i]]
def rwh_primes1(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Returns a list of primes < n """
sieve = [True] * (n/2)
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i/2]:
sieve[i*i/2::i] = [False] * ((n-i*i-1)/(2*i)+1)
return [2] + [2*i+1 for i in xrange(1,n/2) if sieve[i]]
def rwh_primes2(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Input n>=6, Returns a list of primes, 2 <= p < n """
correction = (n%6>1)
n = {0:n,1:n-1,2:n+4,3:n+3,4:n+2,5:n+1}[n%6]
sieve = [True] * (n/3)
sieve[0] = False
for i in xrange(int(n**0.5)/3+1):
if sieve[i]:
k=3*i+1|1
sieve[ ((k*k)/3) ::2*k]=[False]*((n/6-(k*k)/6-1)/k+1)
sieve[(k*k+4*k-2*k*(i&1))/3::2*k]=[False]*((n/6-(k*k+4*k-2*k*(i&1))/6-1)/k+1)
return [2,3] + [3*i+1|1 for i in xrange(1,n/3-correction) if sieve[i]]
def sieve_wheel_30(N):
# http://zerovolt.com/?p=88
''' Returns a list of primes <= N using wheel criterion 2*3*5 = 30
Copyright 2009 by zerovolt.com
This code is free for non-commercial purposes, in which case you can just leave this comment as a credit for my work.
If you need this code for commercial purposes, please contact me by sending an email to: info [at] zerovolt [dot] com.'''
__smallp = ( 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139,
149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311,
313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401,
409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491,
499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683,
691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887,
907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997)
wheel = (2, 3, 5)
const = 30
if N < 2:
return []
if N <= const:
pos = 0
while __smallp[pos] <= N:
pos += 1
return list(__smallp[:pos])
# make the offsets list
offsets = (7, 11, 13, 17, 19, 23, 29, 1)
# prepare the list
p = [2, 3, 5]
dim = 2 + N // const
tk1 = [True] * dim
tk7 = [True] * dim
tk11 = [True] * dim
tk13 = [True] * dim
tk17 = [True] * dim
tk19 = [True] * dim
tk23 = [True] * dim
tk29 = [True] * dim
tk1[0] = False
# help dictionary d
# d[a , b] = c ==> if I want to find the smallest useful multiple of (30*pos)+a
# on tkc, then I need the index given by the product of [(30*pos)+a][(30*pos)+b]
# in general. If b < a, I need [(30*pos)+a][(30*(pos+1))+b]
d = {}
for x in offsets:
for y in offsets:
res = (x*y) % const
if res in offsets:
d[(x, res)] = y
# another help dictionary: gives tkx calling tmptk[x]
tmptk = {1:tk1, 7:tk7, 11:tk11, 13:tk13, 17:tk17, 19:tk19, 23:tk23, 29:tk29}
pos, prime, lastadded, stop = 0, 0, 0, int(ceil(sqrt(N)))
# inner functions definition
def del_mult(tk, start, step):
for k in xrange(start, len(tk), step):
tk[k] = False
# end of inner functions definition
cpos = const * pos
while prime < stop:
# 30k + 7
if tk7[pos]:
prime = cpos + 7
p.append(prime)
lastadded = 7
for off in offsets:
tmp = d[(7, off)]
start = (pos + prime) if off == 7 else (prime * (const * (pos + 1 if tmp < 7 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 11
if tk11[pos]:
prime = cpos + 11
p.append(prime)
lastadded = 11
for off in offsets:
tmp = d[(11, off)]
start = (pos + prime) if off == 11 else (prime * (const * (pos + 1 if tmp < 11 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 13
if tk13[pos]:
prime = cpos + 13
p.append(prime)
lastadded = 13
for off in offsets:
tmp = d[(13, off)]
start = (pos + prime) if off == 13 else (prime * (const * (pos + 1 if tmp < 13 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 17
if tk17[pos]:
prime = cpos + 17
p.append(prime)
lastadded = 17
for off in offsets:
tmp = d[(17, off)]
start = (pos + prime) if off == 17 else (prime * (const * (pos + 1 if tmp < 17 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 19
if tk19[pos]:
prime = cpos + 19
p.append(prime)
lastadded = 19
for off in offsets:
tmp = d[(19, off)]
start = (pos + prime) if off == 19 else (prime * (const * (pos + 1 if tmp < 19 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 23
if tk23[pos]:
prime = cpos + 23
p.append(prime)
lastadded = 23
for off in offsets:
tmp = d[(23, off)]
start = (pos + prime) if off == 23 else (prime * (const * (pos + 1 if tmp < 23 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 29
if tk29[pos]:
prime = cpos + 29
p.append(prime)
lastadded = 29
for off in offsets:
tmp = d[(29, off)]
start = (pos + prime) if off == 29 else (prime * (const * (pos + 1 if tmp < 29 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# now we go back to top tk1, so we need to increase pos by 1
pos += 1
cpos = const * pos
# 30k + 1
if tk1[pos]:
prime = cpos + 1
p.append(prime)
lastadded = 1
for off in offsets:
tmp = d[(1, off)]
start = (pos + prime) if off == 1 else (prime * (const * pos + tmp) )//const
del_mult(tmptk[off], start, prime)
# time to add remaining primes
# if lastadded == 1, remove last element and start adding them from tk1
# this way we don't need an "if" within the last while
if lastadded == 1:
p.pop()
# now complete for every other possible prime
while pos < len(tk1):
cpos = const * pos
if tk1[pos]: p.append(cpos + 1)
if tk7[pos]: p.append(cpos + 7)
if tk11[pos]: p.append(cpos + 11)
if tk13[pos]: p.append(cpos + 13)
if tk17[pos]: p.append(cpos + 17)
if tk19[pos]: p.append(cpos + 19)
if tk23[pos]: p.append(cpos + 23)
if tk29[pos]: p.append(cpos + 29)
pos += 1
# remove exceeding if present
pos = len(p) - 1
while p[pos] > N:
pos -= 1
if pos < len(p) - 1:
del p[pos+1:]
# return p list
return p
def sieveOfEratosthenes(n):
"""sieveOfEratosthenes(n): return the list of the primes < n."""
# Code from: <dickinsm@gmail.com>, Nov 30 2006
# http://groups.google.com/group/comp.lang.python/msg/f1f10ced88c68c2d
if n <= 2:
return []
sieve = range(3, n, 2)
top = len(sieve)
for si in sieve:
if si:
bottom = (si*si - 3) // 2
if bottom >= top:
break
sieve[bottom::si] = [0] * -((bottom - top) // si)
return [2] + [el for el in sieve if el]
def sieveOfAtkin(end):
"""sieveOfAtkin(end): return a list of all the prime numbers <end
using the Sieve of Atkin."""
# Code by Steve Krenzel, <Sgk284@gmail.com>, improved
# Code: https://web.archive.org/web/20080324064651/http://krenzel.info/?p=83
# Info: http://en.wikipedia.org/wiki/Sieve_of_Atkin
assert end > 0
lng = ((end-1) // 2)
sieve = [False] * (lng + 1)
x_max, x2, xd = int(sqrt((end-1)/4.0)), 0, 4
for xd in xrange(4, 8*x_max + 2, 8):
x2 += xd
y_max = int(sqrt(end-x2))
n, n_diff = x2 + y_max*y_max, (y_max << 1) - 1
if not (n & 1):
n -= n_diff
n_diff -= 2
for d in xrange((n_diff - 1) << 1, -1, -8):
m = n % 12
if m == 1 or m == 5:
m = n >> 1
sieve[m] = not sieve[m]
n -= d
x_max, x2, xd = int(sqrt((end-1) / 3.0)), 0, 3
for xd in xrange(3, 6 * x_max + 2, 6):
x2 += xd
y_max = int(sqrt(end-x2))
n, n_diff = x2 + y_max*y_max, (y_max << 1) - 1
if not(n & 1):
n -= n_diff
n_diff -= 2
for d in xrange((n_diff - 1) << 1, -1, -8):
if n % 12 == 7:
m = n >> 1
sieve[m] = not sieve[m]
n -= d
x_max, y_min, x2, xd = int((2 + sqrt(4-8*(1-end)))/4), -1, 0, 3
for x in xrange(1, x_max + 1):
x2 += xd
xd += 6
if x2 >= end: y_min = (((int(ceil(sqrt(x2 - end))) - 1) << 1) - 2) << 1
n, n_diff = ((x*x + x) << 1) - 1, (((x-1) << 1) - 2) << 1
for d in xrange(n_diff, y_min, -8):
if n % 12 == 11:
m = n >> 1
sieve[m] = not sieve[m]
n += d
primes = [2, 3]
if end <= 3:
return primes[:max(0,end-2)]
for n in xrange(5 >> 1, (int(sqrt(end))+1) >> 1):
if sieve[n]:
primes.append((n << 1) + 1)
aux = (n << 1) + 1
aux *= aux
for k in xrange(aux, end, 2 * aux):
sieve[k >> 1] = False
s = int(sqrt(end)) + 1
if s % 2 == 0:
s += 1
primes.extend([i for i in xrange(s, end, 2) if sieve[i >> 1]])
return primes
def ambi_sieve_plain(n):
s = range(3, n, 2)
for m in xrange(3, int(n**0.5)+1, 2):
if s[(m-3)/2]:
for t in xrange((m*m-3)/2,(n>>1)-1,m):
s[t]=0
return [2]+[t for t in s if t>0]
def sundaram3(max_n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/2073279#2073279
numbers = range(3, max_n+1, 2)
half = (max_n)//2
initial = 4
for step in xrange(3, max_n+1, 2):
for i in xrange(initial, half, step):
numbers[i-1] = 0
initial += 2*(step+1)
if initial > half:
return [2] + filter(None, numbers)
################################################################################
# Using Numpy:
def ambi_sieve(n):
# http://tommih.blogspot.com/2009/04/fast-prime-number-generator.html
s = np.arange(3, n, 2)
for m in xrange(3, int(n ** 0.5)+1, 2):
if s[(m-3)/2]:
s[(m*m-3)/2::m]=0
return np.r_[2, s[s>0]]
def primesfrom3to(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Returns a array of primes, p < n """
assert n>=2
sieve = np.ones(n/2, dtype=np.bool)
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i/2]:
sieve[i*i/2::i] = False
return np.r_[2, 2*np.nonzero(sieve)[0][1::]+1]
def primesfrom2to(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Input n>=6, Returns a array of primes, 2 <= p < n """
sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
sieve[0] = False
for i in xrange(int(n**0.5)/3+1):
if sieve[i]:
k=3*i+1|1
sieve[ ((k*k)/3) ::2*k] = False
sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]
if __name__=='__main__':
import itertools
import sys
def test(f1,f2,num):
print('Testing {f1} and {f2} return same results'.format(
f1=f1.func_name,
f2=f2.func_name))
if not all([a==b for a,b in itertools.izip_longest(f1(num),f2(num))]):
sys.exit("Error: %s(%s) != %s(%s)"%(f1.func_name,num,f2.func_name,num))
n=1000000
test(sieveOfAtkin,sieveOfEratosthenes,n)
test(sieveOfAtkin,ambi_sieve,n)
test(sieveOfAtkin,ambi_sieve_plain,n)
test(sieveOfAtkin,sundaram3,n)
test(sieveOfAtkin,sieve_wheel_30,n)
test(sieveOfAtkin,primesfrom3to,n)
test(sieveOfAtkin,primesfrom2to,n)
test(sieveOfAtkin,rwh_primes,n)
test(sieveOfAtkin,rwh_primes1,n)
test(sieveOfAtkin,rwh_primes2,n)
Running the script tests that all implementations give the same result.
Print series of prime numbers in python
You need to check all numbers from 2 to n-1 (to sqrt(n) actually, but ok, let it be n).
If n is divisible by any of the numbers, it is not prime. If a number is prime, print it.
for num in range(2,101):
prime = True
for i in range(2,num):
if (num%i==0):
prime = False
if prime:
print (num)
You can write the same much shorter and more pythonic:
for num in range(2,101):
if all(num%i!=0 for i in range(2,num)):
print (num)
As I've said already, it would be better to check divisors not from 2 to n-1, but from 2 to sqrt(n):
import math
for num in range(2,101):
if all(num%i!=0 for i in range(2,int(math.sqrt(num))+1)):
print (num)
For small numbers like 101 it doesn't matter, but for 10**8 the difference will be really big.
You can improve it a little more by incrementing the range you check by 2, and thereby only checking odd numbers. Like so:
import math
print 2
for num in range(3,101,2):
if all(num%i!=0 for i in range(2,int(math.sqrt(num))+1)):
print (num)
Edited:
As in the first loop odd numbers are selected, in the second loop no
need to check with even numbers, so 'i' value can be start with 3 and
skipped by 2.
import math
print 2
for num in range(3,101,2):
if all(num%i!=0 for i in range(3,int(math.sqrt(num))+1, 2)):
print (num)
Fastest way to find prime in intervals
I wrote an equation based prime finder using MillerRabin, it's not as fast as a next_prime finder that sieves, but it can create large primes and you get the equation it used to do it. Here is an example. Following that is the code:
In [5]: random_powers_of_2_prime_finder(1700)
Out[5]: 'pow_mod_p2(27667926810353357467837030512965232809390030031226210665153053230366733641224969190749433786036367429621811172950201894317760707656743515868441833458231399831181835090133016121983538940390210139495308488162621251038899539040754356082290519897317296011451440743372490592978807226034368488897495284627000283052473128881567140583900869955672587100845212926471955871127908735971483320243645947895142869961737653915035227117609878654364103786076604155505752302208115738401922695154233285466309546195881192879100630465, 2**1700-1, 2**1700) = 39813813626508820802866840332930483915032503127272745949035409006826553224524022617054655998698265075307606470641844262425466307284799062400415723706121978318083341480570093257346685775812379517688088750320304825524129104843315625728552273405257012724890746036676690819264523213918417013254746343166475026521678315406258681897811019418959153156539529686266438553210337341886173951710073382062000738529177807356144889399957163774682298839265163964939160419147731528735814055956971057054406988006642001090729179713'
or use to get the answer directly:
In [6]: random_powers_of_2_prime_finder(1700, withstats=False)
Out[6]: 4294700745548823167814331026002277003506280507463037204789057278997393231742311262730598677178338843033513290622514923311878829768955491790776416394211091580729947152858233850115018443160652214481910152534141980349815095067950295723412327595876094583434338271661005996561619688026571936782640346943257209115949079332605276629723961466102207851395372367417030036395877110498443231648303290010952093560918409759519145163112934517372716658602133001390012193450373443470282242835941058763834226551786349290424923951
The code:
import random
import math
def primes_sieve2(limit):
a = [True] * limit
a[0] = a[1] = False
for (i, isprime) in enumerate(a):
if isprime:
yield i
for n in range(i*i, limit, i):
a[n] = False
def ltrailing(N):
return len(str(bin(N))) - len(str(bin(N)).rstrip('0'))
def pow_mod_p2(x, y, z):
"4-5 times faster than pow for powers of 2"
number = 1
while y:
if y & 1:
number = modular_powerxz(number * x, z)
y >>= 1
x = modular_powerxz(x * x, z)
return number
def modular_powerxz(num, z, bitlength=1, offset=0):
xpowers = 1<<(z.bit_length()-bitlength)
if ((num+1) & (xpowers-1)) == 0:
return ( num & ( xpowers -bitlength)) + 2
elif offset == -1:
return ( num & ( xpowers -bitlength)) + 1
elif offset == 0:
return ( num & ( xpowers -bitlength))
elif offset == 1:
return ( num & ( xpowers -bitlength)) - 1
elif offset == 2:
return ( num & ( xpowers -bitlength)) - 2
def MillerRabin(N, primetest, iterx, powx, withstats=False):
primetest = pow(primetest, powx, N)
if withstats == True:
print("first: ",primetest)
if primetest == 1 or primetest == N - 1:
return True
else:
for x in range(0, iterx-1):
primetest = pow(primetest, 2, N)
if withstats == True:
print("else: ", primetest)
if primetest == N - 1: return True
if primetest == 1: return False
return False
PRIMES=list(primes_sieve2(1000000))
def mr_isprime(N, withstats=False):
if N == 2:
return True
if N % 2 == 0:
return False
if N < 2:
return False
if N in PRIMES:
return True
for xx in PRIMES:
if N % xx == 0:
return False
iterx = ltrailing(N - 1)
k = pow_mod_p2(N, (1<<N.bit_length())-1, 1<<N.bit_length()) - 1
t = N >> iterx
tests = [k+1, k+2, k, k-2, k-1]
for primetest in tests:
if primetest >= N:
primetest %= N
if primetest >= 2:
if MillerRabin(N, primetest, iterx, t, withstats) == False:
return False
return True
def lars_last_modulus_powers_of_two(hm):
return math.gcd(hm, 1<<hm.bit_length())
def random_powers_of_2_prime_finder(powersnumber, primeanswer=False, withstats=True):
while True:
randnum = random.randrange((1<<(powersnumber-1))-1, (1<<powersnumber)-1,2)
while lars_last_modulus_powers_of_two(randnum) == 2 and mr_isprime(randnum//2) == False:
randnum = random.randrange((1<<(powersnumber-1))-1, (1<<powersnumber)-1,2)
answer = randnum//2
# This option makes the finding of a prime much longer, i would suggest not using it as
# the whole point is a prime answer.
if primeanswer == True:
if mr_isprime(answer) == False:
continue
powers2find = pow_mod_p2(answer, (1<<powersnumber)-1, 1<<powersnumber)
if mr_isprime(powers2find) == True:
break
else:
continue
if withstats == False:
return powers2find
elif withstats == True:
return f"pow_mod_p2({answer}, 2**{powersnumber}-1, 2**{powersnumber}) = {powers2find}"
return powers2find
def nextprime(N):
N+=2
while not mr_isprime(N):
N+=2
return N
def get_primes(lower, upper):
lower = lower|1
upper = upper|1
vv = []
if mr_isprime(lower):
vv.append(lower)
else:
vv=[nextprime(lower)]
while vv[-1] < upper:
vv.append(nextprime(vv[-1]))
return vv
Here is an example like yours:
In [1538]: cc = get_primes(1009732533765201, 1009732533767201)
In [1539]: print(cc)
[1009732533765251, 1009732533765281, 1009732533765289, 1009732533765301, 1009732533765341, 1009732533765379, 1009732533765481, 1009732533765493, 1009732533765509, 1009732533765521, 1009732533765539, 1009732533765547, 1009732533765559, 1009732533765589, 1009732533765623, 1009732533765749, 1009732533765751, 1009732533765757, 1009732533765773, 1009732533765821, 1009732533765859, 1009732533765889, 1009732533765899, 1009732533765929, 1009732533765947, 1009732533766063, 1009732533766069, 1009732533766079, 1009732533766093, 1009732533766109, 1009732533766189, 1009732533766211, 1009732533766249, 1009732533766283, 1009732533766337, 1009732533766343, 1009732533766421, 1009732533766427, 1009732533766457, 1009732533766531, 1009732533766631, 1009732533766643, 1009732533766667, 1009732533766703, 1009732533766727, 1009732533766751, 1009732533766763, 1009732533766807, 1009732533766811, 1009732533766829, 1009732533766843, 1009732533766877, 1009732533766909, 1009732533766933, 1009732533766937, 1009732533766973, 1009732533767029, 1009732533767039, 1009732533767093, 1009732533767101, 1009732533767147, 1009732533767159, 1009732533767161, 1009732533767183, 1009732533767197, 1009732533767233]
Fastest way to get the first prime numbers up to a limit
It takes 4.8 s on my machine to calculate the first 1,000,000,000 numbers. It's even faster than reading it from a file!
public static void main(String[] args) {
long startingTime = System.nanoTime();
int limit = 1_000_000_000;
BitSet bitSet = findPrimes(limit);
System.out.println(2);
System.out.println(3);
System.out.println(5);
System.out.println(7);
limit = limit / 2 + 1;
for (int i = 6; i < limit; i++) {
if (!bitSet.get(i)) {
int p = i * 2 - 1;
//System.out.println(p);
}
}
System.out.println("Done in " + (System.nanoTime() - startingTime) / 1_000_000 + " milliseconds");
}
public static BitSet findPrimes(int limit) {
BitSet bitSet = new BitSet(limit / 2 + 1);
int size = (int) Math.sqrt(limit);
size += size % 2;
for (int prime = 3; prime < size; prime += 2) {
if (bitSet.get(prime / 2 + 1)) {
continue;
}
for (int i = prime / 2 + 1; i < bitSet.size(); i += prime) {
bitSet.set(i);
}
}
return bitSet;
}
The most efficient way of computing prime values is to use a Sieve of Eratosthenes, first discovered in the ancient greek millenniums ago.
The problem, however, is that it requires a lot of memory. My Java app is simply crashing when I just declare a boolean array of 1,000,000,000.
So I made two memory optimizations to make it work
- Since all the primes after 2 are odd numbers, I half the array size by mapping the indexes to
int p = i * 2 - 1; - Then instead of using an array of booleans, I used a
BitSetthat works with bit operation on an array of longs.
With those two optimizations using the Sieve of Eratosthenes allows you to compute the btSet in 4.8 seconds. As for printing it out, it's another story ;)
Program to find all the prime numbers up to 1 million in 1 sec or as close to it as possible?
You can use the Sieve of Eratosthenes, implemented below:
EDIT: Added modifications to speed it up a bit, thanks to @rossum.
n = 1000000
is_prime = [False, False] + [True] * (n - 1)
primes = [2]
for j in range(4, n + 1, 2):
is_prime[j] = False
for i in range(3, n + 1, 2):
if is_prime[i]:
primes.append(i)
for j in range(i * i, n + 1, i):
is_prime[j] = False
print len(primes) # ==> 78498, which is really the number of primes below 1 million
print primes[:10] # ==> [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
It takes about 1 to 2 seconds on my computer.
Related Topics
In Python Script, How to Set Pythonpath
Python Multiprocessing: Permission Denied
Django Server Killed Frequently
How to Install Pyodbc on Linux
Simulating Key Press Event Using Python For Linux
Show Matplotlib Plots (And Other Gui) in Ubuntu (Wsl1 & Wsl2)
Run Interactive Bash With Popen and a Dedicated Tty Python
Calling a Python Script from Command Line Without Typing "Python" First
How to Do Sed Like Text Replace With Python
Trouble Installing Scipy in Virtualenv on a Amazon Ec2 Linux Micro Instance
Get the Olson Tz Name For the Local Timezone
How to Open a File With the Standard Application
How to Find and Install the Dependencies For Pygame
How to Install R Packages That Are Not Available in "R-Essentials"
Is Everything an Object in Python Like Ruby
What Does Ruby Have That Python Doesn'T, and Vice Versa