[fpmp] Faster, platform-independent host seed for fp32mp2 cbrt and rcbrt by rmlarsen · Pull Request #12568 · NVIDIA/cccl

rmlarsen · 2026-11-10T05:28:19Z

Results, compared with `main`

On a host with hardware FMA, fp32mp2 cbrt and rcbrt seed their Halley step with a refined polynomial instead of exp2f(±log2f(r)/3). Against current main (07bbb7a):

Performance. ns per call, main → this PR (speedup; in parentheses this PR's share over #12548):

Platform	`cbrt`	`rcbrt`
x86-64: Intel Core i5-11500T, GCC 15.2 `-O3 -march=x86-64-v3`	33.7 → 24.4 (1.38×; 1.30×)	40.3 → 33.4 (1.21×; 1.16×)
aarch64: Raspberry Pi 5 (Cortex-A76), GCC 14 `-O3 -mcpu=cortex-a76`	57.1 → 43.4 (1.32×; 1.24×)	73.9 → 59.6 (1.24×; 1.12×)
aarch64: Apple M4 (Mac mini), Apple clang 21 `-O3 -mcpu=apple-m4`	11.8 → 2.7 (4.3×; 2.0×)	15.4 → 12.5 (1.23×; 1.14×)

One core, 2048 arguments log-uniform in 2^±20 with full lo limbs, the fastest of 25 rounds (M4: 150), three alternating passes that agree within 1%. Two configurations do not change:

Hosts without hardware FMA keep the old seed. With fmaf a library call (x86-64 at -march=x86-64), the polynomial seeds made cbrt 9% and rcbrt 14% slower than [fpmp] Accurate, C99-conforming and faster fp32mp2 pow #12548 (44.8 against 41.1 ns, 63.3 against 55.7 ns), so they are compiled only where FP_FAST_FMAF or __FMA__ is defined.
The device is untouched. It keeps the SFU lg2/ex2 pair; the SASS of a cbrt and an rcbrt kernel is identical to [fpmp] Accurate, C99-conforming and faster fp32mp2 pow #12548's, instruction for instruction, on sm_75 and sm_89.

Accuracy. Unchanged within the Halley step's own error. Maximum relative error on the host against MPFR, in u² (u = 2⁻²⁴):

	`main`	this PR
`cbrt`	4.4	3.6
`rcbrt`	2.5	2.6

Every float in [1, 8) as the hi limb with a random lo limb (25.2M inputs: all significands for each exponent residue mod 3); the same hi limbs with lo = 0 and every 509th positive float, subnormals included, give lower maxima. The mean is 0.27u² (cbrt) and 0.24u² (rcbrt) before and after.

Platform independence. With the new seed, host results are the same bits everywhere: a hash of the results over 29.4M inputs is identical for GCC 15.2 at -O3 and -O0 and Clang 21 on x86-64, GCC 14 on the Pi and Apple clang 21 on the M4. On main the seed comes from the C library, and the M4's results differ from glibc's.

Special values. Unchanged: ±0, ±inf and NaN return before the seed is computed.

Implementation

fpmp_math_impl_pow.h gains __internal_fpmp2_cbrt_seed and __internal_fpmp2_rcbrt_seed, which return r^(1/3) and r^(-1/3) in single precision for the reduced argument r in [1/4, 2). On the device, and on a host without hardware FMA, they are the expression used before. Otherwise:

cbrt: a degree-4 minimax polynomial (relative error ≤ 1.87e-3) and one single-precision Halley step whose residual is fma(s2_hi, s, -r) + s2_lo·s, with s2_hi + s2_lo = s² exact. The seed errs by at most 2⁻²³·⁹⁸.
rcbrt: a degree-3 polynomial (≤ 1.69e-2) and two division-free steps s(1 + u/3 + 2u²/9), u = 1 − r s³. The seed errs by at most 2⁻²³·²⁹.

Both bounds are measured against MPFR over every float in [1/4, 2). The polynomials are Sollya fpminimax fits with supnorm-certified errors; the commands are in the header. Every multiply-add is an explicit fma, so the result does not depend on the compiler's contraction.

The comments of cbrt and rcbrt gave the reduced range as roughly [2⁻¹, 2¹) and [2⁻², 2²); it is [1/4, 2) for both, and they now say so, since the polynomials rely on it.

Testing

New lit test fpmp/math_cbrt.pass.cpp, host and device: 16 inputs covering every exponent residue mod 3, both signs, the ends of the normal range and subnormal hi limbs; cbrt and rcbrt must be within 16u² of MPFR references held as double pairs. rcbrt had no host test before.
Host. math_cbrt, math, math_pow, math_reduction, accuracy and scalar_ops pass with GCC 15.2 and Clang 21 at -march=x86-64-v3 and -march=x86-64 (the old seed), GCC at -O0, GCC 14 on the Pi and Apple clang 21 on the M4.
lit on the GPU. The full cuda/fp/units suite passes on sm_89 (58 passed, 2 unsupported), CUDA 13.4 and GCC 15.2.
No warnings from the new code with GCC or Clang at -Wall -Wextra -Wshadow -Wconversion.

Not run: MSVC (it defines neither macro, so it keeps the old seed) and GPUs other than sm_89.

Merge order

Stacked on #12548 (#12543, #12519 and #12530 have landed): the branch contains its commits, and this PR's change is 395c0e6. It carries #12548's CI override, which is reset before merging.

Checklist

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

🤖 Generated with Claude Code

copy-pr-bot · 2026-11-10T05:28:22Z

This pull request requires additional validation before any workflows can run on NVIDIA's runners.

Pull request vetters can view their responsibilities here.

Contributors can view more details about this message here.

rmlarsen · 2026-11-13T19:11:19Z

/ok to test dc58109

…n kernel, and handle their special values All four find the angle of a point (M, m), 0 <= m <= M, from a table of atan(i/16): atan(m/M) = atan(c) + atan(t), c = i/16, t = (m - c M) / (M + c m), with i chosen so that -1/64 <= t <= 3/64 and the leading words of m and c M cancel exactly. The octant is folded into the table (four rows), and a 4-coefficient kernel replaces the 13- and 19-term compensated Horner evaluations. atan and atan2 take t from one pair division. asin and acos enter as the unit-circle point (x, sqrt((1 - |x|)(1 + |x|))), where the division becomes a product by a tabulated cos(atan(c)) and the kernel is asin; for small x that is x itself, so asin(x) does not depend on the square root there. Special values follow C99 Annex F and are answered on entry. asin(-0) and atan(-0) are -0 (were +0); atan(+-inf) is +-pi/2 (was NaN); atan2 takes the sign of a zero y from its sign bit (atan2(-0, x) was +0 or +pi), accepts subnormal operands (was NaN) and scales its operands by a common power of two, so that tiny or, on the device, huge operands keep full accuracy (the relative error reached 1); asin and acos of a pair just outside [-1, 1] by its low limb are NaN. NaN results carry NaN in both limbs. fp64mp2: asin(+-0), atan(+-0) and atan2(+-0, x) keep the sign of zero, and asin and acos of |x| >= 1/2 form 1 - |x| from the limbs instead of summing them in binary128 first, which lost a low limb beyond its 113 bits: acos((1, -1e-30)) was off by 2e-5 and acos((1, -denorm_min)) was 0; both are now within 2^-106. Against the head of NVIDIA#12568, host, time per call: asin 1.6-2.2x, acos 1.9-2.1x, atan 2.9x and atan2 3.3x faster on x86-64 (2.4x in the geometric mean), 1.6-2.2x on a Cortex-A76 (1.9x). fp32 SASS instructions on sm_89: asin 420 -> 94, acos 432 -> 95, atan 264 -> 76, atan2 287 -> 86. Largest error against MPFR in units of u^2 = 2^-48, host / sm_89 device: asin 6.2 / 6.2 (was 10.9 / 14.1), acos 7.2 / 7.0 (6.2 / 9.6), atan 4.1 / 4.1 (6.4 / 6.9), atan2 5.9 / 5.9 (6.1 on moderate ranges). asin on |x| < 0.5 is less accurate than before (6.2 against 2.2). Host results are the same bits with GCC and Clang, with and without FMA, on x86-64, Cortex-A76 and Apple M4; the host takes the square root from sqrtf and the device from its approximate reciprocal square root, so about 9% of asin and acos results differ between host and device in the last bits. New test fpmp/math_invtrig.pass.cpp; the special-value tables of the four functions in fpmp_spec.rst are regenerated. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>

…set before merging] This PR changes only the fpmp headers and their tests, so limit its CI to the cuda/fp lit suite on a T4 (gcc and clang, CTK 12.X and 13.X) and an H100 (gcc, CTK 13.X), plus an MSVC precompile of the same suite. The override must be reset to empty before the PR can merge. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

