python - Why does my naive implementation of the Sieve of Atkins exclude 5? -

python - Why does my naive implementation of the Sieve of Atkins exclude 5? -

i wrote extremely naive implementation of sieve of atkin, based on wikipedia's inefficient clear pseudocode. wrote algorithm in matlab, , omits 5 prime number. wrote algorithm in python same result.

technically, know why 5 beingness excluded; in step n = 4*x^2 + y^2, n == 5 when x == 1 , y == 1. occurs once, 5 flipped prime nonprime , never flipped back.

why doesn't algorithm match 1 on wikipedia? although have made few superficial adjustments (e.g. calculating x^2 1 time in each iteration, storing value of mod(n, 12) when using in first equation, etc.) shouldn't alter logic of algorithm.

i read on several discussions related sieve of atkin, can't tell differences creating problems in implementations.

python code: def atkin1(limit): res = [0] * (limit + 1) res[2] = 1 res[3] = 1 res[5] = 1 limitsqrt = int(math.sqrt(limit)) x in range(1, limitsqrt+1): y in range(1, limitsqrt+1): x2 = x**2 y2 = y**2 n = 4*x2 + y2 if n == 5: print('debug1') nmod12 = n % 12 if n <= limit , (nmod12 == 1 or nmod12 == 5): res[n] ^= 1 n = 3*x2 + y2 if n == 5: print('debug2') if n <= limit , (n % 12 == 7): res[n] ^= 1 if x > y: n = 3*x2 - y2 if n == 5: print('debug3') if n <= limit , (n % 12 == 11): res[n] ^= 1 ndx = 5 while ndx <= limitsqrt: m = 1 if res[ndx]: ndx2 = ndx**2 ndx2mult =m * ndx2 while ndx2mult < limit: res[ndx2mult] = 0 m += 1 ndx2mult = m * ndx2 ndx += 1 homecoming res matlab code function p = atkin1(limit) % 1. create results list, filled 2, 3, , 5 res = [0, 1, 1, 0, 1]; % 2. create sieve list entry each positive integer; entries of % list should marked nonprime (composite). res = [res, zeros(1, limit-5)]; % 3. each entry number n in sieve list, modulo-sixty remainder r: limitsqrt = floor(sqrt(limit)); x=1:limitsqrt y=1:limitsqrt x2 = x^2; y2 = y^2; % if r 1, 13, 17, 29, 37, 41, 49, or 53, flip entry each % possible solution 4x^2 + y^2 = n. n = 4*x2 + y2; nmod12 = mod(n, 12); if n <= limit && (nmod12 == 1 || nmod12 == 5) res(1, n) = ~res(1, n); end % if r 7, 19, 31, or 43, flip entry each possible solution % 3x^2 + y^2 = n. n = 3*x2 + y2; if n <= limit && mod(n, 12) == 7 res(1, n) = ~res(1, n); end % if r 11, 23, 47, or 59, flip entry each possible solution % 3x^2 - y^2 = n when x > y. if x > y n = 3*x2 - y2; if n <= limit && mod(n, 12) == 11 res(1, n) = ~res(1, n); end end % if r else, ignore completely. end end % 4. start lowest number in sieve list. ndx = 5; while ndx < limitsqrt m = 1; if res(ndx) % 5. take next number in sieve list still marked prime. % 6. include number in results list. % 7. square number , mark multiples of square nonprime. ndx2 = ndx^2; ndx2mult = m * ndx2; while ndx2mult < limit res(ndx2mult) = 0; m = m + 1; ndx2mult = m * ndx2; end end % 8. repeat steps 5 through eight. ndx = ndx + 1; end p = find(res); % find indexes of nonzerogo end wikipedia pseudocode // arbitrary search limit limit ← 1000000 // initialize sieve is_prime(i) ← false, ∀ ∈ [5, limit] // set in candidate primes: // integers have odd number of // representations quadratic forms (x, y) in [1, √limit] × [1, √limit]: n ← 4x²+y² if (n ≤ limit) , (n mod 12 = 1 or n mod 12 = 5): is_prime(n) ← ¬is_prime(n) n ← 3x²+y² if (n ≤ limit) , (n mod 12 = 7): is_prime(n) ← ¬is_prime(n) n ← 3x²-y² if (x > y) , (n ≤ limit) , (n mod 12 = 11): is_prime(n) ← ¬is_prime(n) // eliminate composites sieving n in [5, √limit]: if is_prime(n): // n prime, omit multiples of square; // sufficient because composites managed // on list cannot square-free is_prime(k) ← false, k ∈ {n², 2n², 3n², ..., limit} print 2, 3 n in [5, limit]: if is_prime(n): print n

in wikipedia's text description of algorithm, "results list" , "sieve list" 2 different things. matlab code has 1 vector used both. initial value in sieve list 5 should "non-prime".

python matlab sieve-of-atkin