Pacific-Design.com

    
Home Index

1. Perl5

2. Inline C

Perl5 / Inline C /

#!/usr/bin/perl

use strict;
use warnings;
use feature         qw( say      );
use Readonly;

use Test::More;
use Benchmark       qw( cmpthese );

use Math::Prime::XS qw( primes   );

use Inline C => 'DATA';

#------------------------Subs to test and benchmark -------------------------

Readonly my %bench_subs => (
#   pl_div          =>  \&pl_div,       # Simple Division, pure Perl
    pl_erat         =>  \&pl_erat,      # Sieve of Eratosthenes, pure Perl
#    pl_vec_erat     =>  \&pl_vec_erat,  # Sieve with a bit vector, pure Perl
#   pl_m_erat       =>  \&pl_m_erat,    # Merlyn's bit vector sieve, pure Pl.
#   xs_mpxs_erat    =>  \&xs_mpxs_erat, # XS: Math::Prime::XS::primes()
    ext_cpp_erat    =>  \&ext_cpp_erat, # Extrnl call: primes.exe (C++,sieve)
#   ilc_av_div      =>  \&ilc_av_div,   # Inline::C, division, return AV
    ilc_av_erat     =>  \&ilc_av_erat,  # Inline::C, Sieve, return AV
#   ilc_stk_div     =>  \&ilc_stk_div,  # Inline::C, division, return on Stack
    ilc_stk_erat    =>  \&ilc_stk_erat, # Inline::C, Sieve, return on Stack
#   il_c_gbl_erat   => \&il_c_gbl_erat, # Inline::C, Sieve, direct C call.
);

# ------------------------ Testing Data Sets --------------------------------

Readonly my %known_quantities => (
    -5      => 0,       -1      => 0,       0       => 0,
    1       => 0,       2       => 1,       3       => 2,
    5       => 3,       7       => 4,       10      => 4,
    11      => 5,       13      => 6,       19      => 8,
    3_571   => 500,     100_000 => 9_592,   224_737 => 20_000,
);

Readonly my %known_primes_lists => (
    -1  =>  [],                     0   =>  [],
    1   =>  [],                     2   =>  [2],
    3   =>  [2,3],                  4   =>  [2,3],
    5   =>  [2,3,5],                6   =>  [2,3,5],
    7   =>  [2,3,5,7],              11  =>  [2,3,5,7,11],
    18  =>  [2,3,5,7,11,13,17],     19  =>  [2,3,5,7,11,13,17,19],
    20  =>  [2,3,5,7,11,13,17,19],
);

# ----------------------- Benchmark Data Sets -------------------------------

Readonly my $bench_time   => -5;    # - seconds.

# At 3.3M, my prime sieve and M::P::XS's prime sieve catch up to each other.
# Below 3.3M, my sieve is faster.  Above 3.3M, the M::P::XS version wins.
#Readonly my @bench_inputs => ( 3_300_000 );

Readonly my @bench_inputs => ( 5, 50_000, 250_000, 500_000 );

# ------------------------- Run the tests -----------------------------------

can_ok( 'main', keys %bench_subs ) 
    or BAIL_OUT( "can_ok failed.\n" );
    
while( my ( $name, $sref ) = each %bench_subs ) {
    note "Testing $name.";

    note "\tReturn Value List Sizes";
    while( my ( $n_test, $known_quantity ) = each %known_quantities ) {
        local $Bench::input = $n_test;
        is( 
            scalar @{ $sref->( ) }, 
            $known_quantity, 
            "$name( $n_test ) finds $known_quantity primes."
        ) or BAIL_OUT( "is() failed on $name( $n_test )\n" );
    }

    note"\tReturn Value List Correctness";
    while( my ( $n_test, $listref ) = each %known_primes_lists ) {
        local $Bench::input = $n_test;
        is_deeply( 
            $sref->( ),
            $listref,
            "$name( $n_test ) reports primes of @{$listref}"
        ) or BAIL_OUT( "is_deeply() failed on $name( $n_test ).\n" );
    }

}

done_testing();

# ------------------------- Benchmark the subs ------------------------------

say "\nComparing:\n\t",
    join( "\n\t", sort keys %bench_subs ),
    "\nComparison time: ", -$bench_time, " seconds.";

foreach my $bind_value ( @bench_inputs ) {
    local $Bench::input = $bind_value;
    say "\nInput parameter value of $bind_value";
    cmpthese(
        $bench_time,
        \%bench_subs
    );
}


