[akolesov-nvidia]

Summary

The fp32mp2 math functions are reimplemented in float-float arithmetic, with no FP64 operation on the fast paths. On the normal dataset they reach 45-47 correct bits (of 48) where the previous implementations delivered 41-45, and on the pattern dataset (powers of two, extreme exponents, special structure) the trigonometric functions go from 25-28 bits to 44-45. Functions that were double or binary128 fallbacks become native (erfcx, normcdf, erfinv, erfcinv, hypot, rhypot, 5-30x the throughput), the rounding and remainder families become exact on the whole pair, and the functions return the C99 Annex F values for infinities, NaN and arguments outside their domain.

Highlights. Geomean throughput gain over the 40 changed fp32mp2 cells: 1.58x. normcdf 30x, hypot 13.7x, rhypot 10.6x, erfcinv 10x, erfcx 7.5x, erfinv 4.9x, frexp 2.7x, rcbrt 1.4x, pow 1.4x, exp10 1.35x, sinh 1.2x, atan 1.16x; sin, cos and tan match main's throughput at 5 more bits (41 to 45-46), and the rest of the math is 0.95-1.2x at 1-3 more bits. normcdfinv (0.92x) carries its new domain checks. The exact fpmp2_accuracy::high level of fmod and remainder is new (on main that level ran the default code), so it has no comparison. sqrt gains 2 bits and is 1.13x faster; rsqrt is unchanged.

