[BUG]: fpmp2 math functions mishandle ±0, ±inf, NaN and out-of-domain inputs; spec tables echo the input #12550

Is this a duplicate?

• I confirmed there appear to be no duplicate issues for this bug and that I agree to the Code of Conduct

Type of Bug

Silent Failure

Component

libcu++

Describe the bug

Many cuda::fp32mp2 and cuda::fp64mp2 math functions return wrong results for special inputs: signed zeros, infinities, NaN, arguments outside the domain, subnormals and ±max. Often the result is a plausible finite number where C99 Annex F (and CUDA's own double functions) give NaN or ±inf, so nothing signals the error.

A probe of 50 unary functions (plus sincos and sincospi) at 19 special inputs finds 198 wrong results in 63 (type, function) pairs. It covers both types, on the host and on the device. All of them are on current main (89e1909). The open fpmp PRs #12519, #12530 and #12543 leave every one unchanged. pow is fixed separately in #12548.

Finite results where NaN or ±inf is due. These are the most serious, because nothing looks wrong:

• fp32mp2 log, log2 and log10 ignore the sign, and do not handle 0, inf or NaN:
  • log(x) for x < 0 returns log|x|, e.g. log(−2) = 0.693.
  • log(±0) = −104.7 (ln 2⁻¹⁵¹), log(+inf) = 88.72 and log(NaN) = 89.13.
  • log2(±0) = −151 and log2(+inf) = 128; log10 behaves the same way.
• fp32mp2 normcdfinv returns ±inf instead of NaN outside [0, 1], e.g. normcdfinv(2) = +inf.
• fp64mp2 boys_f0 returns 1 for every x < 1e-15, including x < 0 and −inf (F0(−1) = 1.4627).

NaN where ±inf or a finite value is due:

• fp32mp2 sinh and cosh at ±inf and ±FLT_MAX.
• fp32mp2 atan(±inf), which should be ±π/2.
• fp32mp2 boys_f0 at +inf (should be 0), and at −inf and −FLT_MAX (should be +inf).
• For both types: sqrt(+inf); rsqrt(±0), which should be ±inf; and rsqrt(+inf), which should be +0.
• For both types: ceil, floor, trunc, round, rint and nearbyint of ±inf.

Subnormal inputs, on the device only:

• fp32mp2 sqrt and rsqrt return NaN for subnormal inputs; for example sqrt(2⁻¹⁴⁹) should be 3.74e-23. The host is correct.
• fp64mp2 logb of a subnormal returns INT_MIN instead of −1074 … −1023.

Signed zero. Odd functions return +0 for −0:

• Both types: sin (and sincos), tan, asin, atan, sinh, tanh, erf, expm1, log1p, sqrt and round. On the device this also applies to the placeholders erfinv, j1, cyl_bessel_i1, rint and nearbyint, and tgamma(−0) returns +inf instead of −inf.
• fp32mp2 only: trunc does the same for −0.5 and for negative subnormals. On the device, ceil does too, and floor(−0) returns +0.
• fp64mp2 only: asinh, atanh, cbrt and sinpi (and sincospi) return +0, and rcbrt(−0) returns +inf instead of −inf. On the device, normcdfinv(0.5) returns −0.

Host and device disagree:

• fp32mp2 cospi(±0.5) returns −0 on the host.
• On the host, fp64mp2 sinpi and cospi return sin(πx) rounded in double: sinpi(1) = 1.2e-16, cospi(0.5) = 6.1e-17, and NaN at ±DBL_MAX. The device is exact.
• On the host, fp64mp2 normcdfinv returns ±inf outside [0, 1], and −inf for subnormal inputs.

The specification (docs/libcudacxx/extended_api/fp/fpmp_spec.rst) has 42 unary "Special Values Table"s, and all of them list the input as the result. For example, they give log(−1) = −1, sqrt(−1) = −1, exp(−INF) = −inf and sinh(+maxN) = +maxN, so they document neither the current behavior nor the correct one.

Not covered here:

• two-argument functions other than pow;
• arithmetic at the mid and high levels, where an infinite operand gives NaN (−inf + (−1), −inf · (−1)). That arithmetic is a separate design question.

How to Reproduce

#include <cuda/fpmp>

#include <cfloat>
#include <cmath>
#include <cstdio>

using T = cuda::fp32mp2;

__host__ __device__ void show(const char* expr, const char* expected, T r)
{
  printf("  %-22s = %-12g (expected %s)\n", expr, static_cast<double>(r.hi()), expected);
}

__host__ __device__ void probe()
{
  show("log(+0)", "-inf", cuda::log(T(0.0f)));
  show("log(-2)", "nan", cuda::log(T(-2.0f)));
  show("log(+inf)", "+inf", cuda::log(T(INFINITY)));
  show("log(nan)", "nan", cuda::log(T(NAN)));
  show("log2(-1)", "nan", cuda::log2(T(-1.0f)));
  show("sinh(+inf)", "+inf", cuda::sinh(T(INFINITY)));
  show("cosh(FLT_MAX)", "+inf", cuda::cosh(T(FLT_MAX)));
  show("atan(+inf)", "1.5708", cuda::atan(T(INFINITY)));
  show("boys_f0(+inf)", "0", cuda::boys_f0(T(INFINITY)));
  show("sqrt(+inf)", "+inf", cuda::sqrt(T(INFINITY)));
  show("sqrt(2^-149)", "3.74339e-23", cuda::sqrt(T(0x1p-149f)));
  show("rsqrt(+0)", "+inf", cuda::rsqrt(T(0.0f)));
  show("floor(-inf)", "-inf", cuda::floor(T(-INFINITY)));
  show("normcdfinv(2)", "nan", cuda::normcdfinv(T(2.0f)));
  show("sin(-0)", "-0", cuda::sin(T(-0.0f)));
  show("cospi(0.5)", "0", cuda::cospi(T(0.5f)));
}

__global__ void kernel()
{
  printf("device:\n");
  probe();
}

int main()
{
  printf("host:\n");
  probe();
  kernel<<<1, 1>>>();
  cudaDeviceSynchronize();
}

Built with nvcc -std=c++17 -arch=sm_89 -I <cccl>/libcudacxx/include repro.cu against main, it prints:

host:
  log(+0)                = -104.665     (expected -inf)
  log(-2)                = 0.693147     (expected nan)
  log(+inf)              = 88.7228      (expected +inf)
  log(nan)               = 89.1283      (expected nan)
  log2(-1)               = 0            (expected nan)
  sinh(+inf)             = -nan         (expected +inf)
  cosh(FLT_MAX)          = -nan         (expected +inf)
  atan(+inf)             = -nan         (expected 1.5708)
  boys_f0(+inf)          = -nan         (expected 0)
  sqrt(+inf)             = -nan         (expected +inf)
  sqrt(2^-149)           = 3.74339e-23  (expected 3.74339e-23)
  rsqrt(+0)              = -nan         (expected +inf)
  floor(-inf)            = -nan         (expected -inf)
  normcdfinv(2)          = inf          (expected nan)
  sin(-0)                = 0            (expected -0)
  cospi(0.5)             = -0           (expected 0)
device:
  log(+0)                = -104.665     (expected -inf)
  log(-2)                = 0.693147     (expected nan)
  log(+inf)              = 88.7228      (expected +inf)
  log(nan)               = 89.1283      (expected nan)
  log2(-1)               = 0            (expected nan)
  sinh(+inf)             = nan          (expected +inf)
  cosh(FLT_MAX)          = nan          (expected +inf)
  atan(+inf)             = nan          (expected 1.5708)
  boys_f0(+inf)          = nan          (expected 0)
  sqrt(+inf)             = nan          (expected +inf)
  sqrt(2^-149)           = nan          (expected 3.74339e-23)
  rsqrt(+0)              = nan          (expected +inf)
  floor(-inf)            = nan          (expected -inf)
  normcdfinv(2)          = inf          (expected nan)
  sin(-0)                = 0            (expected -0)
  cospi(0.5)             = 0            (expected 0)

Expected behavior

The C99 Annex F / IEEE 754 results, which are also what CUDA's float and double functions return for these inputs:

• log(x < 0) = NaN, log(±0) = −inf and log(+inf) = +inf;
• sinh(±inf) = ±inf, cosh(±inf) = +inf and atan(±inf) = ±π/2;
• sqrt(+inf) = +inf and rsqrt(±0) = ±inf;
• the rounding functions return ±inf for ±inf;
• odd functions return −0 for −0;
• normcdfinv returns NaN outside [0, 1];
• host and device give the same results.

One way to get there is to handle special inputs explicitly at each function's entry, as #12548 does for pow, instead of relying on the shared kernels to propagate them. The spec's special-value tables should then be regenerated from the implementation.

Reproduction link

No response

Operating System

Ubuntu 26.04.1 LTS (WSL2)

nvidia-smi output

NVIDIA GeForce RTX 4070 Laptop GPU (sm_89), driver 621.14

NVCC version

Cuda compilation tools, release 13.4, V13.4.92; host compiler GCC 15.2