P is for Practical PerlMonks

### Re: Math help: Finding Prime Numbers

by jbert (Priest)
 on Nov 12, 2006 at 23:00 UTC ( #583609=note: print w/replies, xml ) Need Help??

in reply to Math help: Finding Prime Numbers

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;

Replies are listed 'Best First'.
Re^2: Math help: Finding Prime Numbers
by I0 (Priest) on Nov 25, 2006 at 08:25 UTC
```my \$Limit = 1000000;
my \$HighestFactor = sqrt(\$Limit);
my \$is_prime='';                  # Sieve array.
# put in candidate primes: integers which have an odd number of repres
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;

Create A New User
Node Status?
node history
Node Type: note [id://583609]
help
Chatterbox?
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others meditating upon the Monastery: (5)
As of 2018-01-20 19:12 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
How did you see in the new year?

Results (227 votes). Check out past polls.

Notices?