# ------------------- Here's what we're here to see -------------------------

# Find all primes up to '$top' using the simple division test.
# Optimizations are only testing ints, and only testing up to
# sqrt($i).

sub pl_div {
    my $top = $_[0] // $Bench::input;
    my @primes;
    if( $top >= 2 ) { $primes[0] = 2; }
    DIV_PL_OUTER:
    for ( my $i = 3 ; $i <= $top ; $i += 2 ) {
        my $sqrt_i = sqrt($i);
        for ( my $j = 3 ; $j <= $sqrt_i ; $j += 2 ) {
            next DIV_PL_OUTER unless $i % $j;
        }
        push @primes, $i;
    }
    return \@primes;
}


# Sieve of Eratosthenes, in Perl -- Array, no bit-vectors.
sub pl_erat {
    my $top = ( $_[0] // $Bench::input ) + 1;
    return [] if $top < 2;
    my @primes = (1) x $top;
    my $i_times_j;    # Small optimization to eliminate calculating twice.
    # The sieve
    for my $i ( 2 .. sqrt $top ) {
        if ( $primes[$i] ) {
            for ( my $j = $i ; ( $i_times_j = $i * $j ) < $top ; $j++ ) {
                undef $primes[$i_times_j];
            }
        }
    }
    return [ grep { $primes[$_] } 2 .. $#primes ];
}

sub pl_vec_erat {
    my $top = ( $_[0] // $Bench::input ) + 1;
    return [ ] if $top < 2;
    my $primes = '';
    vec( $primes, $top, 1 ) = 0; 
    my $i_times_j;
    # The sieve
    for my $i ( 2 .. sqrt $top ) {
        if ( !vec( $primes, $i, 1 ) ) {
            for ( my $j = $i ; ( $i_times_j = $i * $j ) < $top ; $j++ ) {
                vec( $primes, $i_times_j, 1 ) = 1;
            }
        }
    }
    return [ grep { !vec( $primes, $_, 1 ) } 2 .. $top-1 ];
}


# Merlyn's Unix Review Column 26, June 1999
# http://www.stonehenge.com/merlyn/UnixReview/col26.html
# Modified to store results in an array rather than print.
sub pl_m_erat {
    my $top = ( $_[0] // $Bench::input );
    my $sieve = '';
    my @primes;
    GUESS:
    for ( my $guess = 2 ; $guess <= $top ; $guess++ ) {
        next GUESS if vec( $sieve, $guess, 1 );
        push @primes, $guess;
        for ( my $mults = $guess * $guess ; $mults <= $top ; $mults += $guess )
        {
            vec( $sieve, $mults, 1 ) = 1;
        }
    }
    return \@primes;
}

# A wrapper around the external executable compiled in C++.
# Receives a "big list" of primes via system pipe read.
# Uses Sieve of Eratosthenes implemented in C++.
sub ext_cpp_erat {
    my $top = $_[0] // $Bench::input;
    open my $fh, '-|', 'primes.exe ' . $top;
    local $/ = undef;
    chomp( my @primes = split /\n/, <$fh> );
    return \@primes;
}

# Thin wrapper to bind input params for the benchmark.
# Inline C, Return on stack, division method.
sub ilc_stk_div {
    my $top = $_[0] // $Bench::input;
    return [ il_c_primes_stk($top) ];
}

# Thin wrapper to bind input params for the benchmark.
# Inline C, Return array-ref, division method.
sub ilc_av_div {
    my $top = $_[0] // $Bench::input;
    return il_c_primes_av($top);
}


# Thin wrapper to bind input params for the benchmark.
# Inline C, Return array-ref, Sieve of Eratosthenes method.
sub ilc_av_erat {
    my $top = $_[0] // $Bench::input;
    # This sub already returns a ref.
    return il_c_eratos_primes_av($top);
}



# Using the Math::Prime::XS module (they still beat me!)
sub xs_mpxs_erat {
    my $top = $_[0] // $Bench::input;
    return [  ] if $top < 2;
    return [ primes( $top ) ];
}

# Thin to bind input params for the benchmark.
# Inline C, Return on stack, Sieve of Eratosthenes method.
sub ilc_stk_erat {
    my $top = $_[0] // $Bench::input;
    return [ il_c_eratos_primes_stk($top) ];
}

__DATA__

__C__


#include "math.h"

/* Find all primes up to 'search_to' using the basic division
 * test approach, optimized to avoid testing evens, and to 
 * stop testing a given int when sqrt(int) is reached.
 * Returns a big list.  The wrapper will have to convert it to an array-ref.
 */

void il_c_primes_stk( int search_to ) 
{
    Inline_Stack_Vars;
    Inline_Stack_Reset;
    if( search_to >= 2 ) 
        Inline_Stack_Push( sv_2mortal( newSViv( 2 ) ) );
    else
    {
        Inline_Stack_Done;
        return;
    }
    int i;
    for( i = 3; i <= search_to; i+=2 ) 
    {
        int sqrt_i = sqrt( i );
        int qualifies = 1;
        int j;
        for( j = 3; ( j<=sqrt_i ) && ( qualifies==1 ); j += 2 ) 
        {
            if( i % j == 0 ) qualifies = 0;
        }
        if( qualifies == 1 ) 
            Inline_Stack_Push( sv_2mortal( newSViv( i ) ) );
    }
    Inline_Stack_Done;
}


/* Find all primes up to 'search_to' using the basic division
 * test approach, optimized to avoid testing evens, and to 
 * stop testing a given int when sqrt(int) is reached.
 * Returns aref to big list.
 */

AV * il_c_primes_av( int search_to ) 
{
    AV* av = newAV();
    if( search_to >= 2 ) 
        av_push( av, newSViv( 2 ) );
    else 
        return sv_2mortal( av );
    int i;
    for( i = 3; i <= search_to; i+=2 ) 
    {
        int sqrt_i = sqrt( i );
        int qualifies = 1;
        int j;
        for( j = 3; ( j<=sqrt_i ) && ( qualifies==1 ); j += 2 ) 
        {
            if( i % j == 0 ) qualifies = 0;
        }
        if( qualifies == 1 ) 
            av_push( av, newSViv( i ) );
    }
    return sv_2mortal( av );
}


/* Find all primes up to 'search_to' using the Sieve of Eratosthenes.
 * This function returns a big list on The Stack.
 */

void il_c_eratos_primes_stk ( int search_to ) 
{
    Inline_Stack_Vars;
    Inline_Stack_Reset;
    bool* primes = 0;
    int i;
    if( search_to < 2 ) {
        Inline_Stack_Done;
        return;
    }
    Newxz( primes, search_to + 1 , bool );
    if( ! primes ) croak( "Failed to allocate memory.\n" );
    for( i = 2; i * i <= search_to; i++ ) 
    {
        if( !primes[i] ) 
        {
            int j;
            for( j = i; j * i <= search_to; j++ ) primes[ i * j ] = 1;
        }
    }
    for( i = 2; i <= search_to; i++ ) 
    {
        if( primes[i] == 0 ) Inline_Stack_Push( sv_2mortal( newSViv( i ) ) );
    }
    Safefree( primes );
    Inline_Stack_Done;
}


/* Find all primes up to 'search_to' using the Sieve of Eratosthenes.
 * This function returns an array-ref so that the wrapper doesn't have
 * to make the conversion.
 */

AV * il_c_eratos_primes_av ( int search_to ) 
{
    AV* av = newAV();
    bool* primes = 0;
    int i;
    if( search_to < 2 )
        return sv_2mortal( av );
    Newxz( primes, search_to + 1 , bool );
    if( ! primes ) croak( "Failed to allocate memory.\n" );
    for( i = 2; i * i <= search_to; i++ ) 
    {
        if( !primes[i] ) 
        {
            int j;
            for( j = i; j * i <= search_to; j++ ) primes[ i * j ] = 1;
        }
    }
    for( i = 2; i <= search_to; i++ ) 
    {
        if( primes[i] == 0 ) av_push( av, newSViv( i ) );
    }
    Safefree( primes );
    return sv_2mortal( av );
}

AV * il_c_gbl_erat() 
{
    int search_to = SvIV( get_sv( "Bench::input", 0 ) );
    AV* av = newAV();
    bool* primes = 0;
    int i;
    if( search_to < 2 )
        return sv_2mortal( av );
    Newxz( primes, search_to + 1 , bool );
    if( ! primes ) croak( "Failed to allocate memory.\n" );
    for( i = 2; i * i <= search_to; i++ ) 
    {
        if( !primes[i] ) 
        {
            int j;
            for( j = i; j * i <= search_to; j++ ) primes[ i * j ] = 1;
        }
    }
    for( i = 2; i <= search_to; i++ ) 
    {
        if( primes[i] == 0 ) av_push( av, newSViv( i ) );
    }
    Safefree( primes );
    return sv_2mortal( av );
}