`is_pseudoprime()`: false-negative for some inputs

[trizen]

For some inputs, is_pseudoprime(n, b) seems to return 0, even if n is a Fermat pseudoprime to base b.

Example (OEIS: A020138):

use 5.014;
use Math::Prime::Util::GMP;

say Math::Prime::Util::GMP::is_pseudoprime(4, 9);   #=> 0 (incorrect)
say Math::Prime::Util::GMP::is_pseudoprime(8, 9);   #=> 0 (incorrect)

Notice that Math::Prime::Util works correctly for these inputs:

use 5.014;
use Math::Prime::Util;

say Math::Prime::Util::is_pseudoprime(4, 9);   #=> 1 (correct)
say Math::Prime::Util::is_pseudoprime(8, 9);   #=> 1 (correct)

[danaj]

1. An optimization intended for detecting primes was hard-coding the answer for single digit n. Removed.

2. In cases where base >= n, we didn't check for base mod n == {0,1,n-1} the way the MPU code does. Will investigate.

3. The MPU code has a very early exit for n=2,3 which is not right for n=3 with bases >= 3. Will fix.

[danaj]

Making sure all three (XS, PP, GMP) match. Same for is_euler_pseudoprime.

@trizen I'm trying to do the same with is_strong_pseudoprime(n,base) and there's an interesting question here. The XS code implements the basic definition. But that means things like 367 -- a prime -- will return 0 when given a base = 0 mod 367, e.g. 1101. This seems really wrong as we don't want to return 0 for primes even with screwy bases.

The GMP code made a decision to change this slightly, and return 1 for base = 0 mod n. This means primes will always have a 1 returned regardless of the base. This is important to do for things like large deterministic base sets, though obviously code doing that operation can do the check itself. I like the idea of primes always returning 1, but that doesn't strictly follow the definition.

I think most sources skim over this by insisting the base be less than n. Or in Pari/GP's ispseudoprime it doesn't take a user input base. Wolfram/Mathworld writes:

The algorithm proceeds as follows. Given an odd integer n, let n=2^rs+1 with s odd. Then choose a random integer a with 1<=a<=n-1. If a^s=1 (mod n) or a^(2^js)=-1 (mod n) for some 0<=j<=r-1, then n passes the test. A prime will pass the test for all a.

The restrict the base, and note primes always pass. Which happens because they don't let a = 0 mod n by definition.

Changing either code will make a change to the behavior. Thoughts?

1. Leave the code as is. I don't like how the GMP code has different behavior from XS, and which one should the PP code do?
2. Change GMP to act "normal". Go through other parts of the code, examples, etc. and make sure we use it correctly. Upside is we unconditionally apply the generic algorithm. Downside is of course returning 0 for primes with some bases.
3. Change XS to do the "return 1 (probably prime) for base = 0 mod n" instead of "return 0 (composite)". Has the nice effect of always returning 1 for primes regardless of base, but again it's a change. Could make some cases trying to do related OEIS sequences different.