#!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 nonzero 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:
# RabinMiller aka MillerRabin
# 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 p1.
# 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.
