#!perl -w
use strict;

# This list can be used to eliminate over 80% of odd guesses.
my @smallPrimes = ( 2,   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 ); 

# Use PreCheck to speed things up a little.
# This will eliminate a large percentage of guesses quickly,
# leaving only the more likely guesses for the slower algorithm.
# Note that it returns 0 if the input may be a larger prime
# then covered in @smallPrimes, but a non-zero return may mean
# that the number is a "known" prime in @smallPrimes or is 
# a composite with a prime factor in @smallPrimes.
sub PreCheck
{
    my $p = shift;
    for ( @smallPrimes )
    {
        return $_ if 0 == $p % $_;
    }
    return 0;
}

# Shamelessy stolen from [MrNobo1024]'s node: [id://56906]
sub gcd
{ 
    my($x,$y)=@_; 
    ($x,$y)=($y,$x%$y) while $y; 
    return $x;
}

# Doing $x ** $y can result in very large numbers,
# But if what you want it ( $x ** $y ) % $m, then
# you can pull the modulo into the products and keep
# your numbers smaller.
sub ModuloPwr
{ 
    my($base,$pwr,$mod) = @_;
    my $result = 1;
    while ( $pwr > 0 )
    {
        $result = ( $result * $base ) % $mod;
        --$pwr;
    }
    return $result;
}
 
# And the point of the exercise:
 
# Rabin-Miller aka Miller-Rabin 
# These notes from "Applied Cryptography 2nd ed" (Schneier)
#
# Choose a random number, p, to test. 
# Calculate b, where b is the number of times 2 divides p-1. 
# Then calculate m, such that p = 1 + 2^b * m.
#
# 1. Choose a random number, x, such that x is less than p.
# 2. Set j = 0, and set z = x^m mod p.
# 3. If z = 1, or if z = p - 1, then p passes the test and may be prim
+e.
# 4. If j > 0 and z = 1, then p is not prime.
# 5. Set j = j + 1. If j < b and z <> z^2 mod p go back to step 4.
#    If z = p - 1, then p passes the test and may be prime.
# 6. If j = b and z <> p -1, then p is not prime.
#
sub RabinMiller
{
    my $p = shift;
    my $b = 0;
    my $m = $p - 1;
    my $t = 20; # The probability of error is .25 ** $t
                # Probability of error with $t = 20 is lt 10 ** -12.

    while(  0 == ( 1 & $m )  )
    {
        $m >>= 1;
        ++$b;
    }
    print "B = $b, M = $m\n";

    TESTLOOP:
    for ( 0 .. $t ) 
    {
        my $x;
        do { $x = 1 + int rand( $p - 1 ) } while gcd($x,$p) > 1;
        print "Random X = $x\n";
        my $j = 0;    
        my $z = ModuloPwr( $x, $m, $p );
        next TESTLOOP if 1 == $z;    
        while( $b > ++$j )
        {
            return 0 if 1 == $z;
            next TESTLOOP if $p - 1 == $z;
            $z = ( $z * $z ) % $p;
        }
        if ( $p - 1 != $z )
        {
            # print "Not Prime at Z = $z. B($b) = J($j)\n";
            return 0;
        }
    }
    return 1;
}

__END__
Usage:

Call RabinMiller() with some number, $p, which you want to 
test for primality. If RabinMiller() returns false, then the
number is definatly a composite number. If it returns true, 
then the number may be prime, with a probability of 
error = .25 ** $t.

Bruce Schneier, author of "Applied Cryptography" (Wiley Press),
suggests that in practical applications it is best to weed out 
the obvious composites first. Calling PreCheck will verify that 
$p is non trivial by returning any prime factor less then 300. 
This will eliminate all even guesses and 80% of odd guesses. 
Note that if your number is a prime less then 300, PreCheck will retur
+n
it as a prime factor.