Update June 26, 2015: New results using Math::Prime::Util v0.51 at the end of the post.

Some time ago traveled to an imaginary place filled with prime numbers. In one of the paintings was the following statement. With k initialized to 1, the value alternates from 2 to 1 repeatedly. It was this statement which inspired me to give Algorithm3 a try.

k = 3 - k;
[download]

There are various prime modules on CPAN. Thus, have no plans to create another. All I wanted to do was to compare the Sieve of Eratosthenes with Algorithm3. I created a local copy of Math::Prime::FastSieve by davido and replaced the primes function.

Sieve of Eratosthenes

/* Sieve of Eratosthenes.  Return a reference to an array containing a
+ll
 * prime numbers less than or equal to search_to.  Uses an optimized s
+ieve
 * that requires one bit per odd from 0 .. n.  Evens aren't represente
+d in the
 * sieve.  2 is just handled as a special case.
 */

SV* primes( long search_to )
{
    AV* av = newAV();
    if( search_to < 2 )
        return newRV_noinc( (SV*) av ); // Return an empty list ref.
    av_push( av, newSVuv( 2UL ) );
    // Allocate space for odd numbers (15 bits per 30 values)
    sieve_type primes( search_to/2 + 1, 0 );
    // Sieve over the odd numbers
    for( sieve_size_t i = 3; i * i <= search_to; i+=2 )
        if( ! primes[i/2] )
            for( sieve_size_t k = i*i; k <= search_to; k += 2*i)
                primes[k/2] = 1;
    // Add each prime to the list ref
    for( sieve_size_t i = 3; i <= search_to; i += 2 )
        if( ! primes[i/2] )
          av_push( av, newSVuv( static_cast<unsigned long>( i ) ) );
    return newRV_noinc( (SV*) av );
}
[download]

Sieve of Xuedong Luo (Algorithm3)

/* Sieve of Xuedong Luo (Algorithm3). Return a reference to an array
 * containing all prime numbers less than or equal to search_to.
 *
 *    A practical sieve algorithm for finding prime numbers.
 *    ACM Volume 32 Issue 3, March 1989, Pages 344-346
 *    http://dl.acm.org/citation.cfm?doid=62065.62072
 *
 * Avoid all composites that have 2 or 3 as one of their prime factors
 * where i is odd.
 *
 *    { 0, 5, 7, 11, 13, ... 3i + 2, 3(i + 1) + 1, ..., N }
 *      0, 1, 2,  3,  4, ... list indices (0 is not used)
 */

SV* primes( long search_to )
{
    AV* av = newAV();

    if( search_to < 2 )
        return newRV_noinc( (SV*) av ); // Return an empty list ref.

    sieve_size_t i, j, q = (sieve_size_t) sqrt((double) search_to) / 3
+;
    sieve_size_t M = (sieve_size_t) search_to / 3;
    sieve_size_t c = 0, k = 1, t = 2, ij;

    // Allocate space. Set bits to 1. Unset bit 0.
    sieve_type primes( M + 2, 1 );
    primes[0] = 0;

    // Unset bits greater than search_to.
    if ( 3 * M + 2 > search_to + ((sieve_size_t)search_to & 1) )
        primes[M] = 0;
    if ( 3 * (M + 1) + 1 > search_to + ((sieve_size_t)search_to & 1) )
        primes[M + 1] = 0;

    // Clear composites.
    for ( i = 1; i <= q; i++ ) {
        k  = 3 - k, c = 4 * k * i + c, j = c;
        ij = 2 * i * (3 - k) + 1, t = 4 * k + t;
        if ( primes[i] ) {
            while ( j <= M ) {
                primes[j] = 0;
                j += ij, ij = t - ij;
            }
        }
    }

    // Gather primes.
    if( search_to >= 2 )
        av_push( av, newSVuv( 2UL ) );
    if( search_to >= 3 )
        av_push( av, newSVuv( 3UL ) );

    for ( i = 1; i <= M; i += 2 ) {
        if ( primes[i] )
          av_push( av, newSVuv(static_cast<unsigned long>(3 * i + 2)) 
+);
        if ( primes[i + 1] )
          av_push( av, newSVuv(static_cast<unsigned long>(3 * (i + 1) 
++ 1)) );
    }

    return newRV_noinc( (SV*) av );
}
[download]

