Re: optimizing the miller-rabin algorithm

[rhesa]

The way I read the algorithm at http://en.wikipedia.org/wiki/Miller-Rabin_primality_test, you're looping too much in your middle loop. I believe you're supposed to choose a fixed accuracy k, and then a random set of k bases. You're using *all* numbers less than n as a base, and that is way too strong, especially for very large numbers, for which this algorithm is designed. If you're looping over all numbers less than n, you might as well do trial division :-)

I've tried to implement the algorithm as described on Wikipedia, and for a modest accuracy (say 25), it runs in less than 30 seconds on your entire number range.

Note that I'm using Math::Pari to verify the next prime number. If you don't have it installed, leave it out and do the check in some other way (against your file of primes, for instance).

I also have the GMP library installed, and use that with Math::BigInt. For small numbers (which 91000-93000 are), this may not have a big benefit.

#!/usr/bin/perl
use Math::BigInt lib => 'GMP';
use Math::PariInit qw/ nextprime /;

=pod

The algorithm can be written as follows:

    Inputs: n: a value to test for primality;
            k: a parameter that determines the accuracy of the test
    Output: composite if n is composite, otherwise probably prime
    write n - 1 as 2^s . d by factoring powers of 2 from n - 1
    repeat k times:

        pick a randomly in the range [1, n - 1]
        if a^d mod n != 1 and a^{2^r.d} mod n != -1 for all r
        in the range [0, s - 1] then return composite

    return probably prime 

=cut

$|++;

my $k = shift() || 25; # accuracy


for (my $n=91001; $n<93001; $n+=2) {
    
    my $d = $n - 1;
    my $s = 0;
    while( ! ($d & 1) ) { # using bit manipulation instead of division
+ by 2 (this is a bit of a premature optimization...)
        $d >>=1;
        $s++;
    }

    my $comp = 1;
    my $witn = 0;

  TEST:
    for( 1 .. $k ) { # test in k bases
        my $a = int(rand($n-1))+1; # select a random base between 1 an
+d n - 1
        my $b = Math::BigInt->new($a);

        if( $b->bmodpow($d, $n) == 1 ) { # a is not a witness for comp
+ositeness
            $comp = 0;  
            next TEST ;
        }

        for my $r( 0 .. $s-1 ) {
            if( $b->bmodpow(2**$r * $d, $n) == $n-1 ) { # a is not a w
+itness
                $comp = 0;
                next TEST ;
            }
        }
        $witn++;
    }

    my $p = nextprime($n-1);
    print "$n is probably prime (with $witn witnesses out of $k)\n" un
+less $comp;
    print "$n is not a prime, however: next prime = $p\n" if !$comp an
+d $n != $p;
    print "$n is a prime, but we missed it\n" if $comp and $n == $p;
}

[syphilis]

The way I read the algorithm ..., you're looping too much in your middle loop. I believe you're supposed to choose a fixed accuracy k, and then a random set of k bases. You're using *all* numbers less than n as a base, and that is way too strong, especially for very large numbers, for which this algorithm is designed.

Yes, that's pretty much right. Furthermore, for any number less than 341550071728322, if the Miller-Rabin test doesn't return "composite" for any of the bases 2, 3, 5, 7, 11, 13 and 17, then the number in question is proven to be prime. For numbers less than 1373653 it is sufficient to test only for bases 2 and 3. See http://primes.utm.edu/prove/prove2_3.html

For large numbers Math::BigInt's bmodpow (pure perl) function is just way too slow. You'll be wanting to use something like Math::Pari or Math::GMP (as rhesa has suggested).

Cheers,
Rob