The table compares this branch with main (778da2d, with #12519, #12530 and #12543) on an RTX A6000 (sm_86, CUDA 13.4), measured with the cccl-fpdev suite at the default accuracy level. All rows are fp32mp2 at the default level. Columns: correct bits on the normal and the pattern dataset; throughput in evaluations per clock per SM (one run each on the same GPU in the same order; run-to-run noise on this setup measured at 3.6% or less); SASS instructions of the function's own code (cuobjdump --dump-sass of the bare function compiled with -O2 for sm_86, out-of-line helpers included, CUDA's division and square-root slow-path subroutines excluded); ptxas registers. A dash marks a level that had no implementation of its own on main.

function	level	bits old (normal / pattern)	bits new (normal / pattern)	ev/clk/SM old	ev/clk/SM new	gain	SASS old	SASS new	regs old / new
acos	def	44 / 44	46 / 46	0.46	0.48	1.04x	472	464	21 / 21
acosh	def	43 / 43	45 / 45	0.42	0.44	1.05x	544	504	16 / 17
asin	def	44 / 44	46 / 46	0.43	0.42	0.98x	440	448	21 / 21
asinh	def	43 / 43	45 / 45	0.40	0.43	1.07x	552	512	16 / 17
atan	def	45 / 45	47 / 47	0.37	0.43	1.16x	280	248	22 / 20
atan2	def	45 / 0	46 / 22	0.33	0.38	1.15x	360	312	25 / 22
atanh	def	44 / 44	46 / 46	0.53	0.55	1.04x	496	464	16 / 16
boys_f0	def	43 / 43	44 / 44	0.36	0.44	1.22x	752	824	21 / 21
cbrt	def	45 / 44	47 / 47	1.19	1.32	1.11x	104	96	21 / 20
cos	def	41 / 31	46 / 45	0.48	0.48	1.00x	448	592	25 / 34
cosh	def	43 / 40	46 / 46	0.63	0.64	1.02x	184	176	17 / 18
cospi	def	41 / 41	46 / 45	0.45	0.47	1.04x	320	304	23 / 21
erfc	def	41 / 23	45 / 23	0.18	0.20	1.11x	528	672	32 / 24
erfcinv	def	48 / 48	46 / 46	0.04	0.40	10.00x	408 (+3 slow-path calls)	2928 (+6 slow-path calls)	23 / 33
erfcx	def	48 / 26	45 / 23	0.04	0.30	7.50x	216 (+1 slow-path calls)	1320	18 / 28
erfinv	def	48 / 24	46 / 23	0.13	0.64	4.92x	232 (+1 slow-path calls)	2912 (+6 slow-path calls)	20 / 32
exp	def	43 / 25	46 / 33	0.73	0.82	1.12x	152	144	20 / 18
exp10	def	45 / 30	46 / 30	0.82	1.11	1.35x	200	144	14 / 16
exp2	def	44 / 25	46 / 33	0.75	0.92	1.23x	152	128	16 / 16
expm1	def	46 / 40	46 / 45	0.63	0.70	1.11x	312	272	17 / 20
fmod	def	48 / 48	48 / 48	2.92	2.92	1.00x	288 (+2 slow-path calls)	296 (+2 slow-path calls)	22 / 22
frexp	def	48 / 48	48 / 48	0.99	2.66	2.69x	56	48	10 / 9
hypot	def	48 / 29	46 / 29	0.12	1.64	13.67x	160	88	18 / 14
log	def	44 / 44	45 / 45	0.78	0.83	1.06x	152	152	16 / 18
log10	def	44 / 44	44 / 45	0.72	0.80	1.11x	160	152	16 / 16
log1p	def	43 / 43	45 / 45	0.61	0.66	1.08x	328	280	16 / 16
log2	def	44 / 44	45 / 44	0.73	0.82	1.12x	160	152	16 / 16
normcdf	def	48 / 33	45 / 23	0.02	0.60	30.00x	192	728	34 / 24
normcdfinv	def	41 / 41	45 / 45	0.25	0.23	0.92x	736	840	20 / 19
pow	def	41 / 28	46 / 36	0.32	0.44	1.38x	352	328	17 / 27
rcbrt	def	46 / 46	47 / 47	0.97	1.36	1.40x	136	104	18 / 15
remainder	def	48 / 48	48 / 48	2.25	2.25	1.00x	352 (+2 slow-path calls)	360 (+2 slow-path calls)	26 / 26
rhypot	def	48 / 23	44 / 22	0.15	1.59	10.60x	160	96	18 / 16
round	def	48 / 48	48 / 48	3.33	3.61	1.08x	112	88	10 / 9
sin	def	41 / 33	46 / 45	0.49	0.48	0.98x	440	584	25 / 34
sinh	def	43 / 40	46 / 46	0.63	0.76	1.21x	320	288	17 / 18
sinpi	def	41 / 23	46 / 25	0.43	0.45	1.05x	328	360	26 / 24
sqrt	def	44 / 25	46 / 25	5.22	5.92	1.13x	32	32	12 / 8
tan	def	41 / 31	45 / 44	0.44	0.43	0.98x	456	600	25 / 34
tanh	def	45 / 45	46 / 46	0.60	0.57	0.95x	336	328	17 / 18
fmod	high	-	48 / 48	-	3.10	-	-	1096 (+2 slow-path calls)	- / 32
remainder	high	-	48 / 48	-	1.98	-	-	2152 (+2 slow-path calls)	- / 33

Cells not listed are unchanged in bits and SASS size, with throughput within noise (8%). The pattern dataset is sign-sensitive in places where fpmp2 does not track the sign of zero (documented in this PR), which is what the low pattern bits of atan2, sinpi, erfc and the fallback-replacing functions reflect on both sides; their normal-dataset bits are the meaningful comparison. erfinv and erfcinv call a shared inverse-CDF kernel, which is why their SASS counts grow while the throughput rises 5-10x.

The high rows of fmod and remainder are a new level rather than a replacement: on main the high level ran the default code, so the old-side cells are left empty. On dense pairs the exact path costs remainder 1.98 against the default level's 2.25 evaluations per clock per SM (its wide-reduction block is present whether or not it is taken, 26 to 33 registers); fmod high runs at 3.10 against 2.92. What the level buys shows with sparse pairs in the dataset (one pair in 1024 with its tail 1-36 binades below the head):

cell	main 778da2d	this branch
fmod high, sparse pairs	0 bits, FAIL (the high level ran the default code)	48 bits
remainder high, sparse pairs	0 bits, FAIL	48 bits
fmod / remainder def, sparse pairs	48 bits dense, sparse operands reduced as their 53-bit rounding (by contract, a warning)	same

What changes, per area

• Exponentials and logarithms. exp to 46 bits with gradual underflow; exp2, exp10, expm1, log, log2, log1p to 45-46 bits (log10 stays at 44 on the normal dataset, 45 on the pattern dataset). The exp kernel is split so that erfcx and the other users scale the result themselves.
• Trigonometric and inverse trigonometric. sin, cos, tan, sincos, sinpi, cospi to 45-46 bits; asin, acos, atan, atan2 to 46-47 bits. The argument reduction keeps [fpmp] Faster fp32mp2 math functions; accurate trig near multiples of π/2 and erfc at large x #12530's fast-path accumulation in fp32mp2_high, which is what holds the error at a few u² within 2^-40 of a multiple of pi/2 (math_reduction.pass.cpp); the Payne-Hanek path and the kernels are new.
• Hyperbolic. sinh, tanh, atanh to 46 bits and asinh, acosh to 45, with the overflow and domain fixes; cosh uses the clamped 0.5/e helper from [fpmp] Faster, more accurate fpmp2 arithmetic with single-limb scalars and a shorter pair product #12519.
• Powers, roots, distances. pow is the triple-float implementation of [fpmp] Accurate, C99-conforming and faster fp32mp2 pow #12548 (a triple-float ln|a| and a table-driven exp), taken verbatim with its unit test so that the two branches do not diverge; cbrt, rcbrt to 47 bits; native fp32mp2 hypot (46 bits) and rhypot (44 bits, through rsqrt) instead of the FP64 wrappers.
• Error and probability functions. erfc and normcdfinv to 45 bits, boys_f0 to 44; native erfcx, normcdf, erfinv, erfcinv instead of the double and binary128 fallbacks (5-30x the throughput, at 45-46 bits where the fallbacks read 48).
• Rounding, remainder and conversions. ceil, floor, trunc, round, rint, nearbyint, fdim, remquo, lrint, lround, llrint, llround are exact on the pair (the previous round misrounded half-way pairs such as (9969168, 0.49999997)); fmod, remainder and remquo are exact on sparse pairs (a tail far below the head) in fpmp2_accuracy::high, and compute on the pair rounded to 53 bits in the other levels, as before; the default level is unchanged in code path, accuracy and throughput.
• Classification. isfinite, isinf, isnan, signbit and the fpmp_ spellings look at both limbs.
• sqrt. A new form for fp32mp2, 46 bits at 1.13x the throughput; rsqrt is unchanged. Both take no branch for special arguments, which is now documented.
• Special values. log, log2, log10, atan, the six rounding functions, normcdfinv and boys_f0 return the C99 Annex F values at zero, negative arguments, the infinities and NaN, through branch-free selects (cost measured at 0-7% on the affected functions). The sign of a zero result is not tracked by fpmp2, documented next to the sqrt/rsqrt note.
• Kernels take scalars as scalars (x * 2.0f rather than x * fp32mp2(2.0f)), so the single-limb operator overloads of [fpmp] Faster, more accurate fpmp2 arithmetic with single-limb scalars and a shorter pair product #12519 apply inside the math as well.

Tests

• cuda/fp/units/fpmp/math.pass.cpp: normcdfinv at the ends of (0, 1), the remainder family, more functions on both host and device.
• cuda/fp/units/fpmp/math_special_values.pass.cpp (new): class and sign of the result at zero, negative arguments, the largest finite values, the infinities and NaN for the functions above.
• cuda/fp/units/fpmp/math_pow.pass.cpp (new, from [fpmp] Accurate, C99-conforming and faster fp32mp2 pow #12548).
• cuda/fp/units/fpmp/math_reduction.pass.cpp (from [fpmp] Faster fp32mp2 math functions; accurate trig near multiples of π/2 and erfc at large x #12530, on main): passes against the new kernels, worst case 5.8 u² of the 32 u² bound.

Validation

With the cccl-fpdev suite against this branch: unit matrix (host, device, NVRTC, each at C++17 and C++20) all passing; the CI accuracy set (171 cells, pattern dataset, all levels) passes with 0 FAIL and 103 OK / 68 WARN cells, against 86 OK / 85 WARN on main (the warnings are the by-design special-input classes); the special-value probe (fp32mp2, 41 functions x 19 arguments, host and device) reports no failure, only the documented warnings (sign of zero, sqrt/rsqrt, subnormal results).

[copy-pr-bot]

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.