#!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 prime. # 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 return it as a prime factor.