use strict;
$|++;

my $generator = do {
  my $num = 1;
  sub { ++$num }
};

while (1) {
  my $prime = $generator->();
  print "$prime\n";
  ## don't let any multiple of these escape:
  my $oldgenerator = $generator;
  $generator = sub {
    {
      my $next_prime = $oldgenerator->();
      redo if $next_prime % $prime == 0;
      return $next_prime;
    }
  };
}
[download]

-- Randal L. Schwartz, Perl hacker
Be sure to read my standard disclaimer if this is a reply.

The key insight that makes the Sieve of Eratosthenes so nice is that for each prime you only have to do an operation on multiples of that prime. This causes the running time to be just slightly supralinear (IIRC it is O(n*log(log(n))) to get all the primes up to n) rather than being close to O(n**1.5) (yours should be O(n**1.5/log(n)).

Here is an actual Sieve of Eratosthenes implemented with closures. It isn't that fast, but it scales right. Memory usage should scale O(sqrt(n)).

#! /usr/bin/perl -w
use strict;

sub build_sieve {
  my $n = 0;
  my @upcoming_factors = ();
  my ($next_small_p, $sub_iter);
  my $gives_next_square = 5;
  
  return sub {
    LOOP: {
      if (not $n++) {
        return 2; # Special case
      }
      if (not defined $upcoming_factors[0]) {
        if ($n == $gives_next_square) {
          if (not defined ($sub_iter)) {
            # Be lazy to avoid an infinite loop...
            $sub_iter = build_sieve();
            $sub_iter->(); # Throw away 2
            $next_small_p = $sub_iter->();
          }
  
          push @{$upcoming_factors[$next_small_p]}, $next_small_p;
          $next_small_p = $sub_iter->();
          my $next_p2 = $next_small_p * $next_small_p;
          $gives_next_square = ($next_p2 + 1)/2;
          shift @upcoming_factors;
          redo LOOP;
        }
        else {
          shift @upcoming_factors;
          return 2*$n-1;
        }
      }
      else {
        foreach my $i (@{$upcoming_factors[0]}) {
          push @{$upcoming_factors[$i]}, $i;
        }
        shift @upcoming_factors;
        redo LOOP;
      }
    }
  }
}

# Produce as many primes as are asked for, or 100.
my $sieve = build_sieve();
print $sieve->(), "\n" for 1..(shift || 100);
[download]

Wow, this looks like some complex stuff. How exactly would this be benneficial in a program? Or is it basically just concept and application?

Closures allow you to organize code differently. The resulting style is as distinct from OO and procedural as they are from each other.

Right now Dominus is writing a book on the subject. You can find some information in this article on iterators, a number of my nodes talk about closures a bit, see Why I like functional programming for one of the better ones. Or you can search for more information.

Of course in this case the technique is just for fun. A somewhat shorter implementation of a sieve of Eratosthenes is:

sub sieve {
grep{@_[map$a*$_,$_..@_/($a=$_)]=0if$_[$_]>1}@_=0..pop
}

print join " ", sieve(100);
[download]

(Implementation due to tye at (tye)Re3: (Golf): Sieve of Eratosthenes.)

tilly's is a very nice closure implementation. It isn't the fastest pure Perl SoE though. For example, for the first 5761455 primes (primes up to 10^8) tilly's code takes 35.7s on my Macbook. The vector sieve on RosettaCode here takes 37.9s, and the string sieve here takes only 10.2 seconds, albeit with more memory and not an iterator. Using the pure-Perl segmented sieve in Math::Prime::Util::PP takes 7.1 seconds. All of these are terrifically slow compared to C.

If you're looking for even faster SoEs, I recommend Prime Sieve Benchmarks for a list of benchmarks. The C code from AnonymousMonk back in 2003 works fine, but is about 2x slower than the "Trivial 4-line SoE" shown on that page, because it doesn't get the inner loop start correct. Of course they're all in C, though one is built into a Perl module. In general primesieve is the fastest and easiest solution for most projects.

To compare using XS modules, printing primes to 10^8. At this point a majority of our time is actually spent printing, which is why the specialized "print_primes()" routine stomps on everything else, since it has optimized C code (even C's printf is a major bottleneck once the sieving is fast enough).

• 3.4s Math::Prime::XS
• 2.8s Math::Prime::FastSieve
• 2.0s Math::Prime::Util "say for @{primes(1e8)}"
• 1.9s Math::Prime::Util "forprimes { say } 1e8"
• 0.15s Math::Prime::Util "print_primes(1e8);"

Slightly better, here we're naively counting, except for the last one which uses a fast exact prime count method.

• 1.7s Math::Prime::XS
• 0.89s Math::Prime::FastSieve
• 0.39s Math::Prime::Util "$n++ for @{primes(1e8)}"
• 0.18s Math::Prime::Util "forprimes { $n++ } 1e8;"
• 0.00026s Math::Prime::Util "say prime_count(1e8)"

Just for fun (not particularly useful because it's slow): here is a SoE implementation in SQL (postgres). (I didn't compare with the above; it's undoubtedly (very) much slower)

This is a great idea. I always wanted to implement a sieve that can just go on finding new primes, i.e. without specifying a maximum up to which the sieving will be conducted. This would require not only to store the primes up to a point but also kind of the largest multiple of each prime up to that point. Also, I think using the modulus operator % violates the sieving philosophy. I never was able to make it work as I got confused keeping the primes and the multiples.

Using closures is ideal for that. Every time one finds a new prime, one generates a new sifter that is responsible for that prime and to store the progress to date. Here is my formulation:

use strict;
use warnings;

sub sifter {
    my $p = shift;
    my $c = $p;
    return sub { $c += $p while $c < $_[0]; return $c-$_[0]; }
}

my @sieves;
my $n = 1;

loop: while( $n++ ){
    $_->( $n ) or next loop for @sieves;
    push @sieves, sifter( $n );
    print "$n\n";
}
[download]

Clearly this is not a fast algorithm to find primes but I liked the idea...

Of course if you want fast, say 1 million primes in a second you would never use Perl. Perl is great for high level stuff. For raw algorithms.....

#include <stdio.h>
#include <malloc.h>
#include <time.h>

#define TEST(f,x)    (*(f+(x)/16)&(1<<(((x)%16L)/2)))
#define SET(f,x)    *(f+(x)/16)|=1<<(((x)%16L)/2)

int
main(int argc, char *argv[])
{
    unsigned char *ptr=NULL;
    unsigned long max, mom, s=0, e=1;
    register unsigned long test=1;
    register unsigned long count=1;
    time_t begin;
    
    max = ( argc>1 ) ? atol(argv[1]) : 0xFFFFFFL;
    ptr = (unsigned char *) calloc( ((max>>4)+1L), sizeof(char) );
    if ( ! ptr ) {
        printf( "Can't allocate enough memory!\n" );
        system("pause");
      return 1;
    }

    printf( "Searching prime numbers to : %ld\n", max );
    
    begin = time (NULL);
    while ( (test+=2) < max) if (!TEST(ptr, test) ) {
        if  ( ++count%0xFFL==0 ) {
            printf ( "%ld prime number\x0d", count );
            fflush(stdout);
        }
        for ( mom=3L*test; mom<max; mom+=test<<1 ) SET (ptr, mom);
    }

    printf( " %ld prime numbers found in %ld secs.\n\n", count, time(N
+ULL)-begin);

    printf( "Show prime numbers ( Enter end < start to exit )" );
    while ( s<e ) {
        printf ("\n\nStart of Area : "); 
        fflush (stdout); scanf ("%ld", &s);
        printf ("End of Area   : ");     
        fflush (stdout); scanf ("%ld", &e);
        count = s-2;
        if ( s%2==0 ) count++;
        while ((count+=2)<e) if (!TEST(ptr,count)) printf ("%ld\t", co
+unt);
    }

    free(ptr);
  return 0;
}

/*

Searching prime numbers to : 16777215
 1077871 prime numbers found in 1 secs.

Show prime numbers ( Enter end < start to exit )

Start of Area : 1
End of Area   : 1000
1       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

Start of Area :

*/
[download]