Fast Polynomial Transcendentals for LLMs

• The paper targets low-degree polynomial transcendentals in attention, dense activations, and routed experts.Source: REPORT FOCUS · p. 1
• Mixture of experts is also among the paper’s stated application areas.Source: REPORT FOCUS · pp. 1–2
• Its stated technical scope includes GPU kernels, polynomial approximation, transcendental functions, FlashAttention, and BF16.Source: REPORT FOCUS · p. 2

• Blackwell’s faster Tensor Core scaling relative to exponential and shared-memory resources exposed non-matrix work inside optimized attention [24].Source: §1 · p. 2
• The study traces polynomial candidates from target-precision error through component runtime and complete-step throughput to task and training behavior.Source: §1 · p. 2
• The four interventions use degree-3 or degree-4 polynomials, with B3 replacing sigmoid and its derivative but not softmax exponential.Source: Table 1 · p. 3
• “Native” means the unchanged implementation of the same transcendental for the same integration task.Source: §1 · p. 3
• The fitting workflow combines analytical range reduction, structural constraints, and replay of coefficient and intermediate rounding in deployed order.Source: §1 · p. 3
• Complete-step throughput improves 2.7%, 2.9%, and 8.0% for B1, B2, and B4; B3 improves complete-attention forward by 7.4%.The paper also releases fitting tools, packed kernels, integration and evaluation code, and supporting artifacts.Source: §1 · p. 3
• Near 100 billion tokens, final smoothed training-loss differences (polynomial minus native) are +0.079, +0.003, −0.107, and +0.003 for B1–B4.Source: §1 · p. 3

2.1 Asymmetric accelerator scaling

• Blackwell’s faster matrix-throughput growth made non-matrix work more visible in attention, while later vendor-reported hardware changes illustrate that bottlenecks can move across generations.All measurements in the paper use GB200.Source: §2.1 · p. 3

2.2 Target-precision numerical approximation

• Unlike RLIBM’s correctly rounded library target, this work targets lower-latency programs inside specific BF16/FP16 neural-network kernels.The paper also contrasts its fixed symmetry, saturation, and FMA-budget constraints with AutoNumerics-Zero’s search over limited-precision programs.Source: §2.2 · p. 4

2.3 Nonlinear operators in LLM components

• Prior work spans low-precision tanh-derived activations and CPU polynomial operators, whereas this work targets packed BF16 CUDA/CuTe programs in Blackwell attention and expert kernels.The cited operator settings include sigmoid self-attention and routed SwiGLU experts [19].Source: §2.2 · p. 4, §2.3 · p. 4

3 Target-Precision Polynomial Programs

• The optimization target includes range reduction, coefficients, instruction order, intermediate rounding, tail behavior, and the backward rule in BF16 or FP16.The method replays the exact target-precision schedule before kernel integration.Source: §3 · p. 4

3.1 Analytical reduction

• Analytical identities restrict fitting to a half-domain residual and recover known offsets, linear terms, or reflected values at runtime.For sigmoid and tanh, odd symmetry supports sign reconstruction; subtracting x/2 from SiLU exposes an even residual.Source: §3.1 · pp. 5–6, Figure 1 · p. 5

3.2 Constrained fitting and rounded-program replay

• The fitting procedure searches the polynomial interval Li separately from the runtime clamp Lc, while endpoint-vanishing basis terms preserve fixed values.Candidates are rounded to BF16 and scored after target-format Horner replay; final selection replays the deployed packed-FMA order.Source: §3.2 · pp. 6–7

3.3 Forward and backward programs

• Forward and backward programs use parity identities, while direct derivative fitting can control gradient error independently of the forward fit.B3 pairs a direct D3 sigmoid value with a separately targeted, factored D4 derivative surrogate.Source: §3.3 · p. 8

3.4 Packed CUDA and inline PTX programs

