// arbitrary search limit
limit ← 1000000         

// initialize the sieve
is_prime(i) ← false, i ∈ [5, limit] 

// put in candidate primes: 
// integers which have an odd number of
// representations by certain quadratic forms
for (x, y) in [1, √limit] × [1, √limit]:
    n ← 4x²+y²
    if (n ≤ limit) ∧ (n mod 12 = 1 ∨ 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 by sieving
for n in [5, √limit]:
    if is_prime(n):
        // n is prime, omit multiples of its square; this is
        // sufficient because composites which managed to get
        // on the list cannot be square-free
        is_prime(k) ← false, k ∈ {n², 2n², 3n², ..., limit} 

print 2, 3
for n in [5, limit]:
    if is_prime(n): print n