Bug Report: Off-by-One Indexing in cufinufft GPU Kernel

[DiamonDinoia]

Issue

A long-standing bug in the interval function caused out-of-bound access in CUDA kernel arrays during spreading/interpolation. It was very hard to trigger — it required a non-uniform point to fall on a grid point to floating-point accuracy, i.e., x + ns/2 exactly an integer. In that case, the computed support window became ns + 1 wide instead of ns, resulting in illegal access to ker[ns].

Before summer 2024, this was masked by kernel arrays padded to MAX_NSPREAD = 16, so out-of-bound reads rarely affected anything unless ns = 16. After kernel arrays were resized to exactly ns, the bug became real, with invalid memory reads introducing small but nonzero garbage into the output. This mainly affected GPU type 3 transforms at tighter tolerances (eps = 1e-8), especially with upsampfac = 1.25.

Cause

• interval() returns a window of length ns + 1 when x + ns/2 is an integer (float precision)
• Results in access to ker[ns], out of bounds for T ker[ns]
• CUDA memory reuse can alias into adjacent arrays (e.g., ker3[0])
• Sanitizers missed it due to valid-but-wrong memory accesses

Fix

• Fixed interval() logic to enforce window size = ns
• Corrected all kernel index arithmetic to prevent ker[ns] access
• Reapplied ceil() fix for nf computation on GPU

Validation

• type 3 transform errors resolved
• With kerevalmeth=1:
  • abs error: 1e-15 to 1e-10
  • rel L2 error: ~2.31e-10

PR

• Fix in PR #662