• Packed two-lane arithmetic is used for FP16 sweeps and BF16 integrations, with CUDA vectorization and fused SwiGLU choices reducing scalar work and intermediate materialization.The SwiGLU program combines activation and gating in forward and computes both projected-input gradients in one backward kernel.Source: §3.4 · p. 8
• B3 embeds a branch-free, two-lane BF16 sigmoid evaluator in the FA4 instruction schedule using a packed mask, immediate coefficients, and three packed FMAs.Its forward block has 17 PTX instructions; the coefficients are fitted for sequence length 4096 and pre-scaled by 1/n.Source: §3.4 · p. 9
• B3 backward reuses the forward conversion, clamp, mask, and Horner value to return a factored D4 surrogate with the D3 probability.The coupled block has 22 PTX instructions and four packed FMAs, versus 17 plus 22 instructions for separate forward and derivative blocks.Source: §3.4 · p. 9

4 Experimental Setup

• All performance measurements use NVIDIA GB200, with every operator and model setting outside the named substitution held fixed.Source: §4 · p. 9

4.1 Function error and coefficient control

• Target-precision plots replay BF16 coefficient and intermediate rounding, while a same-form Sollya minimax control is evaluated separately in real arithmetic.The control uses the same degree, sparse monomials, clamp, and reconstruction, over 20,001 points on [−Lc, Lc].Source: §4.1 · p. 9

4.2 Isolated function timing

• The isolated FP16 sweep compares packed polynomial kernels with native PyTorch at L2-resident and HBM-streaming sizes, excluding forward work from backward timings.Native sigmoid and tanh backward use precomputed outputs and algebraic rules, whereas polynomial comparators evaluate separately fitted D4 derivative programs.Source: §4.2 · pp. 9–10

4.3 Component and training-step timing

• B1–B3 full-model timings bracket forward, backward, and complete GPU optimizer steps; B4 runs the complete DeepSeek-V3 27A4B architecture, changing only 29 routed grouped-expert SwiGLU modules.B1–B3 use Llama 3 8B with 100 optimizer steps and medians over steps 20–100. B4 uses four GB200s with EP=4, sequence length 4096, local batch one, and global batch four.Source: §4.3 · p. 10

4.4 Same-checkpoint evaluation

• Same-checkpoint evaluation compares native, constrained-fit, and same-form Sollya-fit degree-4 substitutions across open-weight models using accuracy tasks and WikiText-2 perplexity.The reported task mean is the unweighted mean of seven accuracy tasks; throughput estimates are single-run point estimates.Source: §4.4 · p. 10

4.5 Paired pre-training

• Each production comparison pairs one native and one polynomial BF16 trajectory, with sampled training loss smoothed using a one-billion-token half-life and compared on a shared token grid.Runs use 32 GB200 GPUs and sequence length 4096; the logged training loss is the sole long-run loss metric.Source: §4.5 · p. 11

5.1 Function-level approximation

• The constrained D3 sigmoid fit has maximum absolute error 5.44 × 10−3 versus 5.47 × 10−3 for its Sollya control; Sollya has lower error in the other seven rows.Target-precision sweeps show that in each of six sweeps at least one higher-degree program has greater maximum error than a lower-degree program after BF16 replay; B3 pairs a D3 forward fit with a separately fitted D4 derivative.Source: §5.1 · p. 11, Table 2 · p. 11

5.1.1 Isolated FP16 function speed

• The cache-resident isolated FP16 endpoint ranges from 1.187x to 2.193x; at the HBM endpoint, sigmoid and tanh backward measure 1.004x and 1.003x, while GELU backward measures 1.615x.Seven of eight rows have smaller speedups at the HBM endpoint; GELU forward measures 1.673x in L2 and 1.702x in HBM.Source: §5.1.1 · p. 12, Table 3 · p. 12

5.2 Integration into LLM components

