Quick search for modest numbers of prime numbers

[GrandFather]

Having asked if there were any quick ways of finding the first n prime numbers and receiving the useful, but unsatisfying, answer - here's a place to download them, I was stimulated to action.

The code below keeps a list of all the primes it has found so far and checks new numbers for primality against that list.

This code generates the first 100,000 primes in about 10 seconds and the first 1,000,000 in about 262 seconds.

use strict;
use warnings;
use Time::HiRes;

my $startTime = [Time::HiRes::gettimeofday ()];
my $toFind = 201;
my @primes;
my $scan = 3;

outer: while (@primes < $toFind - 1) {
    my $root = sqrt $scan;
    for (@primes) {
      next outer if $scan % $_ == 0;
      last if $_ >= $root;
    }
    push @primes, $scan;
  } continue {
  $scan += 2; # Don't need to check even number
  }

print "Found $toFind primes in " . Time::HiRes::tv_interval ($startTim
+e) . " seconds\n";
print "2 ";
print "@primes[$_*10..$_*10+9]\n" for 0..($toFind / 10)-1;

Found 201 primes in 0.002236 seconds
2 3 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 401 409 419
421 431 433 439 443 449 457 461 463 467
479 487 491 499 503 509 521 523 541 547
557 563 569 571 577 587 593 599 601 607
613 617 619 631 641 643 647 653 659 661
673 677 683 691 701 709 719 727 733 739
743 751 757 761 769 773 787 797 809 811
821 823 827 829 839 853 857 859 863 877
881 883 887 907 911 919 929 937 941 947
953 967 971 977 983 991 997 1009 1013 1019
1021 1031 1033 1039 1049 1051 1061 1063 1069 1087
1091 1093 1097 1103 1109 1117 1123 1129 1151 1153
1163 1171 1181 1187 1193 1201 1213 1217 1223 1229

[syphilis]

On my 1GHz machine that code takes 22 seconds to find the first 100,000 primes. Even if you're looking for primes in terms of the "the first x primes", as opposed to terms of "all primes less than y", you can still use a sieve of Eratosthenes - which is siginificantly faster - though you probably need to settle for a little wastefulness. The following code, which is wasteful to the extent that it actually finds the first 105807 primes, does so in less than 2 seconds on the same 1GHz machine. (The wastefulness could be reduced by perhaps lowering $allowance a little closer to 0.5.) If $to_find is set to 1,000,000 it takes about 28 seconds to do the calculation ... and a helluvalot longer to print them out :-)

use warnings;
use strict;

my $t0 = time;

my $to_find = 100000;
my $allowance = 0.6;
my $max = int($to_find * log($to_find) * $allowance);
my $count = 1;

my $imax = sqrt($max * 2 + 1) + 1;
$imax -= 1;
$imax /= 2;

my @bits = (1) x $max;

$bits[0] = 0;

for(my $i = 0; $i < $imax; $i++) {
   if($bits[$i]) {
     my $k = $i * 2 + 1;
     for(my $j = $k + $i; $j < $max; $j += $k) {
        $bits[$j] = 0;
        }
     }
   }

my $t1 = time;

for(my $i = 0; $i < $max; $i++) {
   if($bits[$i]) {
     $count++;
     print $i * 2 + 1, " ";
     if($count == $to_find){last}
     }
   }

print "\nFound at least $count primes in ", $t1 - $t0, " seconds\n";

I haven't pushed the primes into their own array. In the end @bits encodes the primes, and the creation of a separate array of primes amounts to needless replication. (Of course you can still create the array of primes on the fly if you so desire - it's not massively expensive in terms of time.)

You could also use Bit::Vector's implementation (in C) of the Eratosthenes sieve in much the same way. It's faster and consumes less memory despite the fact that it needlessly encodes even numbers. (@bits in the above code encodes only only odd numbers.)

I don't claim that the above script is thoroughly debugged - nor do I claim that it can't be improved upon .... nor do I claim that it is simpler than the script posted by Grandfather. It's just a lot faster ... and that's about all it has going for it :-)

Cheers,
Rob

[johngg]

I had a play years ago with primes and thought I had a pretty good algorithm until someone told me about Eratosthenes. Using his method was much faster. I have tried one variation since which was to use a vector to hold the primes rather than an array/list. The script will print primes as it goes until it reaches the threshhold of sqrt(limit) at which point it knows it has filled the vector up to the limit we have given. After that it prints the rest as fast as the frame buffer can handle.

#!/usr/local/bin/perl -w
#
$limit = shift or die "No limit supplied\n";
die "Limit is not integer\n" unless $limit =~ /^\d+$/;
$sqrtLimit = sqrt $limit;
$sieve = "";
vec($sieve, 0, 1) = 1;
vec($sieve, 1, 1) = 1;
vec($sieve, $limit, 1) = 0;
$reached = 1;
while($reached < $sqrtLimit)
{
    $prime = $reached + 1;
    ++ $prime while vec($sieve, $prime, 1);
    print "$prime is a prime\n";
    $fill = 2 * $prime;
    while($fill <= $limit)
    {
        vec($sieve, $fill, 1) = 1;
    $fill += $prime;
    }
    $reached = $prime;
}
foreach $value ($reached + 1 .. $limit)
{
    print "$value is a prime\n" unless vec($sieve, $value, 1);
}

On an ancient SPARCstation 20 with 150MHz processor it will calculate and print primes up to 100,000 in under 10 seconds.

Cheers

JohnGG

[szbalint]

Oh this one brings back memories.

Around the age of ten/twelve when I first started to take my baby steps with programming languages I was fascinated by primes.

After a lot of trial and error and basically no help from any mathematical textbooks or anything I've implemented more or less the same way of searching for primes back then - in Turbo Pascal. :)

That code was twenty times longer than this, but oh boy I was proud.

I've been in love with programming and computers ever since, but you know how it is, you never forget the first...I thought I'd share.

[ambrus]

See Re: Simple primality testing and thread.