[rmlarsen]

Results, compared with fpmp-hyp-compare

fp32mp2 and fp64mp2 sqrt and rsqrt take one Newton step plus its second-order term; hypot, rhypot, norm3d, norm4d, rnorm3d, rnorm4d and fdim are computed in pair arithmetic instead of in double; and these functions, with fmax, fmin, copysign and signbit, get the C99 Annex F special values. Every comparison below is against the parent commit, fpmp-hyp-compare (6960b2b), not against main: the earlier branches of the stack also change fpmp_impl_divsqrt.h, and no run against main was made.

Special values. A probe of 30,249 special operands per type over 13 functions (sqrt, rsqrt, fabs, copysign, fdim, fmax, fmin, hypot, rhypot and the four norms), against MPFR, on the host and on the device (sm_75):

Failures of 30,249 checks	parent	this branch
fp32mp2, host	9,235	0
fp32mp2, device	5,465	0
fp64mp2, host	16,347	0
fp64mp2, device	8,023	0

• 5,450 of the failures in each fp32mp2 row (5,449 for fp64mp2) were a result of the right class with the wrong limbs, mostly a NaN result with a zero low limb, in fdim, hypot, rhypot and the norms.
• The others were results of the wrong class (fp32mp2 host: 3,781, of which 3,774 in the norms; device: 13) or a finite value beyond the tolerance (fp32mp2: 4 and 2, all in sqrt; fp64mp2: 1,141 on the host and 2,466 on the device, in sqrt, fdim, hypot, rhypot and the norms).
• fmax(+0, -0) is +0 and fmin(-0, +0) is -0; an infinite operand of hypot or a norm wins over a NaN; an infinite or NaN result has canonical limbs.

Accuracy. Largest relative error against MPFR in u² (u = 2⁻²⁴ for fp32mp2, 2⁻⁵³ for fp64mp2), host / device (sm_75), over 6 ranges of operands for sqrt and rsqrt and 8 for the others (on this branch 1,048,576 operands per range for fp32mp2 and 262,144 for fp64mp2):

	fp32mp2 parent	fp32mp2 this branch	fp64mp2 parent	fp64mp2 this branch
sqrt, operands ≥ 2⁻⁷⁸ (2⁻⁹¹⁶)	5.5 / 8.4	0.50 / 0.50	5.5 / 2.8	0.50 / 1.8
rsqrt	7.3 / 10.6	2.7 / 3.9	7.4 / 4.8	2.5 / 2.4
fdim	0.53 / 0.53	2.4 / 2.4	4.5e20 / 4.5e20	2.2 / 2.2
hypot	0.54 / 0.55	3.2 / 3.2	2.2e16 / 2.3e16	3.2 / 3.2
rhypot	0.55 / 0.54	3.9 / 4.4	2.4e16 / 1.7e16	4.5 / 3.8
norm3d	0.55 / 0.55	3.5 / 3.5	not finite / 2.5e16	4.3 / 4.3
rnorm3d	0.56 / 0.54	4.3 / 4.3	not finite / 1.7e16	3.9 / 3.4
norm4d	0.55 / 0.56	4.5 / 4.5	not finite / 2.3e16	3.9 / 4.4
rnorm4d	0.55 / 0.54	5.4 / 4.4	not finite / 1.8e16	4.5 / 4.2

• sqrt is within 0.5u² of the exact root over the whole range on the host, subnormal operands included (mean 0.12u²; the parent's mean is 0.60u²). On the parent, fp32mp2 operands in the lowest two ranges (subnormal, and normal ones below about 2⁻⁷⁸) gave results off by up to 3.5e7u² on the host, and on the device none of the 262,144 subnormal operands gave a finite result, for sqrt or for rsqrt (the parent's device figures in the table leave that range out); the fp64mp2 sqrt was off by up to 1.8e16u² there.
• The fp64mp2 fdim, hypot and norms carried 53 bits (errors of 1e16u² and more), and on the host the norms returned a non-finite result for 140,175 to 148,318 of the sampled operands; they are now within 4.5u².
• Less accurate than before: the fp32mp2 fdim, hypot and norms were computed in double and came out correctly rounded to within 0.56u². In pair arithmetic they err by up to 2.4u² (fdim), 3.2u² (hypot) and 3.5u² to 5.4u² (the others), with a mean of 0.14u² to 0.37u² in the worst range (0.10u² to 0.25u² before). In exchange the device code has no FP64 instruction.
• asinh and acosh add the square root to x with an error-free sum that does not assume which operand is larger: the high limb of the new root can round across that of x, and acosh(0x1.fffffep24 + 1) was off by 2⁻²⁹ without it. Their largest errors on this branch are 6.3u² and 7.4u² (1,048,576 pairs in each of 10 and 8 ranges).