The override paired CTK 12.X with the bare 'gcc', 'clang' and 'msvc', which expand to GCC 16, Clang 22 and MSVC 14.50. CTK 12.9 supports none of them, so those three jobs failed before building anything. CTK 12.X now runs with gcc14, clang19 and msvc2022, as the weekly matrix does. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>

…en exp pow(a, b) = exp(b ln|a|), and exp turns an absolute error in its argument into the same relative error in its result. ln|a| and b ln|a| were rounded to fp32mp2's own precision, so the error grew with |b ln|a||: up to 767u^2 (u = 2^-24) at 88. The special values also departed from C99 Annex F: -0 lost its sign, -inf with a non-integer exponent and +-FLT_MAX^+-FLT_MAX gave NaN, integer exponents above 2^24 were not recognized when the integer lived in the lo limb, and results below FLT_MIN flushed to zero. ln|a| is now formed to about 2^-58 as three floats: |a| = 2^k z, a 128-entry table of 1/c (at most 9 bits, so that fma(z_hi, 1/c, -1) is exact) and ln c, with c = 1 for the bucket around 1, and ln(1 + r) with r exact as three floats. b ln|a| comes from three exact products, is reduced by ln2/128, and exp uses a 128-entry table of 2^(j/128) and a degree-5 polynomial. The result is scaled in two steps so that it underflows gradually. Special values are handled explicitly before the main path, and overflow and underflow are decided from b_hi * ln|a| before the products, so no intermediate is infinite or NaN. Against binary128 the error is now at most 2.9u^2 (mean 0.29u^2) at every magnitude of b ln|a|, and all 685 special-value cases match glibc's pow. Results are bit-identical on the host and the device. pow is 2.4-3.8x faster on the tested CPUs and 1.5x faster on an RTX 4070 Laptop GPU when the inputs are read coalesced; where every input is itself a scattered gather, the two table lookups compete for L1 and it is about 15% slower. The pow special-values table of fpmp_spec.rst is regenerated; its accuracy and performance rows still describe the previous implementation. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…l bases The pow lookup tables were inline variables with default visibility, so a shared library that calls cuda::pow exported both as unique global symbols, and two libraries built against different versions of the tables would read one copy. They are now hidden, like the functions that read them. math_pow.pass.cpp gains two accuracy cases on the scaled path of ln|a|, a subnormal base and a tiny base whose lo limb is subnormal, and includes <cuda/std/limits> for numeric_limits. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