• B1’s unary D3 SiLU probe measures 1.055x forward and 1.007x backward; B2’s complete causal attention measures 2.053x forward and 1.046x backward.B1’s probe uses a direct D3 derivative fit, while deployed B1 recomputes a D3 sigmoid approximation and forms the SiLU derivative algebraically; B2 differentiates its selected D4 tanh coefficients for backward.Source: §5.2.1 · p. 13, §5.2.2 · p. 13

5.2.3 B3: sigmoid attention

• At sequence length 4096, B3 uses a D3 sigmoid forward and coupled, factored D4 derivative surrogate; complete-attention forward measures 1.074x and backward 1.001x.A runtime guard enforces the fitted sequence length.Source: §5.2.3 · p. 13

5.2.4 B4: routed expert SwiGLU

• B4 changes only the routed-expert SwiGLU activation, leaving the router native; its one-kernel activation benchmark measures 1.003x forward and 1.005x backward.Backward recomputes the D3 sigmoid approximation, forms SiLU and its derivative algebraically, and computes gradients for both projected inputs in one pass. An auxiliary fractional-exponential simplification measures 1.898x for its isolated core and 1.0008x at the model step.Source: §5.2.4 · p. 13

5.3.1 Complete-model timing

• B1 and B2 retain complete-step gains in the Llama 3 probes, B3’s attention-forward change is smaller at the complete-model step, and B4 improves both full-model phases.Table 5 reports B1–B3 on matched 100-step Llama 3 8B probes on one GB200 and B4 on the complete DeepSeek-V3 27A4B architecture on four GB200s.Source: §5.3.1 · p. 14, Table 5 · p. 14

5.3.2 Same-checkpoint open-weight evaluation

• Qwen3 and Kimi prefill improves by 5.1–5.6% and 8.4%, respectively, while Qwen2.5 prefill improves by 3.5–3.6% as its word perplexity rises from 8.734 to 197.677 and 298,700.735.GLM prefill ratios are 0.969x and 0.966x; GPT-OSS ratios are 1.001x and 0.999x. Table 6 compares degree-4 substitutions with same-form Sollya controls.Source: §5.3.2 · p. 14, Table 6 · p. 14

5.3.3 Paired pre-training to 100 billion tokens

• At common horizons near 100 billion tokens, B2 and B4 finish within 0.003 final EWM training-loss units of native, while B1 is 0.0786 higher and B3 is 0.1073 lower.Figure 3 shows paired trajectories; Table 7 reports the final EWM at each pair’s common observed horizon.Source: Table 7 · p. 15, §5.3.3 · p. 16, Figure 3 · p. 16

6 Discussion and Conclusion

• Phase-level gains differ by integration site: dense SiLU gains more in backward, tanh softcap in forward, and routed-expert SwiGLU in both.Sigmoid attention improves its forward kernel while the complete model step changes much less.Source: §6 · p. 16
• The measurements span GB200 BF16/FP16 programs, four principal integration sites, same-checkpoint substitutions, and paired pre-training trajectories measured with training loss.The evaluation covers symmetry, coefficient rounding, backward rules, packed schedules, and consuming operators as one integrated kernel program.Source: §6 · pp. 16–17
• The workflow locates transcendental work that remains exposed in generative modelling as matrix, special-function, and memory pipelines change across hardware generations.Source: §6 · p. 17

Generative-AI Disclosure

• Generative AI tools assisted with code development, experiment orchestration, analysis, figure generation, and drafting and editing the report.Source: Generative-AI Disclosure · p. 17
• The author reviewed and verified the code, experiments, numerical results, figures, citations, and final text.Source: Generative-AI Disclosure · p. 17
• The author takes responsibility for the report's content.Source: Generative-AI Disclosure · p. 17

A Extended Numerical Results

• The D3–D6 sweeps replay BF16-rounded coefficients in runtime evaluation order at BF16-representable inputs.Forward and derivative fits are evaluated separately.Source: A Extended Numerical Results · p. 19
• Figure 4 reports forward absolute error against the FP32 reference for D3–D6 BF16 target-precision programs at BF16-representable inputs.Source: Figure 4 · p. 19
• Figure 5 reports backward absolute error against the FP32 reference for D3–D6 BF16 target-precision programs at BF16-representable inputs.Source: Figure 5 · p. 20