Below is the time taken to find all prime numbers smaller than 1 billion.

my $primes = primes( 1_000_000_000 );

Sieve of Eratosthenes   4.879 seconds
Sieve of Xuedong Luo    3.751 seconds

There are 50,847,534 prime numbers between 1 and 1 billion.
[download]

I was so fascinated by Algorithm3 that I decided to parallelize it. It took me 5 weekends just to get the math to work. But I wanted faster and placed it aside. Two years later tried again and created Sandboxing with Perl + MCE + Inline::C. Also, tried Math::Prime::Util by Dana Jacobsen and primesieve by Kim Walisch.

Testing was done on a Haswell Core i7 Macbook Pro running at 2.6 GHz configured with 1600 MHz memory.

#
# Count primes
#

perl algorithm3.pl 1_000_000_000
Prime numbers : 50847534
Compute time  : 0.146 sec

perl primesieve.pl 1_000_000_000
Prime numbers : 50847534
Compute time  : 0.064 sec

perl primeutil.pl  1_000_000_000
Prime numbers : 50847534
Compute time  : 0.024 sec

#
# Sum primes ( 203_280_221 prime numbers )
#

perl algorithm3.pl 4_294_967_296 --sum
Sum of primes : 425649736193687430
Compute time  : 1.082 sec

perl primesieve.pl 4_294_967_296 --sum
Sum of primes : 425649736193687430
Compute time  : 0.369 sec

perl primeutil.pl  4_294_967_296 --sum
Sum of primes : 425649736193687430
Compute time  : 2.207 sec

#
# Print primes ( 2.0 GB, beware... )
#

perl algorithm3.pl 4_294_967_296 --print >/dev/null
Compute time  :  2.071 sec

perl primesieve.pl 4_294_967_296 --print >/dev/null
Compute time  :  1.395 sec

perl primeutil.pl  4_294_967_296 --print >/dev/null
Compute time  : 13.470 sec
[download]

Fast is possible in Perl. Thus, Perl is fun. One is not likely to print that many prime numbers. Math::Prime::Util is powerful with many features. Algorithm3 was mainly an exercise exploring Perl + MCE + Inline::C possibilities.

Update June 26, 2015 using Math::Prime::Util v0.51

The mce-sandbox was updated to call the new sum_primes/print_primes functions in Math::Prime::Util v0.51 for the primeutil.pl example.

Count primes

   $ perl algorithm3.pl 4294967296
   Prime numbers : 203280221
   Compute time  : 0.623 sec

   $ perl primesieve.pl 4294967296
   Prime numbers : 203280221
   Compute time  : 0.252 sec

   $ perl primeutil.pl 4294967296
   Prime numbers : 203280221
   Compute time  : 0.210 sec

Sum of primes

   $ perl algorithm3.pl 4294967296 --sum
   Sum of primes : 425649736193687430
   Compute time  : 1.090 sec

   $ perl primesieve.pl 4294967296 --sum
   Sum of primes : 425649736193687430
   Compute time  : 0.367 sec

   $ perl primeutil.pl 4294967296 --sum
   Sum of primes : 425649736193687430
   Compute time  : 0.768 sec

Print primes ( outputs 2 GB containing 2032802221 prime numbers )

   $ perl algorithm3.pl 4294967296 --print >/dev/null
   Compute time  : 2.086 sec

   $ perl primesieve.pl 4294967296 --print >/dev/null
   Compute time  : 1.397 sec

   $ perl primeutil.pl 4294967296 --print >/dev/null
   Compute time  : 1.925 sec
[download]

Kind regards, Mario