#!/usr/bin/perl -w
# primes - generate primes
# Written for the PPT initiative by Jonathan Feinberg.
# The algorithm was substantially modified by Benjamin Tilly.
# See docs for license.
use strict;
#use integer; # faster, but cuts the maxint down
$|++;
my @primes = (2, 3, 5, 7, 11); # None have been tested
my @next_primes = (); # Avoid redoing work
my $VERSION = '1.001';
END {
close STDOUT || die "$0: can't close stdout: $!\n";
$? = 1 if $? == 255; # from die
}
chomp(my $start = @ARGV ? $ARGV[0] : );
my $end = $ARGV[1] || 2**32 - 1;
for ($start, $end) {
s/^\s*\+?(\d{1,10}).*/$1/ || die "$0: $_: illegal numeric format\n";
$_ > 2**32 - 1 && die "$0: $_: Numerical result out of range\n";
}
primes ($start, $end);
exit 0;
sub primes {
my ($start, $end) = @_;
return if $end <= $start;
if ($start <= 2 and 2 < $end) {
print "2\n";
}
$end--; # Reindex
# Initialize the list of primes
&more_primes($primes[-1]+1, int(sqrt($end)));
while (scalar @next_primes) {
push @primes, @next_primes;
# Careful, we need to ensure that
# we get a prime past sqrt($end)...
&more_primes($primes[-1]+1, int(2*sqrt($end)));
}
my $from = $start-1;
my $to = $from;
until ($to == $end) {
$from = $to + 1;
$to = $from + 99999; # By default do 100,000
$to = $end if $end < $to; # Unless I can finish in one pass
&more_primes($from, $to);
print map {"$_\n"} @next_primes; # Print primes
}
}
sub more_primes {
# This adds to the list of primes until it reaches $max
# or the square of the largest current prime (assumed odd)
my $base = shift;
my $max = shift;
my $square = $primes[-1] * $primes[-1];
$max = $square if $square < $max; # Determine what to find primes to
$base++ unless $base % 2; # Make the base odd
$max-- if $max %2; # Make the max odd
$max = ($max - $base)/2; # Make $max into a count of odds
return @next_primes = () if $max < 0; # Sanity check
my @more = map {0} 0..$max; # Initialize array of 0's for the
# odd numbers in our range
shift @primes; # Remove 2
foreach my $p (@primes) {
my $start;
if ($base < $p * $p) {
$start = ($p * $p - $base)/2; # Start at the square
if ($max < $start) { # Rest of primes don't matter!
last;
}
}
else { # Start at first odd it divides
$start = $base % $p; # Find remainder
$start = $p - $start if $start; # Distance to first thing it divides
$start += $p if $start %2; # Distance to first odd it divides
$start = $start/2; # Reindex for counting over odd!
}
for (my $i = $start; $i <= $max; $i += $p) {
$more[$i] = 1;
}
}
unshift @primes, 2; # Replace 2
# Read off list of primes
@next_primes = map {$_ + $_ + $base} grep {$more[$_] == 0} 0..$max;
}
__END__
=head1 NAME
B - generate primes
=head1 SYNOPSIS
B [I [I]]
=head1 DESCRIPTION
The B utility prints primes in ascending order, one per line,
starting at or above I and continuing until, but not including
I. The I value must be at least 0 and not greater than
stop. The I value must not be greater than 4294967295. The
default value of I is 4294967295.
When the primes utility is invoked with no arguments, I is read
from standard input. I is taken to be 4294967295. The I
value may be preceded by a single C<+>. The I value is
terminated by a non-digit character (such as a newline).
=head1 BUGS
B has no known bugs. The algorithm is non-optimal.
=head1 AUTHOR
The Perl implementation of I was originally written by Jonathan
Feinberg, I and modified by Benjamin Tilly,
I.
=head1 COPYRIGHT and LICENSE
This program is copyright (c) Jonathan Feinberg and Benjamin Tilly (1999).
This program is free and open software. You may use, modify, distribute,
and sell this program (and any modified variants) in any way you wish,
provided you do not restrict others from doing the same.