You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
{{ message }}
Repository navigation
[fpmp] Faster, more accurate fpmp2 arithmetic with single-limb scalars and a shorter pair product - #12519
Arithmetic between an fpmp2 and a scalar that one limb holds exactly no longer converts the scalar to a pair first. The pair product is also shorter, and its high level more accurate. All figures are against current main (89e1909). This PR is stacked on #12543, which force-inlines the fpmp2 operators and should merge first (see Merge order). u² is 2⁻⁴⁸ for fp32mp2 and 2⁻¹⁰⁶ for fp64mp2.
x * c is 1.65–1.76× faster at low and 1.42–1.51× at mid (the default).
x + c is 1.9–2.0× faster at high and 1.1× at mid; x += c is 1.27–1.34× faster at high, x / c 1.04–1.19×, and x * y 1.10–1.16× at mid.
Two operations get slower, both at high and both in exchange for accuracy (below): x * c runs at 0.86–0.90× and x * y at 0.97–1.02×.
CPU (one x86-64 core): a dependent fp32mp2x * y takes 3.17 ns instead of 3.78 (1.19×), and x * c 2.77 ns instead of 3.77 (1.36×). At high, x * c runs at 0.68–0.76× and x * y at 0.81–0.83×. On a DGX Spark's Cortex-X925 (GCC 13), x * y takes 3.13 ns instead of 3.37 (1.07×) and x * c 2.58 ns instead of 3.36 (1.30×); at high they run at 0.87× and 0.73×.
fp32mp2 math functions, which pick up the kernels without changes of their own: the geometric mean over 30 functions is 1.028–1.061× faster, depending on the platform, GPUs from sm_75 to sm_121 included. The slowest function is 0.97× (asin on the Intel host), and the largest gains are 1.1–1.3× in asinh, acosh, atanh, pow, sinh and cosh.
Accuracy. Largest relative error, main → this PR, fp32mp2 / fp64mp2:
Operation
Level
main
This PR
First-order bound
x * c
low, mid
2.97 / 2.95u²
1.98 / 1.99u²
2u²
x * c
high
2.97 / 2.95u²
1.50 / 1.50u²
2u²
x * y
mid
5.00 / 4.92u²
5.00 / 4.99u²
6u²
x * y
high
5.00 / 4.92u²
3.99 / 4.00u²
5u²
x + c, x - c, c - x and x / c return the same values as on main.
Of the 31 fp32mp2 math functions, only asin, acos and pow return different values on the device, and only in the last bits. Their largest error against binary128 is unchanged: 9.9, 9.1 and 329u².
Host and device now agree more often. With GCC's default contraction, every product returns the same bits on the host as on the device; main's mid product disagreed on 20% of operands. Over the math functions (31 functions × 2¹⁸ inputs, GCC 15.2 -O2 -march=x86-64-v3), host results that differ from the device fall from 534,891 to 494,113, with cosh going from 21,135 to 4,344 and pow from 19,653 to 4,675. With -ffp-contract=off the count is unchanged (300,525 → 300,617).
The same bits on every GPU. Every operator result is bit-identical on sm_89, sm_90, sm_100, sm_103, sm_107, sm_120 and sm_121 (2²⁰ random operand pairs for each of x * c, x * y, x + c, x / c and x * x, at every level, for both types), on main and here, so the figures above hold on all of them.
On a DGX Spark's Cortex-X925 (GCC 13.3, default contraction), every product now returns the device's bits; on main, x * y at mid and high differed in 22% of results. The fp32mp2low quotient still differs in 13%, as on main.
Special values. No result moves away from IEEE 754.
Arithmetic. The test grid is the spec's 15 × 15 special operands (signed zeros, subnormals, ±min, ±1, ±max, infinities and NaN). It covers every scalar operator (x ± c, c - x, x * c, x / c, c / x and the compound forms), plus x * y and x / y, at every accuracy level, for both types, on the host and on the device.
The high limb changes in 9 cells of x * c (and of x *= c) at mid, for both types, and each change is to the IEEE result. In those cells x is a positive operand of magnitude at most min and c a negative one, and the product now underflows to −0 instead of +0.
The mirrored cells (x negative, c positive) still give +0, as x * y does on main and here.
Every other change is to a low limb only: the sign of an exactly zero low limb, ±inf instead of NaN next to an infinite high limb in an overflowing low quotient, and NaN payloads.
Math functions. On the 17 special inputs of each function, results are bit-identical to main's. The one exception, boys_f0(NaN), returns a subnormal instead of NaN on both main and this PR; only its bits differ. That value is undefined: a NaN reaches boys_f0's last branch with an uninitialized variable ([fpmp] Faster fp32mp2 math functions; accurate trig near multiples of π/2 and erfc at large x #12530).
Implementation
Single-limb kernels. Some scalars c fit exactly in one limb: the limb type, float for fp64mp2, and integers up to the limb's significand width. Such scalars now take dedicated kernels in the mixed operators, in the compound +=, -=, *= and /=, and in add/sub/mul/div<Acc>. Any other scalar (a double with fp32mp2, say) is still converted to a pair.
x + c, x - c, c - x: DWPlusFP, the existing accumulate. It gives the same value as the pair addition, in 10 instead of 20 instructions at high and 10 instead of 11 at mid.
x * c: TwoProd, one FMA and Fast2Sum: 6 instead of 9 instructions at mid, 3 instead of 5 at low. With x_hi and c scaled into [1, 2), the FMA rounds a value below 4u, so the error is below 2u² to first order. high uses Algorithm 7 of Joldes, Muller & Popescu (2017), in 10 instructions.
x / c: the existing division without the FMA on the divisor's zero low limb. It gives the same value in one instruction fewer.
Pair product.
mid sums the cross terms with one FMA (fma(x_lo, y_hi, x_hi * y_lo)), off the path of the exact residual of x_hi * y_hi, then adds that sum to the residual. In a chain of products x_lo arrives last and enters only that FMA. The product takes 8 instructions instead of 9, with the measured error of Algorithm 10.
high adds x_lo * y_lo and the cross terms, smallest first, through two FMAs off the residual's path, in 9 instructions.
Chaining two FMAs onto the residual would save one more instruction, but it lengthens the dependent path and measured 9% slower on the CPU.
__fpmp2_high_acc drops a second Fast2Sum that never changed its already normalized input: 13 → 10 instructions, bit-identical (0 differences over 5·10⁷ operand pairs, cancellation included).
Host/device agreement. No sequence feeds a separately rounded product into an add, so a host compiler has nothing to contract.
Where a product must stay unfused, it goes through the new __fpmp_mul_rn_nocontract, which computes x * y as an FMA with a +0 addend. Clang folds a −0 addend back into a multiply, so the addend has to be +0. The low quotients and the high scalar product use it.
On the device this changes only the sign of an exactly zero low limb, which is now +0.
sinh and cosh. Formed unfused, the low limb of their 0.5/e term is subnormal for |x| > 70, and Intel cores take a microcode assist for every product with a subnormal result.
The new __internal_fpmp2_half_inv_exp therefore returns 0 for that term once e = exp(|x|) exceeds 2³², where the term is under 2⁻⁶³ of the result. It also clamps the divisor there, so no subnormal is formed.
On sm_89, sinh and cosh return the same bits as on main over 2¹⁸ inputs each.
fpmp2_stat forwards its mixed and compound scalar operators to the wrapped type's. It accepts exactly the scalars fpmp2 accepts, and stays bit-identical to it.
Its deleted scalar compound operators are _CCCL_HOST_DEVICE_API, like its other deleted members. Under clang-CUDA, a host-only deleted overload is not viable in device code, so a call would otherwise fall through to the narrowing single-limb constructor that the deletion exists to block.
The new kernels are inline compositions of primitives and of entry points that exist in both the inline and the library mode, so CCCL_FPMP_LIB needs no new libcufp symbols. The changed pair products and __fpmp2_high_acc reach library users when libcufp is rebuilt from these headers.
Measurements
Accuracy method. Errors were measured on the host against an exact reference: the error accumulated from error-free transforms. The inputs were 5·10⁷ normalized operand pairs from random and structured regimes, plus a hill-climbing search for the worst case, for general operands and for squares. Squares have a maximum of 4.00u² at mid and 3.50u² at high, against 4.00u² on main. Math-function errors are against binary128 (libquadmath) over the 2¹⁸ inputs per function described in #12530.
FP-pipe SASS instructions on sm_100, main → this PR:
low
mid
high
x * c
5 → 3
9 → 6
9 → 10
x + c, x - c
8 → 7
11 → 10
20 → 10
x += c
7
10
13 → 10
x / c
−1
−1
unchanged
x * y
5
9 → 8
9
mad
13
16
29 → 28
GPU throughput, this PR over main, fp32mp2 / fp64mp2. Each thread runs 8 independent chains, with 8 blocks of 256 threads per SM; the figure is the best of two alternated runs, with CUDA 13.4. The operations that did not change (x += c at low and mid, x * y at low) measured 0.99–1.05.
Operation
Level
sm_89 (Ada)
sm_120 (Blackwell)
sm_75 (Turing)
x * c
low
1.76 / 1.67
1.69 / 1.66
1.65 / 1.67
x * c
mid
1.46 / 1.50
1.50 / 1.50
1.49 / 1.50
x * c
high
0.90 / 0.90
0.89 / 0.90
0.89 / 0.90
x + c
mid
1.08 / 1.10
1.10 / 1.10
1.09 / 1.10
x + c
high
2.02 / 2.00
1.98 / 2.00
1.99 / 2.00
x += c
high
1.33 / 1.30
1.33 / 1.30
1.30 / 1.30
x / c
mid
1.12 / 1.06
1.07 / 1.09
1.08 / 1.09
x * y
mid
1.16 / 1.12
1.10 / 1.12
1.13 / 1.12
x * y
high
0.97 / 1.00
0.99 / 1.00
0.99 / 1.00
CPU latency on an AVX-512 x86-64 workstation CPU (AMD Ryzen Threadripper PRO 7965WX, one core, GCC 15.2 -O3 -march=native). The figures are the latency of one dependent product, as the Google Benchmark median of 10 repetitions; the coefficient of variation is below 1.7%.
Product
float
double
x * y on main (Algorithm 10 as GCC compiles it)
3.78 ns
4.57 ns
x * y here (mid)
3.17 ns
4.02 ns
x * c on main (converted to a pair)
3.77 ns
4.59 ns
x * c here (mid)
2.77 ns
3.45 ns
CPU throughput, this PR over main, for the operators in a loop of 1 dependent chain and of 8 or 16 chains, which GCC vectorizes into one or two dependent vector chains. The figures are for fp32mp2; the fp64mp2 ratios are within 0.15 of them.
Operation
Level
1 chain
8 chains
16 chains
x * c
low
2.40
1.48
1.38
x * c
mid
1.36
0.97
0.95
x * c
high
0.76
0.69
0.68
x + c
high
1.55
1.35
1.31
x * y
mid
1.19
1.04
1.03
x * y
high
0.83
0.81
0.82
More GPUs, this PR over main, fp32mp2 / fp64mp2, the same benchmark: H100 (sm_90), B200 (sm_100), B300 (sm_103) and Rubin (sm_107), built with CUDA 13.4 and GCC 13, and the GB10 GPU of a DGX Spark (sm_121), built with its CUDA 13.0. Clocks could not be locked; the figure is the best of three alternated runs. The unchanged operations measured 0.99–1.00, and an H200 and a second Rubin system give the same ratios within 0.01.
Operation
Level
H100 (sm_90)
B200 (sm_100)
B300 (sm_103)
Rubin (sm_107)
GB10 (sm_121)
x * c
low
1.76 / 1.66
1.72 / 1.65
1.72 / 1.67
1.68 / 1.67
1.66 / 1.67
x * c
mid
1.42 / 1.46
1.50 / 1.46
1.50 / 1.50
1.51 / 1.50
1.49 / 1.50
x * c
high
0.87 / 0.89
0.86 / 0.89
0.86 / 0.90
0.89 / 0.90
0.90 / 0.90
x + c
mid
1.09 / 1.09
1.10 / 1.09
1.10 / 1.10
1.10 / 1.10
1.09 / 1.10
x + c
high
1.99 / 1.97
1.99 / 1.97
1.99 / 2.00
1.94 / 2.00
1.99 / 2.00
x += c
high
1.34 / 1.31
1.29 / 1.29
1.29 / 1.30
1.27 / 1.30
1.33 / 1.30
x / c
mid
1.10 / 1.19
1.16 / 1.09
1.16 / 1.09
1.04 / 1.09
1.06 / 1.09
x * y
mid
1.10 / 1.12
1.12 / 1.13
1.12 / 1.13
1.11 / 1.13
1.11 / 1.12
x * y
high
0.98 / 0.99
1.02 / 1.00
1.02 / 1.00
0.99 / 1.00
1.00 / 1.00
CPU latency on a DGX Spark: one dependent fp32mp2 operation on one Cortex-X925 core, GCC 13.3 -O3 -mcpu=neoverse-v2 (GCC 13 has no -mcpu=cortex-x925), best of three alternated runs; fp64mp2 takes the same time. With 8 or 16 chains, GCC 13 vectorizes main's and this PR's loops differently, and their ratios swing between 0.55× and 3.5× with the build options, so they are not shown.
Operation
Level
main
This PR
Speedup
x * c
low
2.06 ns
1.03 ns
2.00
x * c
mid
3.36 ns
2.58 ns
1.30
x * c
high
3.36 ns
4.63 ns
0.73
x + c
mid
3.73 ns
3.60 ns
1.04
x + c
high
5.78 ns
3.60 ns
1.61
x += c
high
5.14 ns
3.60 ns
1.43
x / c
mid
8.49 ns
7.71 ns
1.10
x * y
mid
3.37 ns
3.13 ns
1.07
x * y
high
3.37 ns
3.88 ns
0.87
Both slowdowns are at high, and both are accuracy trades. x * c there is Algorithm 7, which halves the error of the pair product (2.97 → 1.50u²). x * y there is the 9-instruction product, which takes the error from 5.0 to 4.0u².
fp32mp2 math functions, speedup over main (old time / new time). The figure is the geometric mean over 30 functions, from the harness and platforms described in #12530; the "with #12530" column adds that PR, which is stacked on this one.
aarch64: DGX Spark Cortex-X925, GCC 13.3, over #12543, inlined / called out of line
1.028 / 1.030
1.096 / 1.093
Per function, the largest gains are in asinh and acosh (1.11–1.31), atanh (1.04–1.21) and pow (1.02–1.18). On the i5-11500T, sinh and cosh are also 1.18× faster each.
The slowest results are asin 0.97× and acos 0.97× on the i5-11500T, acos 0.97× on the Pi, and cos 0.98× on the M4.
New cuda/fp/units/fpmp/scalar_ops.pass.cpp, run on the host and the device for all six types, checks:
x + c and x * c within 2u², and x * y within its first-order bound, of an exact expansion of error-free transforms;
x + c and x / c equal to the pair operation;
the compound forms, both operand orders and the integer scalars;
pinned operands on which main's product errs by 2.9u². These fail on main and pass here.
accuracy_mixed.pass.cpp: its stated contract now says that mul with a single-limb scalar is not bit-identical to the wrapped call. Its operands give exact products, so the checks themselves are unchanged.
fptool/stat.pass.cpp checks that all compound forms match fpmp2's, for a single-limb scalar and for one that needs both limbs (16777217), and that both types accept the same scalars. The new stat_implicit_casts.pass.cpp checks a double under CCCL_FPMP_EXPLICIT_CASTS=0.
Bit-exact comparisons against main:
the special-operand grid above, for both types and all levels, on the host and on sm_89;
the 31 fp32mp2 math functions × 2¹⁸ inputs, on sm_89 and on the host (GCC 15.2, -march=x86-64-v3, with and without -ffp-contract=off).
The same bits on more GPUs. The operator results above and the 33 fp32mp2 math functions × 2¹⁸ inputs are bit-identical on sm_89, sm_90, sm_100, sm_103, sm_107, sm_120 and sm_121, on main and here (boys_f0(NaN) aside).
It also passed on sm_75 and sm_120 without the low-quotient and sinh/cosh parts of items 4 and 5.
On sm_120, fpemu/charbool.pass.cpp fails the same way on main.
Documentation
Updated:
the error-bounds rows;
the mixed-arithmetic and compound-assignment text in fpmp.rst;
the optimization guide's cost table;
the spec's SASS counts for mul, acc and mad.
Not regenerated: the spec's measured accuracy and performance tables for multiplication (B300/B200). They need the spec harness. The PR's own benchmarks on H100, B200, B300, Rubin and DGX Spark are under Measurements.
Merge order
Stacked on #12543, which force-inlines the fpmp2 operators and math functions and should merge first. This branch merges it (278c4ae, 9bc91f2), so its diff includes #12543's commits until #12543 lands. With it, the new scalar operators of this PR are force-inlined like the others, with _CCCL_HOST_DEVICE_API _CCCL_FORCEINLINE (b8b4f8b). #12543's CI override is reset, so this PR runs the full CI.
Checklist
New or existing tests cover these changes.
The documentation is up to date with these changes.
…ir product
Arithmetic between an fpmp2 and a scalar that one limb holds exactly (the limb
type, float for fp64mp2, integers up to the limb's significand width) no longer
promotes the scalar to a pair:
- x + c, x - c, c - x use DWPlusFP (the existing accumulate): the same value as
the pair addition, at 10 instead of 20 instructions at high.
- x * c is TwoProd plus one fma and a Fast2Sum, 6 instead of 9 instructions at
mid, with error below 2u^2 to first order where the promoted product reaches
3u^2; high uses JMP17 Algorithm 7 (measured 1.5u^2).
- x / c drops the term on the divisor's zero low limb (same value).
- The compound forms +=, -=, *=, /= and the accuracy-selecting add/sub/mul/div
take the same kernels; fpmp2_stat forwards to them to stay bit-identical.
The default pair product folds both cross terms into the exact residual with
two fmas (7 instead of 9 instructions, first-order bound 5u^2, the bound of
Algorithm 10); high is the low product plus renormalization (8 instructions,
measured 3.9u^2 against 4.9u^2). __fpmp2_high_acc drops a second Fast2Sum that
never changed its already normalized result.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
It looks like this branch is under active development. To avoid overwhelming you with review comments due to an influx of new commits, CodeRabbit has automatically paused this review. You can configure this behavior by changing the reviews.auto_review.auto_pause_after_reviewed_commits setting.
Mixed arithmetic uses a direct scalar path when a value is exactly representable in one limb; other values continue to use pair-based arithmetic.
Compound addition, subtraction, multiplication, and division support exactly representable one-limb scalars. High-accuracy multiplication and scalar division are supported.
fp32mp2 hyperbolic functions return a zero half-inverse-exponent result for sufficiently large exponent values.
Documentation
Updated performance, instruction-count, and error-bound guidance for scalar and pair operations.
Tests
Expanded coverage of scalar arithmetic across accuracy levels, operand orders, compound operations, and statistical parity.
fpmp2 now uses single-limb kernels for eligible scalar arithmetic and adds scalar compound assignments. The change updates multiplication and division kernels, extends fpmp2_stat forwarding and tracing, changes fp32mp2 hyperbolic calculations, and revises tests, documentation, math wrapper annotations, and CI overrides.
Changes: fpmp2 Arithmetic
Layer / File(s)
Summary
Arithmetic kernels and bounds libcudacxx/include/cuda/__fp/fpmp_impl.h, libcudacxx/include/cuda/__fp/fpmp_impl_muladd.h, libcudacxx/include/cuda/__fp/fpmp_impl_divsqrt.h, libcudacxx/test/libcudacxx/cuda/fp/units/fpmp/scalar_ops.pass.cpp, docs/libcudacxx/extended_api/fp/fpmp_error_bounds.rst, docs/libcudacxx/extended_api/fp/fpmp_spec.rst
Adds scalar-limb add, multiply, and divide kernels. Pair multiplication and accumulation implementations change. Low-accuracy division uses a non-contracting residual product. Tests and documentation give accuracy bounds and instruction counts.
Mixed and compound operations use scalar kernels when one limb holds the scalar. Other scalar inputs use pair arithmetic. Tests cover scalar operations across accuracy levels. Documentation describes overloads, instruction counts, and error bounds.
Scalar operations in fpmp2_stat libcudacxx/include/cuda/__fp/fptool_stat.h, libcudacxx/test/libcudacxx/cuda/fp/units/fptool/stat*.pass.cpp
Adds scalar compound-assignment forwarding and support checks to fpmp2_stat. Mixed arithmetic calls wrapped-type overloads and traces operations. Tests compare assignment support, arithmetic results, implicit-cast behavior, and operation counts.
Adds a helper that returns zero for large exponential high limbs and otherwise computes 0.5/e. The fp32mp2 sinh and cosh implementations use the helper.
Math wrapper annotations libcudacxx/include/cuda/__fp/fpmp_math.h, libcudacxx/include/cuda/__fp/fpmp_math_impl.h
Math wrapper declarations and polynomial helpers switch to _CCCL_TRIVIAL_HOST_DEVICE_API. Their signatures and forwarding behavior remain unchanged.
FP unit-test CI lanes ci/matrix.yaml
Adds four GPU and CPU target-project overrides that run tests or precompile tests for cuda/fp/units.
The temporary CI override restricts PR validation to four targeted jobs and must be cleared before merge. The arithmetic changes carry no remaining blocking concern in the supplied evidence, and the two earlier test-coverage concerns are resolved.
Security architecture risk: Low · up to 278c4
The arithmetic changes stay inside the floating-point library. They change how scalar operands are sent to existing kernels and add no new interfaces, dependencies, trust boundaries or stored state. The main risk is the CI change: it swaps the normal pull-request checks for a narrow set of floating-point unit-test jobs. Its own comment says to undo it before merging, but nothing found in the pipeline enforces that. If it is merged as is, later pull requests could merge after running only these floating-point jobs.
Retained concerns
Low · reliability · inferred: The PR fills in the override list in ci/matrix.yaml. The pull-request workflow always allows overrides, so a non-empty list replaces the whole normal set of PR jobs with cuda/fp unit-test runs. No enforcement step was found that stops this from being merged. If it lands on main, every later PR would be checked with this narrow set instead of the full matrix, until someone resets it.
Security review details
Security Blast Radius
inferred — The worst case for the arithmetic changes is a numerical-accuracy regression in code that uses fpmp2 or the fp32mp2 math functions, not exposure of identities, tenants or data. If the CI override is merged, every later pull request into main would run the narrow cuda/fp job set instead of the full matrix.
Trust Boundaries and Controls
observed — The only thing keeping the override from reaching main is a "Reset before merging" comment. No workflow or action check was found that rejects a non-empty override at merge time.
Hardening Proposals
proposed — Add a required check that fails when the override list in ci/matrix.yaml is non-empty on the merge target, or that blocks merging while an override is active. That would make the existing "reset before merging" practice automatic.
- [ ] Autofix · Keep fixing CodeRabbit findings and required CI, and resolving merge conflicts
Add an integer case such as 16777217. With CCCL_FPMP_EXPLICIT_CASTS=0, also test a double that does not fit in a float. Compare limbs after each operation. The current test cannot detect narrowing before stat_t forwards the scalar.
Included review availability: This review used your included allowance. Your plan provides up to 12 included reviews per hour; 11 remain after this review.
…accepts
fpmp2_stat's +=, -=, *= and /= took a plain _FpType, so a scalar that one
limb cannot hold was narrowed first: fp32mp2_stat *= 16777217 multiplied by
16777216, where fp32mp2 converts the integer to an exact pair. With explicit
casts, a double also reached the fpmp2_stat overload through the implicit
single-limb constructor, which fp32mp2 rejects.
The four operators now take any arithmetic scalar the wrapped type's own
compound operator accepts and forward it unconverted, and are deleted for the
others. stat.pass.cpp checks the integer case and that both types accept the
same scalars; stat_implicit_casts.pass.cpp checks a double under
CCCL_FPMP_EXPLICIT_CASTS=0.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…function templates
Use _FpType{1} and _FpType{0} for the limb constants in the increment and
decrement operators (fpmp2 and fpmp2_stat) and in the single-limb division,
drop the redundant inline from the new __fpmp2_*_scalar function templates,
and mark the element-type aliases in the accuracy-selecting mixed functions
_CCCL_NODEBUG.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
suggestion: Check scalar multiplication and division counts.
The parity checks cannot detect a missing trace. Add device-record assertions for scalar assignments such as s *= int and s /= float. The existing increment and decrement checks already cover scalar += float and -= float indirectly.
suggestion: Test scalar-first subtraction against independent values.
The current assertion shares negation and scalar addition with sub, so it cannot detect an error in those operations. Add finite test cases with independently calculated results, including a non-normalized low pair.
Included review availability: This review used your included allowance. Your plan provides up to 12 included reviews per hour; 10 remain after this review.
…duct more accurate
The two-fma default product saved two instructions but made the residual of
x_hi * y_hi feed two dependent fmas. In a chain of products on the CPU that
made it 9% slower than the old code (whose p1 + p2 GCC fuses). The default
product now sums the cross terms with one fma off the residual's path, so
x_lo, which arrives last, enters only that fma: 8 instructions instead of 9,
1.19x lower latency on the host than before, 1.10-1.13x GPU throughput, and
the same measured error as Algorithm 10 (5.0u^2, 4.0u^2 for squares).
The high product (low product plus renormalization) turned out to reach
5.0u^2 on float squares, no better than mid. It is now x_lo * y_lo and the
cross terms, smallest first, summed by two fmas off the residual's path:
9 instructions, measured 4.0u^2 (3.5u^2 for squares).
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…vice
In the high x * c (Algorithm 7) the separately rounded x_lo * c fed an add,
which GCC contracts on the host by default: host and device disagreed on 4%
of products. Rounding it through an fma with a zero addend gives the same
value and leaves nothing to contract; on the device it is still one FFMA.
Every product now agrees bit for bit between host and device on 2^20 random
operands, where the old default product disagreed on 20%.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
cd5c542: the two-FMA mid product saved two instructions, but it lengthened the dependent path, and measured 9% slower than main on the CPU in chains of products. mid now sums the cross terms with one FMA off the residual's path, so x_lo enters last. That gives 8 instead of 9 instructions, 1.19× lower latency on the host, and 1.10–1.16× GPU throughput, with Algorithm 10's measured error. The earlier high product (low plus renormalization) turned out to reach 5.0u² on float squares, no better than mid, so the "3.9u²" in the first description was wrong. high is now a 9-instruction product measuring 4.0u² (3.5u² for squares).
f88376b: the highx * c rounded x_lo * c with a plain multiply that GCC contracts on the host, so host and device disagreed on 4% of products. With the fix, every product agrees bit for bit on 2²⁰ random operands. The old default product disagreed on 20%.
98047c8 / e5f7b61: the review fixes for fpmp2_stat and the style suggestions.
…ently
The fpmp2_stat counting kernel now also runs s *= 2 (a scalar that takes the
pair path), s /= 2.0f and s * 0.5f, and checks that each is traced once.
scalar_ops checks x - c and c - x against an exact expansion instead of
against the addition they share code with, and checks c - x on fixed
operands with hand-computed results, including a pair whose low limb is not
normalized.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
I addressed the two nitpicks from CodeRabbit's second review in 319912f:
Trace counts for scalar *= and /= (stat.pass.cpp): the counting kernel now also runs s *= 2, s /= 2.0f and s * 0.5f, and the expected mul and div counts rise to 4 and 3. The int takes the wrapped type's pair path, and the two float scalars take the single-limb path. A missing trace on any of them now fails the count assertions.
Scalar-first subtraction checked independently (scalar_ops.pass.cpp): x - c and c - x are now checked against an exact expansion of error-free transforms, the same reference as x + c, rather than against the addition they share code with. The test also checks c - x on fixed operands with hand-computed results:
3 - (1, u/2) gives the limbs (2, -u/2), which the normalizing levels must return exactly;
3 - (1, 4u), whose low limb is not normalized, must equal 2 - 4u exactly.
Both tests pass on host and device on sm_75, sm_89 and sm_120.
…er promotes
The mixed operator+ now adds the float product with the single-limb kernel
instead of converting it to a pair first.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…cannot contract it
__fpmp2_low_div and the low branch of __fpmp2_div_limb returned their lo word as a
plain product. GCC, under its default -ffp-contract=fast on FMA hardware, fuses such a
product into whatever adds the caller does next, which gives different bits from the
device, where __fmul_rn is never contracted.
Add __fpmp_mul_rn_nocontract, x * y as an fma with a zero addend, and use it for those
lo words and for the x_lo * c term of the high scalar product, which spelled the same
construct out by hand.
Device results are unchanged, except that an exactly zero lo word is now +0. Over the
fp32mp2 math functions (2^18 inputs each), results are bit-identical on sm_89 and on
the host with -ffp-contract=off; with GCC 15.2 -march=x86-64-v3 and default
contraction, cosh's results that differ from the device drop from 21135 to 4344.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… fpmp2_stat
- fptool_stat.h: mark the four deleted scalar compound operators _CCCL_HOST_DEVICE_API,
as every other deleted member is. Under clang-CUDA a host-only deleted overload is not
viable in device code, so the call would fall through to operator+=(const fpmp2_stat&)
via the narrowing single-limb constructor that the deletion exists to block. nvcc and
host compilers are unaffected.
- fpmp.h: the mixed operators and the mixed add/sub/mul/div<Acc> convert a scalar
dividend to a pair; say so in both comments. Rewrap the compound-assignment comment.
- fpmp_impl_muladd.h: the default product's dependent path ends in a Fast2Sum, not a
TwoSum, and in the scalar product |x_lo * c + residual| is below 4u, not 3u, when
x_hi * c lies in [2, 4); the 2u^2 rounding bound is unchanged.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…mb is subnormal
dd2262e made the low quotient's lo word a separate product so that host compilers
cannot fuse it into the next add. In the 0.5/e of fp32mp2 sinh and cosh, that product
is subnormal once e > 2^101 (|x| > 70), and Intel cores take a microcode assist for
each product with a subnormal result; fused into the add, it was normal. On an Intel
Core i5-11500T, where 18% of the benchmark's cosh inputs have |x| > 70, cosh went
from 44.5 to 51.8 ns per call. With FTZ/DAZ set, both versions run alike, and an AMD
Ryzen 5 3500U shows no difference.
For e > 2^32, 0.5/e is under 2^-63 of either result. The new
__internal_fpmp2_half_inv_exp returns 0 there, after dividing by 2^32 instead of e so
that no subnormal is formed, without a branch. On the i5-11500T, cosh now takes 39.6
ns and sinh 82.6 ns, against 51.8 and 94.9 before and 44.5 and 91.6 on main; no other
function changes.
Over 2^18 inputs each, sinh and cosh are bit-identical to before on sm_89 and on the
host with GCC 15.2, with and without -ffp-contract=off. Results can differ only where
the dropped term would have reached the low limb, by under 2^-63 relative. fp64mp2
sinh and cosh are binary128 wrappers and are unchanged.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Twenty expression sites in exp, the hyperbolic functions, the special functions and
asin/acos promoted a float that one limb holds exactly to a pair first (__ffloat(0.5f),
__ffloat(1.0f), __ffloat(__y0), ...) or doubled a pair as x + x. They now pass the scalar
to the mixed-type operators, or multiply by 2.0f. Today the operators promote the scalar
themselves, so values and code are unchanged (verified: identical maximum errors on the
accuracy suite); once NVIDIA/cccl PR NVIDIA#12519 lands, the same sites take its single-limb
kernels (6 operations and a 2u^2 bound for the product), measured at 1.053x throughput
over the fp32mp2 math set. Coefficient tables and pairs built from two limbs are unchanged.
(cherry picked from commit 47c0752)
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…-picked there) into dev/akolesov/fp32mp2_math_update_scalar
The rewrite now lives on math_update itself; this branch keeps the PR NVIDIA#12519 merge on top.
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
# Conflicts:
# libcudacxx/include/cuda/__fp/fpmp_math_impl_hyperbolic.h
__fpmp2_add_limb, __fpmp2_mul_limb, __fpmp2_high_mul, __fpmp2_div_limb and
__fpmp_mul_rn_nocontract were declared _CCCL_TRIVIAL_HOST_DEVICE_API, which also
marks them nodebug and artificial. They are now
_CCCL_HOST_DEVICE_API _CCCL_FORCEINLINE, like the operators that call them, so
they stay force-inlined and a debugger can step into them.
Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
…e CPU (#12543)
* [fpmp] Force-inline the fpmp2 arithmetic operators and the polynomial evaluators
The fpmp2 arithmetic and comparison operators, add/sub/mul/div, the __fpmp2_*_scalar
helpers and the three polynomial evaluators only dispatch to kernels that are already
force-inlined, but they were plain inline functions. GCC 14 for aarch64 calls several
of them out of line in the fp32mp2 math functions: the low-accuracy pair product and
quotient, the pair-by-float product and three __fpmp_poly_horner_mixed instances.
tanh, atanh, asinh, acosh and pow pay for those calls, and tanh became 0.72x on a
Raspberry Pi 5 when the single-limb kernels changed which operator it calls.
Declare those 44 functions in fpmp.h and the evaluators in fpmp_math_impl.h
_CCCL_TRIVIAL_HOST_DEVICE_API. The constructors, conversion and assignment operators,
volatile members and the named math functions (sqrt, fma, ...) are left as they are.
Timed with each function behind one out-of-line call, so that whether the compiler
inlines the math function into the loop does not enter, the geometric mean over 30
fp32mp2 math functions is 1.092x faster on a Raspberry Pi 5 (GCC 14; tanh 1.41x,
atanh 1.65x, pow 1.76x), 1.044x on an Intel Core i5-11500T (GCC 15.2) and 1.071x on
an Apple M4 (Apple clang 21), and no function is slower than 0.98x. When the math
function may be inlined into the loop, the larger sin and cosh no longer are, so the
plain benchmark shows sin 0.85x and cosh 0.84x on the Pi; its geometric means are
1.078x, 1.027x and 1.060x.
The device SASS is unchanged. Over 2^18 inputs per function, results are
bit-identical on sm_89; on the hosts (x86-64 with GCC 15.2, with and without
-ffp-contract=off, and the Pi) only the sign of a NaN result differs, for a few
special inputs.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] Force-inline the fpmp2 math function wrappers too
With the operators force-inlined, the bodies of the public math functions look larger
to GCC's inliner, and it stops inlining exp, log, log2, atan and cosh into their
callers. A loop such as `for (...) y[i] = cuda::log(x[i]);`, which GCC inlined and
vectorized on main, then calls the function for every element.
Their implementations are already force-inlined (_CCCL_FPMP_CORE_API); the 90 public
wrappers in fpmp_math.h and the 13 cuda::std forwarders were the only levels left to
the compiler. Declare them _CCCL_TRIVIAL_HOST_DEVICE_API as well.
ns per element of y[i] = f(x[i]) over 4096 fp32mp2 values, measured when this change sat on
top of #12519 and #12530: main / #12530 / #12530 with the operators force-inlined / with
this commit too:
AMD Ryzen Threadripper PRO 7965WX, GCC 15.2 -O3 -march=native
exp 1.93 / 1.84 / 23.6 / 1.82 atan 3.77 / 3.77 / 35.7 / 3.81
log 2.19 / 2.08 / 22.3 / 2.08 cosh 2.64 / 2.50 / 29.2 / 2.51
log2 27.2 / 2.25 / 24.4 / 2.23
Cortex-X925 (DGX Spark), GCC 13.3 -O3 -mcpu=neoverse-v2
exp 5.87 / 4.97 / 17.8 / 5.03 cosh 7.87 / 7.66 / 22.6 / 7.28
log 17.8 / 16.4 / 16.3 / 5.31 log2 19.8 / 18.4 / 18.5 / 6.69
sin, tanh, expm1, atanh and pow are equal or faster in the same loops. Called through
cuda::std, the wrappers alone are not enough (exp 22.2 ns, atan 33.6 ns on the
Threadripper); with the forwarders forced too, exp takes 1.83 ns and atan 3.77 ns.
Behind one out-of-line call, the geometric mean over 30 functions became 1.139x faster than
#12530 instead of 1.124x on the Threadripper, and 2.22x instead of 2.18x on the X925.
The device SASS is unchanged (sm_89, sm_90, sm_100, sm_120). On the host, results keep
their values; on x86-64 with GCC 15.2, one special-input result of erfc changes its
bits.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] CI override: run the cuda/fp unit tests on sm_75 and sm_90, precompile on arm64 [reset before merging]
This PR changes only fpmp headers, 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 a precompile of
the same suite with MSVC and on an arm64 runner with gcc, the compiler this change is
for. The override must be reset to empty before the PR can merge.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] Reset the CI override
Empties the override matrix again, so this PR runs the full CI before merging.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] Force-inline with _CCCL_FORCEINLINE so debuggers can still step in
_CCCL_TRIVIAL_HOST_DEVICE_API also marks a function nodebug and artificial, so a
debugger cannot step into it. The operators, polynomial evaluators and math functions
that this PR force-inlines are user-facing, so declare them
`_CCCL_HOST_DEVICE_API _CCCL_FORCEINLINE` instead, as libcudacxx does elsewhere
(cuda/__stream/invalid_stream.h). The fpmp kernels behind them keep
_CCCL_FPMP_CORE_API.
The generated host code (GCC 15.2) and the sm_89 SASS are unchanged.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] Review: say why the operators and math functions are force-inlined
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---------
Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
Twenty expression sites in exp, the hyperbolic functions, the special functions and
asin/acos promoted a float that one limb holds exactly to a pair first (__ffloat(0.5f),
__ffloat(1.0f), __ffloat(__y0), ...) or doubled a pair as x + x. They now pass the scalar
to the mixed-type operators, or multiply by 2.0f. Today the operators promote the scalar
themselves, so values and code are unchanged (verified: identical maximum errors on the
accuracy suite); once NVIDIA/cccl PR NVIDIA#12519 lands, the same sites take its single-limb
kernels (6 operations and a 2u^2 bound for the product), measured at 1.053x throughput
over the fp32mp2 math set. Coefficient tables and pairs built from two limbs are unchanged.
(cherry picked from commit 47c0752)
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
NVIDIA#12543 and NVIDIA#12519 have landed and this branch now sits directly on main, so the
override that limited CI to the cuda/fp unit tests is emptied and the full PR matrix
runs.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The single-limb `(s, e) - pl` and `+ x_lo` of the Cody-Waite reduction return the same
values as the pair adds they replaced, in fewer operations, and make sin, cos, sincos
and tan 5-10% faster on sm_100 and sm_120. On sm_90, nvcc schedules them worse, with
the same register use and occupancy: on an H100, sin and cos ran at 0.96x and 0.97x of
#12519, which is 0.97x and 0.98x of main.
sm_90 now keeps the pair adds (NV_IF_ELSE_TARGET(NV_IS_EXACTLY_SM_90, ...)); every other
target, the host included, compiles the same code as before. On an H100 (8 chains per
thread, median of 15 interleaved launches), over #12519:
sin 0.963 -> 0.996 sincos 0.989 -> 1.011
cos 0.966 -> 1.001 tan 1.043 -> 1.003
Against main, sin and cos are now 1.006x and 1.013x, and tan stays 1.029x faster. An
H200 measures the same. Results are bit-identical to before on sm_89, sm_90 and sm_120
over 2^18 inputs per function.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… π/2 and erfc at large x (#12530)
* [fpmp] Use the single-limb kernels in the fp32mp2 math functions
Where a math kernel promoted a scalar that one limb holds exactly to a pair before an
add, subtract or multiply (__ffloat(c), __afloat(c), __ffloat(c, 0.0f), a double literal,
or a pair built from a float variable), pass the float scalar instead, so the operation
takes the single-limb kernel. Doublings x + x become x * 2.0f; cbrt and rcbrt square
their float seed with one TwoProd; erf forms 1 - 2^n with one Fast2Sum.
The products that now reach a single-limb add directly, n * C3 in the sin/cos/tan and
exp10 reductions and f * pi2, f * pi3 in sincospi, go through __fpmp_mul_rn_nocontract,
so a host compiler that contracts gives the same bits as before.
Results are bit-identical on the device (31 fp32mp2 functions, 2^18 inputs each, sm_89)
and on the host with -ffp-contract=off, and unchanged on the host under GCC's default
contraction. FP operations saved per call, estimated from the kernels' operation counts:
30 in the sin/cos/sincos/tan Cody-Waite reduction, 30 in exp10's reduction, 20 in
sincospi, 17 in atanh, 15 in erfc, about 13 in cbrt and 8 in rcbrt, 10 in log1p, expm1
and boys_f0, 9 in erf, 7 in log and tanh, and 1 to 6 in exp, exp2, asin, acos, sinh,
cosh and normcdfinv.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] Extend the trig Cody-Waite reduction to a four-piece pi/2 with an exact n*C3
Near a multiple of pi/2 the reduced argument r is far smaller than n*C3, and two errors
that are negligible elsewhere dominate it: rounding n*C3 to one float, up to
2^-24 * n*C3, and the 9.5e-23 by which C1 + C2 + C3 misses pi/2, which adds
n * 9.5e-23. For cos(28781.701171875 - 7.6e-5), with n = 18323, each contributes about
1.7e-18 to r, a relative error of 3.8e-7 in all.
Compute n*C3 exactly with two_mult_fma, add n*C4 (C4 = -0x1.cceba4p-74; the four pieces
miss pi/2 by 8.5e-32) to its low word with one fma, and subtract the pair. Relative
error against binary128 over 2^18 inputs per function, 30% of them near a multiple of
pi/2, on sm_89 (the host gives the same figures), in u^2 = 2^-48:
function median 99.9th percentile max
sin 0.149 -> 0.115 8.3e8 -> 52 2.3e11 -> 9.3e3
cos 0.100 -> 0.077 8.5e8 -> 56 1.7e11 -> 8.5e3
tan 0.608 -> 0.31 1.6e9 -> 64 2.8e11 -> 6.7e3
The last correction costs 12 more FP operations than the single-limb version; the
reduction path as a whole is 73 operations, against 91 before the single-limb kernels.
sincospi keeps the three-piece pi, since pi * f does not cancel.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] Correct erfc's third ln 2 term
erfc reduces -x^2 by i * ln 2 with ln 2 as the pair __ln2_hi plus
__LN2_LO = ln 2 - 6.9314718055994529e-1. The pair holds that double only to 48 bits and
misses ln 2 by +8.78e-17, so adding __LN2_LO (+2.32e-17) left the term's error at
1.11e-16 * |i| rather than reducing it. Use ln 2 minus the pair, -0x1.950d88p-54 (within
3e-24), as a float, multiplied by i once and subtracted with the single-limb kernel:
8 FP operations instead of 11.
The accuracy gain is below noise. erfc's error is dominated by the pair rounding of -x^2
and by the low-accuracy reduction at magnitude x^2, both about u^2 * x^2. Against
binary128 over 2^18 inputs (sm_89), the median relative error goes from 5.54 to 5.53 u^2,
the 99.9th percentile from 80 to 79.6 u^2, and the maximum stays at 212 u^2.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] Reduce erfc's exp(-x^2) from an exact x^2
erfc(x) = erfcx(a) * exp(-a^2) with a = |x|. The exponent was rounded to a pair,
-a^2 ~ xx, and reduced by i * ln 2 in low-accuracy pair arithmetic at magnitude a^2;
each step erred by about u^2 * a^2 in the exponent, which is the relative error of the
result. The correction meant to restore -a^2 - xx recomputed -a * a as the same pair,
so it was always zero.
Take a^2 exactly instead: a_hi^2 and 2 a_hi a_lo by two_mult_fma, plus a_lo^2. Split
ln 2 as L1 + L2 + L3 with L1 to 12 bits, so that i * L1 is exact for |i| <= 1091 and
-a_hi^2 - i * L1 is exact by Sterbenz's lemma, sum the remaining small terms as a pair,
and add the leading term last. The reduced argument then errs by about u^2, and the
always-zero correction is dropped, so the exponent path also takes fewer operations.
Relative error against binary128 over 2^18 inputs (sm_89), in u^2 = 2^-48: median
5.53 -> 1.33, 99.9th percentile 79.6 -> 57.9, maximum 212 -> 71.1. For x >= 2 the
largest error in each interval of width 0.5 drops from 17..212, growing with x^2, to
5..8. The remaining maximum, about 70, is at |x| < 0.5 and comes from the erfcx kernel;
it is unchanged.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] Correct boys_f0's expansion point in its comments
The 4 <= a < 11.46 branch evaluates its minimax polynomial in 3.46 - a (0x1.baf1a8p1);
its leading coefficient is F0(3.4605). The comments said 6.92 - a.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] Make the polynomial evaluators round alike on the host and the device
__fpmp_poly_horner_mixed's float phase and __fpmp_poly_horner_comp's Phase 0 computed
v * x_hi + c as a plain expression, which nvcc contracts into an FMA by default and a
host compiler may or may not; they now call __fpmp_fma_rn. horner_mixed's transition
step added the plain product v * x_hi to a pair through the single-limb kernel, and GCC
under its default -ffp-contract=fast fused it into that kernel's TwoSum; the product
now goes through __fpmp_mul_rn_nocontract.
Device results are bit-identical (31 fp32mp2 math functions, 2^18 inputs each, sm_89),
and no function differs from the device on the host more often than before. With GCC
15.2 -march=x86-64-v3 and default contraction, exp, exp2, exp10, expm1 and erf now agree
with the device bit for bit (exp10 differed in 25991 of 262144 results), erfc's
differences drop from 46035 to 6996, and normcdfinv's from 76727 to 21881. With
-ffp-contract=off, exp, exp2, exp10 and erf agree bit for bit, and erfc's differences
drop from 11947 to 4863.
Also correct horner_mixed's comments: its float phase takes M - 1 steps, not M.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] Test the trig and erfc reductions where the reduced argument cancels
New math_reduction.pass.cpp checks fp32mp2 sin, cos, sincos and tan near a multiple of
pi/2, and erfc at large x, against binary128 references held as double pairs, to within
32u^2 (u = 2^-24), on the host and the device.
The inputs come from a sweep of 2^18 arguments per function. On them the earlier
reductions erred by 1.8e9 to 2.1e10 u^2 (trig) and 94 to 113 u^2 (erfc); the current
ones err by at most 6u^2 on sm_89 and on the host, with and without -ffp-contract=off.
The test aborts on the earlier headers. One extra cos case is the example from the
reduction commit, cos(28781.701171875 - 7.6e-5).
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] CI override: run the cuda/fp unit tests on sm_75 and sm_90 [reset 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>
* [fpmp] Keep the trig reduction's pair adds on sm_90
The single-limb `(s, e) - pl` and `+ x_lo` of the Cody-Waite reduction return the same
values as the pair adds they replaced, in fewer operations, and make sin, cos, sincos
and tan 5-10% faster on sm_100 and sm_120. On sm_90, nvcc schedules them worse, with
the same register use and occupancy: on an H100, sin and cos ran at 0.96x and 0.97x of
#12519, which is 0.97x and 0.98x of main.
sm_90 now keeps the pair adds (NV_IF_ELSE_TARGET(NV_IS_EXACTLY_SM_90, ...)); every other
target, the host included, compiles the same code as before. On an H100 (8 chains per
thread, median of 15 interleaved launches), over #12519:
sin 0.963 -> 0.996 sincos 0.989 -> 1.011
cos 0.966 -> 1.001 tan 1.043 -> 1.003
Against main, sin and cos are now 1.006x and 1.013x, and tan stays 1.029x faster. An
H200 measures the same. Results are bit-identical to before on sm_89, sm_90 and sm_120
over 2^18 inputs per function.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* [fpmp] CI override: name CTK 12.X's compilers
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>
* [fpmp] Reset the CI override before merging
#12543 and #12519 have landed and this branch now sits directly on main, so the
override that limited CI to the cuda/fp unit tests is emptied and the full PR matrix
runs.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---------
Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
Twenty expression sites in exp, the hyperbolic functions, the special functions and
asin/acos promoted a float that one limb holds exactly to a pair first (__ffloat(0.5f),
__ffloat(1.0f), __ffloat(__y0), ...) or doubled a pair as x + x. They now pass the scalar
to the mixed-type operators, or multiply by 2.0f. Today the operators promote the scalar
themselves, so values and code are unchanged (verified: identical maximum errors on the
accuracy suite); once NVIDIA/cccl PR NVIDIA#12519 lands, the same sites take its single-limb
kernels (6 operations and a 2u^2 bound for the product), measured at 1.053x throughput
over the fp32mp2 math set. Coefficient tables and pairs built from two limbs are unchanged.
(cherry picked from commit 47c0752)
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
Results, compared with
mainArithmetic between an
fpmp2and a scalar that one limb holds exactly no longer converts the scalar to a pair first. The pair product is also shorter, and itshighlevel more accurate. All figures are against currentmain(89e1909). This PR is stacked on #12543, which force-inlines thefpmp2operators and should merge first (see Merge order). u² is 2⁻⁴⁸ forfp32mp2and 2⁻¹⁰⁶ forfp64mp2.Performance.
x * cis 1.65–1.76× faster atlowand 1.42–1.51× atmid(the default).x + cis 1.9–2.0× faster athighand 1.1× atmid;x += cis 1.27–1.34× faster athigh,x / c1.04–1.19×, andx * y1.10–1.16× atmid.highand both in exchange for accuracy (below):x * cruns at 0.86–0.90× andx * yat 0.97–1.02×.fp32mp2x * ytakes 3.17 ns instead of 3.78 (1.19×), andx * c2.77 ns instead of 3.77 (1.36×). Athigh,x * cruns at 0.68–0.76× andx * yat 0.81–0.83×. On a DGX Spark's Cortex-X925 (GCC 13),x * ytakes 3.13 ns instead of 3.37 (1.07×) andx * c2.58 ns instead of 3.36 (1.30×); athighthey run at 0.87× and 0.73×.fp32mp2math functions, which pick up the kernels without changes of their own: the geometric mean over 30 functions is 1.028–1.061× faster, depending on the platform, GPUs from sm_75 to sm_121 included. The slowest function is 0.97× (asin on the Intel host), and the largest gains are 1.1–1.3× in asinh, acosh, atanh, pow, sinh and cosh.main, which is why [fpmp] Force-inline the fpmp2 operators for faster fp32mp2 math on the CPU #12543 lands first.Accuracy. Largest relative error,
main→ this PR,fp32mp2/fp64mp2:mainx * clow,midx * chighx * ymidx * yhighx + c,x - c,c - xandx / creturn the same values as onmain.fp32mp2math functions, only asin, acos and pow return different values on the device, and only in the last bits. Their largest error against binary128 is unchanged: 9.9, 9.1 and 329u².main'smidproduct disagreed on 20% of operands. Over the math functions (31 functions × 2¹⁸ inputs, GCC 15.2-O2 -march=x86-64-v3), host results that differ from the device fall from 534,891 to 494,113, with cosh going from 21,135 to 4,344 and pow from 19,653 to 4,675. With-ffp-contract=offthe count is unchanged (300,525 → 300,617).x * c,x * y,x + c,x / candx * x, at every level, for both types), onmainand here, so the figures above hold on all of them.main,x * yatmidandhighdiffered in 22% of results. Thefp32mp2lowquotient still differs in 13%, as onmain.Special values. No result moves away from IEEE 754.
min, ±1, ±max, infinities and NaN). It covers every scalar operator (x ± c,c - x,x * c,x / c,c / xand the compound forms), plusx * yandx / y, at every accuracy level, for both types, on the host and on the device.x * c(and ofx *= c) atmid, for both types, and each change is to the IEEE result. In those cells x is a positive operand of magnitude at mostminand c a negative one, and the product now underflows to −0 instead of +0.x * ydoes onmainand here.lowquotient, and NaN payloads.main's. The one exception, boys_f0(NaN), returns a subnormal instead of NaN on bothmainand this PR; only its bits differ. That value is undefined: a NaN reaches boys_f0's last branch with an uninitialized variable ([fpmp] Faster fp32mp2 math functions; accurate trig near multiples of π/2 and erfc at large x #12530).Implementation
cfit exactly in one limb: the limb type,floatforfp64mp2, and integers up to the limb's significand width. Such scalars now take dedicated kernels in the mixed operators, in the compound+=,-=,*=and/=, and inadd/sub/mul/div<Acc>. Any other scalar (adoublewithfp32mp2, say) is still converted to a pair.x + c,x - c,c - x: DWPlusFP, the existing accumulate. It gives the same value as the pair addition, in 10 instead of 20 instructions athighand 10 instead of 11 atmid.x * c: TwoProd, one FMA and Fast2Sum: 6 instead of 9 instructions atmid, 3 instead of 5 atlow. Withx_hiandcscaled into [1, 2), the FMA rounds a value below 4u, so the error is below 2u² to first order.highuses Algorithm 7 of Joldes, Muller & Popescu (2017), in 10 instructions.x / c: the existing division without the FMA on the divisor's zero low limb. It gives the same value in one instruction fewer.midsums the cross terms with one FMA (fma(x_lo, y_hi, x_hi * y_lo)), off the path of the exact residual ofx_hi * y_hi, then adds that sum to the residual. In a chain of productsx_loarrives last and enters only that FMA. The product takes 8 instructions instead of 9, with the measured error of Algorithm 10.highaddsx_lo * y_loand the cross terms, smallest first, through two FMAs off the residual's path, in 9 instructions.__fpmp2_high_accdrops a second Fast2Sum that never changed its already normalized input: 13 → 10 instructions, bit-identical (0 differences over 5·10⁷ operand pairs, cancellation included).__fpmp_mul_rn_nocontract, which computesx * yas an FMA with a +0 addend. Clang folds a −0 addend back into a multiply, so the addend has to be +0. Thelowquotients and thehighscalar product use it.__internal_fpmp2_half_inv_exptherefore returns 0 for that term once e = exp(|x|) exceeds 2³², where the term is under 2⁻⁶³ of the result. It also clamps the divisor there, so no subnormal is formed.mainover 2¹⁸ inputs each.fpmp2_statforwards its mixed and compound scalar operators to the wrapped type's. It accepts exactly the scalarsfpmp2accepts, and stays bit-identical to it._CCCL_HOST_DEVICE_API, like its other deleted members. Under clang-CUDA, a host-only deleted overload is not viable in device code, so a call would otherwise fall through to the narrowing single-limb constructor that the deletion exists to block.The new kernels are inline compositions of primitives and of entry points that exist in both the inline and the library mode, so
CCCL_FPMP_LIBneeds no new libcufp symbols. The changed pair products and__fpmp2_high_accreach library users when libcufp is rebuilt from these headers.Measurements
Accuracy method. Errors were measured on the host against an exact reference: the error accumulated from error-free transforms. The inputs were 5·10⁷ normalized operand pairs from random and structured regimes, plus a hill-climbing search for the worst case, for general operands and for squares. Squares have a maximum of 4.00u² at
midand 3.50u² athigh, against 4.00u² onmain. Math-function errors are against binary128 (libquadmath) over the 2¹⁸ inputs per function described in #12530.FP-pipe SASS instructions on sm_100,
main→ this PR:lowmidhighx * cx + c,x - cx += cx / cx * ymadGPU throughput, this PR over
main,fp32mp2/fp64mp2. Each thread runs 8 independent chains, with 8 blocks of 256 threads per SM; the figure is the best of two alternated runs, with CUDA 13.4. The operations that did not change (x += catlowandmid,x * yatlow) measured 0.99–1.05.x * clowx * cmidx * chighx + cmidx + chighx += chighx / cmidx * ymidx * yhighCPU latency on an AVX-512 x86-64 workstation CPU (AMD Ryzen Threadripper PRO 7965WX, one core, GCC 15.2
-O3 -march=native). The figures are the latency of one dependent product, as the Google Benchmark median of 10 repetitions; the coefficient of variation is below 1.7%.floatdoublex * yonmain(Algorithm 10 as GCC compiles it)x * yhere (mid)x * conmain(converted to a pair)x * chere (mid)CPU throughput, this PR over
main, for the operators in a loop of 1 dependent chain and of 8 or 16 chains, which GCC vectorizes into one or two dependent vector chains. The figures are forfp32mp2; thefp64mp2ratios are within 0.15 of them.x * clowx * cmidx * chighx + chighx * ymidx * yhighMore GPUs, this PR over
main,fp32mp2/fp64mp2, the same benchmark: H100 (sm_90), B200 (sm_100), B300 (sm_103) and Rubin (sm_107), built with CUDA 13.4 and GCC 13, and the GB10 GPU of a DGX Spark (sm_121), built with its CUDA 13.0. Clocks could not be locked; the figure is the best of three alternated runs. The unchanged operations measured 0.99–1.00, and an H200 and a second Rubin system give the same ratios within 0.01.x * clowx * cmidx * chighx + cmidx + chighx += chighx / cmidx * ymidx * yhighCPU latency on a DGX Spark: one dependent
fp32mp2operation on one Cortex-X925 core, GCC 13.3-O3 -mcpu=neoverse-v2(GCC 13 has no-mcpu=cortex-x925), best of three alternated runs;fp64mp2takes the same time. With 8 or 16 chains, GCC 13 vectorizesmain's and this PR's loops differently, and their ratios swing between 0.55× and 3.5× with the build options, so they are not shown.mainx * clowx * cmidx * chighx + cmidx + chighx += chighx / cmidx * ymidx * yhighBoth slowdowns are at
high, and both are accuracy trades.x * cthere is Algorithm 7, which halves the error of the pair product (2.97 → 1.50u²).x * ythere is the 9-instruction product, which takes the error from 5.0 to 4.0u².fp32mp2math functions, speedup overmain(old time / new time). The figure is the geometric mean over 30 functions, from the harness and platforms described in #12530; the "with #12530" column adds that PR, which is stacked on this one.__internal_fpmp2_half_inv_exp(item 5), which changes only sinh and cosh. [fpmp] Faster fp32mp2 math functions; accurate trig near multiples of π/2 and erfc at large x #12530 has the per-function table.Testing
cuda/fp/units/fpmp/scalar_ops.pass.cpp, run on the host and the device for all six types, checks:x + candx * cwithin 2u², andx * ywithin its first-order bound, of an exact expansion of error-free transforms;x + candx / cequal to the pair operation;main's product errs by 2.9u². These fail onmainand pass here.accuracy_mixed.pass.cpp: its stated contract now says thatmulwith a single-limb scalar is not bit-identical to the wrapped call. Its operands give exact products, so the checks themselves are unchanged.fptool/stat.pass.cppchecks that all compound forms matchfpmp2's, for a single-limb scalar and for one that needs both limbs (16777217), and that both types accept the same scalars. The newstat_implicit_casts.pass.cppchecks adoubleunderCCCL_FPMP_EXPLICIT_CASTS=0.main:fp32mp2math functions × 2¹⁸ inputs, on sm_89 and on the host (GCC 15.2,-march=x86-64-v3, with and without-ffp-contract=off).fp32mp2math functions × 2¹⁸ inputs are bit-identical on sm_89, sm_90, sm_100, sm_103, sm_107, sm_120 and sm_121, onmainand here (boys_f0(NaN) aside).cuda/fp/unitssuite passes on sm_89 (55 passed). Before that, it passed on sm_89 (55 passed;fpemu/float64andfpmp/atomic_ddunsupported; CUDA 13.4, GCC 15.2).fpemu/charbool.pass.cppfails the same way onmain.Documentation
Updated:
fpmp.rst;mul,accandmad.Not regenerated: the spec's measured accuracy and performance tables for multiplication (B300/B200). They need the spec harness. The PR's own benchmarks on H100, B200, B300, Rubin and DGX Spark are under Measurements.
Merge order
Stacked on #12543, which force-inlines the
fpmp2operators and math functions and should merge first. This branch merges it (278c4ae, 9bc91f2), so its diff includes #12543's commits until #12543 lands. With it, the new scalar operators of this PR are force-inlined like the others, with_CCCL_HOST_DEVICE_API _CCCL_FORCEINLINE(b8b4f8b). #12543's CI override is reset, so this PR runs the full CI.Checklist
🤖 Generated with Claude Code