Performance. ns per call on the host, parent → this branch (ratio above 1 is faster):

fp32mp2	x86-64	Pi 5	M4
sqrt	3.65 → 3.27 (1.12×)	9.10 → 8.43 (1.08×)	1.42 → 1.57 (0.90×)
rsqrt	3.33 → 3.73 (0.89×)	8.78 → 9.99 (0.88×)	1.31 → 1.47 (0.89×)
fdim	14.77 → 3.15 (4.69×)	17.89 → 6.49 (2.76×)	3.73 → 2.02 (1.85×)
hypot	17.50 → 9.18 (1.91×)	25.57 → 20.88 (1.22×)	4.88 → 4.35 (1.12×)
rhypot	19.08 → 9.90 (1.93×)	26.88 → 22.67 (1.19×)	5.08 → 4.21 (1.20×)
norm3d	20.41 → 11.09 (1.84×)	24.75 → 24.88 (0.99×)	4.70 → 5.55 (0.85×)
rnorm3d	21.36 → 11.69 (1.83×)	27.09 → 26.34 (1.03×)	4.87 → 5.30 (0.92×)
norm4d	26.08 → 13.58 (1.92×)	32.16 → 28.43 (1.13×)	5.67 → 6.67 (0.85×)
rnorm4d	26.62 → 14.30 (1.86×)	33.58 → 29.77 (1.13×)	5.85 → 6.62 (0.88×)

fp64mp2	x86-64	Pi 5	M4
sqrt	3.97 → 3.67 (1.08×)	11.56 → 10.22 (1.13×)	1.43 → 1.57 (0.91×)
rsqrt	3.71 → 4.08 (0.91×)	10.79 → 12.07 (0.89×)	1.43 → 1.46 (0.98×)
fdim	1.58 → 3.23 (0.49×)	5.00 → 6.62 (0.76×)	2.02 → 2.02 (1.00×)
hypot	2.60 → 9.72 (0.27×)	8.57 → 22.27 (0.38×)	2.02 → 4.48 (0.45×)
rhypot	2.89 → 10.53 (0.27×)	9.09 → 24.05 (0.38×)	2.04 → 4.28 (0.48×)
norm3d	1.54 → 11.59 (0.13×)	4.11 → 26.23 (0.16×)	1.34 → 5.51 (0.24×)
rnorm3d	2.57 → 12.37 (0.21×)	5.27 → 27.27 (0.19×)	1.34 → 5.37 (0.25×)
norm4d	1.72 → 14.22 (0.12×)	4.39 → 29.85 (0.15×)	1.36 → 6.69 (0.20×)
rnorm4d	1.81 → 15.07 (0.12×)	5.13 → 31.70 (0.16×)	1.39 → 6.50 (0.21×)

• x86-64: 11th Gen Intel(R) Core(TM) i5-11500T @ 1.50GHz, GCC -O3 -march=x86-64-v3. Pi 5: Cortex-A76, GCC 14 -O3 -mcpu=cortex-a76. M4: Apple M4, Apple clang 21.
• One pinned core, operands in [1/4, 4), the fastest round of two alternating passes on an idle machine.
• The Pi run of the committed tree was not kept. Its sqrt and rsqrt rows come from the run of the kernel layout that was committed, and its other rows from a run made before the last change to the sqrt and rsqrt wrappers; on x86-64 that change moved the same rows by at most 4% (hypot from 9.54 to 9.18 ns).
• Slower than before: the fp32mp2 rsqrt takes 12% to 14% longer on all three hosts (0.88× to 0.89×; the fp64mp2 one 0.89× to 0.98×), sqrt 10% longer on the M4, and the fp32mp2 norms 9% to 18% longer on the M4 (norm3d and norm4d 0.85×).
• The fp64mp2 fdim, hypot and norms take 2 to 8 times as long (up to 0.12× on x86-64): on the parent they were a single double operation, which is also why they carried only 53 bits.

