Fast Polynomial Transcendentals for LLMs

Robert Hu

cs.LGSep 3, 2026

Fast Polynomial Transcendentals for LLMs

Authors: Robert Hu

Organizations: GRAPHCORE RESEARCH

Abstract

Graphics processing unit (GPU) generations scale matrix, special-function, and memory pipelines at different rates, so kernel bottlenecks move as hardware evolves. FlashAttention-4 exposed this imbalance inside attention on NVIDIA Blackwell. We test whether short polynomial programs can accelerate other special-function-unit (SFU) operations in large language models (LLMs). We first compare native PyTorch evaluation with packed fused multiply--add (FMA) programs in an isolated IEEE binary16 (FP16) sweep spanning L2-resident and high-bandwidth-memory (HBM)-resident working sets. We then replace native sigmoid, tanh, and sigmoid linear unit (SiLU) with degree-3 or degree-4 bfloat16 (BF16) programs in four GB200 integration tasks: dense SiLU, tanh-softcapped attention, sigmoid attention, and routed-expert Swish-gated linear unit (SwiGLU). The programs combine analytical symmetry, target-format rounding, and packed arithmetic inside consuming kernels. The isolated paths improve by 1.19--2.19x in L2 and 1.00--1.70x in HBM. The dense-SiLU, tanh-softcapped-attention, and routed-expert substitutions improve complete training-step throughput by 2.7%, 2.9%, and 8.0%, respectively. The sigmoid-attention substitution improves complete-attention forward by 7.4% and the complete GPU step by 0.3%. Same-checkpoint open-weight ablations and one paired pre-training comparison per task extend the evaluation to model behavior. At common horizons near 100 billion tokens, the final smoothed training-loss differences (polynomial minus native) range from $-0.107$ to $+0.079$ across the four tasks.

Figures & tables

Case	Baseline	Polynomial intervention	Question
B1	Native SiLU in dense MLPs	Packed D3 SiLU	Does packed SiLU increase dense-model training-step throughput?
B2	Tanh-softcapped attention	D4 tanh inside FA4	Does packed tanh reduce complete softcapped-attention runtime?
B3	Native sigmoid attention	Direct D3 sigmoid and D4 derivative	Does polynomial sigmoid reduce complete sigmoid-attention runtime?
B4	Native SiLU in routed experts	Fused D3 polynomial Swish-gated linear unit (SwiGLU)	Does packed SwiGLU increase routed-expert training-step throughput?

Table 1: The four principal interventions. D3 and D4 denote degree-3 and degree-4 polynomials. B3 replaces sigmoid and its derivative; softmax exponential is outside its scope.

Figure 1: Analytic reduction before fitting. The lightly shaded positive half-domain is the only side that needs a polynomial. (a) Centering sigmoid exposes an odd residual, so the negative half is recovered by a sign-bit operation. (b) Subtracting

x/2

from SiLU exposes an even residual despite SiLU itself being asymmetric. (c) Centering the SiLU derivative exposes an odd residual and enforces

s^{\prime}(x)+s^{\prime}(-x)=1

by construction.

Figure 2: Constrained fitting and target-program replay, illustrated with the selected D4 centered-sigmoid candidate. (a) In normalized coordinates, every free basis term $u^{j}-u$ vanishes at both endpoints, so optimization cannot violate the two fixed values. (b) Fit interval $L_{i}$ and runtime clamp $L_{c}$ are searched independently; gray cells violate the $L_{c}\leq L_{i}+1$ search constraint and the star marks the selected $(5.0,5.28125)$ configuration. (c) Absolute-error curves separate the FP64 fit, BF16 coefficient rounding, the fitter’s split multiply/add replay, and a packed-FMA-style replay that rounds once per multiply-add. The execution schedule changes both error magnitude and where the worst error occurs.

Direction	Function	Reduced program	Endpoint-constrained least squares	Sollya minimax
			Maximum absolute error
Forward	Sigmoid	D3, $\tfrac{1}{2}+$ odd	5.44	5.47
Forward	Tanh	D4, odd	9.61	7.22
Forward	SiLU	D3, $x/2+$ even	21.74	19.95
Forward	GELU	D5, $x/2+$ even	1.06	0.50
Backward	Sigmoid	D4, even	4.46	1.97

Table 2: Maximum absolute error (

\times 10^{-3}

) for representative constrained fits and same-form Sollya fpminimax controls. Both use BF16-rounded coefficients under the same real-arithmetic replay; intermediate target-precision rounding is evaluated in appendix A .

Direction	Function	Polynomial program	4M elements (L2 cache)	268M elements (HBM)
Forward	Sigmoid	D3 polynomial	1.702x	1.249x
Forward	Tanh	D4 polynomial	1.187x	1.143x
Forward	SiLU	D3 polynomial	1.536x	1.267x
Forward	GELU	D5 polynomial	1.673x	1.702x
Backward	Sigmoid	D4 direct derivative	2.135x	1.004x
Backward	Tanh	D4 direct derivative	2.184x	1.003x

