[cufinufft] fold_rescale/interval fails in for float32 inputs on large 1D (~1e8) uniform grids

[blackwer]

While working on the review for PR #663, I came into a segfault bug occurring in index calculation from the combination of fold_rescale with interval. This happens roughly when N has more digits than a floating point number can accurately represent (~7 digits for float32, ~15 digits for float64). So the segfault is really only practically valid in 1D when using float32, though I think the issue is more broad than that. This bug is not only present in the new "output driven" method, though doesn't produce a segfault in the other methods, just spreads to the wrong(?) points. Looking at the spread_1d_nuptsdriven code (all comments are mine)

    const auto x_rescaled     = fold_rescale(x[idxnupts[i]], nf1); // float32 of O(nf1), possible loss of precision)
    const auto cnow           = c[idxnupts[i]];
    const auto [xstart, xend] = interval(ns, x_rescaled); // last digits of xstart possibly wrong
    const T x1                = (T)xstart - x_rescaled;
    if constexpr (KEREVALMETH == 1)
      eval_kernel_vec_horner<T, ns>(ker1, x1, sigma);
    else
      eval_kernel_vec<T, ns>(ker1, x1, es_c, es_beta);

    for (auto xx = xstart; xx <= xend; xx++) {
      // segfault averted by wrapping all values of xx, though density spreads to the wrong gridpoints
      auto ix          = xx < 0 ? xx + nf1 : (xx > nf1 - 1 ? xx - nf1 : xx); 
      const T kervalue = ker1[xx - xstart];
      atomicAdd(&fw[ix].x, cnow.x * kervalue);
      atomicAdd(&fw[ix].y, cnow.y * kervalue);
    }

Similar safety wrappers are written on every other fold_rescale operation to prevent illegal memory address access and properly handle periodic boundary issues. I think it's reasonable to ask if this is actually a bug or a feature, since ~100 million grid points is... a lot, and I'm unsure if spreading to slightly the wrong(?) grid points will do anything actually noticeable. That said, as written, what is normally the infrequent off-by-one errors in spreading can be considerably amplified to off-by-100 errors if the user requests large uniform grids for whatever reason.

Remediation:

• upcast to double?
  • I think this is almost equally wrong for low N -- it's just adding random digits to our float, but the added intermediate precision should make index-finding slightly more accurate. If my thinking is clear, it should at least fix the issue where large N leads to very wrong index calculation.
  • I have verified that somewhat frequently, even for low N, upcasting to double, not surprisingly, does give different index calculations
• Tell the user "don't do that" -- probably mathematically justified, but I'm not the expert here
• I don't know. Just need to bring this to discussion since it will give negative indices in the output-driven method as currently written, and it's worth considering how to handle these issues

[DiamonDinoia]

Maybe we could use some bitwise trickery if N is too big:

inline std::uint64_t fold_rescale_bitwise(float x, std::uint64_t N) noexcept
{
    float f = std::fmaf(x, INV_2PI_F, 0.5f);
    f -= std::floor(f);
    union { float f; std::uint32_t u; } v{f};
    // 24-bit mantissa (implicit leading 1)
    std::uint32_t mant = (v.u & 0x7FFFFFu) | 0x800000u;
    // unbiased exponent
    std::uint32_t exp  = (v.u >> 23) & 0xFFu;
    // exact floor(f * N) == floor( mant * N / 2^(23 + 127 - exp) )
    // so do one 64-bit multiply and one shift:
    const unsigned shift = 23 + (127u - exp);  // == 150 - exp
    return (std::uint64_t(mant) * N) >> shift;
}

[mreineck]

From my own experience I heartily recommend doing all location-dependent calculations in double precision; they are very subdominant in comparison to the actual kernel evaluations and grid accesses. In my library I even use long double for this purpose in double precision NUFFTs, and haven't noticed any slowdown.

[mreineck]

Independent of this, a 1D single-precision NUFFT wih 1e8 grid points and single -precision location calculations is just an expensive, but very bad random number generator ...

[DiamonDinoia]

Hi @mreineck,

Yes, doing them in double on the CPU is totally fine. On consumer GPU there might be slowdowns as sometimes double is 64x slower than single. I recommend some benchmarking to measure if it makes a difference.

Thanks,
Marco

[DiamonDinoia]

Hi @mreineck,

Is it necessary to use long double? Unless |x| is huge extracting the fractional part should not go too wrong.

[mreineck]

Yes, doing them in double on the CPU is totally fine. On consumer GPU there might be slowdowns as sometimes double is 64x slower than single. I recommend some benchmarking to measure if it makes a difference.

I agree that benchmarking this is important. On the GPU, and (unfortunately) especially in1D, it could be noticeable.

Is it necessary to use long double? Unless |x| is huge extracting the fractional part should not go too wrong.

That's more a philosophical question than a mathematical one, I fear, and it's one where @ahbarnett and I agree to disagree, so please take my opinion with a big grain of salt :)

Assuming that the input coordinates are single precision, and all calculations on them are done in single precision as well, the best accuracy (i.e. eps) you can have in a NUFFT, is 1e-7*<largest number of grid points in any dimension>.
(That's also why I wrote that under those circumstances a 1D NUFFT with a grid size of 1e8 makes no real sense). If you carry out the coordinate calculations in double precision instead (and you assume that you have no discretization errors in your input coordinates), this error goes down to about 1e-16*<largest number of grid points in any dimension>. You can see this phenomenon (for double precision) in the first plot of FIg. 4.1 of the original Finufft paper.
I want to avoid this effect as much as I can, so I even go to long double. @ahbarnett says that if we have to assume that the input coordinates are inaccurate to within a machine epsilon anyway, then that's just how things are and there is no need to make things look better.