Subject: | Math::Primality and semi-large numbers |

Multiplication in the Miller Rabin test is not done mod n, so will (1)
take far too long and consume huge amounts of memory, and (2) eventually
crash unable to continue as the number of digits in residue become
absurdly large.
I suggest replacing all the code inside the map { } of
is_strong_pseudoprime() with:
Rmpz_powm($residue, $residue, GMP->new(2), $n);
if (Rmpz_cmp($residue, $m) == 0) {
debug "$_:$residue == $m => spsp!" if DEBUG;
return 1;
}
Alternately one can do a multiply and mod, or pull the GMP '2' out of
the loop. For example:
time perl -Iblib/lib -Iblib/arch -MMath::Primality=:all -Mbigint -E '$n
= (10**25)->bstr; do { $n = next_prime($n); } for (1..10); say $n;'
Takes 9.5s with the old code, 0.04s with the new code, and of course
produces the same result.
10**48 runs in 0.06s with the new code. With the old code, my GMP
crashes after a bit over a minute, perhaps unable to get more than 32GB
of memory to store the product of two 3.2 million digit numbers.
The results with the new code run with 10 ** 100 and 10k next_primes
match those of Math::Prime::Util::GMP, which uses different code for
next_prime and the strong MR and strong Lucas tests. It also matches
Pari/GP 'n=10^100; for (x=1,10000,n=nextprime(n+1)); print(n)' which
uses completely different code.