Prime Sieves in Perl

prime_sieves.pl

#!/usr/bin/env perl
package Acme::Math::ExampleSieves;
use warnings;
use strict;

# A simple example, showing some new prime sieves in Perl (all SoE).
# Copyright 2012 by Dana Jacobsen (dana@acm.org)
# This program is free software using the Perl 5 artistic license v2.

# Note, the string sieve will be 2x faster using Perl 5.16.0 or later.
# On my machine it can generate and print the first 10 million primes
# in under 1 second.
#
# The modules Math::Prime::Util and Math::Prime::FastSieve will generate
# them much faster.

__PACKAGE__->run unless caller;

# Using a vector.  1 bit per odd number.
sub dj_vector {
  my($end) = @_;
  return () if $end < 2;  return (2) if $end < 3;
  $end-- if ($end & 1) == 0;   # Ensure end is odd

  my ($sieve, $n, $limit, $s_end) = ( '', 3, int(sqrt($end)), $end >> 1 );
  while ( $n <= $limit ) {
    for (my $s = ($n*$n) >> 1; $s <= $s_end; $s += $n) {
      vec($sieve, $s, 1) = 1;
    }
    do { $n += 2 } while vec($sieve, $n >> 1, 1) != 0;
  }
  my @primes = (2);
  do { push @primes, 2*$_+1 if !vec($sieve,$_,1) } for (1..int(($end-1)/2));
  @primes;
}

# Using an array.  1 array entry per odd number.
sub dj_array {
  my($end) = @_;
  return () if $end < 2;  return (2) if $end < 3;

  my @sieve;
  my ($n, $limit) = ( 3, int(sqrt($end)) );
  while ( $n <= $limit ) {
    for (my $s = $n*$n; $s <= $end; $s += 2*$n) {
      $sieve[$s>>1] = 1;
    }
    do { $n += 2 } while $sieve[$n>>1];
  }
  my @primes = (2);
  do { push @primes, 2*$_+1 if !$sieve[$_] } for (1 .. int(($end-1)/2));
  @primes;
}

# Using a string.  1 byte per odd number.
# Using Perl 5.16.0 or later, this is twice as fast as the others.
# With earlier versions, it is comparable in speed.
sub dj_string {
  my($end) = @_;
  return () if $end < 2;  return (2) if $end < 3;
  $end-- if ($end & 1) == 0;
  my $s_end = $end >> 1;

  my $whole = int( ($end>>1) / 15);    # prefill with 3 and 5 marked
  my $sieve = '100010010010110' . '011010010010110' x $whole;
  substr($sieve, ($end>>1)+1) = '';
  my ($n, $limit) = ( 7, int(sqrt($end)) );
  while ( $n <= $limit ) {
    for (my $s = ($n*$n) >> 1; $s <= $s_end; $s += $n) {
      substr($sieve, $s, 1) = '1';
    }
    do { $n += 2 } while substr($sieve, $n>>1, 1);
  }
  # If you just want the count, it's very fast:
  #       my $count = 1 + $sieve =~ tr/0//;
  my @primes = (2, 3, 5);
  $n = 7-2;
  foreach my $s (split("0", substr($sieve, 3), -1)) {
    $n += 2 + 2 * length($s);
    push @primes, $n if $n <= $end;
  }
  @primes;
}

sub run {
  my $end = $ARGV[0] || 100;
  print join("\n", dj_string($end)), "\n";;
}