Table 3: Isolated FP16 function speed on GB200. Speedup is native PyTorch time divided by packed polynomial CUDA time, rounded to three decimals. The two columns report 4,194,304 elements (L2-cache resident) and 268,435,456 elements (HBM resident). Backward polynomial rows are direct derivative fits evaluated from

x

; native sigmoid and tanh backward use precomputed forward outputs.

Case	Timed component	Program	Forward speedup	Backward speedup
B1	One-kernel SiLU activation	D3 direct-fit probe	1.055x	1.007x
B2	Complete causal GQA attention	D4 tanh softcap	2.053x	1.046x
B3	Complete causal GQA attention	D3/D4 sigmoid	1.074x	1.001x
B4	One-kernel expert activation	D3 SwiGLU	1.003x	1.005x

Table 4: B1–B4 component performance. Speedup is native time divided by polynomial time; ratios are rounded to three decimals. Each ratio uses the medians of five alternating, 100-iteration round means. B1 evaluates 58,720,256 values; B2 and B3 use causal grouped-query attention (GQA) with batch 1, sequence length 4096, 32 query heads, 8 key/value heads, and head dimension 128; B4 evaluates 24,576 rows with hidden size 1280. The B1 row is a unary direct-derivative probe; the complete-step B1 measurements use the deployed packed SwiGLU path described below.

Case	Model / intervention	Forward speedup	Backward speedup	GPU-step speedup	Training-throughput ratio
B1	Llama 3 8B, dense D3 SiLU	1.005x	1.038x	1.027x	1.027x
B2	Llama 3 8B, D4 tanh softcap	1.107x	1.004x	1.029x	1.029x
B3	Llama 3 8B, D3/D4 sigmoid attention	1.011x	0.997x	1.003x	1.002x
B4	DeepSeek-V3 27A4B, expert D3 SwiGLU	1.087x	1.077x	1.080x	1.080x

Table 5: Complete-model performance. B1–B3 use matched 100-step Llama 3 8B probes on one GB200. B4 uses the complete DeepSeek-V3 27A4B architecture on four GB200s. Timing speedup is native time divided by polynomial time; training-throughput ratio is polynomial tokens/s divided by native tokens/s. Ratios use medians over steps 20–100, are rounded to three decimals, and favor the polynomial above one.

Model	Operator	Fit method	Mean accuracy change (percentage points)	WikiText-2 word perplexity	Prefill throughput ratio	Decode throughput ratio
Qwen2.5 7B	Dense SwiGLU	Endpoint LS	-0.51	$8.734\!\rightarrow\!197.677$	1.035x	1.013x
		Sollya minimax	-39.80	$8.734\!\rightarrow\!298{,}700.735$	1.036x	1.015x
Qwen3 30B-A3B	Expert SwiGLU	Endpoint LS	-0.10	$7.648\!\rightarrow\!7.655$	1.056x	1.025x
		Sollya minimax	+0.04	$7.648\!\rightarrow\!7.651$	1.051x	1.025x
GLM-4.7 Flash	Expert SwiGLU + router sigmoid	Endpoint LS	-0.30	$11.406\!\rightarrow\!11.400$	0.969x	1.004x
		Sollya minimax	-0.18	$11.406\!\rightarrow\!11.614$	0.966x	1.003x

Table 6: Same-checkpoint degree-4 ablations. “Endpoint LS” is this paper’s endpoint-constrained least-squares fit; “Sollya minimax” is the same-form control. Mean accuracy is the unweighted mean of seven tasks; its change is variant minus native and is rounded to 0.01 percentage points. WikiText-2 word perplexity is native

\rightarrow

variant and is rounded to three decimals. Throughput ratios are variant divided by native, are rounded to three decimals, and favor the polynomial above one.

Case	Intervention	Shared horizon (billion tokens)	Native EWM training loss	Polynomial EWM training loss	Training-loss change
B1	Dense D3 SiLU	100.2701	2.6264	2.7050	$+0.0786$
B2	D4 tanh softcap	100.2701	2.5415	2.5439	$+0.0025$
B3	D3/D4 sigmoid attention	100.2701	4.9471	4.8397	$-0.1073$
B4	Routed-expert D3 SwiGLU	100.2439	2.1532	2.1560	$+0.0028$

Table 7: Paired pre-training at the common observed horizon. Loss is the exponentially weighted mean (EWM) of sampled training loss with a half-life of one billion tokens. Training-loss change is polynomial minus native. Losses and changes are rounded to four decimals. Each row is one paired trajectory.

Figure 3: Paired B1–B4 pre-training loss. Gray is the native/SFU control and coral is the polynomial candidate. Thin lines are uniformly thinned sampled training loss; thick lines are token-domain exponentially weighted means (EWMs) with a half-life of one billion tokens, interpolated after smoothing onto a shared grid with 50 million tokens between points. Each panel ends at the common observed horizon. Insets show the final 25 billion tokens.

Appendix figures & tables5 assets

Supplementary material from the paper’s appendix.

Appendix