SASS instructions per kernel on sm_89, fp32 + fp64 (total):

fp32mp2	parent	this branch
sqrt	18 + 0 (48)	55 + 0 (96)
rsqrt	17 + 0 (40)	45 + 0 (88)
fdim	22 + 9 (144)	26 + 0 (72)
hypot	23 + 22 (176)	61 + 0 (112)
rhypot	23 + 23 (176)	56 + 0 (104)
norm3d	31 + 31 (256)	78 + 0 (128)
rnorm3d	31 + 31 (248)	73 + 0 (120)
norm4d	39 + 39 (320)	95 + 0 (152)
rnorm4d	39 + 39 (312)	90 + 0 (144)
copysign	0 + 2 (40)	0 + 0 (40)

sqrt and rsqrt grow threefold: the second-order term is new, and the kernel is inlined a second time in the branch for scaled (tiny) operands. The norms, fdim and copysign lose every FP64 instruction and are shorter in total. The fp64mp2 kernels grow throughout (sqrt 33 → 70 fp64 instructions, hypot 17 → 72, norm4d 30 → 106). GPU timings on the B200/B300 are not in this branch.

Platform independence. A hash of the host results over every range of the nine functions is identical for GCC on x86-64, GCC 14 on Cortex-A76 and Apple clang 21 on the M4, for both types.

Implementation

• sqrt, rsqrt (fpmp_impl_divsqrt.h). From a single-limb seed y ≈ 1/√a (__fpmp_sqrt_seed in fpmp_impl.h: on the host the correctly rounded sqrtf and 0.5/y, on the device the approximate reciprocal square root after one single-limb step) the kernels take one Newton step and add its second-order term, with the residual formed exactly by fma.
• Checked wrappers. The public sqrt and rsqrt go through __fpmp2_sqrt_checked and __fpmp2_rsqrt_checked, which answer the IEEE special values on entry and scale operands below 2⁻⁷⁸ (fp32mp2) and 2⁻⁹¹⁶ (fp64mp2) by a power of four, where the kernel's residual would underflow. Internal callers keep the bare kernels.
• Norms. hypot, rhypot, norm3d, norm4d, rnorm3d and rnorm4d share one pair kernel, __internal_fpmp2_norm: the squares and their rounding errors are summed as a pair, scaled by a power of two where they could overflow or underflow, and the root or reciprocal root is taken with the kernels above.
• fdim is a pair subtraction (__fpmp2_high_sub).
• fmax, fmin, copysign, signbit. The zeros are ordered by sign, and the sign is tested with the float signbit, without a conversion to double on the device.
• asinh, acosh add the root to x with __fpmp_two_sum.

Testing

• New lit test fpmp/math_sqrt.pass.cpp, host and device, both types: the nine functions against MPFR references over the whole range of each (subnormal operands and results, and operands whose squares overflow or underflow, included), the results that must be exact, and the IEEE 754 special values, with those of fabs, copysign, fmax and fmin.
• lit on the GPU. cuda/fp/units/fpmp: 30 passed, 1 unsupported on sm_89 and on sm_75.
• Host. math_sqrt, math_hyperbolic and math pass built with Apple clang 21 on the M4.
• The special-value tables of these functions in fpmp_spec.rst are regenerated, and the error bounds and the list of dedicated functions in the documentation are updated.

Not run: MSVC, and the B200/B300 benchmarks.

Merge order

On top of fpmp-hyp-compare (f425f30), which is on top of fpmp-explog-compare (624b214) and NVIDIA#12568 (PR numbers are those of the upstream repository). This branch's changes are 00e56da and e71c6b0; the second declares the two __fpmp_sqrt_seed overloads _CCCL_HOST_DEVICE_API _CCCL_FORCEINLINE, so that a debugger can step into them.

Checklist

• New or existing tests cover these changes.
• The documentation is up to date with these changes.

🤖 Generated with Claude Code