Fix wrong residues for FP32+M61 and M31+M61 FFTs at low BPW (i96 sloppy carry word)…

[wadeghcp]

Symptom

Master at 4d0e759, make with g++ 13.3 (Ubuntu 24.04):

prpll -device 2 -ll 3021377        # M37, a known prime
...
3021377 FFT: 256K 2:256:2:256:101 (11.53 bpw)
3021377     20000 d7979bd162116bee
...
3021377   3021375 a0ec9abad5afd56e
{"status":"C", ..., "res64":"a0ec9abad5afd56e", "fft-type":"FP32+M61", "fft-length":262144, "shift-count":0}

NVIDIA RTX 3070 Ti, driver 595.84. The same residue on every rerun, so it is deterministic. -prp 3021377 on the same config fails its Gerbicz check at the first block twice and stops with "consistent error".

The interim residue after 20,000 squarings is checkable against other implementations:

implementation	res64 after 20,000 squarings
prime95 v30.19 b20, InterimResidues=20000 (its "iteration 20002" line; it counts the seed as iteration 2)	CE811E3772129824

Bits-per-word sweep, type 2 at 256K on the 3070 Ti, fp64 as reference: wrong at 11.53 and 12.59 bpw, correct from 13.35 up (11.53 to 19.07 tested). The boundary is not sharp: with -use STATS=1 the 13.35 case fails as well (both -ll and -prp).

Cause

-carry long (the expanded carryA/carryB kernels) is correct; the fused kernel is not. The fused kernel differs by using weightAndCarryPairSloppy, i.e. carryStepUnsignedSloppy(i96) / carryStepSignedSloppy(i96). In carryStepSignedSloppy(i96 x, i64 *outCarry, ...) the 32-bit "sloppy word" branch is taken whenever SLOPPY_MAXBPW >= 3200 (MAXBPW is 3457 for 2:256:2:256) with no lower bound on BPW, so for these FFT types the second word of every pair is a full 32-bit word whenever |x| >= 2^31. The transform itself copes with 32-bit words (that is what MAXBPW ~34 promises), but the convolution digits then come out around 2^60-2^63 instead of ~2^(2BPW+9), and the pair carry is correspondingly larger. That carry is i64 through the shuttle, fine, and is then added to the next line's first word in carryFinal, whose carry is a hard-coded i32 tmpCarry computed by carryStep(i64, i32*) as xtract32(x, nBits): it needs about 63 - 2nBits bits, which exceeds 31 once nBits <= 15. The high bits are silently dropped, i.e. a multiple of 2^(32+nBits) is lost at that word, and the residue is wrong from then on.

Measured with a probe in the fused kernel (max log2 |(u.x + inCarry) >> nBits| over the sampled iterations, reported through the ROE stat), unpatched:

| BPW | bits tmpCarry must hold | fits i32? | result | |---|---|---|---|
| 11.53 | 39.5 | no | wrong (Gerbicz EE) |
| 12.59 | 36.2 | no | wrong |
| 13.35 | 33.9 | no | wrong |
| 14.50 | 29.5 | yes | correct |
| 17.93 | 23.4 | yes | correct |

Same probe with the guard below: 9.3-9.5 bits at 11.53, 13.35 and 17.93 BPW.

The i64 x, i32 *outCarry variant of carryStepSignedSloppy already guards its equivalent branch with EXP / NWORDS >= 23 and carries the comment "the sloppy case fails if BPW is too low ... Probably ... would require CARRY_LONG to work properly"; this is that failure, in the i96 variant, with the mechanism identified.

Fix

Give the i96 variant the same EXP / NWORDS >= 23 guard as the i64 variant. Below that the branch falls through to carryStep, as it already does when SLOPPY_MAXBPW < 3200. The requirement derived above is nBits >= 16; 23 keeps the same margin the i64 variant uses. One condition change plus a comment; no other file touched. (The alternative would be widening tmpCarry in carryFinal to i64, which would keep the 32-bit sloppy word at low BPW; I left that to you since it touches the common path.)

With the fix, same device, same commands:

test	result
type 2 at 11.53, 12.59, 13.35 bpw, 20,000 iterations	ce811e3772129824 / matches fp64 / matches fp64
type 2 at 13.35 bpw with -use STATS=1	matches fp64

Timing at 11.53 bpw on the 3070 Ti: 60 us/iteration before, 59 after.

Why it is worth fixing even though production never hits it

The code is new and unreleased, wavefront exponents sit near 19 bpw where the i32 temporary fits, and PRP with the Gerbicz check stops on it, so no result was ever at risk. What breaks is the small case: the quick known-prime run (M37 in three minutes) that people use to decide whether to trust a new GPU, a new driver, or a long test on the same machine. A self-check that fails only at the small end is the worst place for a weak link, because that is the case everything else gets trusted on. With the guard, LL and PRP on the smallest exponents the program accepts are correct, and the fast path is untouched where it was already correct.

[wadeghcp]

@gwoltman this is in the sloppy-carry paths from #340; the probe table is in the description.