Panel of Figure 4 — Figure 4: Forward numerical accuracy for D3–D6 BF16 target-precision programs. Absolute error is measured against the FP32 reference at BF16-representable inputs.

Panel of Figure 5 — Figure 5: Backward numerical accuracy for D3–D6 BF16 target-precision programs. Absolute error is measured against the FP32 reference at BF16-representable inputs.

Family	Objective	Fractional core	Maximum relative error	Root-mean-square relative error
D1	Endpoint constrained	1 packed FMA	$6.148\times 10^{-2}$	$4.466\times 10^{-2}$
PWL2	Least squares	2 packed FMAs + hinge	$1.526\times 10^{-2}$	$6.754\times 10^{-3}$
PWL2	Sampled minimax	2 packed FMAs + hinge	$1.074\times 10^{-2}$	$7.174\times 10^{-3}$
D2	Least squares	2 packed FMAs	$3.204\times 10^{-3}$	$1.987\times 10^{-3}$
D2	Sampled minimax	2 packed FMAs	$2.679\times 10^{-3}$	$1.921\times 10^{-3}$
FA4 D3	Sollya/minimax	3 packed FMAs; $p(1)=1.99983$	$8.763\times 10^{-5}$	$6.029\times 10^{-5}$

Appendix

Table 8: Fractional

2^{f}

fits on

f\in[0,1]

. Errors are measured in FP64 on a dense grid before coefficient quantization. “Core” excludes common range reduction and exponent reconstruction. D1, PWL2, and D2 preserve both endpoints exactly.

Measurement scope	FA4 D3	PWL2 FP16	D2 FP32
Fractional core, repeated compute	–	1.898x	–
Full $2^{x}$ , repeated compute	0.884x	0.863x	–
Full $2^{x}$ , 4 MiB input + 4 MiB output (L2-sized)	0.995x	1.002x	–
Full $2^{x}$ , 256 MiB input + 256 MiB output (HBM-sized)	0.995x	0.990x	–
Causal FA4 forward	1.012x	1.010x	1.012x
Causal FA4 forward + backward	1.001x	1.005x	1.006x

Appendix

Table 9: The exp2 performance ladder on GB200. Entries are speedups over the matching native/all-SFU path. Full-exp rows include unchanged FP32 reduction and reconstruction; FA4 rows use identical causal routing. Values are rounded to three decimals, except the 8B optimizer-step probes, which are rounded to four decimals.

Measured quantity	FA4 D3	PWL2 FP16	D2 FP32
Maximum over output and input-gradient tensors	$8.641\times 10^{-4}$	$3.458\times 10^{-3}$	$2.019\times 10^{-3}$

Appendix

Table 10: Maximum tensor-relative

L_{2}

error against all-SFU FA4 on GB200.

GPU Kernel Optimization

GPU kernel optimization improves the speed and memory efficiency of computations by tailoring code to GPU hardware and workloads.

See all 202 papers

Hardware-Aware FP4 FlashAttention-4
Sep 3, 2026Robert HuMixed-Precision Quantization GPU Kernel Optimization
GRAPHCORE RESEARCH
LampAttention: Look-Ahead Mixed-Precision FlashAttention for Dedicated Accelerators
Sep 30, 2026Stanislav Budzinskiy, Marian Gloser, Tolunay Yilmaz +5AI Accelerator Inference LLM Inference Acceleration
Faculty of Mathematics University of Vienna, Austria · Huawei Heisenberg Research Center, Munich, Germany · Huawei Technologies Co. Ltd
Tile-Level Activation Overlap for Efficient LLM Inference
May 5, 2026Abhinav Jangda, Tyler Sorensen, Sebastian Burckhardt +3LLM Inference Efficiency GPU Kernel Optimization
Microsoft Research · Microsoft

Sep 3, 2026J/K move · Enter open · S save

Hardware-Aware FP4 FlashAttention-4

Robert Hu

GRAPHCORE RESEARCH

Mixed-Precision Quantization GPU Kernel Optimization Transformer Attention Efficient Attention Attention Quantization Attention Mechanisms

Open paper AskPDF arXiv

Blackwell's 4-bit floating-point (FP4) tensor cores do not automatically make attention faster because softmax conversion and on-chip dependencies dominate once its matrix products shrink. We address this with Direct-P for noncausal inference and a causal path that passes the forward quantization directly into backward. Direct-P maps scores directly to FP4 probabilities and reaches up to 2.13

\times

the bfloat16 (BF16) forward throughput on an NVIDIA GB200. The causal path reconstructs probabilities from saved quantized queries and keys and uses 8-bit floating-point (FP8) gradient operands, accelerating a complete single-GPU 8-billion-parameter update by up to 1.14

\times

. Matched distributed training retains FP8 probabilities and values; every tested MXFP4 probability/value training trajectory diverges.

Abstract

GPU Kernel Optimization

Explore similar work

Hardware-Aware FP4 FlashAttention-4

LampAttention: Look-Ahead Mixed-Precision FlashAttention for Dedicated Accelerators

Tile-Level Activation Overlap for Efficient LLM Inference