CUDA 13.2 Double-Precision Math Bugs on sm_80: Device Implementation Defects in jn/yn and Constant Folding Errors in llrint/lrint/copysign

Accelerated Computing GPU-Accelerated Libraries
1

Hi CUDA Team and Community,

We identified five correctness failures when testing CUDA double-precision math intrinsics on sm_80 (NVIDIA A30 / A100) using CUDA 13.2 (nvcc -O3).

These failures stem from two distinct root causes:

  1. Core Device Math Algorithm Defects (jn, yn): The underlying device math implementation contains logic bugs for boundary inputs, causing both dynamic runtime execution (arg) and compile-time evaluation (const) to return incorrect values.

  2. Compiler Constant-Folding / Optimization Discrepancies (llrint, lrint, copysign): Host-side constant folding during compilation yields results that diverge from correct dynamic runtime outputs.

Environment

  • CUDA Toolkit: CUDA 13.2 (nvcc -O3)

  • Target Architecture: sm_80 (NVIDIA A30 / A100)

Detailed Analysis of Observed Failures

Category 1: Device Implementation Defects (jn, yn)

Both dynamic runtime execution (arg) and compile-time evaluation (const) fail for these functions because the bug is located directly inside the device math algorithm (e.g., misidentifying near-zero or signed zero as pole/domain branch errors).

  • yn(2, -0.0) $\rightarrow$ Returns +inf instead of -inf

    • CUDA Specification Contract Violation: The CUDA Math API documentation explicitly states: yn(n, ±0) returns -inf.

    • Standard Compliance: C Standard (Annex F.10.8.5) and glibc both specify -inf (0xfff0000000000000).

    • Actual Output: Dynamic and static calls return +inf (0x7ff0000000000000), directly violating CUDA’s written API contract.

  • jn(2, +0.0 / -0.0 / 2^-1074) $\rightarrow$ Returns -qNaN instead of +0.0

    • Mathematical Definition & Standards: For order $n=2$, $J_2(0) = 0$, $J_2(-0) = J_2(+0) = +0$, and $J_2(2^{-1074}) \approx x^2/8$, which underflows to +0. C Standard (Annex F.10.8.4) and glibc specify +0 for $n>0$.

    • CUDA Specification: CUDA documentation explicitly lists all NaN-returning conditions ($x=\text{NaN}$, $n<0$, or $x=+\infty$). Inputs $x \in \{+0, -0, 2^{-1074}\}$ for $n=2$ do not match any specified NaN condition.

    • Actual Output: Dynamic and static calls return -qNaN, indicating that zero/underflow inputs incorrectly trigger a domain/pole branch check error inside the device implementation.

Category 2: Compiler Constant-Folding Discrepancies (llrint, lrint, copysign)

  • llrint / lrint (x = 0x1.0000000000001p-1):

    • Input value 0x3fe0000000000001ULL is slightly larger than 0.5.

    • Dynamic execution (arg) correctly rounds half-up to 1.

    • Compile-time constant folding (const) incorrectly truncates/folds to 0.

  • copysign (x = qNaN or sNaN, y = -1.0):

    • Dynamic execution with runtime arguments correctly sets the sign bit to negative (0xfff8... / 0xfff0...).

    • Mixed execution (x dynamic, y = -1.0 literal) loses the negative sign bit and returns a positive NaN (0x7ff8... / 0x7ff0...).

Reproducer & Details

Questions for NVIDIA

  1. Bug Confirmation: Are these five cases confirmed as bugs in CUDA 13.2 (either inside the device math algorithm or host-side constant folding)?

  2. Normative Reference Standard: In cases where host constant folding diverges from device execution, should dynamic device execution always be treated as the normative runtime baseline?

2

Thanks for this report. I can confirm the (1) yn/jn are library implementation bugs. We shall improve the implementations in the future releases.

The llrint/lrint and copysign are the effects of const-folding, and copysign is a new report.

>In cases where host constant folding diverges from device execution, should dynamic device execution always be treated as the normative runtime baseline?

It may serve as a good starting point, but in general it cannot be treated as normative. Runtime implementation may have bugs (like in above yn/jn). Additionally, user program may run into undefined or unspecified behavior cases for which runtime will (or will not) return some result, but that outcome cannot be relied upon in general.

Some math functions in CUDA are not specified to deliver one and only possible result, in current documentation we allow a range of possibilities limited by expected tolerances (ulp-level errors). For such functions I would suggest that if the compile-time result is within the specified tolerance, then it can be different from the runtime (which also complies with the tolerance) and the difference should not be treated as a bug.

3

I observe that llrint, lrint and copysign have a documented maximum ULP error of zero. Given that the compile-time constant-folded results lie outside this error tolerance, does this indicate that bugs exist in these three operators?

4

Yes, these are bugs and they were reported to the NVCC compiler, as the const-folding optimization is guilty here.

5

As an additional diagnostic method we recommend that users always check the effects of compiler optimizations by compiling with optimizations and with all device optimizations disabled (e.g. -G compiler option) and running tests in both cases. Please refer to: 5.5. Floating-Point Computation — CUDA Programming Guide , specifically, Compiler Flags and Optimizations section.

6

OK, thank you for your reply.