面向 LLM 的快速多项式超越函数

# Fast Polynomial Transcendentals for LLMs
Source: [https://arxiv.org/html/2610.00049](https://arxiv.org/html/2610.00049)
Sept 2026

###### 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\-0\.107to\+0\.079\+0\.079across the four tasks\.

Keywords:GPU kernels, polynomial approximation, transcendental functions, FlashAttention, BF16, mixture of experts

## 1Introduction

Graphics processing unit \(GPU\) generations scale matrix, special\-function\-unit \(SFU\), and memory pipelines at different rates\. FlashAttention\-4 \(FA4\) provides a concrete example: on Blackwell, Tensor Core throughput increased more quickly than the resources used for exponentials and shared\-memory traffic, exposing non\-matrix work inside an optimized attention kernel\[[24](https://arxiv.org/html/2610.00049#bib.bib1)\]\. FA4 recovered performance by redesigning the algorithm and its instruction pipeline around that imbalance\.

Large language models \(LLMs\) expose the same resource imbalance beyond softmax\. Dense multilayer perceptrons \(MLPs\), attention variants, and routed experts interleave matrix multiplication with sigmoid, tanh, and sigmoid linear unit \(SiLU\)\. The native forward paths studied here evaluate these transcendentals through SFU\-backed instruction sequences\. Native backward paths either reevaluate the required transcendental or apply an algebraic rule to a forward result, depending on the operator\. A short polynomial can instead use packed bfloat16 \(BF16\) or IEEE binary16 \(FP16\) arithmetic on CUDA\-core pipelines\. We trace each candidate from target\-precision error through consuming\-component runtime and complete\-step throughput, then measure same\-checkpoint task accuracy, word perplexity, prefill and decode throughput, and paired pre\-training loss\.

We investigate the four integration tasks summarized in[table1](https://arxiv.org/html/2610.00049#S1.T1)\.

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\.Throughout, “native” denotes the unchanged implementation of the same transcendental for the same integration task\.

This paper makes three contributions:

- •We construct target\-precision programs by combining analytical range reduction, structurally constrained fitting, and replay of coefficient and intermediate rounding in deployed order\. Same\-form Sollya minimax fits provide coefficient controls\.
- •We lower the accepted programs to packed CUDA and inline Parallel Thread Execution \(PTX\) code in four consuming LLM components\. B1, B2, and B4 improve complete\-step throughput by 2\.7%, 2\.9%, and 8\.0%; B3 improves complete attention forward by 7\.4%\. We release the fitting tools, packed kernels, portable integration and evaluation code, and supporting artifacts; repository and reproduction details appear in[appendixE](https://arxiv.org/html/2610.00049#A5)\.
- •We evaluate degree\-4 substitutions in five open\-weight models and compare native and polynomial pre\-training trajectories for each case\. At common horizons near 100 billion tokens, final smoothed training\-loss differences are\+0\.079\+0\.079,\+0\.003\+0\.003,−0\.107\-0\.107, and\+0\.003\+0\.003for B1–B4\.

## 2Background and Related Work

### 2\.1Asymmetric accelerator scaling

Tensor Cores execute the matrix operations that dominate Transformer layers, whereas CUDA\-core and SFU pipelines execute scalar and packed arithmetic around those matrices\. The L2 cache and high\-bandwidth memory \(HBM\) supply different working\-set regimes\. Accelerating one resource changes which of the others is visible in an integrated kernel\.

FA4 quantified this effect on Blackwell: matrix throughput increased much more quickly than the exponential units used by softmax, allowing non\-matrix work to occupy a substantial share of the attention schedule\[[24](https://arxiv.org/html/2610.00049#bib.bib1)\]\. Blackwell Ultra subsequently doubled SFU throughput relative to Blackwell for key attention instructions; NVIDIA identifies the hardware exponential asMUFU\.EX2in NVIDIA machine code \(SASS\)\[[3](https://arxiv.org/html/2610.00049#bib.bib14),[14](https://arxiv.org/html/2610.00049#bib.bib13)\]\. NVIDIA reports that Rubin, the GPU in the Vera Rubin platform, provides up to 22 TB/s of HBM4 bandwidth, doubles Tensor Core throughput per clock relative to Blackwell, and raises exponential throughput by 2x for 32\-bit floating\-point \(FP32\) arithmetic and 4x for BF16/FP16\[[2](https://arxiv.org/html/2610.00049#bib.bib15)\]\. These vendor\-reported changes illustrate why the active bottleneck can move between hardware generations; all measurements in this paper use GB200\.

### 2\.2Target\-precision numerical approximation

Sollya provides high\-precision function analysis, minimax construction, coefficient rounding, and code\-generation support\[[5](https://arxiv.org/html/2610.00049#bib.bib7)\]\. We use it as the reference for unconstrained minimax fits\. RLIBM instead constructs polynomials around rounding intervals, demonstrating that output representation and rounding behavior belong inside the approximation problem\[[15](https://arxiv.org/html/2610.00049#bib.bib6)\]\. RLIBM targets correctly rounded libraries; our target is a lower\-latency program inside a specific BF16/FP16 neural\-network kernel\.

AutoNumerics\-Zero searches directly over limited\-precision programs and finds approximations that trigger favorable compiler paths\[[20](https://arxiv.org/html/2610.00049#bib.bib2)\]\. Our search instead fixes analytical symmetry, saturation, and a fused\-multiply–add \(FMA\) budget before optimizing the rounded program that executes\.

K\-TanH approximates tanh using low\-precision integer operations, derives sigmoid, Swish, and Gaussian error linear unit \(GELU\) variants, and reports both CPU speed and Google Neural Machine Translation \(GNMT\) model convergence\[[13](https://arxiv.org/html/2610.00049#bib.bib4)\]\. Tensor Processing Primitives use a range\-limited tanh polynomial evaluated with four Horner FMAs in a portable CPU tensor library\[[11](https://arxiv.org/html/2610.00049#bib.bib5)\]\. Our work targets packed BF16 using CUDA and CuTe, NVIDIA’s C\+\+ CUDA template abstractions for hierarchical tensor layouts, inside Blackwell attention and expert kernels rather than portable CPU operators\.

Non\-uniform linear interpolation \(NLI\) selects interpolation cut points by dynamic programming and applies the resulting piecewise\-linear representation to nonlinear operators including SiLU, root mean square layer normalization \(RMSNorm\), and Softmax\[[23](https://arxiv.org/html/2610.00049#bib.bib3)\]\. NLI proposes a universal hardware unit; our programs remain register\-resident and use packed arithmetic on an unmodified GPU\.

### 2\.3Nonlinear operators in LLM components

Low\-Cost FlashAttention proposes fused exponential\-and\-multiplication hardware\[[1](https://arxiv.org/html/2610.00049#bib.bib12)\]; we instead change the software program mapped to existing hardware\. Sigmoid self\-attention can match softmax attention when normalized and stabilized appropriately, and FlashSigmoid provides a memory\-efficient implementation\[[19](https://arxiv.org/html/2610.00049#bib.bib8)\]\. B3 adopts that algorithmic setting and compares an SFU sigmoid with a packed polynomial sigmoid\.

Gemma 2 applies

softcapα⁡\(x\)=α​tanh⁡\(x/α\)\\operatorname\{softcap\}\_\{\\alpha\}\(x\)=\\alpha\\tanh\(x/\\alpha\)\(1\)to attention logits\[[10](https://arxiv.org/html/2610.00049#bib.bib9)\], giving B2 a direct open\-weight analogue\. SwiGLU was established as a strong Transformer feed\-forward choice\[[22](https://arxiv.org/html/2610.00049#bib.bib10)\], and DeepSeek\-V3 makes routed SwiGLU experts a representative production\-scale mixture\-of\-experts \(MoE\) workload\[[9](https://arxiv.org/html/2610.00049#bib.bib11)\]\. These architectures motivate B1 and B4\.

## 3Target\-Precision Polynomial Programs

Our optimization object is the complete low\-precision program in BF16 or FP16: analytical range reduction, polynomial coefficients, instruction order, intermediate rounding, tail behavior, and the backward rule\. A real\-arithmetic polynomial can change materially after its coefficients and Horner stages are rounded\. We therefore choose the mathematical form first, fit its residual behavior, and replay the exact target\-precision schedule before integrating it into a kernel\. We write Dkkfor a degree\-kkpolynomial\.

### 3\.1Analytical reduction

Analytical identities reduce the fitted domain and remove coefficients for known offsets, linear terms, and reflected values\. Let

x^=clip⁡\(x,−Lc,Lc\),t=\|x^\|\.\\hat\{x\}=\\operatorname\{clip\}\(x,\-L\_\{c\},L\_\{c\}\),\\qquad t=\|\\hat\{x\}\|\.\(2\)For sigmoid, the centered functiongσ​\(x\)=σ⁡\(x\)−12g\_\{\\sigma\}\(x\)=\\sigma\(x\)\-\\tfrac\{1\}\{2\}is odd becauseσ⁡\(−x\)=1−σ⁡\(x\)\\sigma\(\-x\)=1\-\\sigma\(x\)\. We therefore fit a polynomialpDp\_\{D\}only on the positive half\-domain and reconstruct

σ~​\(x\)=12\+sign⁡\(x\)​pD​\(t\),pD​\(0\)=0\.\\widetilde\{\\sigma\}\(x\)=\\frac\{1\}\{2\}\+\\operatorname\{sign\}\(x\)p\_\{D\}\(t\),\\qquad p\_\{D\}\(0\)=0\.\(3\)The same construction applies to the odd function tanh\. Sigmoid and tanh derivatives are even, so their direct backward programs evaluate a polynomial ofttwithout sign reconstruction\. Saturating targets use a clamp or an exact tail outside the active interval\. This both prevents polynomial blow\-up on outliers and lets a short interior polynomial spend its error budget where the function has useful curvature\.

Althoughs⁡\(x\)=x​σ​\(x\)s\(x\)=x\\sigma\(x\)is neither odd nor even, its asymmetric part is known exactly:

s⁡\(x\)−s⁡\(−x\)\\displaystyle s\(x\)\-s\(\-x\)=x,\\displaystyle=x,\(4\)s⁡\(x\)−x2\\displaystyle s\(x\)\-\\frac\{x\}\{2\}=x⁡\(σ⁡\(x\)−12\)\.\\displaystyle=x\\left\(\\sigma\(x\)\-\\frac\{1\}\{2\}\\right\)\.\(5\)The residual in the second line is even, because it is the product of two odd functions\. Consequently the forward program can be written as

s~​\(x\)=x2\+t​pD​\(t\),\\widetilde\{s\}\(x\)=\\frac\{x\}\{2\}\+t\\,p\_\{D\}\(t\),\(6\)without constructing a full sigmoid tensor\. Exact GELU,x​Φ​\(x\)x\\Phi\(x\), has the same decomposition sinceΦ⁡\(−x\)=1−Φ⁡\(x\)\\Phi\(\-x\)=1\-\\Phi\(x\)\. These identities are exact; the approximation is confined to one half\-domain residual\.

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\) Subtractingx/2x/2from SiLU exposes an even residual despite SiLU itself being asymmetric\. \(c\) Centering the SiLU derivative exposes an odd residual and enforcess′​\(x\)\+s′​\(−x\)=1s^\{\\prime\}\(x\)\+s^\{\\prime\}\(\-x\)=1by construction\.[Figure1](https://arxiv.org/html/2610.00049#S3.F1)shows why these transformations reduce more than the sample domain\. They remove coefficients that would otherwise relearn an exact offset, linear term, or reflected half of the function\. The resulting runtime reconstruction uses absolute value, sign manipulation, and at most one exact affine term around a short half\-domain Horner chain\.

### 3\.2Constrained fitting and rounded\-program replay

Classical minimax construction, target\-format rounding, and finite\-precision program search answer different questions\. Sollya supplies the real\-arithmetic minimax control\[[5](https://arxiv.org/html/2610.00049#bib.bib7)\]; RLIBM shows why rounding intervals belong in the specification\[[15](https://arxiv.org/html/2610.00049#bib.bib6)\]; and AutoNumerics\-Zero demonstrates direct search over executable low\-precision programs\[[20](https://arxiv.org/html/2610.00049#bib.bib2)\]\. Our workflow combines these ideas under a fixed instruction budget\.

The fitting procedure searches two boundaries\.LiL\_\{i\}is the interval used to construct the polynomial, whileLcL\_\{c\}is the runtime clamp\. Separating them is useful because the least\-squares endpoint and the best post\-quantization tail transition need not coincide\. For each degree andLiL\_\{i\}, we solve in 64\-bit floating point \(FP64\) a dual endpoint\-constrained least\-squares problem

P⁡\(0\)=f⁡\(0\),P⁡\(Li\)=f⁡\(Li\)\.P\(0\)=f\(0\),\\qquad P\(L\_\{i\}\)=f\(L\_\{i\}\)\.\(7\)Eliminating the first two coefficients gives the numerically convenient basis

P⁡\(x\)=f⁡\(0\)\+x​f⁡\(Li\)−f⁡\(0\)Li\+∑j=2Daj\(xj−xLij−1\),\\begin\{split\}P\(x\)=\{\}&f\(0\)\+x\\frac\{f\(L\_\{i\}\)\-f\(0\)\}\{L\_\{i\}\}\\\\ &\+\\sum\_\{j=2\}^\{D\}a\_\{j\}\\left\(x^\{j\}\-xL\_\{i\}^\{j\-1\}\\right\),\\end\{split\}\(8\)whose free basis terms vanish at both endpoints\. The remaining coefficients are solved from 2,000 Chebyshev nodes\. Odd, centered\-odd, and even transformations are applied before this solve, so every candidate preserves the corresponding global identity by construction rather than approximately\.

![Three panels showing endpoint-constrained basis functions, a search over fit and clamp intervals, and error curves for exact and rounded program evaluations.](https://arxiv.org/html/2610.00049v1/constrained_fit_replay.png)Figure 2:Constrained fitting and target\-program replay, illustrated with the selected D4 centered\-sigmoid candidate\. \(a\) In normalized coordinates, every free basis termuj−uu^\{j\}\-uvanishes at both endpoints, so optimization cannot violate the two fixed values\. \(b\) Fit intervalLiL\_\{i\}and runtime clampLcL\_\{c\}are searched independently; gray cells violate theLc≤Li\+1L\_\{c\}\\leq L\_\{i\}\+1search constraint and the star marks the selected\(5\.0,5\.28125\)\(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\.The three panels in[figure2](https://arxiv.org/html/2610.00049#S3.F2)correspond to three design choices: enforce global structure algebraically, search the active interval jointly with the tail transition, and score the program that executes rather than only its real\-arithmetic polynomial\. The BF16 search applies these choices as follows:

1. 1\.enumerate degrees D3–D6 and the function\-specificLi,LcL\_\{i\},L\_\{c\}ranges at 0\.25\-unit spacing;
2. 2\.solve[equation8](https://arxiv.org/html/2610.00049#S3.E8)in FP64 on the reduced domain;
3. 3\.round every coefficient to BF16;
4. 4\.evaluate the configured Horner chain with target\-BF16 rounding at each multiply and add;
5. 5\.apply the same absolute\-value reduction, clamp, sign restoration, and offset as the CUDA implementation;
6. 6\.score maximum absolute error on 8,000 points over\[0,12\]\[0,12\], rejectingLc\>Li\+1L\_\{c\}\>L\_\{i\}\+1and any shape or boundary violation;
7. 7\.retain distribution\-weighted alternatives only after the uniform error and monotonicity gates pass\.

The search emulator rounds multiplication and addition as separate BF16 operations, whereas the deployed\_\_hfma2or PTX instruction rounds the fused result once\. Final selection replays the deployed packed\-FMA order on device\. The last panel in[figure2](https://arxiv.org/html/2610.00049#S3.F2)also explains why a coefficient set can improve after BF16 rounding in real Horner evaluation yet regress once intermediate round points are introduced\. Final acceptance therefore uses the deployed packed\-FMA order\.

### 3\.3Forward and backward programs

We reduce forward and backward separately\. Centered sigmoid is odd and its derivative is even\. Tanh is odd and its derivative is even\. Subtractingx/2x/2from SiLU or exact GELU produces an even residual; subtracting1/21/2from either derivative produces an odd residual\. We impose these identities as algebraic fitting constraints\.

For SiLU,

s′​\(x\)=σ⁡\(x\)\+x​σ′​\(x\),s′​\(x\)\+s′​\(−x\)=1,s^\{\\prime\}\(x\)=\\sigma\(x\)\+x\\sigma^\{\\prime\}\(x\),\\qquad s^\{\\prime\}\(x\)\+s^\{\\prime\}\(\-x\)=1,\(9\)sos′​\(x\)−12s^\{\\prime\}\(x\)\-\\tfrac\{1\}\{2\}is odd even though the forward function is asymmetric\. Exact GELU obeys the same derivative identity\. Our direct backward fit targets this centered odd residual and preserves the complementary\-gradient relationship\.

We use three backward constructions\. Algebraic identities express a derivative through a function value:

σ′​\(x\)\\displaystyle\\sigma^\{\\prime\}\(x\)=σ​\(x\)​\(1−σ​\(x\)\),\\displaystyle=\\sigma\(x\)\(1\-\\sigma\(x\)\),\(10\)tanh′⁡\(x\)\\displaystyle\\tanh^\{\\prime\}\(x\)=1−tanh2⁡\(x\)\.\\displaystyle=1\-\\tanh^\{2\}\(x\)\.\(11\)For softcap, the scale cancels:

dd​x​\[α​tanh⁡\(x/α\)\]=1−tanh2⁡\(x/α\)\.\\frac\{\\mathrm\{d\}\}\{\\mathrm\{d\}x\}\\left\[\\alpha\\tanh\(x/\\alpha\)\\right\]=1\-\\tanh^\{2\}\(x/\\alpha\)\.\(12\)The value can be saved from forward or recomputed in backward\. B1 and B4 recompute a D3 sigmoid approximationp⁡\(x\)p\(x\), formSiLU~​\(x\)=x​p​\(x\)\\widetilde\{\\operatorname\{SiLU\}\}\(x\)=xp\(x\), and evaluate its backward rule asp⁡\(x\)\+x​p​\(x\)​\(1−p⁡\(x\)\)p\(x\)\+xp\(x\)\(1\-p\(x\)\)\. A differentiated\-polynomial backward instead derives its coefficients from the selected forward polynomial; B2 differentiates its D4 tanh program and evaluates the resulting D3 program at the recomputed score\. A direct fit tof′f^\{\\prime\}controls gradient error independently and can exploit the simpler backward parity\. This direct program supplies a surrogate gradient rather than the analytic derivative of the selected forward polynomial\. The accepted B3 design couples a direct D3 sigmoid value with a separately targeted, factored D4 derivative surrogate\. We measure gradient error and training loss for the deployed construction in each case\.

### 3\.4Packed CUDA and inline PTX programs

We focus on 16\-bit floating\-point execution because NVIDIA exposes packed two\-lane arithmetic for both IEEE binary16 \(FP16\) and BF16\. The isolated function sweep uses FP16: each 32\-bit\_\_half2operand holds two values, and\_\_hfma2advances two Horner chains with one packed FMA\. The B1–B4 integrations use the corresponding BF16 path, with\_\_nv\_bfloat162and, for B3, explicitfma\.rn\.bf16x2PTX\.

The implementation removes scalar work from the instruction schedule in three steps\. First, B1, B2, and B4 process two BF16 values per arithmetic instruction\. On aligned workloads divisible by eight and above the vectorization threshold, one thread loads a 16\-byteint4chunk, reinterprets it as four\_\_nv\_bfloat162pairs, and evaluates each Horner stage with\_\_hfma2\. Second, the fused SwiGLU program combines activation and gating in forward and computes both projected\-input gradients in one backward kernel\. These choices amortize address calculation and avoid materializing an intermediate activation tensor\.

Third, B3 places the evaluator directly inside the FA4 instruction schedule\. A separate device routine would introduce an interface boundary and hide the register\-level schedule from the surrounding CuTe kernel, so B3 uses onellvm\.inline\_asmblock\. PTX instruction names and lane semantics follow NVIDIA’s PTX instruction set architecture \(ISA\)\[[18](https://arxiv.org/html/2610.00049#bib.bib16)\]\.

The B3 forward block is a branch\-free, two\-lane BF16 evaluator\. Onecvt\.rn\.bf16x2\.f32packs two FP32 scores into a 32\-bit register;min\.bf16x2andmax\.bf16x2clamp both lanes\. A packedset\.ge\.u32\.bf16x2produces the piecewise mask, andlop3\.b32with truth table0xcaselects embedded BF16 coefficient bit patterns independently in each lane\. The mask makes region selection branch\-free, and immediate coefficient bits remove coefficient loads\. Threefma\.rn\.bf16x2instructions then evaluate the D3 Horner chain\. The coefficients are fitted for sequence length 4096 and pre\-scaled by1/n1/n, so the evaluator does not materialize−log⁡n\-\\log nor apply a final normalization multiply\. The complete forward block contains 17 PTX instructions, including three packed FMAs, and is reproduced in\.

Backward reuses the forward conversion, clamp, region mask, and Horner value\. We express the factored piecewise\-D4 surrogate asP​\(x\)​\(a\+b​x\)P\(x\)\(a\+bx\)and return it from the same inline block as the D3 probability\. The derivative adds one packed FMA and one packed multiply\. The coupled value/derivative block therefore contains 22 PTX instructions and four packed FMAs, compared with 17 instructions for forward plus 22 for a standalone D4 derivative\. The complete shared schedule appears in\.

## 4Experimental Setup

All performance measurements use NVIDIA GB200\. Every operator and model setting outside the named substitution is held fixed\.

### 4\.1Function error and coefficient control

The target\-precision plots evaluate BF16\-rounded coefficients at BF16\-representable inputs and replay the runtime Horner order\. The reference value is computed in FP32\. We also compare the constrained coefficients with Sollya’sfpminimax\[[5](https://arxiv.org/html/2610.00049#bib.bib7)\]\. For each row, the control uses the same degree, sparse monomials, clamp, and algebraic reconstruction\. Sollya constructs an absolute minimax polynomial with each coefficient constrained to 8\-bit floating\-point precision; the parsed coefficients are then cast to BF16\. The comparison in[table2](https://arxiv.org/html/2610.00049#S5.T2)evaluates both coefficient sets in real arithmetic on 20,001 uniformly spaced points over\[−Lc,Lc\]\[\-L\_\{c\},L\_\{c\}\]\. This coefficient control omits intermediate rounding and is separate from the target\-precision replay\.

### 4\.2Isolated function timing

The standalone size sweep compares packed FP16 CUDA polynomial kernels with the corresponding native PyTorch paths on GB200\. Each path receives 1,000 warm\-up calls\. We then alternate native\-first and polynomial\-first order over nine rounds of 1,000 CUDA\-event\-timed calls and report the ratio of the median per\-round times\. The retained endpoints contain 4,194,304 elements, whose forward and backward working sets fit in L2 cache, and 268,435,456 elements, which stream through HBM\. Forward work is excluded from every backward timing\. Native sigmoid and tanh backward consume precomputed forward outputs and applyg​y​\(1−y\)g\\,y\(1\-y\)andg⁡\(1−y2\)g\(1\-y^\{2\}\), respectively\. Their polynomial comparators consumexxand evaluate separately fitted D4 derivative programs\. Those two rows therefore compare direct derivative fits with native algebraic backward, rather than comparing SFU instructions with FMAs\. This sweep exercises the FP16\_\_half2implementation; the B1–B4 component and model measurements use BF16\.

### 4\.3Component and training\-step timing

Speedup is

S=TnativeTpoly,S=\\frac\{T\_\{\\mathrm\{native\}\}\}\{T\_\{\\mathrm\{poly\}\}\},\(13\)so values above one favor the polynomial\. The retained component benchmark warms up for 25 iterations, then alternates native\-first and polynomial\-first order over five rounds of 100 CUDA\-event\-timed iterations\. We report the ratio of the median per\-round mean times\. B1 evaluates 58,720,256 BF16 activation elements\. 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\.

For B1–B3 full\-model timing, CUDA events bracket model\-plus\-loss forward, loss backward, and the complete GPU optimizer step\. Each Llama 3 8B run performs 100 optimizer steps with compilation enabled; medians use steps 20–100\. Native and polynomial runs use seed 1234, batch size 4, sequence length 4096, and the same fixed sequence of tokenized examples\. Token throughput is the training harness’s wall\-clock tokens per second\.

B4 uses the complete DeepSeek\-V3 27A4B architecture on four GB200s with expert parallelism \(EP\) of four, sequence length 4096, local batch one, and global batch four\. Here 27A4B denotes approximately 27 billion total parameters and 4 billion active parameters per token\. The same 100\-step protocol and steps 20–100 aggregation apply\. Only the 29 routed grouped\-expert SwiGLU modules change; routers, dense MLPs, and shared MLPs remain native\.

### 4\.4Same\-checkpoint evaluation

The auxiliary open\-weight study evaluates degree\-4 substitutions, which are separate from the D3 B1 and B4 programs\. Native, constrained\-fit, and same\-form Sollya\-fit variants run in the same process over identical requests\. The Qwen3 and Kimi checkpoint labels use A3B for approximately 3 billion parameters active per token\. The suite records Massive Multitask Language Understanding \(MMLU; 5\-shot\)\[[12](https://arxiv.org/html/2610.00049#bib.bib18)\], HellaSwag \(10\-shot\)\[[25](https://arxiv.org/html/2610.00049#bib.bib19)\], the Challenge set of the AI2 Reasoning Challenge \(ARC\-Challenge; 25\-shot\)\[[6](https://arxiv.org/html/2610.00049#bib.bib20)\], WinoGrande \(5\-shot\)\[[21](https://arxiv.org/html/2610.00049#bib.bib21)\], Physical Interaction: Question Answering \(PIQA; 0\-shot\)\[[4](https://arxiv.org/html/2610.00049#bib.bib22)\], Grade School Math 8K \(GSM8K; 5\-shot\)\[[7](https://arxiv.org/html/2610.00049#bib.bib23)\], TruthfulQA\-MC2 \(0\-shot\)\[[16](https://arxiv.org/html/2610.00049#bib.bib24)\], and WikiText\-2 word perplexity \(0\-shot\)\[[17](https://arxiv.org/html/2610.00049#bib.bib25)\]\. The reported task mean is the unweighted mean of the seven accuracy tasks\. The historical campaign launch source sets synthetic throughput to batch 1 and sequence length 2048 for Qwen2\.5; Qwen3, GLM\-4\.7 Flash, GPT\-OSS, and Kimi use sequence length 512\. Its defaults request 20 prefill measurements after five warm\-ups and ten 64\-token decode measurements after two warm\-ups\. Results are single\-run point estimates\. The retained GPT\-OSS helpers use clamps of 5\.25 for the constrained fit and 5\.28125 for the Sollya fit\.

### 4\.5Paired pre\-training

Each production comparison contains one native and one polynomial trajectory\. All runs use BF16 on 32 GB200 GPUs, sequence length 4096, and 2,000 learning\-rate warm\-up steps\. B1–B3 use 32\-way fully sharded data parallelism \(FSDP\) and global batch 32; B4 uses global batch 128 and EP=8\. Run fragments are ordered and joined by token count\.

The long\-run metric is sampled training loss\. We smooth each arm causally in token space with an exponentially weighted mean \(EWM\) whose half\-life is one billion tokens, then interpolate the smoothed curves onto a shared grid spaced every 50 million tokens through the last token count observed by both arms\. Logged training loss is the sole long\-run loss metric\.

## 5Results

### 5\.1Function\-level approximation

How does maximum absolute error after BF16 coefficient rounding compare with same\-form Sollya minimax controls?[Table2](https://arxiv.org/html/2610.00049#S5.T2)evaluates both coefficient sets with the same real\-arithmetic replay\.

Table 2:Maximum absolute error \(×10−3\\times 10^\{\-3\}\) for representative constrained fits and same\-form Sollyafpminimaxcontrols\. Both use BF16\-rounded coefficients under the same real\-arithmetic replay; intermediate target\-precision rounding is evaluated in[appendixA](https://arxiv.org/html/2610.00049#A1)\.The constrained D3 sigmoid fit has maximum absolute error5\.44×10−35\.44\\times 10^\{\-3\}, compared with5\.47×10−35\.47\\times 10^\{\-3\}for its Sollya control\. Sollya has lower error in the other seven rows\. The target\-precision sweeps measure intermediate BF16 rounding in runtime order\. GELU is an auxiliary function\-level result outside the B1–B4 integrations\.

The extended sweeps in[figures4](https://arxiv.org/html/2610.00049#A1.F4)and[5](https://arxiv.org/html/2610.00049#A1.F5)replay target\-precision coefficients in the runtime evaluation order\. Their lower panels use logarithmic absolute\-error axes\. In each of the six sweeps, at least one higher\-degree program has greater maximum error than a lower\-degree program after BF16 replay\.

The backward plots include cases where a low\-error value fit produces a larger derivative error after coefficient quantization\. B3 therefore pairs its D3 forward fit with a separately fitted D4 derivative\. The following measurements move from standalone functions to their consuming components\.

#### 5\.1\.1Isolated FP16 function speed

How does the packed polynomial path compare with the native PyTorch path before integration, and how does that comparison change with the memory regime?[Table3](https://arxiv.org/html/2610.00049#S5.T3)reports the cache\-resident and HBM\-resident endpoints of the standalone FP16 size sweep\.

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 fromxx; native sigmoid and tanh backward use precomputed forward outputs\.The cache\-resident endpoint ranges from 1\.187x to 2\.193x\. At the HBM\-resident endpoint, sigmoid and tanh backward measure 1\.004x and 1\.003x, while GELU backward measures 1\.615x\. The sigmoid and tanh backward paths each read two FP16 arrays and write one; their median times differ by less than 0\.4% in HBM\. Seven of the eight rows have a smaller speedup at the HBM\-resident endpoint; GELU forward measures 1\.673x in L2 and 1\.702x in HBM\.

### 5\.2Integration into LLM components

Which function\-level programs reduce component runtime after integration?[Table4](https://arxiv.org/html/2610.00049#S5.T4)reports matched component measurements before complete\-model timing\. B1 is a unary SiLU probe, B4 isolates the routed\-expert activation, and B2 and B3 time complete causal attention\.

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\.#### 5\.2\.1B1: dense SiLU

The unary D3 SiLU probe measures 1\.055x forward speedup and 1\.007x backward speedup on 58\.7M BF16 values; its backward is the standalone direct D3 derivative fit\. The deployed B1 packed SwiGLU backward instead recomputes the D3 sigmoid approximation from the gate input and forms the SiLU derivative algebraically\. The complete\-step measurements below exercise that deployed path\.

#### 5\.2\.2B2: tanh softcap

B2 replaces tanh in[equation1](https://arxiv.org/html/2610.00049#S2.E1)inside FA4\. Forward uses packed D4\. Backward differentiates the selected D4 coefficients and evaluates the resulting D3 polynomial at the recomputed, scaled attention score\. Complete causal attention measures 2\.053x forward speedup and 1\.046x backward speedup\.

#### 5\.2\.3B3: sigmoid attention

The direct sigmoid program accelerates complete attention forward\. For sequence lengthnn, attention scoresi​js\_\{ij\}is mapped to probabilityPi​jP\_\{ij\}, which weights the value matrixVVto produce outputOO:

Pi​j=σ⁡\(si​j−log⁡n\),O=P​V\.P\_\{ij\}=\\sigma\(s\_\{ij\}\-\\log n\),\\qquad O=PV\.\(14\)At sequence length 4096, the evaluator computes this normalized sigmoid with D3 forward and a coupled, factored D4 derivative surrogate; a runtime guard enforces that fitted length\. Complete\-attention forward measures 1\.074x and backward measures 1\.001x\.

#### 5\.2\.4B4: routed expert SwiGLU

B4 applies D3 SiLU inside routed\-expert SwiGLU\. The router implementation remains native; only the activation inside the routed\-expert operation changes\. The fused expert operation is

Y=SiLU~​\(X​Wg\)⊙\(X​Wu\)\.Y=\\widetilde\{\\operatorname\{SiLU\}\}\(XW\_\{g\}\)\\odot\(XW\_\{u\}\)\.\(15\)Its backward recomputes the D3 sigmoid approximation from the gate input, forms SiLU and its derivative algebraically, and computes gradients for both projected inputs in one pass\. The one\-kernel activation benchmark measures 1\.003x in forward and 1\.005x in backward\.

##### Auxiliary integrations\.

Further simplification of FA4’s fractional exponential measures 1\.898x for its isolated fractional core and 1\.0008x at the model step; complete fits and timings appear in[sectionB\.1](https://arxiv.org/html/2610.00049#A2.SS1)\. Appendix[B\.2](https://arxiv.org/html/2610.00049#A2.SS2)describes the packed sine/cosine implementation and its integration with rotary position embedding \(RoPE\)\.

### 5\.3Model\-scale results

#### 5\.3\.1Complete\-model timing

How do phase speedups and training\-throughput ratios change at complete\-model scale?[Table5](https://arxiv.org/html/2610.00049#S5.T5)separates forward, backward, GPU\-step, and wall\-clock throughput ratios\.

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\.B1 and B2 retain complete\-step gains in the Llama 3 probes\. B3’s attention\-forward change is smaller at the complete\-model step\. B4 improves both full\-model phases in the 27A4B probe\.

#### 5\.3\.2Same\-checkpoint open\-weight evaluation

How do fixed\-checkpoint task accuracy, word perplexity, and inference throughput change under the degree\-4 substitutions?[Table6](https://arxiv.org/html/2610.00049#S5.T6)compares constrained and same\-form Sollya coefficients through the same degree and integration site\.

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→\\rightarrowvariant and is rounded to three decimals\. Throughput ratios are variant divided by native, are rounded to three decimals, and favor the polynomial above one\.Qwen3 and Kimi improve prefill by 5\.1–5\.6% and 8\.4%, respectively; their task\-mean changes span−0\.10\-0\.10to\+0\.04\+0\.04percentage points and their WikiText\-2 perplexity changes by at most 0\.0071\. Qwen2\.5 improves prefill by 3\.5–3\.6% while its word perplexity rises from 8\.734 to 197\.677 and 298,700\.735\. GLM prefill measures 0\.969x and 0\.966x, and GPT\-OSS measures 1\.001x and 0\.999x\. The Qwen2\.5 rows also show that similar kernel throughput across coefficient sources can accompany sharply different model outputs\.

#### 5\.3\.3Paired pre\-training to 100 billion tokens

How do the training\-loss trajectories change when the polynomial is used throughout pre\-training?[Table7](https://arxiv.org/html/2610.00049#S5.T7)reports the final EWM value at the last token count shared by each native/polynomial pair, and[figure3](https://arxiv.org/html/2610.00049#S5.F3)shows the complete trajectories\.

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\.At the common horizons near 100 billion tokens, B2 and B4 finish within 0\.003 final EWM training\-loss units of their native pairs\. B1 finishes 0\.0786 higher and B3 finishes 0\.1073 lower\.

## 6Discussion and Conclusion

Across B1–B4, the phase\-level measurements differ by integration site\. Dense SiLU measures its larger gain 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\. The auxiliary fractional exponential measures 1\.898x in isolation and 1\.0008x at the model step after range reduction, data movement, and fusion are included\.

The measurements cover GB200, BF16/FP16 programs, four principal integration sites, model\-dependent same\-checkpoint substitutions, and paired pre\-training trajectories measured with training loss\. Function, component, and model measurements evaluate symmetry, coefficient rounding, backward rules, packed instruction schedules, and the consuming operator as one integrated kernel program\. Additional engineering observations are collected in[appendixD](https://arxiv.org/html/2610.00049#A4), and the auxiliary integrations appear in[appendixB](https://arxiv.org/html/2610.00049#A2)\.

As matrix, special\-function, and memory pipelines change across hardware generations, this workflow locates transcendental work that remains exposed in modern generative modelling\.

## Acknowledgements

The author thanks Tom Cashman, Paul Balança, Douglas Orr, Carlo Luschi, and Andrew Fitzgibbon for discussions and feedback around this work\.

## Generative\-AI Disclosure

Generative AI tools assisted with code development, experiment orchestration, analysis, figure generation, and drafting and editing this report\. The author reviewed and verified the code, experiments, numerical results, figures, citations, and final text and takes responsibility for the content\.

## References

- \[1\]K\. Alexandridis, V\. Titopoulos, and G\. Dimitrakopoulos\(2025\)Low\-cost flashattention with fused exponential and multiplication hardware operators\.External Links:2505\.14314,[Link](https://arxiv.org/abs/2505.14314)Cited by:[§2\.3](https://arxiv.org/html/2610.00049#S2.SS3.p1.1)\.
- \[2\]E\. Alvarez, V\. Mehta, and F\. Ghodsian\(2026\)Inside NVIDIA Rubin GPU architecture: powering the era of agentic AI\.Note:NVIDIA Technical BlogExternal Links:[Link](https://developer.nvidia.com/blog/inside-nvidia-rubin-gpu-architecture-powering-the-era-of-agentic-ai/)Cited by:[§2\.1](https://arxiv.org/html/2610.00049#S2.SS1.p2.1)\.
- \[3\]K\. Aubrey and N\. Stam\(2025\)Inside NVIDIA Blackwell Ultra: the chip powering the AI factory era\.Note:NVIDIA Technical BlogExternal Links:[Link](https://developer.nvidia.com/blog/inside-nvidia-blackwell-ultra-the-chip-powering-the-ai-factory-era/)Cited by:[§B\.1](https://arxiv.org/html/2610.00049#A2.SS1.SSS0.Px5.p2.1),[§2\.1](https://arxiv.org/html/2610.00049#S2.SS1.p2.1)\.
- \[4\]Y\. Bisk, R\. Zellers, R\. L\. Bras, J\. Gao, and Y\. Choi\(2020\)PIQA: reasoning about physical commonsense in natural language\.InProceedings of the AAAI Conference on Artificial Intelligence,External Links:1911\.11641,[Link](https://arxiv.org/abs/1911.11641)Cited by:[§4\.4](https://arxiv.org/html/2610.00049#S4.SS4.p1.1)\.
- \[5\]S\. Chevillard, M\. Joldes, and C\. Lauter\(2010\)Sollya: an environment for the development of numerical codes\.InMathematical Software – ICMS 2010,pp\.28–31\.External Links:[Link](https://www.sollya.org/)Cited by:[§2\.2](https://arxiv.org/html/2610.00049#S2.SS2.p1.1),[§3\.2](https://arxiv.org/html/2610.00049#S3.SS2.p1.1),[§4\.1](https://arxiv.org/html/2610.00049#S4.SS1.p1.1)\.
- \[6\]P\. Clark, I\. Cowhey, O\. Etzioni, T\. Khot, A\. Sabharwal, C\. Schoenick, and O\. Tafjord\(2018\)Think you have solved question answering? try ARC, the AI2 reasoning challenge\.External Links:1803\.05457,[Link](https://arxiv.org/abs/1803.05457)Cited by:[§4\.4](https://arxiv.org/html/2610.00049#S4.SS4.p1.1)\.
- \[7\]K\. Cobbe, V\. Kosaraju, M\. Bavarian, M\. Chen, H\. Jun, L\. Kaiser, M\. Plappert, J\. Tworek, J\. Hilton, R\. Nakano, C\. Hesse, and J\. Schulman\(2021\)Training verifiers to solve math word problems\.External Links:2110\.14168,[Link](https://arxiv.org/abs/2110.14168)Cited by:[§4\.4](https://arxiv.org/html/2610.00049#S4.SS4.p1.1)\.
- \[8\]\(2026\)CUDA Binary Utilities 13\.0\.NVIDIA Corporation\.External Links:[Link](https://docs.nvidia.com/cuda/archive/13.0.2/cuda-binary-utilities/index.html)Cited by:[§B\.2](https://arxiv.org/html/2610.00049#A2.SS2.p3.1)\.
- \[9\]DeepSeek\-AI, A\. Liu, B\. Feng, B\. Xue, B\. Wang, B\. Wu, C\. Lu, C\. Zhao, C\. Deng, C\. Zhang, C\. Ruan, D\. Dai,et al\.\(2024\)DeepSeek\-v3 technical report\.External Links:2412\.19437,[Link](https://arxiv.org/abs/2412.19437)Cited by:[§2\.3](https://arxiv.org/html/2610.00049#S2.SS3.p2.2)\.
- \[10\]Gemma Team, M\. Riviere, S\. Pathak, P\. G\. Sessa, C\. Hardin, S\. Bhupatiraju, L\. Hussenot, T\. Mesnard, B\. Shahriari, A\. Rame, J\. Ferret, P\. Liu,et al\.\(2024\)Gemma 2: improving open language models at a practical size\.External Links:2408\.00118,[Link](https://arxiv.org/abs/2408.00118)Cited by:[§2\.3](https://arxiv.org/html/2610.00049#S2.SS3.p2.2)\.
- \[11\]E\. Georganas, D\. Kalamkar, S\. Avancha, M\. Adelman, D\. Aggarwal, C\. Anderson, A\. Breuer, J\. Bruestle, N\. Chaudhary, A\. Kundu, D\. Kutnick, F\. Laub, V\. Md, S\. Misra, R\. Mohanty, H\. Pabst, B\. Retford, B\. Ziv, and A\. Heinecke\(2022\)Tensor processing primitives: a programming abstraction for efficiency and portability in deep learning and hpc workloads\.Frontiers in Applied Mathematics and Statistics8\.External Links:[Document](https://dx.doi.org/10.3389/fams.2022.826269),[Link](https://doi.org/10.3389/fams.2022.826269)Cited by:[§2\.2](https://arxiv.org/html/2610.00049#S2.SS2.p3.1)\.
- \[12\]D\. Hendrycks, C\. Burns, S\. Basart, A\. Zou, M\. Mazeika, D\. Song, and J\. Steinhardt\(2021\)Measuring massive multitask language understanding\.InInternational Conference on Learning Representations,External Links:2009\.03300,[Link](https://arxiv.org/abs/2009.03300)Cited by:[§4\.4](https://arxiv.org/html/2610.00049#S4.SS4.p1.1)\.
- \[13\]A\. Kundu, A\. Heinecke, D\. Kalamkar, S\. Srinivasan, E\. C\. Qin, N\. K\. Mellempudi, D\. Das, K\. Banerjee, B\. Kaul, and P\. Dubey\(2019\)K\-TanH: efficient tanh for deep learning\.External Links:1909\.07729,[Link](https://arxiv.org/abs/1909.07729)Cited by:[§2\.2](https://arxiv.org/html/2610.00049#S2.SS2.p3.1)\.
- \[14\]J\. Li, A\. Zhurkevich, V\. Agarwalla, S\. Lee, and R\. Anders\(2026\)Making softmax more efficient with NVIDIA Blackwell Ultra\.Note:NVIDIA Technical BlogExternal Links:[Link](https://developer.nvidia.com/blog/making-softmax-more-efficient-with-nvidia-blackwell-ultra/)Cited by:[§B\.1](https://arxiv.org/html/2610.00049#A2.SS1.SSS0.Px5.p2.1),[§2\.1](https://arxiv.org/html/2610.00049#S2.SS1.p2.1)\.
- \[15\]J\. P\. Lim and S\. Nagarakatte\(2021\)RLIBM\-32: high performance correctly rounded math libraries for 32\-bit floating point representations\.External Links:2104\.04043,[Link](https://arxiv.org/abs/2104.04043)Cited by:[§2\.2](https://arxiv.org/html/2610.00049#S2.SS2.p1.1),[§3\.2](https://arxiv.org/html/2610.00049#S3.SS2.p1.1)\.
- \[16\]S\. Lin, J\. Hilton, and O\. Evans\(2022\)TruthfulQA: measuring how models mimic human falsehoods\.InProceedings of the 60th Annual Meeting of the Association for Computational Linguistics,External Links:2109\.07958,[Link](https://arxiv.org/abs/2109.07958)Cited by:[§4\.4](https://arxiv.org/html/2610.00049#S4.SS4.p1.1)\.
- \[17\]S\. Merity, C\. Xiong, J\. Bradbury, and R\. Socher\(2017\)Pointer sentinel mixture models\.InInternational Conference on Learning Representations,External Links:1609\.07843,[Link](https://arxiv.org/abs/1609.07843)Cited by:[§4\.4](https://arxiv.org/html/2610.00049#S4.SS4.p1.1)\.
- \[18\]\(2026\)PTX ISA 9\.0\.NVIDIA Corporation\.External Links:[Link](https://docs.nvidia.com/cuda/archive/13.0.2/parallel-thread-execution/index.html)Cited by:[§3\.4](https://arxiv.org/html/2610.00049#S3.SS4.p3.1)\.
- \[19\]J\. Ramapuram, F\. Danieli, E\. Dhekane, F\. Weers, D\. Busbridge, P\. Ablin, T\. Likhomanenko, J\. Digani, Z\. Gu, A\. Shidani, and R\. Webb\(2024\)Theory, analysis, and best practices for sigmoid self\-attention\.External Links:2409\.04431,[Link](https://arxiv.org/abs/2409.04431)Cited by:[§2\.3](https://arxiv.org/html/2610.00049#S2.SS3.p1.1)\.
- \[20\]E\. Real, Y\. Chen, M\. Rossini, C\. de Souza, M\. Garg, A\. Verghese, M\. Firsching, Q\. V\. Le, E\. D\. Cubuk, and D\. H\. Park\(2023\)AutoNumerics\-zero: automated discovery of state\-of\-the\-art mathematical functions\.External Links:2312\.08472,[Link](https://arxiv.org/abs/2312.08472)Cited by:[§2\.2](https://arxiv.org/html/2610.00049#S2.SS2.p2.1),[§3\.2](https://arxiv.org/html/2610.00049#S3.SS2.p1.1)\.
- \[21\]K\. Sakaguchi, R\. L\. Bras, C\. Bhagavatula, and Y\. Choi\(2020\)WinoGrande: an adversarial winograd schema challenge at scale\.InProceedings of the AAAI Conference on Artificial Intelligence,External Links:1907\.10641,[Link](https://arxiv.org/abs/1907.10641)Cited by:[§4\.4](https://arxiv.org/html/2610.00049#S4.SS4.p1.1)\.
- \[22\]N\. Shazeer\(2020\)GLU variants improve transformer\.External Links:2002\.05202,[Link](https://arxiv.org/abs/2002.05202)Cited by:[§2\.3](https://arxiv.org/html/2610.00049#S2.SS3.p2.2)\.
- \[23\]J\. Yu, X\. Han, X\. Hu, C\. Xu, Z\. Jiang, and D\. Yang\(2026\)NLI: non\-uniform linear interpolation approximation of nonlinear operations for efficient LLM inference\.External Links:2602\.02988,[Document](https://dx.doi.org/10.48550/arXiv.2602.02988),[Link](https://arxiv.org/abs/2602.02988)Cited by:[§2\.2](https://arxiv.org/html/2610.00049#S2.SS2.p4.1)\.
- \[24\]T\. Zadouri, M\. Hoehnerbach, J\. Shah, T\. Liu, V\. Thakkar, and T\. Dao\(2026\)FlashAttention\-4: algorithm and kernel pipelining co\-design for asymmetric hardware scaling\.External Links:2603\.05451,[Document](https://dx.doi.org/10.48550/arXiv.2603.05451),[Link](https://arxiv.org/abs/2603.05451)Cited by:[§1](https://arxiv.org/html/2610.00049#S1.p1.1),[§2\.1](https://arxiv.org/html/2610.00049#S2.SS1.p2.1)\.
- \[25\]R\. Zellers, A\. Holtzman, Y\. Bisk, A\. Farhadi, and Y\. Choi\(2019\)HellaSwag: can a machine really finish your sentence?\.InProceedings of the 57th Annual Meeting of the Association for Computational Linguistics,External Links:1905\.07830,[Link](https://arxiv.org/abs/1905.07830)Cited by:[§4\.4](https://arxiv.org/html/2610.00049#S4.SS4.p1.1)\.

## Appendix AExtended Numerical Results

The complete D3–D6 sweeps below replay BF16\-rounded coefficients in the runtime evaluation order at BF16\-representable inputs\. Forward and derivative fits are evaluated separately\.

\(a\)Sigmoid forward\.
\(b\)Tanh forward\.
\(c\)SiLU forward\.

Figure 4:Forward numerical accuracy for D3–D6 BF16 target\-precision programs\. Absolute error is measured against the FP32 reference at BF16\-representable inputs\.\(a\)Sigmoid derivative\.
\(b\)Tanh derivative\.
\(c\)SiLU derivative\.

Figure 5:Backward numerical accuracy for D3–D6 BF16 target\-precision programs\. Absolute error is measured against the FP32 reference at BF16\-representable inputs\.
## Appendix BAuxiliary Integration Results

Two auxiliary experiments pair inner\-arithmetic measurements with the surrounding FA4 and rotary\-position\-embedding operations\.

### B\.1Furtherexp2simplification in FlashAttention\-4

FA4 decomposesxxinto an integer and a fractional component\. Its software path writes

x=n\+f,n=⌊x⌋,f∈\[0,1\),2x=2n​2f,x=n\+f,\\qquad n=\\lfloor x\\rfloor,\\quad f\\in\[0,1\),\\qquad 2^\{x\}=2^\{n\}2^\{f\},\(16\)then approximates only2f2^\{f\}\. The integer scale2n2^\{n\}is reconstructed from the exponent bits\. All comparisons below retain FA4’s FP32 lower clamp, range reduction, exponent reconstruction, masking rules, and mixed SFU/software routing\. Only the fractional evaluator changes within the existing attention contract\.

##### Fits and invariants\.

Online softmax repeatedly rescales its running row state when the row maximum changes\. A fractional approximation must therefore satisfyp⁡\(0\)=1p\(0\)=1,p⁡\(1\)=2p\(1\)=2, and continuity between adjacent exponent intervals\. These conditions leave onlyp1​\(f\)=1\+fp\_\{1\}\(f\)=1\+ffor a single line\. For a two\-piece piecewise\-linear \(PWL2\) approximation, we use the selection\-free hinge form

pPWL2​\(f\)=1\+a​f\+2​\(1−a\)​max⁡\(f−12,0\),a=0\.7980513010,p\_\{\\mathrm\{PWL2\}\}\(f\)=1\+af\+2\(1\-a\)\\max\(f\-\\tfrac\{1\}\{2\},0\),\\qquad a=0\.7980513010,\(17\)and for the quadratic control we use

pD2​\(f\)=1\+a​f\+\(1−a\)​f2,a=0\.6602339719\.p\_\{\\mathrm\{D2\}\}\(f\)=1\+af\+\(1\-a\)f^\{2\},\\qquad a=0\.6602339719\.\(18\)The parameters in[equations17](https://arxiv.org/html/2610.00049#A2.E17)and[18](https://arxiv.org/html/2610.00049#A2.E18)minimize sampled maximum relative error on a1,000,0011\{,\}000\{,\}001\-point grid\. We also solve the same one\-parameter families by least squares to expose the objective trade\-off; least squares givesa=0\.7852775264a=0\.7852775264for PWL2 anda=0\.6563658611a=0\.6563658611for D2\. For comparison, FA4’s existing D3 is

pD3​\(f\)=1\+0\.6951461434​f\+0\.2275643945​f2\+0\.0771190897​f3\.p\_\{\\mathrm\{D3\}\}\(f\)=1\+0\.6951461434f\+0\.2275643945f^\{2\}\+0\.0771190897f^\{3\}\.\(19\)[Table8](https://arxiv.org/html/2610.00049#A2.T8)compares both objectives with FA4’s existing Sollya/minimax D3 polynomial\. The D2 sampled\-minimax fit is more accurate than PWL2 despite having the same two\-FMA Horner depth because it needs no hinge\.

Table 8:Fractional2f2^\{f\}fits onf∈\[0,1\]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\.The packed FP16 PWL2 coefficients area=0\.7978515625a=0\.7978515625and2​\(1−a\)=0\.403808593752\(1\-a\)=0\.40380859375\. They preserve both endpoints exactly after rounding and have maximum relative error1\.132×10−21\.132\\times 10^\{\-2\}\. The handwritten core is 11 inline PTX instructions, including twofma\.rn\.f16x2instructions\. The D1 core needs only one packed FMA, but its fit error is much larger and it is no faster than D2 in causal FA4\. The BF16 PWL2 core is also slower and less accurate than FP16\.

For online softmax, endpoint constraints are part of the approximation contract:p⁡\(0\)=1p\(0\)=1preserves the identity update andp⁡\(1\)=2p\(1\)=2preserves an integer exponent step\. The accepted PWL2 comparison therefore includes only fits that satisfy both constraints after target\-precision rounding\. Pointwise error and endpoint invariants are both carried through the recurrent state transition\.

##### Isolated evaluator\.

[Table9](https://arxiv.org/html/2610.00049#A2.T9)separates fractional\-core timing from the complete exponential\. The PWL2 fractional timing uses the reported coefficient set; the constrained fit has the same instruction shape and latency\. Repeated register\-resident arithmetic exposes the shorter fractional core, while L2\-sized and HBM\-sized working sets and the shared FP32 reduction and reconstruction change the measured balance\. PWL2 replaces the saved FMA with hinge subtraction and maximum instructions\.

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\.For the FA4 accuracy comparison, tensor\-relativeL2L\_\{2\}error is∥yvariant−ySFU∥2/max⁡\(∥ySFU∥2,10−30\)\\lVert y\_\{\\mathrm\{variant\}\}\-y\_\{\\mathrm\{SFU\}\}\\rVert\_\{2\}/\\max\(\\lVert y\_\{\\mathrm\{SFU\}\}\\rVert\_\{2\},10^\{\-30\}\)\. We report the maximum over the attention output and the query, key, and value input gradients\.

Table 10:Maximum tensor\-relativeL2L\_\{2\}error against all\-SFU FA4 on GB200\.
##### FA4 forward and backward\.

The integrated target is FA4, which deliberately interleaves a minority of software evaluations with SFU work to balance two execution pipelines\. The integrated rows of[table9](https://arxiv.org/html/2610.00049#A2.T9)compare D3, PWL2, and D2 under one routing contract\. The causal Llama\-8B GQA geometry isB=1,S=4096,Hq=32,Hk​v=8,D=128B=1,S=4096,H\_\{q\}=32,H\_\{kv\}=8,D=128, whereBBis batch size,SSis sequence length,HqH\_\{q\}is the number of query heads,Hk​vH\_\{kv\}is the number of key/value heads, andDDis head dimension\. Timings use a requested 1300 MHz streaming\-multiprocessor clock \(1305 MHz observed\), and speedups are medians of 15 paired per\-round ratios with 100 iterations per round\. The matched rows route 33\.3% of eligible forward and 12\.5% of backward lane pairs through software\. The all\-SFU reference measures 0\.2749 ms forward and 0\.9174 ms forward plus backward; matched D3, PWL2, and D2 measure 0\.2715/0\.9169, 0\.2723/0\.9127, and 0\.2714/0\.9113 ms, respectively\.

The compared paths differ in four relevant ways\. FA4’s D3 core uses three packed FP32 FMAs and runs concurrently with SFU work rather than replacing a serial native exponential\. PWL2 replaces one FMA with a hinge subtraction and maximum, and its FP16 path adds conversions\. Every software variant retains clamp, floor emulation, fractional extraction, and exponent\-bit reconstruction\. Only a minority of eligible lane pairs use software: final fragments and causal masked tiles remain on SFU, while query–key \(QK\) and probability–value \(PV\) matrix multiplication, online reductions, and synchronization are unchanged\.

##### 8B model probes\.

Opposite\-order 80\-step Llama\-8B probes use the same sequence length and clock\. PWL2 forward\-plus\-backward measures 1\.0025x when the SFU control runs first and 0\.9990x in reverse order, for a 1\.0008x geometric mean\. D2 measures 0\.9904x and 1\.0158x\. The opposed orders span a wider range than their geometric means in[table9](https://arxiv.org/html/2610.00049#A2.T9)\. Final candidate\-minus\-control loss deltas range from \-0\.323 to \-0\.050 after 80 optimizer steps\.

##### Exposure analysis for further exp2 degree reduction\.

The bottleneck can be quantified with Amdahl’s law\. If the isolated fractional core accelerates bys=1\.898s=1\.898but matched FA4 forward accelerates by onlyS=1\.010S=1\.010, the effective serial share exposed to that optimization is

p=1−S−11−s−1≈0\.021\.p=\\frac\{1\-S^\{\-1\}\}\{1\-s^\{\-1\}\}\\approx 0\.021\.\(20\)The measured saving is 2\.6μ\\mus in FA4 forward and 4\.7μ\\mus in forward plus backward\. Even multiplying the latter by all 32 Llama layers gives approximately 0\.15 ms against a roughly 275 ms optimizer step, about 0\.05% of step time, consistent with the complete\-step timing\.

Blackwell Ultra doublesMUFU\.EX2throughput for key attention instructions and advertises up to 2x attention\-layer compute over Blackwell\[[14](https://arxiv.org/html/2610.00049#bib.bib13),[3](https://arxiv.org/html/2610.00049#bib.bib14)\]\.

### B\.2Paired sine/cosine and fused RoPE

Rotary position embedding \(RoPE\) supplies a paired\-function experiment because sine is odd, cosine is even, and both share one periodic reduction\. For angular frequencyωi\\omega\_\{i\}, the host stores the unsigned Q0\.32 turn incrementΔi=round⁡\(ωi​232/2​π\)\\Delta\_\{i\}=\\operatorname\{round\}\(\\omega\_\{i\}2^\{32\}/2\\pi\)in a 32\-bit word\. At positionpp, the kernel forms\(p​Δi\)mod232\(p\\Delta\_\{i\}\)\\bmod 2^\{32\}, reduces it to a signed Q15 half\-turn coordinate, and evaluates a packed FP16 D5 sine and D6 cosine with shared squared input\. Each thread processes four frequencies as two packed pairs\. The fitting, numerical\-check, repeated\-evaluator, and fused\-RoPE drivers are undersrc/sfu\_repro/rope/in the standalone reproduction repository\. The retained timing summaries do not identify a final quantitative comparison, so this appendix reports the implementation and integration boundary\.

The fused experiment reads contiguous FP16 Q in\[B,S,Hq,D\]\[B,S,H\_\{q\},D\]and K in\[B,S,Hk​v,D\]\[B,S,H\_\{kv\},D\], whereBBis batch size,SSis sequence length,HqH\_\{q\}is the number of query heads,Hk​vH\_\{kv\}is the number of key/value heads, andDDis head dimension\. It computes four sine/cosine pairs once per thread, reuses them across the corresponding heads, applies the packed complex rotation, and writes the transposed layouts\[B,Hq,S,D\]\[B,H\_\{q\},S,D\]and\[B,Hk​v,S,D\]\[B,H\_\{kv\},S,D\]\. Cached\-table, on\-the\-fly native\-SFU, and on\-the\-fly polynomial variants share the same rotation and memory instructions\.

Inspection with NVIDIA’s binary utilities can distinguish the polynomial instruction path from SFU, FP32\-FMA, or packed\-BF16 fallbacks\[[8](https://arxiv.org/html/2610.00049#bib.bib17)\]\. Query/key movement, packed complex rotation, transpose, and store remain shared across variants\.

## Appendix CPacked PTX Listings

The following listings are emitted by the B3 inline\-PTX generator inflash\-attention/flash\_attn/cute/handwritten\_spline\_ptx\.py\. These are generator\-emitted source\-PTX listings; reported counts cover executable instructions and exclude register declarations\.

1\{

2\.reg\.b32%fa4\_input,%fa4\_mask\_lo,%fa4\_mask\_hi,

3%fa4\_mask,%fa4\_clamped,%fa4\_value,%fa4\_active;

4cvt\.rn\.bf16x2\.f32%fa4\_input,$2,$1;

5mov\.b32%fa4\_mask\_lo,0x3fa03fa0;

6set\.ge\.u32\.bf16x2%fa4\_mask,%fa4\_input,%fa4\_mask\_lo;

7mov\.b32%fa4\_mask\_lo,0x40c040c0;

8min\.bf16x2%fa4\_clamped,%fa4\_input,%fa4\_mask\_lo;

9mov\.b32%fa4\_mask\_lo,0xc0c0c0c0;

10max\.bf16x2%fa4\_clamped,%fa4\_clamped,%fa4\_mask\_lo;

11lop3\.b32%fa4\_value,%fa4\_mask,

120x3a093a09,0x37e837e8,0xca;

13lop3\.b32%fa4\_mask\_hi,%fa4\_mask,

140xbb11bb11,0x39093909,0xca;

15fma\.rn\.bf16x2%fa4\_value,%fa4\_clamped,

16%fa4\_value,%fa4\_mask\_hi;

17lop3\.b32%fa4\_mask\_hi,%fa4\_mask,

180x3b8c3b8c,0x39873987,0xca;

19fma\.rn\.bf16x2%fa4\_value,%fa4\_clamped,

20%fa4\_value,%fa4\_mask\_hi;

21lop3\.b32%fa4\_mask\_hi,%fa4\_mask,

220xbb09bb09,0x39803980,0xca;

23fma\.rn\.bf16x2%fa4\_value,%fa4\_clamped,

24%fa4\_value,%fa4\_mask\_hi;

25mov\.b32%fa4\_mask\_lo,0;

26max\.bf16x2%fa4\_value,%fa4\_value,%fa4\_mask\_lo;

27mov\.b32$0,%fa4\_value;

28\}

Listing 1:Source PTX for the fixed\-length \(n=4096n=4096\) B3 direct\-D3 FlashSigmoid forward evaluator\. Two FP32 scores are rounded into one BF16x2 register; packed comparison andlop3\.b32select one of two coefficient rows, and three packed FMAs evaluate the cubic\. Hex immediates contain duplicated BF16 values after1/40961/4096coefficient scaling\.1\{

2\.reg\.b32%fa4\_input,%fa4\_mask\_lo,%fa4\_mask\_hi,

3%fa4\_mask,%fa4\_clamped,%fa4\_value,%fa4\_grad;

4cvt\.rn\.bf16x2\.f32%fa4\_input,$3,$2;

5mov\.b32%fa4\_mask\_lo,0x3fa03fa0;

6set\.ge\.u32\.bf16x2%fa4\_mask,%fa4\_input,%fa4\_mask\_lo;

7mov\.b32%fa4\_mask\_lo,0x40c040c0;

8min\.bf16x2%fa4\_clamped,%fa4\_input,%fa4\_mask\_lo;

9mov\.b32%fa4\_mask\_lo,0xc0c0c0c0;

10max\.bf16x2%fa4\_clamped,%fa4\_clamped,%fa4\_mask\_lo;

11lop3\.b32%fa4\_value,%fa4\_mask,

120x3a093a09,0x37e837e8,0xca;

13lop3\.b32%fa4\_mask\_hi,%fa4\_mask,

140xbb11bb11,0x39093909,0xca;

15fma\.rn\.bf16x2%fa4\_value,%fa4\_clamped,

16%fa4\_value,%fa4\_mask\_hi;

17lop3\.b32%fa4\_mask\_hi,%fa4\_mask,

180x3b8c3b8c,0x39873987,0xca;

19fma\.rn\.bf16x2%fa4\_value,%fa4\_clamped,

20%fa4\_value,%fa4\_mask\_hi;

21lop3\.b32%fa4\_mask\_hi,%fa4\_mask,

220xbb09bb09,0x39803980,0xca;

23fma\.rn\.bf16x2%fa4\_value,%fa4\_clamped,

24%fa4\_value,%fa4\_mask\_hi;

25mov\.b32%fa4\_mask\_lo,0;

26max\.bf16x2%fa4\_value,%fa4\_value,%fa4\_mask\_lo;

27lop3\.b32%fa4\_mask\_hi,%fa4\_mask,

280xbb38bb38,0xbc10bc10,0xca;

29lop3\.b32%fa4\_mask\_lo,%fa4\_mask,

300x3f813f81,0x3f803f80,0xca;

31fma\.rn\.bf16x2%fa4\_grad,%fa4\_clamped,

32%fa4\_mask\_hi,%fa4\_mask\_lo;

33mul\.rn\.bf16x2%fa4\_grad,%fa4\_value,%fa4\_grad;

34mov\.b32$0,%fa4\_value;

35mov\.b32$1,%fa4\_grad;

36\}

Listing 2:Source PTX for the coupled B3 direct\-D3 probability and factored\-D4 derivative evaluator\. Packing, region selection, clamping, and the complete three\-FMA D3 chain are shared\. The derivative tail selects a region\-specific linear factor, evaluates it with one packed FMA, and multiplies it by the D3 value\.
## Appendix DLessons and Boundaries

1. 1\.Generated instruction schedules determine latency\.High\-level CuTe expressions sometimes lower to long instruction sequences; the workflow inspects the generated PTX\.
2. 2\.Backward has its own error profile\.Direct backward fitting is used when the differentiated forward approximation exceeds the training error budget\.
3. 3\.Algorithmic invariants accompany pointwise error\.Online\-softmax approximations preserve interval endpoints alongside the scalar error bound\.
4. 4\.Integration determines exposure\.The fractional exponential core measures 1\.898x before its surrounding FA4 work and 1\.0008x at the model step\. B1 measures 1\.005x in full forward and 1\.038x in backward, producing a 1\.027x complete\-step result\.

## Appendix ECode and Reproducibility

The report builds withmake \-C report; the artifact map and regeneration commands are inreport/README\.md\. The standalone reproduction repository is[https://github\.com/MrHuff/fast\-polynomial\-transcendentals](https://github.com/MrHuff/fast-polynomial-transcendentals)\. It contains the paper\-specific integration layer without the original model\-training framework, plus the fitting, packed\-kernel, open\-weight, and materialized\-evidence paths\. Filenames below without a directory are underreport/data/\.

Function errorautonumerics\_zero/cuda\_benchmarks/analysis\_results/sollya\_device\_bf16\.jsoncontains the constrained and same\-form Sollya comparison\.

Isolated function timingisolated\_function\_speedups\_fp16\.csvcontains the FP16 endpoint ratios in[table3](https://arxiv.org/html/2610.00049#S5.T3), andisolated\_function\_speedups\_fp16\.provenance\.jsonrecords the benchmark protocol and manuscript lineage\. Raw alternating\-round timings and environment metadata are inisolated\_function\_speedups\_fp16\_gb200\.json\. The runnable driver isautonumerics\_zero/experiments/benchmark\_report\_function\_speed\.py\.

Integrated timingb1\_b4\_100\_iteration\_phase\_probes\_gb200\.json,b1\_b4\_full\_model\_100step\_gb200\.json, andb4\_27a4b\_local\_100step\_gb200\.jsoncontain the component and complete\-model probes\.

Paired pre\-trainingb1\_b4\_paired\_loss\_curves\.csvcontains the sampled histories used in[figure3](https://arxiv.org/html/2610.00049#S5.F3);b1\_b4\_paired\_loss\_curves\.provenance\.jsonrecords the retained run fragments\. Complete launch manifests for initialization and data order are not part of the artifact set\.

Open\-weight evaluationsollya\_eval\_results\.csvcontains[table6](https://arxiv.org/html/2610.00049#S5.T6);open\_weight\_eval\_protocol\.jsonrecords the source\-derived launch configuration and runtime\-attestation status\.

ImplementationsIn the standalone repository, fitting code is underautonumerics\_zero/evolution/, packed kernels are underautonumerics\_zero/spline\_ops/, portable B1–B4 plumbing is undersrc/sfu\_repro/, and B3’s inline\-PTX source is in the pinnedflash\-attentionsubmodule\.