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;
}
-
Are you posting in the right place? Check out Where do I post X? to know for sure.
-
Posts may use any of the Perl Monks Approved HTML tags. Currently these include the following:
<code> <a> <b> <big>
<blockquote> <br /> <dd>
<dl> <dt> <em> <font>
<h1> <h2> <h3> <h4>
<h5> <h6> <hr /> <i>
<li> <nbsp> <ol> <p>
<small> <strike> <strong>
<sub> <sup> <table>
<td> <th> <tr> <tt>
<u> <ul>
-
Snippets of code should be wrapped in
<code> tags not
<pre> tags. In fact, <pre>
tags should generally be avoided. If they must
be used, extreme care should be
taken to ensure that their contents do not
have long lines (<70 chars), in order to prevent
horizontal scrolling (and possible janitor
intervention).
-
Want more info? How to link
or How to display code and escape characters
are good places to start.
|