B Auxiliary Integration Results

• Two auxiliary experiments pair inner-arithmetic measurements with the surrounding FA4 and rotary-position-embedding operations.Source: B Auxiliary Integration Results · p. 20
• The auxiliary experiments consider measurements alongside surrounding operations, rather than inner arithmetic alone.Source: B Auxiliary Integration Results · p. 20
• The surrounding operations named for these auxiliary measurements are FA4 and rotary-position embedding.Source: B Auxiliary Integration Results · p. 20

B.1 Further exp2 simplification in FlashAttention-4

• Online-softmax fits preserve p(0) = 1, p(1) = 2, and interval continuity; D2 is more accurate than PWL2 at equal two-FMA Horner depth.The sampled maximum-relative-error objective uses a 1,000,001-point grid.Source: §B.1 · p. 21
• FP16 PWL2 has maximum relative error 1.132 × 10^-2 and an 11-instruction core; its hinge subtraction and maximum replace a saved FMA.Repeated register-resident arithmetic exposes the shorter fractional core, while working-set size and shared reduction and reconstruction change the measured balance.Source: §B.1 · pp. 21–22
• FA4 keeps range reduction, exponent reconstruction, masking, and most eligible lane pairs on the existing path; matched variants differ only in the fractional evaluator.The matched rows route 33.3% of eligible forward and 12.5% of backward lane pairs through software.Source: §B.1 · p. 22
• 1.0008x geometric-mean PWL2 forward-plus-backward speedup in 80-step Llama-8B probes; candidate-minus-control loss deltas after 80 steps range from -0.323 to -0.050.The probes report PWL2 ratios of 1.0025x and 0.9990x in opposite execution orders.Source: §B.1 · p. 23
• 2.6 µs saved in FA4 forward scales to about 0.15 ms across 32 Llama layers, roughly 0.05% of a 275 ms optimizer step.The estimate is consistent with the complete-step timing.Source: §B.1 · p. 23

B.2 Paired sine/cosine and fused RoPE

• RoPE forms a shared periodic reduction for sine and cosine, evaluates packed FP16 D5/D6 polynomials, and reuses pairs across heads; retained timings do not establish a final quantitative comparison.The fused kernel applies the packed complex rotation and writes transposed Q and K layouts.Source: §B.2 · p. 23

C Packed PTX Listings

• The B3 listing selects coefficient rows with packed comparisons and lop3, then uses three packed FMAs; its derivative evaluator adds a factored linear tail.The probability and derivative paths share packing, region selection, clamping, and the cubic chain.Source: C Packed PTX Listings · pp. 24–25

D Lessons and Boundaries

• Generated instruction schedules, backward error, and algorithmic invariants constrain optimization; the isolated 1.898x fractional-core speedup becomes 1.0008x at the model step.Direct backward fitting is used when differentiated forward error exceeds the training budget, and online-softmax approximations preserve interval endpoints.Source: D Lessons and Boundaries · p. 25

E Code and Reproducibility

• The standalone repository contains the paper-specific integration layer and fitting, packed-kernel, open-weight, and materialized-evidence paths.The report build command and artifact map are documented in the report README.Source: E Code and Reproducibility · p. 25
• Retained artifacts cover isolated and integrated timing, paired pre-training loss curves, and open-weight evaluation results.The section names the corresponding data files and protocol records.Source: E Code and Reproducibility · pp. 25–26
• Complete launch manifests for initialization and data order are not part of the artifact set.The paired-loss file contains sampled histories and retained run fragments.Source: E Code and Reproducibility · p. 25
• Implementation directories cover fitting, packed kernels, portable B1–B4 plumbing, and B3 inline PTX.The inline-PTX source resides in a pinned FlashAttention submodule.Source: E Code and Reproducibility · p. 26