my \$Limit = 1000000; my \$HighestFactor = sqrt(\$Limit); my \$is_prime=''; # Sieve array. # put in candidate primes: integers which have an odd number of representations by certain quadratic forms. for \$x ( 1..\$HighestFactor){ my \$x2 = \$x*\$x; last if \$x2*2 >= \$Limit; for \$y ( 1..\$HighestFactor ){ my \$y2 = \$y*\$y; next if (\$n = 3*\$x2 - \$y2) > \$Limit; vec(\$is_prime,\$n,1) ^= 1 if \$x > \$y && \$n % 12 == 11; next if ( \$n = 3*\$x2 + \$y2 ) > \$Limit; vec(\$is_prime,\$n,1) ^= 1 if \$n % 12 == 7; next if ( \$n = 4*\$x2 + \$y2 ) > \$Limit; vec(\$is_prime,\$n,1) ^= 1 if \$n % 12 == 1 || \$n % 12 == 5; } } # eliminate composites by sieving # if n is prime, omit all multiples of its square; this is sufficient because # composites which managed to get on the list cannot be square-free for \$n (5..\$HighestFactor ){ next unless vec(\$is_prime,\$n,1); for( \$k=\$n2=\$n*\$n; \$k <= \$Limit; \$k += \$n2 ){ vec(\$is_prime,\$k,1) = 0 }; } # Present the results. \$\="\n"; print for 2,3, grep vec(\$is_prime,\$_,1), 5..\$Limit;