Suyama PRP call doesn’t have correctly calculated Selfridge-Hurwitz residues

[xanthe-cat]

A quick example to demonstrate, with $M_{130439}$ which is semiprime with a small factor $260879$. The worktodo line is:

PRP=N/A,1,2,130439,-1,99,0,3,5,"260879"

My customised Mlucas test build prints S–H residues before and after the final modular division step:

Before mod div: Selfridge-Hurwitz residues Res64,Res35m1,Res36m1 = 0xFF122E035D795FEB, 5596574508,36217851704.
After mod div:  Selfridge-Hurwitz residues Res64,Res35m1,Res36m1 = 0x553AE8AB0A62D1C4, 8242891979,65108185009.

The Res64 values are correct, but the “After” Res35m1 and Res36m1 are not. For comparison, the output from cofact 0.74c:

                                        |      Selfridge - Hurwitz residues
                       mod 2^64 (hex)   |   mod 2^36    mod 2^36-1   mod 2^35-1
Calculated A residue matches proof file residue;
 Undivided  residue: 0xFF122E035D795FEB |  14453137387  36217851704   5596574508
  Divided   residue: 0x553AE8AB0A62D1C4 |  47418888644  65108186657   8674905291

Mlucas’ output from the Suyama call reports the Suyama (A) value but again only the Res64 is correctly calculated:

Suyama-PRP on cofactors of M130439: using FFT length 8K = 8192 8-byte floats.
Fermat-PRP residue (A)     = 0x553AE8AB0A62D1C4,          0,          0
Processed 17 bits in binary modpow; MaxErr = 0.000244141
3^(F-1) residue (B)        = 0x553AE8AB0A62D1C4, 8674905291,65108186657
(A - B) Res64 = 0x0000000000000000, C Res64 = 0xED5CFA1BF90E6C11
(A - B)/C: Quotient = 0, Remainder Res64 = 0x0000000000000000
This cofactor is PROBABLE PRIME [PRP39261].

Again, to compare with cofact:

Testing the cofactor for primality using the following 1 supplied factor: 260879 
Cofactor is 39261 digits long
                                        |      Selfridge - Hurwitz residues
                       mod 2^64 (hex)   |   mod 2^36    mod 2^36-1   mod 2^35-1
Suyama    A residue: 0x553AE8AB0A62D1C4 |  47418888644  65108186657   8674905291
Suyama    B residue: 0x553AE8AB0A62D1C4 |  47418888644  65108186657   8674905291
(A-B) mod C residue: 0x0000000000000000 |            0            0            0
M130439 cofactor is (probably) prime!

I tried looking at how the (A) residue is passed from the main Mlucas code to see why Mlucas reports zeroes for the Res35m1 and Res36m1 values but couldn’t make head or tail of fixing it.

[tdulcet]

Maybe try saving the residue to a savefile and then checking it with the --check option of my GIMPS status script. That would tell you if the full residue is somehow getting corrupted or if the Selfridge-Hurwitz residues are just being calculated incorrectly by Mlucas.

If it is the latter, the modular divide by 9 could change the size of the residue, so it could be that the res_SH() function calculating the Selfridge-Hurwitz residues is not correctly taking that into consideration. Maybe check if the len value being passed to the res_SH() function is correct: https://github.com/primesearch/Mlucas/blob/18398583da19f270eed22a036c27d9b6beb9973d/src/Mlucas.c#L5968-L5975

[xanthe-cat]

I think I have solved it: At the earlier point just after the modular division for Mersenne exponents, the arrtmp array is j in size, not i. At the later point during the call to the Suyama_CF_PRP() we add a limb of the initial if loop for Mersennes,

res_SH(ci,n,&itmp64,&Res35m1,&Res36m1);

The modular division ends up in ci (rather than arrtmp) in this part of the code, and since the Res64 is correct we use the itmp64 value instead.