From here I found a fast method to make prime list, the python version works well, but my Mathematica version does not. Obviously ,169
is not a prime number. What's wrong with my program?
'''Python code'''
n = 200
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
print ((3 * np.nonzero(sieve)[0] + 1) | 1)
(*Mathematica code*)
Clear["`*"];
n = 200;
p = ConstantArray[1, Quotient[n, 3] + Boole[Mod[n, 6] == 2]];
p[[1]] = 0;
Do[
If[p[[i]] != 0,
k = BitOr[3 (i - 1) + 1, 1];
p[[Quotient[k^2, 3] + 1 ;; -1 ;; 2 k]] = 0;
p[[Quotient[(k^2 + 4 k - 2 k BitAnd[i - 1, 1]), 3] + 1 ;; -1 ;; 2 k]] = 0;],
{i, 1, Floor[n^0.5]/3}];
res = BitOr[3 (Flatten@SparseArray[p]["NonzeroPositions"] - 1) + 1, 1];
Pick[res, PrimeQ @ res, False]
Prime[Range[50]]
gives the first 50 primes. $\endgroup$Clear["`*"];
. $\endgroup$Prime@Range@PrimePi[10^7] // Length
on my laptop takes about3s
, that method only need0.4s
. $\endgroup${i, 1, Floor[n^0.5]/3+2}
$\endgroup$