What is the best, most performant algorithm to find all primes up to a given number?
Sieve of Eratosthenes. This algorithm can generate all prime numbers up to
n. Time complexity -O(nlog(n)), memory complexity -O(n)BPSW primality test. This algorithm can check if
nis pseudoprime. It was tested on first 10^15 numbers. Time complexity -O(log(n)).
UPDATE:
I did some research and wrote simple implementation of generating prime numbers in c#. Main idea when we check number N for primality - we just need to check if it divisible by any prime number that less than sqrt(N).
First implementation:
public static List<int> GeneratePrimes(int n)
{
var primes = new List<int>();
for(var i = 2; i <= n; i++)
{
var ok = true;
foreach(var prime in primes)
{
if (prime * prime > i)
break;
if (i % prime == 0)
{
ok = false;
break;
}
}
if(ok)
primes.Add(i);
}
return primes;
}
Test results:
10^6 - 0.297s
10^7 - 6.202s
10^8 - 141.860s
Second implementation using parallel computing:
1. Generate all primes up to sqrt(N)
2. Generate all primes from sqrt(N) + 1 to N using primes up to sqrt(N) using parallel computing.
public static List<int> GeneratePrimesParallel(int n)
{
var sqrt = (int) Math.Sqrt(n);
var lowestPrimes = GeneratePrimes(sqrt);
var highestPrimes = (Enumerable.Range(sqrt + 1, n - sqrt)
.AsParallel()
.Where(i => lowestPrimes.All(prime => i % prime != 0)));
return lowestPrimes.Concat(highestPrimes).ToList();
}
Test results:
10^6 - 0.276s
10^7 - 4.082s
10^8 - 78.624
Finding the list of prime numbers in shortest time
The good thing in your primality test is that you only divide by primes.
private static boolean checkMod( int num) {
for (int i : primes){
if( num % i == 0){
return false;
}
}
return true;
}
The bad thing is that you divide by all primes found so far, that is, all primes smaller than the candidate. That means that for the largest prime below one million, 999983, you divide by 78497 primes to find out that this number is a prime. That's a lot of work. So much, in fact, that the work spent on primes in this algorithm accounts for about 99.9% of all work when going to one million, a larger part for higher limits. And that algorithm is nearly quadratic, to find the primes to n in this way, you need to perform about
n² / (2*(log n)²)
divisions.
A simple improvement is to stop the division earlier. Let n be a composite number (i.e. a number greter than 1 that has divisors other than 1 and n), and let d be a divisor of n.
Now, d being a divisor of n means that n/d is an integer, and also a divisor of n: n/(n/d) = d.
So we can naturally group the divisors of n into pairs, each divisor d gives rise to the pair (d, n/d).
For such a pair, there are two possibilities:
d = n/d, which meansn = d², ord = √n.- The two are different, then one of them is smaller than the other, say
d < n/d. But that immediately translates tod² < nord < √n.
So, either way, each pair of divisors contains (at least) one not exceeding √n, hence, if n is a composite number, its smallest divisor (other than 1) does not exceed √n.
So we can stop the trial division when we've reached √n:
private static boolean checkMod( int num) {
for (int i : primes){
if (i*i > n){
// We have not found a divisor less than √n, so it's a prime
return true;
}
if( num % i == 0){
return false;
}
}
return true;
}
Note: That depends on the list of primes being iterated in ascending order. If that is not guaranteed by the language, you have to use a different method, iterate by index through an ArrayList or something like that.
Stopping the trial division at the square root of the candidate, for the largest prime below one million, 999983, we now only need to divide it by the 168 primes below 1000. That's a lot less work than previously. Stopping the trial division at the square root, and dividing only by primes, is as good as trial division can possibly get and requires about
2*n^1.5 / (3*(log n)²)
divisions, for n = 1000000, that's a factor of about 750, not bad, is it?
But that's still not very efficient, the most efficient methods to find all primes below n are sieves. Simple to implement is the classical Sieve of Eratosthenes. That finds the primes below n in O(n*log log n) operations, with some enhancements (eliminating multiples of several small primes from consideration in advance), its complexity can be reduced to O(n) operations. A relatively new sieve with better asymptotic behaviour is the Sieve of Atkin, which finds the primes to n in O(n) operations, or with the enhancement of eliminating the multiples of some small primes, in O(n/log log n) operations.
The Sieve of Atkin is more complicated to implement, so it's likely that a good implementation of a Sieve of Eratosthenes performs better than a naive implementation of a Sieve of Atkin. For implementations of like optimisation levels, the performance difference is small unless the limit becomes large (larger than 1010; and it's not uncommon that in practice, a Sieve of Eratosthenes scales better than a Sieve of Atkin beyond that, due to better memory access patterns). So I would recommend beginning with a Sieve of Eratosthenes, and only when its performance isn't satisfactory despite honest efforts at optimisation, delve into the Sieve of Atkin. Or, if you don't want to implement it yourself, find a good implementation somebody else has already seriously tuned.
I have gone into a bit more detail in an answer with a slightly different setting, where the problem was finding the n-th prime. Some implementations of more-or-less efficient methods are linked from that answer, in particular one or two usable (though not much optimised) implementations of a Sieve of Eratosthenes.
What is the fastest way of finding prime between n and 2n
The fastest way would probably to pre-compute and store a 1-dimensional array of size 2^32, where the value for index n is the desired prime number between n and 2n. This would be an outrageous use of memory, of course, but it probably is fastest.
A slightly slower way that uses much less memory is to pre-compute and store a list of all "Bertrand primes", where the first element is the first
prime number and each element after the first is the greatest prime number
less than double the previous element. You can use a binary search of that list to quickly find the desired prime number. If you want 1 < n < 2^32 you need the last prime in that list to be above 2^32 to catch all such n. That would need a list of just 34 prime numbers, very doable. By the way, if you want to do this up to 2^64 you need only 66 prime numbers.
Here is Python 3.5 code to implement that algorithm. It uses a binary search function in the standard library. The list of Bertrand primes was found with another simple Python routine, though it is also available in The Online Encyclopedia of Integer Sequences, sequence A006992.
from bisect import bisect_right
_bertrand_primes = [
2, 3, 5, 7, 13, 23,
43, 83, 163, 317, 631, 1259,
2503, 5003, 9973, 19937, 39869, 79699,
159389, 318751, 637499, 1274989, 2549951, 5099893,
10199767, 20399531, 40799041, 81598067, 163196129, 326392249,
652784471, 1305568919, 2611137817, 5222275627]
def prime_between_n_and_2n(n):
"""Find a prime number p such that n < p < 2n. The returned value
will be the first 'Bertrand prime' <https://oeis.org/A006992>
greater than n. n is limited to 1 < n < 2**32 but need not be an
integer. Outside those limits, None is returned.
"""
if 1 < n < 2**32:
return _bertrand_primes[bisect_right(_bertrand_primes, n)]
Fastest way to prime factorise a number up to 10^18
Approach
Lets say you have a number n that goes up to 1018 and you want to prime factorise it. Since this number can be as small as unity and as big as 1018, all along it can be prime as well as composite, this would be my approach -
- Using miller rabin primality testing, make sure that the number is composite.
- Factorise
nusing primes up to 106, which can be calculated using sieve of Eratosthenes.
Now the updated value of n is such that it has prime factors only above 106 and since the value of n can still be as big as 1018, we conclude that the number is either prime or it has exactly two prime factors (not necessarily distinct).
- Run Miller Rabin again to ensure the number isn't prime.
- Use Pollard rho algorithm to get one prime factor.
You have the complete factorisation now.
Lets look at the time-complexity of the above approach:
- Miller Rabin takes
O(log n) - Sieve of Eratosthenes takes
O(n*log n) - The implementation of Pollard rho I shared takes
O(n^0.25)
Time Complexity
Step 2 takes maximum time which is equal to O(10^7), which is in turn the complexity of the above algorithm. This means you can find the factorisation within a second for almost all programming languages.
Space Complexity
Space is used only in the step 2 where sieve is implemented and is equal to O(10^6). Again, very practical for the purpose.
Implementation
Complete Code implemented in C++14. The code has a hidden bug. You can either reveal it in the next section, or skip towards the challenge ;)
Bug in the code
In
line 105, iterate tilli<=np. Otherwise, you may miss the cases whereprime[np]=999983is a prime factor
Challenge
Give me a value of n, if any, where the shared code results in wrong prime factorisation.
Bonus
How many such values of n exist ?
Solution
For such value of n, assertion in
Line 119may fail.
Bonus Solution
Lets call
P=999983. All numbers of the formn = p*q*rwhere p, q, r are primes >=Psuch that at least one of them is equal toPwill result in wrong prime factorisation.
There are exactly four such numbers: {P03, P02P1, P02P2, P0P12}, where P0 = P = 999983, P1 = next_prime(P0) = 1000003, P2 = next_prime(P1) = 1000033.
fastest method for finding number of prime numbers between two large numbers x and y
Use a segmented sieve of Eratosthenes.
That is, use a bit set to store the numbers between x and y, represented by x as an offset and a bit set for [0,y-x). Then sieve (eliminate multiples) for all the primes less or equal to the square root of y. Those numbers that remain in the set are prime.
With y at most 1012 you have to sieve with primes up to at most 106, which will take less than a second in a proper implementation.
Better algorithm to find nth prime number?
For sufficiently large N (e.g. more than a million or so), the best algorithm is to use an approximation (e.g. Logarithmic Integral or Riemann's R function), then use a fast prime count method such as LMO, then sieve the small remainder. This is many orders of magnitude faster than sieving.
See https://math.stackexchange.com/questions/507178/most-efficient-algorithm-for-nth-prime-deterministic-and-probabilistic
There are at least two open source implementations.
- https://metacpan.org/pod/ntheory
- https://github.com/kimwalisch/primecount
The latter has progressed past the first, and is also multithreaded.
Add: Will Ness has also pointed out a nice post from Daniel Fischer that provides a walkthrough of different ways to solve this: Calculating and printing the nth prime number
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.
Related Topics
How to Use Boost in Visual Studio 2010
Unicode Encoding For String Literals in C++11
Does C++ Support 'Finally' Blocks? (And What's This 'Raii' I Keep Hearing About)
C99 Stdint.H Header and Ms Visual Studio
When Do I Use a Dot, Arrow, or Double Colon to Refer to Members of a Class in C++
Are There Benefits of Passing by Pointer Over Passing by Reference in C++
How to Enable C++17 Compiling in Visual Studio
How to Make a Http Request With C++
Adding External Library into Qt Creator Project
Simple Example of Threading in C++
What Is an 'Undeclared Identifier' Error and How to Fix It
How to Print (Using Cout) a Number in Binary Form