# What's wrong with my code for finding primes?

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]

• The simple way: Prime[Range[50]] gives the first 50 primes. Aug 26, 2013 at 17:32
• Might I suggest you set your "Notebook's Default Context" to "Unique to this Notebook"? It eliminates some of the need to run Clear["*"];. Aug 26, 2013 at 17:32
• @bills,Prime@Range@PrimePi[10^7] // Length on my laptop takes about 3s, that method only need 0.4s. Aug 26, 2013 at 17:41
• {i, 1, Floor[n^0.5]/3+2} Aug 27, 2013 at 2:27
• explorer, I believe "code-review" is best reserved for critique of working code, though it doesn't say that in the tag wiki. I have replaced the tag with "broken-code" which I think is more accurate. Aug 27, 2013 at 6:43

n = 400;
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];
x = Quotient[k^2, 3] + 1;
y = Quotient[(k^2 + 4 k - 2 k BitAnd[i - 1, 1]), 3] + 1;
If[x <= Length@p, p[[x ;; -1 ;; 2 k]] = 0];
If[y <= Length@p, p[[y ;; -1 ;; 2 k]] = 0];],
{i, 1, Floor[Sqrt[n]]/3 + 2}
];
res = BitOr[3 (Flatten@SparseArray[p]["NonzeroPositions"] - 1) + 1, 1]
Pick[res, PrimeQ@res, False]

(* {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} *)
(* {} *)