The Rabin-Miller test is a fast method to determine either a) if a number is composite or b) a (arbitrarily high) probability that a number is prime. "Introduction to Algorithms" by Cormen, et al has a good section on this (they call it Miller-Rabin). Here's my (very lightly documented) implementation--

use warnings;
use strict;

use Getopt::Long;
use Math::BigInt;

my %pars = ( verify => 0 );
GetOptions ( \%pars, 'verify', 'help', 'usage' )
    or die "Bad GetOptions";

die "primes number [ more numbers ] \nreturns the number's primality s
+tatus"
    if $pars{help} || $pars{usage};

my @inputs = @ARGV 
    or die "no numbers given for primality test\n";

foreach ( @inputs ) {

    my $number = Math::BigInt->new ( $_ );
    my ( $primality_prob, $witness ) = miller_rabin_prime ( $number );

    print "$number is ",
          ( $primality_prob )
              ? "prime with prob $primality_prob\n"
              : "not prime, witness $witness\n";
}

sub miller_rabin_prime {

    my $n       = shift () or die "no n in MillerRabin";
    my $s       = shift () || 25;

    my $a = Math::BigInt->new ( 0 );

    for ( my $j = 1; $j <= $s; ++$j ) {

        # $a is a random int in [1,n-1]
        #
        my $ranlim = $n-1;
        do {
            $a = Math::BigInt->new ( sprintf "%d", rand ( $ranlim ) + 
+1 );
            $ranlim /= 2;
        } while ( 0 > $a );

        if ( 1 == Witness ( $a, $n ) ) {
            return 0, $a;
        }
    }
    my $prob = 1 - (3/4)**$s;
    return $prob, undef;
}

sub Witness {

    my $a = shift () or die "no a in Witness";
    my $n = shift () or die "no n in Witness";

    # n-1 = ( 2**t ) * u
    #
    my ( $t, $u ) = get_t_u ( $n - 1 );

    my @x;
    $x[0] = mod_exp ( $a, $u, $n );

    for ( my $i = 1; $i <= $t; ++$i ) {
        $x[$i] = mod_exp ( $x[$i-1], 2, $n );
        return 1
            if (1 == $x[$i]) && (1 != $x[$i-1]) && (( $n - 1 ) != $x[$
+i-1]);
    }
    return 1 
        if 1 != $x[$t];

    return 0;
}

sub get_t_u {

    my $m = shift () or die "no m in get_t_u";

    my ( $t, $u ) = ( 0, 0 );

    while ( 0 == $m % 2 ) {
        ++$t;
        $m /= 2;
    }
    $u = $m;
    return ( $t, $u );
}

sub mod_exp {

    # i**j mod n
    #
    # works with Math::BigInt in addition to regular scalars
    #
    my $i = shift () or die "no i in mod_exp";
    my $j = shift () or die "no j in mod_exp";
    my $n = shift () or die "no n in mod_exp";

    # return 1 for zero exponents
    #
    my $result = $i - $i + 1;
    return $result unless $j;

    my $pow2 = $i;

    while (1) {
        if ( $j % 2 ) {
            $result = ( $pow2 * $result ) % $n;
            return $result unless --$j;
        }
        $j /= 2;
        $pow2 = ( $pow2 * $pow2 ) % $n;
    }
}
[download]

---

---

Posts are HTML formatted. Put <p> </p> tags around your paragraphs. Put <code> </code> tags around your code and data!

Titles consisting of a single word are discouraged, and in most cases are disallowed outright.

Read Where should I post X? if you're not absolutely sure you're posting in the right place.

Please read these before you post! —

• How do I compose an effective node title?
• How do I post a question effectively?
• Markup in the Monastery

Posts may use any of the Perl Monks Approved HTML tags:

a, abbr, b, big, blockquote, br, caption, center, col, colgroup, dd, del, details, div, dl, dt, em, font, h1, h2, h3, h4, h5, h6, hr, i, ins, li, ol, p, pre, readmore, small, span, spoiler, strike, strong, sub, summary, sup, table, tbody, td, tfoot, th, thead, tr, tt, u, ul, wbr

You may need to use entities for some characters, as follows. (Exception: Within code tags, you can put the characters literally.)

	For:		Use:
	&		&
	<		&lt;
	>		&gt;
	[		&#91;
	]		&#93;

Link using PerlMonks shortcuts! What shortcuts can I use for linking?

See Writeup Formatting Tips and other pages linked from there for more info.