ADR-1453: float_moment_cuda adds the float squares the CPU adds, and is bit-identical while the CPU's own sum is exact¶

• Status: Accepted
• Date: 2026-10-02
• Deciders: lusoris
• Tags: cuda, gpu-parity, numerics, float-moment, testing, ci, rc3, fork-local

Context¶

float_moment reports the mean and the mean square of the reference and the distorted luma plane. moment.c adds the samples, and their squares, into one double per output and divides by the pixel count. picture_copy() has divided each sample by 1, 4, 16 or 256 before, and compute_2nd_moment() forms each square in float (const float term = pic_ * pic_).

float_moment_cuda accumulates four uint64 sums on the device and divides them on the host. Its second sums were the exact integer squares of the raw samples. Up to 12 bits a square has at most 24 significant bits, the float square is the integer square, and the twin returned the CPU's bits. At 16 bits the float square is the integer square rounded to 24 bits and the twin did not. ADR-1447 found and fixed the defect in the HIP twin, measured the CUDA twin on two fixtures and left it in T-GPU-FLOAT-MOMENT-16BIT-SQUARES-2026-10-02; ADR-1449 fixed the SYCL twin.

Measured on an RTX 4090 at --precision max against --backend cpu on master 2096bd1bb, second moments identical to the CPU's and their largest difference:

The repository's 16-bit Netflix fixture is 8-bit content shifted left, whose squares have few significant bits, which is why the parity tests never saw it. The first moments were identical everywhere.

Decision¶

We will make the 16-bit kernel add the CPU's term, as ADR-1447 and ADR-1449 did. moment_float_square() in integer_moment/moment_score.cu converts the raw sample to float, multiplies it by itself with __fmul_rn() (one fp32 product, rounded to nearest even, as on the CPU) and converts the result, an integer below 2^32, to unsigned long long. Dividing a sample by a power of two before squaring does not change which bits the rounding drops, so this integer is the CPU's term in units of 1 / scaler^2. sample_square<T>() selects it for uint16_t samples and keeps the integer square for uint8_t, where the two are the same number. The reduction, the readback and the host's two divisions stay as they are.

Every term is a multiple of the unit. The CPU's running double sum is therefore exact, and equal to the device's integer sum, while it is below 2^53 units. A term is below 2^32 units, so that holds for every frame of up to 2^21 = 2 097 152 pixels at any content (1920x1080 has 2 073 600) and for every frame at 8, 10 and 12 bits. scripts/ci/exact_twins.d/float_moment.cuda declares the twin exact for that range.

On a 16-bit frame of more than 2^21 pixels whose second moment times the pixel count reaches 2^37 the CPU's sum passes 2^53 units. From there the CPU rounds every add of a term that is not a multiple of the sum's last place, and the twin, which holds the exact sum and rounds once, can differ from it by at most the bound ADR-1447 derives,

(pixels - 2^21 + 1) / pixels * 2^(e - 69) + 2^-37

with e the binade of the sum in units. test_cuda_float_moment_parity asserts equality in the exact range and this bound past it, through the cases of core/test/float_moment_twin_parity.h that the SYCL test already uses.

Alternatives considered¶

Consequences¶

• Positive: measured on an RTX 4090 at --precision max, the four outputs of every frame equal --backend cpu on all 262 frames measured: the 250 frames of the table above and 12 frames of 40x40 to 64x64 noise (1048 outputs; before, the second moments of all 77 16-bit frames with real low bits differed). 17 of the 32 16-bit 3840x2160 frames are past 2^53 units and identical too.
• Positive: no measurable cost; the numbers are in the CUDA backend page.
• Negative: stored 16-bit float_moment_cuda second moments change by up to 1.0e-4.
• Neutral / follow-ups:
• Past 2^53 the twin is within the bound above and not bit-identical in general: 2.7e-7 on the test's 2560x1440 frame (bound 6.6e-6). The parity gate script compares an exact cell at 0 for every frame; it has no per-frame range, so a 16-bit fixture in that range would have to be added to it with this bound. T-HIP-FLOAT-MOMENT-PAST-2-53-2026-10-02 covers every twin.
• float_moment_metal stays in T-GPU-FLOAT-MOMENT-16BIT-SQUARES-2026-10-02.
• test_hip_float_moment_parity.c carries its own copy of the fixtures and the bound that float_moment_twin_parity.h holds for the SYCL and CUDA tests; moving the HIP test onto the header is left to the RC5 dedupe pass.
• Guards: test_cuda_float_moment_parity and _large (noise at 8, 10, 12 and 16 bits and a bright 16-bit 1920x1080 frame with ==; the 16-bit cases fail on the old twin; a 2560x1440 frame past 2^53 against the bound) and test_cuda_float_moment_exact_contract.py (five planted regressions, no device).

References¶

• req (coordinator brief for the CUDA lane, 2026-10-02): "Already known from the HIP lane's cross-check on the 4090: float_moment_cuda is 2.8e-5 / 1.0e-4 off at 16 bit (integer squares instead of the CPU's float squares; fixed for HIP in #1789 / ADR-1447 and SYCL in #1793 / ADR-1449, incl. the derived bound past 2^53: reuse their header and bound)".
• ADR-1447, ADR-1449, ADR-1212, ADR-1392, ADR-1403, ADR-1421, ADR-1428, ADR-0214.
• docs/state.md: T-CUDA-FLOAT-MOMENT-16BIT-SQUARES-2026-10-02, T-GPU-FLOAT-MOMENT-16BIT-SQUARES-2026-10-02, T-HIP-FLOAT-MOMENT-PAST-2-53-2026-10-02.