Re: Math help: Finding Prime Numbers

[jbert]

The sieve of Eratosthenes is a good algorithm. It requires N bits of storage for primes to N.

Here is a pure-perl implementation, using 'vec' for the bit array and done as an OO perl module built around a blessed scalar reference (to the bit array). It isn't heavily tested, I'm afraid, but you use it as follows:

use Sieve; # Rename as appropriate...

my $sieve = Sieve->new;
foreach my $number (qw/13 20 35 3/) {
    print "$number is ",
    $sieve->is_prime($number) ? "prime" : "composite", "\n";
}
print join(", ", $sieve->primes_to(300)), "\n";

And here is the code:

package Sieve;
use strict;
use warnings;

my $BITS_PER_BYTE = 8;
my $INITIAL_SIZE = $BITS_PER_BYTE ** 2;

sub new {

    # A sieve is a bit array, where 'true' => composite
    my $sieve = '';
    # 0 isn't prime
    vec($sieve, 0, 1) = 1; 
    # 1 isn't prime
    vec($sieve, 1, 1) = 1; 

    my $s = \$sieve;
    bless $s, __PACKAGE__;
    # Pre-extend the array
    vec($sieve, $INITIAL_SIZE, 1) = 0; 

    # And fill it in
    $s->_run;
    return $s;
}

sub is_prime {
    my $s = shift;
    my $n = shift;

    $s->_extend($n);

    return !vec($$s, $n, 1);
}

sub primes_to {
    my $s = shift;
    my $n = shift;
    $s->_extend($n);
    return grep { $s->is_prime($_) } 1..$n;
}

sub _run {
    my $s = shift;

    my $i;
    my $limit = sqrt ($s->_size);
    for ($i = 2; $i < $limit; ++$i) {
        next unless $s->is_prime($i);
        $s->_mark_multiples($i);
    }
}

sub _extend {
    my $s = shift;
    my $to = shift;

    return 0 if $to <= $s->_size;
    vec($$s, $to, 1) = 0;

    return $s->_run;
}

sub _mark_multiples {
    my $s = shift;
    my $p = shift;

    my $i;
    my $limit = $s->_size;
    for ($i = 2 * $p; $i < $limit; $i += $p) {
        vec($$s, $i, 1) = 1;
    }
}

sub _size {
    my $s = shift;
    return (length $$s) * $BITS_PER_BYTE;
}

1;

[I0]

my $Limit = 1000000;
my $HighestFactor = sqrt($Limit);
my $is_prime='';                  # Sieve array.
# put in candidate primes: integers which have an odd number of repres
+entations 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;