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] Single-limb kernels for wide constants that fit one limb; 6-operation single-limb product at high - #12640
Arithmetic between an fpmp2 and a constant of a wider type whose value one limb holds, such as the 2.0 in y * 2.0 or the 2 in y + 2, now takes the single-limb kernels of #12519 instead of converting the constant to a pair. Code ported from double is full of such constants. x * 2.5 and x * 3 drop from 8 to 6 FP32 instructions (from 9 to 6 at high), and x + 2 at high from 20 to 10. The single-limb product at high now runs the 6-operation kernel that mid uses, so x * 2.5f at high drops from 10 to 6 as well. Scalars known only at run time keep their current path.
Expression (fp32mp2, mid unless noted)
main
This PR
x * 2.5, x * 3, x *= 3
8
6
x + 2.0, x + 2, 1 - x, x += 2
11
10
x + 2 at high
20
10
x * 2.5, x * 3.0 at high
9
6
x * 3 at high
9 (+3)
6
x * 2.5f at high (a float already)
10
6
x / 2.0, x / 3
12–13
11
fp64mp2: x * 3LL, x + 2LL
8, 11
6, 10
fp64mp2 at high: x * 2.5, x * 3LL
10, 9 (+3)
6, 6
FP-pipe SASS instructions on sm_89 (CUDA 13.4). On main the integer forms there also spend three conversion instructions, which are gone. sm_120 gives the same counts. With g++ 15.2 on x86-64 the new counts are the same, and main's are at most one lower. Constants that need both limbs (x * 0.1, x * 16777217) are unchanged.
The single-limb product at high. Since #12519 it ran Algorithm 7 of Joldes, Muller & Popescu (2017): 10 operations, one more than the pair product at high (9), so with this PR's dispatch x * 2.5 at high would have gone from 9 to 10. It now runs the kernel mid uses (one fma on the low terms and one renormalization, 6 operations; 966b254), and mid and high return the same bits.
Cost at high, time per product against pair × pair: 32% less on sm_89 and sm_120 (0.157 against 0.230 ps and 0.181 against 0.266 ps, 4M threads × 1024 products), 16% less on an x86-64 host (0.145 against 0.173 ns, g++ 15.2 -O3, AVX2, one core). With Algorithm 7 the same product took about 11% more than pair × pair on all three.
Accuracy at high: the largest error over 2·10⁸ random products is 1.91u², where Algorithm 7 gave 1.50u². The documented bound for x * c, 2u², is the same at every level and still holds; the pair product at high measures 4.0u².
A product of a high pair with a one-limb scalar therefore changes in its last bits against main. No math function uses that product: the code of all 88 is identical with and without this commit (below).
The change
A wide scalar is a double next to a float limb, or an integer type with more value bits than one significand (__fpmp2_is_wide_scalar_v).
It fits one limb if it is an integer with $|c| \le 2^{24}$ ($2^{53}$ for fp64mp2), or a double that is zero or a normal float: its low 29 significand bits are zero and $2^{-126} \le |c| < 2^{128}$ (__fpmp2_value_fits_limb, integer operations only).
The mixed operators +, -, * and /, in both operand orders, and add/sub/mul/div<Acc> check _CCCL_IS_CONSTANT_P(c) && __fpmp2_value_fits_limb(c) before converting c to a pair. Inside the force-inlined operators __builtin_constant_p is true for a literal, so the check folds away and only the chosen kernel is emitted. A scalar known only at run time never reaches the value test. A scalar dividend has no single-limb kernel and is converted as before.
New compound operators +=, -=, *= and /= take the wide scalars the pair already converts from implicitly (x += 2); these went through operator+=(const fpmp2&). x += 1.5 on an fp32mp2 is still rejected under the default CCCL_FPMP_EXPLICIT_CASTS=1.
Why only constants
Checking the value at run time would choose the same kernel on every compiler, but on the host it made per-element double scalars 1.6× slower when none fit and about 3× slower when half did. A run-time double in a pair expression already costs a full conversion, which the optimization guide says to hoist out of hot loops.
What can differ between builds
Where the compiler does not see the constant, the expression takes the pair path: in an unoptimized build, and with compilers that lack __builtin_constant_p (MSVC, and nvcc device code before CUDA 13.4).
Sums, differences and quotients return the same values on both paths. A product by a constant that is not a power of two can differ in its last bits. Both results are within their bounds, 2u² for the single-limb product and 6u² for the pair product (5u² at high).
Where the single-limb kernel is taken, x op 2.0 returns exactly what x op 2.0f returns. On a 15 × 15 grid of special operands (signed zeros, subnormals, ±min, ±1, ±max, infinities and NaN) this changes no high limb. (With Algorithm 7 at high, a tiny positive x times c = -FLT_MIN underflowed to +0 where IEEE 754 gives −0; the 6-operation kernel returns −0, as the pair path does.) The changes are to low limbs only (the sign of a zero, NaN payloads, and ±inf instead of NaN next to an infinite high limb in a low quotient), again matching the float constant.
Validation
cuda/fp/units, every fpmp and fptool test, compiled directly against the headers (CUDA 13.4, g++ 15.2), on this branch and on main:
host, g++ with the CI warning set, -Wconversion and -Werror, C++17 and C++20: 32 of 32 pass;
device, nvcc running the body on the host and in a kernel, sm_89, C++17 and C++20: 31 of 31 pass (atomic_dd needs sm_90);
8 tests that use mixed arithmetic (scalar_ops, stat, stat_implicit_casts, accuracy_mixed, api, implicit_conv, int, math) on sm_120 and sm_75, and through NVRTC on sm_89 (stat is unsupported there).
No test passes on main and fails here.
scalar_ops.pass.cpp now checks wide scalars that fit one limb against the limb-typed scalar, and ones that need both limbs against the pair, for every form and both operand orders. It also checks literal constants. Products accept either result, since a test cannot tell whether the compiler saw the constant.
scalar_ops passes with g++ 15.2 and clang 21 at -O0, -O1 and -O3, and with nvcc -G -O0 on sm_89, where constants are not folded.
Run-time scalars: on sm_120 the SASS of x * d and x + d with a double variable is identical to main's. On sm_89 it has one more integer instruction in the existing double to pair conversion. With g++ the instructions are reordered, and the int path has one more move.
After the kernel change (966b254): every cuda/fp/units test, built with the lit wrapper's force_include.h, passes on the host (g++ 15.2, -Wall -Wextra -Wconversion -Werror; 57 of 57) and on the host and in a kernel on sm_89 (nvcc 13.4; 55 of 55, atomic_dd and one test unsupported with nvcc aside), in C++17 and C++20.
Math functions: 88 functions (every public fpmp math function, sqrt, rsqrt, fma, mad) × the 8 pair types, compiled for the host and for sm_89, have the same instruction sequence on main, with this PR's dispatch, and with the kernel change. The same holds on top of [fpmp] Faster, platform-independent host seed for fp32mp2 cbrt and rcbrt #12568. So this PR changes code that mixes a pair with a constant, and nothing inside the library.
Not run: MSVC and lit itself.
Merge order
Branched from main at ae772fa, which contains #12519. It merges cleanly with current main (89b4901, which adds #12530), and every cuda/fp/units test passes on that merge: 58 on the host and 56 on the host and in a kernel on sm_89, in C++17 and C++20.
Checklist
New or existing tests cover these changes.
The documentation is up to date with these changes.
A constant of a type wider than one limb whose value one limb holds now takes the
single-limb kernels in x + c, x - c, c - x, x * c, x / c, the compound forms and
add/sub/mul/div<Acc>: a double that is zero or a normal float next to an fp32mp2, and
an integer of magnitude at most 2^24 (2^53 for fp64mp2). This covers the constants of
code ported from double, as in y * 2.0 and y + 2. Constants that need both limbs, and
scalars known only at run time, are converted to the pair as before.
- The choice is made where the compiler knows the value (_CCCL_IS_CONSTANT_P), so a
run-time scalar never runs the test. Where it does not know a constant (no
__builtin_constant_p, or an unoptimized build), a product by a constant that is not a
power of two can differ in its last bits; sums, differences and quotients are the same
either way.
- New compound operators take the wide scalars the pair converts from implicitly
(x += 2, x *= 3), which otherwise went through operator+=(const fpmp2&).
- scalar_ops.pass.cpp checks wide scalars that fit one limb and ones that need both,
for every form and both operand orders, and literal constants.
- fpmp.rst and the optimization guide describe the new cases.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…an the pair's
Review: x * 2.5 on an fp32mp2_high goes from 9 to 10 FP32 instructions, since it
now takes the single-limb product, the code x * 2.5f gets on main. fpmp.rst gave only
the product's saving at mid, and the optimization guide said a one-limb scalar is
always cheaper than a pair; at high the product costs 10 against the pair's 9 and
has the smaller error. An integer constant still comes out cheaper there, as the pair
path converts it first (12 against 10 on sm_89, CUDA 13.4).
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…high
The note added in 582b8d3 said an integer constant comes out cheaper in a product
at high. That holds for a constant of integer type (x * 3: 12 against 10 on sm_89,
CUDA 13.4), which the pair path converts with three instructions. An integer-valued
double such as the 3.0 in x * 3.0 costs 10 against 9, like x * 2.5.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The product of a pair with a scalar that one limb holds ran Algorithm 7 of Joldes,
Muller & Popescu 2017 at high: 10 operations, one more than the pair product there.
With this PR's dispatch, x * 2.5 at high would have gone from 9 to 10 instructions.
high now runs the kernel mid uses (6 operations), and the two return the same bits.
- x * 2.5, x * 3 and x * 2.5f at high take 6 FP instructions on sm_89 (CUDA 13.4),
where they took 10 and the pair product takes 9; fp64mp2 likewise.
- The largest error over 2e8 random products at high is 1.91u^2, where Algorithm 7
gave 1.50u^2. The documented bound, 2u^2, is the same at every level.
- No math function uses this product: their code is unchanged.
- The cost table, fpmp.rst and the error-bounds table follow; the two notes on the
one-instruction penalty at high are gone again.
Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
Mixed-precision arithmetic now uses the single-limb path for compile-time-known wider scalar values that fit in one limb. Other values, including runtime scalars, continue through pair conversion.
Compound assignments accept wider scalar types that can be converted to the pair type.
Documentation
Updated guidance on scalar handling, instruction counts, and accuracy bounds, including notes on possible last-bit differences when a constant’s value is not visible to the compiler.
The change routes compile-time-known wide scalar values that fit one limb through single-limb fpmp operations and adds wide-scalar compound assignments. It updates one-limb multiplication, adds scalar-operation tests, and revises the related documentation. Runtime scalars and values requiring both limbs retain pair conversion.
Changes: fpmp scalar arithmetic
Layer / File(s)
Summary
Wide-scalar classification and dispatch libcudacxx/include/cuda/__fp/fpmp_impl.h, libcudacxx/include/cuda/__fp/fpmp.h, libcudacxx/test/libcudacxx/cuda/fp/units/fpmp/scalar_ops.pass.cpp, docs/libcudacxx/extended_api/fp/fpmp.rst, docs/libcudacxx/extended_api/fp/fp_optimization_guide.rst
The implementation detects wide scalar constants that fit one limb and routes mixed arithmetic and accuracy-explicit operations through single-limb helpers. It adds compound assignments for convertible wide scalars. Tests cover fitting and promoted values, literals, and compound operators.
One-limb multiplication and documented bounds libcudacxx/include/cuda/__fp/fpmp_impl_muladd.h, docs/libcudacxx/extended_api/fp/fpmp_error_bounds.rst, docs/libcudacxx/extended_api/fp/fpmp.rst, docs/libcudacxx/extended_api/fp/fp_optimization_guide.rst
One-limb multiplication uses an FMA to calculate the low product limb and residual. The documentation updates scalar multiplication instruction counts and describes the bound across accuracy levels.
Suggested reviewers:miscco
Change: Feature
Merge Risk: Low · up to 966b2
Mixed arithmetic with constants that fit one limb now uses cheaper single-limb operations, and runtime scalars keep the old behavior. One documentation sentence wrongly calls the literal 1.5 a float. Fix that wording; the change is otherwise low risk to merge.
Security architecture risk: Minimal · up to 966b2
The inspected changes remain within numerical arithmetic and caller-owned values. No introduced or worsened risk to security boundaries, failure containment, rollback, or ownership guarantees was established. Multiplication results can change in their last bits, but no security-sensitive downstream consequence was identified.
Retained concerns
No architecture-level concerns identified.
Security review details
Trust Boundaries and Controls
inferred — The new dispatch predicate selects a numerical representation, not an identity or authority level. For wide operands not known to the compiler, the added branch does not replace pair conversion. These inspected paths do not establish expanded privileged reachability; security-sensitive downstream use remains outside this assessment.
Resilience and Maintainability Implications
inferred — The new compound operators compute an arithmetic result and assign it back to the same object, matching the existing pair compound operators. Their bodies introduce no reservation, persistent intermediate state, retry, cleanup, or ownership-transfer protocol. This does not establish atomicity or safety for concurrent unsynchronized access.
- [ ] Autofix · Keep fixing CodeRabbit findings and required CI, and resolving merge conflicts
Included review availability: This review used your included allowance. Your plan provides up to 12 included reviews per hour; 11 remain after this review.
The reason will be displayed to describe this comment to others. Learn more.
Minor · Functional Correctness · Quick win
suggestion: Correct the type of 1.5.
1.5 has type double, not float. (eel.is) Say “since the constant 1.5 is exactly representable in a float” to distinguish value-based dispatch from type-based dispatch.
This branch has not been deployed
No deployments
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}
Arithmetic between an
fpmp2and a constant of a wider type whose value one limb holds, such as the2.0iny * 2.0or the2iny + 2, now takes the single-limb kernels of #12519 instead of converting the constant to a pair. Code ported fromdoubleis full of such constants.x * 2.5andx * 3drop from 8 to 6 FP32 instructions (from 9 to 6 athigh), andx + 2athighfrom 20 to 10. The single-limb product athighnow runs the 6-operation kernel thatmiduses, sox * 2.5fathighdrops from 10 to 6 as well. Scalars known only at run time keep their current path.fp32mp2,midunless noted)mainx * 2.5,x * 3,x *= 3x + 2.0,x + 2,1 - x,x += 2x + 2athighx * 2.5,x * 3.0athighx * 3athighx * 2.5fathigh(afloatalready)x / 2.0,x / 3fp64mp2:x * 3LL,x + 2LLfp64mp2athigh:x * 2.5,x * 3LLFP-pipe SASS instructions on sm_89 (CUDA 13.4). On
mainthe integer forms there also spend three conversion instructions, which are gone. sm_120 gives the same counts. With g++ 15.2 on x86-64 the new counts are the same, andmain's are at most one lower. Constants that need both limbs (x * 0.1,x * 16777217) are unchanged.The single-limb product at
high. Since #12519 it ran Algorithm 7 of Joldes, Muller & Popescu (2017): 10 operations, one more than the pair product athigh(9), so with this PR's dispatchx * 2.5athighwould have gone from 9 to 10. It now runs the kernelmiduses (one fma on the low terms and one renormalization, 6 operations; 966b254), andmidandhighreturn the same bits.high, time per product against pair × pair: 32% less on sm_89 and sm_120 (0.157 against 0.230 ps and 0.181 against 0.266 ps, 4M threads × 1024 products), 16% less on an x86-64 host (0.145 against 0.173 ns, g++ 15.2-O3, AVX2, one core). With Algorithm 7 the same product took about 11% more than pair × pair on all three.high: the largest error over 2·10⁸ random products is 1.91u², where Algorithm 7 gave 1.50u². The documented bound forx * c, 2u², is the same at every level and still holds; the pair product athighmeasures 4.0u².highpair with a one-limb scalar therefore changes in its last bits againstmain. No math function uses that product: the code of all 88 is identical with and without this commit (below).The change
doublenext to afloatlimb, or an integer type with more value bits than one significand (__fpmp2_is_wide_scalar_v).fp64mp2), or adoublethat is zero or a normalfloat: its low 29 significand bits are zero and__fpmp2_value_fits_limb, integer operations only).+,-,*and/, in both operand orders, andadd/sub/mul/div<Acc>check_CCCL_IS_CONSTANT_P(c) && __fpmp2_value_fits_limb(c)before convertingcto a pair. Inside the force-inlined operators__builtin_constant_pis true for a literal, so the check folds away and only the chosen kernel is emitted. A scalar known only at run time never reaches the value test. A scalar dividend has no single-limb kernel and is converted as before.+=,-=,*=and/=take the wide scalars the pair already converts from implicitly (x += 2); these went throughoperator+=(const fpmp2&).x += 1.5on anfp32mp2is still rejected under the defaultCCCL_FPMP_EXPLICIT_CASTS=1.Why only constants
Checking the value at run time would choose the same kernel on every compiler, but on the host it made per-element
doublescalars 1.6× slower when none fit and about 3× slower when half did. A run-timedoublein a pair expression already costs a full conversion, which the optimization guide says to hoist out of hot loops.What can differ between builds
__builtin_constant_p(MSVC, and nvcc device code before CUDA 13.4).high).x op 2.0returns exactly whatx op 2.0freturns. On a 15 × 15 grid of special operands (signed zeros, subnormals, ±min, ±1, ±max, infinities and NaN) this changes no high limb. (With Algorithm 7 athigh, a tiny positivextimesc = -FLT_MINunderflowed to +0 where IEEE 754 gives −0; the 6-operation kernel returns −0, as the pair path does.) The changes are to low limbs only (the sign of a zero, NaN payloads, and ±inf instead of NaN next to an infinite high limb in alowquotient), again matching thefloatconstant.Validation
cuda/fp/units, everyfpmpandfptooltest, compiled directly against the headers (CUDA 13.4, g++ 15.2), on this branch and onmain:-Wconversionand-Werror, C++17 and C++20: 32 of 32 pass;atomic_ddneeds sm_90);scalar_ops,stat,stat_implicit_casts,accuracy_mixed,api,implicit_conv,int,math) on sm_120 and sm_75, and through NVRTC on sm_89 (statis unsupported there).mainand fails here.scalar_ops.pass.cppnow checks wide scalars that fit one limb against the limb-typed scalar, and ones that need both limbs against the pair, for every form and both operand orders. It also checks literal constants. Products accept either result, since a test cannot tell whether the compiler saw the constant.scalar_opspasses with g++ 15.2 and clang 21 at-O0,-O1and-O3, and with nvcc-G -O0on sm_89, where constants are not folded.Run-time scalars: on sm_120 the SASS of
x * dandx + dwith adoublevariable is identical tomain's. On sm_89 it has one more integer instruction in the existingdoubleto pair conversion. With g++ the instructions are reordered, and theintpath has one more move.After the kernel change (966b254): every
cuda/fp/unitstest, built with the lit wrapper'sforce_include.h, passes on the host (g++ 15.2,-Wall -Wextra -Wconversion -Werror; 57 of 57) and on the host and in a kernel on sm_89 (nvcc 13.4; 55 of 55,atomic_ddand one test unsupported with nvcc aside), in C++17 and C++20.Math functions: 88 functions (every public
fpmpmath function,sqrt,rsqrt,fma,mad) × the 8 pair types, compiled for the host and for sm_89, have the same instruction sequence onmain, with this PR's dispatch, and with the kernel change. The same holds on top of [fpmp] Faster, platform-independent host seed for fp32mp2 cbrt and rcbrt #12568. So this PR changes code that mixes a pair with a constant, and nothing inside the library.Not run: MSVC and lit itself.
Merge order
Branched from
mainat ae772fa, which contains #12519. It merges cleanly with currentmain(89b4901, which adds #12530), and everycuda/fp/unitstest passes on that merge: 58 on the host and 56 on the host and in a kernel on sm_89, in C++17 and C++20.Checklist
🤖 Generated with Claude Code