http://qs321.pair.com?node_id=82388
Category: Cryptography
Author/Contact Info Adam - Feel free to msg me.
Description: These functions provide a fast approach to testing the primality of numbers. Very useful in public key cryptography. This is based on the Rabin-Miller algorithm.

Many thanks to MrNobo1024 (and Euclid) for gcd and to Tye for helping me debug.

#!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.
Replies are listed 'Best First'.
Re (tilly) 1: Rabin Miller
by tilly (Archbishop) on May 23, 2001 at 02:19 UTC
    The following will speed up your calculation significantly:
    # 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) = @_; return 1 if $pwr < 1; # I probably don't need the 0.1, but doesn't hurt. my $result = ModuloPwr($base, int(0.1 + $pwr/2), $mod); $result = ($result * $result) % $mod; if ($pwr%2) { $result = ( $result * $base ) % $mod; } return $result; }
    However there appears to be a bug somewhere. It works fine for small primes, but for instance 97397 is reported not a prime, while it is. (I checked and double-checked. My change does not cause the issue.)

    UPDATE
    I am a moron. The bug is simple and obvious. Square then reduce modulo the prime works fine until the prime is large enough to go outside of native integer arithmetic. At that point things go haywire.

      But Perl's native integer range goes to about 2**56 so I think I introduced the bug when I handed Adam some code that uses &1 instead of %2. I suspect that if you fix that then you'll be able to handle past 8-digit primes.

              - tye (but my friends call me "Tye")
Rabin Miller may be deterministic
by tilly (Archbishop) on Jul 14, 2001 at 12:40 UTC
    I randomly ran across an interesting note that if the Extended Riemann Hypothesis is true then with the exception of the easily recognized Carmichael numbers there is always a witness of size at most (log2n)2. This is a polynomial primality test.

    Another interesting point is that Rabin-Miller raises some philosophical questions about mathematics. The fact is that no mathematical operation can be carried out without error. Indeed this is why mathematics, often popularly supposed to be a bastion of pure truth, is in fact a subject like any other where people have to make informed decisions about what results they will and will not believe. (You can think of mathematical proof as being like trying to write and verify a computer program without ever having the opportunity to run it and see if it works. This is not far wrong.)

    So the question is at what point does showing that a result is likely true an acceptable substitute for showing that it is true? This question is not hypothetical. Suppose that you want to use RSA in practice. Well first you need to have 2 big primes. Well from the prime number theorem, primes are not hard to find. The density of the primes near n is 1/log(n). But it is not hard to just look at numbers that you know are congruent to 1 or 5 mod 6, which drops 2/3 of the composites from your net. If you are, say, interested in 100 digit numbers this means that about 1% of the numbers you look at are prime. Grab enough numbers and you probably have quite a few primes in there, the issue is figuring out which ones they are. Well the answer to that is to apply something like Rabin-Miller. If a number survives a dozen rounds of Rabin-Miller then the conditional odds of it being prime pass 99.999%. At some point most people will trust their credit card to it...

Re: Rabin Miller
by Adam (Vicar) on May 29, 2001 at 20:36 UTC
    I've just noticed that there is something on this algorithm in O'Reilly's Mastering Algorithms with Perl. (Miller-Rabin test, pgs 520-522). I don't actually have a copy of that book, but I am curious as to if it has an implementation, and if so then what differences there are with mine.
Re: Rabin Miller
by ambrus (Abbot) on Feb 09, 2009 at 16:53 UTC
Re: Rabin Miller
by Anonymous Monk on Mar 17, 2011 at 19:51 UTC
    This code will pass 4294967296 (2**32) as a prime!

    Quick fix:

    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. return 0 if ($p % 2 == 0 && $p > 2); .....