The special-value test compared only the hi limb where pow returns exactly 1 or -1, such as pow(NaN, 0), pow(1, NaN) and pow(-1, +-inf), so a nonzero lo limb would have changed the value hi + lo and still passed. Those results must now have a zero lo limb, like the zero, infinite and subnormal ones. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>

In a device pass without RDC, _CCCL_GLOBAL_CONSTANT has internal linkage, and nvcc with a Clang host rejects a visibility attribute there (error NVIDIA#1407-D), so every cuda/fp test failed to compile with CTK 13.4 and Clang 22. The tables now take _CCCL_VISIBILITY_HIDDEN outside the device pass only, which is where a shared library would export them. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>

On a host whose fmaf is a hardware instruction, the single-precision seed of the cbrt / rcbrt Halley step is now a minimax polynomial in the reduced argument r in [1/4, 2), refined in single precision: one Halley step with an fma residual for cbrt, two division-free steps for rcbrt. It replaces exp2f(+-log2f(r) / 3), which stays on hosts where fmaf is a library call (the polynomial seeds are slower there). The device keeps the SFU lg2/ex2 pair, and its generated code is unchanged. The seeds err by at most 2^-23.98 (cbrt) and 2^-23.29 (rcbrt) over every float in [1/4, 2). Every multiply-add in them is an explicit fma, so those host results no longer depend on the C library or on the compiler's contraction. math_cbrt.pass.cpp checks cbrt and rcbrt against MPFR references across the argument reduction, on the host and on the device. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>

…n kernel, and handle their special values All four find the angle of a point (M, m), 0 <= m <= M, from a table of atan(i/16): atan(m/M) = atan(c) + atan(t), c = i/16, t = (m - c M) / (M + c m), with i chosen so that -1/64 <= t <= 3/64 and the leading words of m and c M cancel exactly. The octant is folded into the table (four rows), and a 4-coefficient kernel replaces the 13- and 19-term compensated Horner evaluations. atan and atan2 take t from one pair division. asin and acos enter as the unit-circle point (x, sqrt((1 - |x|)(1 + |x|))), where the division becomes a product by a tabulated cos(atan(c)) and the kernel is asin; for small x that is x itself, so asin(x) does not depend on the square root there. Special values follow C99 Annex F and are answered on entry. asin(-0) and atan(-0) are -0 (were +0); atan(+-inf) is +-pi/2 (was NaN); atan2 takes the sign of a zero y from its sign bit (atan2(-0, x) was +0 or +pi), accepts subnormal operands (was NaN) and scales its operands by a common power of two, so that tiny or, on the device, huge operands keep full accuracy (the relative error reached 1); asin and acos of a pair just outside [-1, 1] by its low limb are NaN. NaN results carry NaN in both limbs. fp64mp2: asin(+-0), atan(+-0) and atan2(+-0, x) keep the sign of zero, and asin and acos of |x| >= 1/2 form 1 - |x| from the limbs instead of summing them in binary128 first, which lost a low limb beyond its 113 bits: acos((1, -1e-30)) was off by 2e-5 and acos((1, -denorm_min)) was 0; both are now within 2^-106. Against the head of NVIDIA#12568, host, time per call: asin 1.6-2.2x, acos 1.9-2.1x, atan 2.9x and atan2 3.3x faster on x86-64 (2.4x in the geometric mean), 1.6-2.2x on a Cortex-A76 (1.9x). fp32 SASS instructions on sm_89: asin 420 -> 94, acos 432 -> 95, atan 264 -> 76, atan2 287 -> 86. Largest error against MPFR in units of u^2 = 2^-48, host / sm_89 device: asin 6.2 / 6.2 (was 10.9 / 14.1), acos 7.2 / 7.0 (6.2 / 9.6), atan 4.1 / 4.1 (6.4 / 6.9), atan2 5.9 / 5.9 (6.1 on moderate ranges). asin on |x| < 0.5 is less accurate than before (6.2 against 2.2). Host results are the same bits with GCC and Clang, with and without FMA, on x86-64, Cortex-A76 and Apple M4; the host takes the square root from sqrtf and the device from its approximate reciprocal square root, so about 9% of asin and acos results differ between host and device in the last bits. New test fpmp/math_invtrig.pass.cpp; the special-value tables of the four functions in fpmp_spec.rst are regenerated. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>

github-project-automation Bot added this to CCCL Nov 10, 2026

github-project-automation Bot moved this to Todo in CCCL Nov 10, 2026

cccl-authenticator-app Bot moved this from Todo to In Progress in CCCL Nov 10, 2026

rmlarsen force-pushed the fpmp-cbrt-host-seed branch from 653ded7 to dc58109 Compare November 13, 2026 18:27

rmlarsen mentioned this pull request Nov 13, 2026

[fpmp] Faster fp32mp2 math functions; accurate trig near multiples of π/2 and erfc at large x #12530

Merged

1 of 2 tasks

This was referenced Nov 13, 2026

[fpmp] Single-limb kernels for wide constants that fit one limb rmlarsen/cccl#7

Closed

[fpmp] Single-limb kernels for wide constants that fit one limb; 6-operation single-limb product at high #12640

Open

rmlarsen and others added 7 commits November 13, 2026 15:29

rmlarsen force-pushed the fpmp-cbrt-host-seed branch from dc58109 to 395c0e6 Compare November 13, 2026 23:31

rmlarsen mentioned this pull request Nov 13, 2026

[fpmp] Accurate, C99-conforming and faster fp32mp2 pow #12548

Open

2 tasks done

Navigation Menu

Repository navigation

[fpmp] Faster, platform-independent host seed for fp32mp2 cbrt and rcbrt - #12568

[fpmp] Faster, platform-independent host seed for fp32mp2 cbrt and rcbrt#12568
rmlarsen wants to merge 7 commits into
NVIDIA:mainfrom
rmlarsen:fpmp-cbrt-host-seed

rmlarsen commented Nov 10, 2026 •

edited

Loading

Uh oh!

copy-pr-bot Bot commented Nov 10, 2026

Uh oh!

rmlarsen commented Nov 13, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Conversation

rmlarsen commented Nov 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Results, compared with main

Implementation

Testing

Merge order

Checklist

Uh oh!

copy-pr-bot Bot commented Nov 10, 2026

Uh oh!

rmlarsen commented Nov 13, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

rmlarsen commented Nov 10, 2026 •

edited

Loading

Results, compared with `main`