The fix seems to be to intervene just at the last minute to write 3 as the PRP base to disk on the final iteration, with an if statement inserted at about L2223:
In theory, you may run the Pépin test with any basea that has the correct Jacobi symbol (a/Fm) = –1, however there is no ability within Mlucas’ code or the worktodo format to modify the Pépin test to use a different base, and the test format only allows the exponent m of the Fermat number Fm to be chosen, e.g.:
Fermat,Test=30
This is not the case with the PRP worktodo format, so it is quite possible to run Mersenne cofactor testing using a different base; but the results will obviously not match the default values, typically base 3.
This is a weird necessity of the Fermat modular squaring code: at about
Mlucas.cL1212 Ernst set the PRP_BASE to 2, even though the starting base for the Pépin test is 3. When writing residues to disk though, the wrong base gets written and subsequently causes a problem for any ensuing cofactor testing. https://github.com/tdulcet/Mlucas/blob/ba49d68b8e62f272bf682f51f3ca77c1bb0e89fa/src/Mlucas.c#L1212-L1215The fix seems to be to intervene just at the last minute to write 3 as the PRP base to disk on the final iteration, with an
ifstatement inserted at about L2223:In theory, you may run the Pépin test with any base a that has the correct Jacobi symbol (a/Fm) = –1, however there is no ability within Mlucas’ code or the
worktodoformat to modify the Pépin test to use a different base, and the test format only allows the exponent m of the Fermat number Fm to be chosen, e.g.:This is not the case with the PRP
worktodoformat, so it is quite possible to run Mersenne cofactor testing using a different base; but the results will obviously not match the default values, typically base 3.