Fast math, FMA and precision, one flag at a time

Take x = 1 + 2^-23 and y = -(1 + 2^-22). The exact value of x*x + y is 2^-46. A CPU computing in plain single precision says zero.

A GPU with default compiler settings says 2^-46, exactly right, and it uses one instruction where the CPU uses two.

People meet this as a bug report: "my GPU results don't match my CPU results, single precision, same code." Then they find -fmad=false, the results match, and they conclude that the GPU was wrong. The GPU fused the multiply and the add into one instruction with one rounding, and one rounding is more accurate than two.

This optimization improves accuracy. The rest of the fast-math bundle trades accuracy for speed, so test each setting before you enable the full bundle.

What -use_fast_math actually is

-use_fast_math enables several settings. The nvcc manual says it implies --ftz=true --prec-div=false --prec-sqrt=false --fmad=true (https://docs.nvidia.com/cuda/cuda-compiler-driver-nvcc/index.html , checked 2026-08-30), and on top of those four it swaps a list of CUDA Math API calls for their intrinsic counterparts.

The defaults it overrides are -ftz=false, -prec-div=true and -prec-sqrt=true. -fmad=true is already the default, so that member of the bundle changes nothing.

Piece by piece:

• --ftz=true: single-precision operations flush subnormal results to zero instead of keeping them. Subnormals are the values between the smallest normal float, 2^-126, and zero; they exist so that underflow is gradual.
• --prec-div=false: x / y becomes __fdividef(x, y), 2 ULP instead of correctly rounded, and 1 / x becomes a 1 ULP reciprocal.
• --prec-sqrt=false: sqrtf(x) goes from correctly rounded, 0 ULP, to 1 ULP.
• The substitution list: expf becomes __expf, logf becomes __logf, sinf becomes __sinf, powf becomes __powf, and so on through twelve entries. The full table is in the programming guide's "--use_fast_math Effect" section (https://docs.nvidia.com/cuda/cuda-programming-guide/05-appendices/mathematical-functions.html , checked 2026-08-30).

The intrinsics are cheap because they run on the special function units, which the glossary files under CUDA core, and because their error bounds are looser and grow with the argument. expf is documented at 2 ULP everywhere. __expf is documented at 2 + floor(|1.173 * x|) ULP, so at x = -16 you are allowed about 20 ULP and at x = -80 about 95.

This bound helps predict how the softmax below behaves.

Three of the four changes have a source-level spelling. You can write __expf or __fdividef yourself in one kernel and leave every other kernel precise, and you can write __fmul_rn and __fadd_rn to forbid fusing at one expression. Flush-to-zero is the exception: there is no per-callsite spelling because it is a build setting.

FMA: faster and more accurate, which surprises people

A fused multiply-add computes a * b + c as one instruction with one rounding at the end. The unfused form rounds twice: once after the multiply, once after the add. One rounding step fewer means less error, and one instruction instead of two means less time.

One source of cross-platform mismatch remains: the CPU compiler may not fuse the same expression.

The x*x + y example at the top of the page is the programming guide's own, from its Floating Point section: the product x*x rounds to float and the 2^-46 tail is exactly what rounding throws away, so the unfused form cancels to zero while the fused form keeps it.

Day 47's program measures this with a pair of kernels that differ only in whether the compiler is allowed to fuse:

__global__ void residualFused(const float* __restrict__ x,
                              const float* __restrict__ y,
                              float* __restrict__ out, size_t n) {
    const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
    if (i < n) {
        out[i] = x[i] * x[i] + y[i];
    }
}

// The same value with the contraction taken away, and taken away in a form
// no command-line flag can put back. The programming guide states that
// __fadd_rn and __fmul_rn "map to addition and multiplication operations
// that the compiler never merges into the FFMA or DFMA instructions", so
// this kernel is two roundings under every build in the README, including
// the one with --use_fast_math.
__global__ void residualSeparate(const float* __restrict__ x,
                                 const float* __restrict__ y,
                                 float* __restrict__ out, size_t n) {
    const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
    if (i < n) {
        out[i] = __fadd_rn(__fmul_rn(x[i], x[i]), y[i]);
    }
}

The inputs are built so the two cannot agree: x = 1 + 2^-k and y = -(1 + 2^(1-k)) for k from 12 to 23, so the exact answer is 2^-2k and the unfused form returns zero every time. Twelve bytes and two flops per element makes this kernel memory bound on the measured GPU. Timing should barely separate the two versions.

Here, FMA changes rounding rather than run time. The instruction counts in the shipped Nsight Compute report are what prove the instruction was removed even when the clock cannot see it.

Check what the compiler built

Because -use_fast_math rewrites calls, the program must identify which settings the compiler used. It cannot read its own compile command, so it tests the binary's behavior:

__global__ void buildFingerprint(const float* __restrict__ probe,
                                 int* __restrict__ same,
                                 float* __restrict__ tiny) {
    int expDiff = 0;
    int divDiff = 0;
    for (int p = 0; p + 1 < static_cast<int>(kProbeCount); ++p) {
        const float u = probe[p];
        const float v = probe[p + 1];
        expDiff |= (expf(u) == __expf(u)) ? 0 : 1;
        divDiff |= ((u / v) == __fdividef(u, v)) ? 0 : 1;
    }
    const float x = probe[kProbeCount];
    const float y = probe[kProbeCount + 1];
    const int madDiff = ((x * x + y) == __fadd_rn(__fmul_rn(x, x), y)) ? 0 : 1;

    same[0] = 1 - expDiff;
    same[1] = 1 - divDiff;
    same[2] = 1 - madDiff;
    tiny[0] = expf(probe[kProbeCount + 2]);
}

If expf(u) and __expf(u) return the same bits over five probe values, the compiler already substituted the intrinsic for you. If x*x + y matches the never-fused form, contraction is off. And expf(-90), whose true value is about 8.19e-40, comes back either subnormal or as zero, which identifies the flush-to-zero setting.

Read this block first in every transcript because the other row labels depend on it.

The softmax, four ways

The measurement the curriculum asks for: accuracy and speed on a softmax, 2048 rows of 1024 floats, one block per row, the block-wide max and sum reductions from day 24 in shared memory. One template kernel, four instantiations: expf or __expf, the IEEE divide or __fdividef. No instantiation changes an address, so all four move identical bytes and the times are comparable. Three lines contain the math:

    for (size_t c = tid; c < cols; c += blockDim.x) {
        const float shifted = in[rowBase + c] - rowMax;
        const float e = kFastExp ? __expf(shifted) : expf(shifted);
        out[rowBase + c] = kFastDiv ? __fdividef(e, denom) : (e / denom);
    }

Each variant reports its maximum relative error against a double-precision reference with a compensated sum, next to its mean time over 10 runs after 3 warm-ups, timed with CUDA events the way day 9 established.

Subtracting the row maximum also limits the intrinsic's error. __expf's error bound grows with |x|, and the shift guarantees every argument lies in [-spread, 0], where the spread here is about 16. The same subtraction that keeps the softmax from overflowing is also what keeps __expf cheap to use.

A softmax over wider logits can have more error under the same flag.

The fourth part of the program tests --ftz=true near underflow. It computes 1024 exponentials from exp(-87) down to exp(-107.46), a range that crosses the smallest normal float near -87.34 and the smallest subnormal one near -103.28.

The host classifies every result as normal, subnormal or zero. A build with flush-to-zero has an empty subnormal row.

--prec-sqrt=false goes unmeasured, deliberately. A softmax has no square root, and its effect is documented rather than dramatic: sqrtf moves from 0 ULP to 1. Measuring it wants a kernel built around a norm, which is a different day.

Results

Re-verified without material numerical or behavioral drift on a Tesla T4 with driver 580.173.02 and CUDA 13.0 V13.0.88 on 2026-09-02. All fingerprints, errors, normal/subnormal/zero counts and compiler substitutions reproduced.

The earlier CUDA 13.2 compatibility addendum remains in evidence. The default softmax rows still differed by only about three percent under 13.0, with a different fastest row within that small spread. Existing NCU reports remain the instruction-count evidence; the CUDA 13 PTX, SASS and Nsight Systems artifacts are listed in front matter.

Originally measured on a Tesla T4, driver 595.84, CUDA 12.6 (V12.6.85), four builds of one source file with nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo plus the flag each row names. Captured 2026-09-01 on the project's verification node; all four transcripts are in the page's evidence file.

That GPU's ridge point was 33 flop/byte, which classified the residual kernels as memory bound for this run.

The default build's part 2 provides the main comparison:

exp      divide       max rel err   time (ms)   vs row 1
-------  ----------  ------------  ----------  ---------
expf     IEEE /         7.406e-07      0.1073      1.000
__expf   IEEE /         1.494e-06      0.1110      1.035
expf     __fdividef     7.518e-07      0.1093      1.018
__expf   __fdividef     1.512e-06      0.1097      1.022

And the -use_fast_math build's, where the source spellings stop mattering because the compiler has already made every substitution:

exp      divide       max rel err   time (ms)   vs row 1
-------  ----------  ------------  ----------  ---------
expf     IEEE /         1.512e-06      0.1167      1.000
__expf   IEEE /         1.512e-06      0.1175      1.007
expf     __fdividef     1.512e-06      0.1168      1.001
__expf   __fdividef     1.512e-06      0.1168      1.000

The four errors collapsing to the same 1.512e-06 is the fingerprint argument made by the data: under -use_fast_math all four variants are one kernel.

One held, exactly. In the default build residualFused reports a maximum relative error of exactly 0 and residualSeparate exactly 1.0, at 0.0511 against 0.0510 ms. In the -fmad=false build the fused kernel's error is 1.0 too: same source, one flag, the answer gone. The instruction counts below say why.

The error bound held, but the ratio did not. The expf softmax lands at 7.406e-07 and swapping in __expf moves it to 1.494e-06, a factor of 2.0, not the "roughly ten" implied by the documented ULP bound. The bound 2 + floor(|1.173 * x|) is a ceiling, and on these inputs the intrinsic sits far under its ceiling. The divide swap moves the error by a factor of 1.015.

Both stay under 1e-5.

The speed prediction did not hold. No variant beats another by tens of percent: the default build's spread is 3.5 percent with plain expf fastest, and the fast build's is under 1 percent. This softmax at this shape is memory bound, the flag bought nothing on the clock, and the result does not support a speed gain. The flags changed the error column without improving time.

Four held. Default build census over the 1024-point tail ramp: 17 normal, 832 subnormal, 175 exactly zero, smallest nonzero 1.401298e-45. Under -ftz=true: 17 normal, 0 subnormal, 1007 zero, smallest nonzero 1.195105e-38, the smallest normal float. The 832 subnormals moved into the zero row to the last count, and the fingerprint's expf(-90) line flipped from subnormal to zero in both the -ftz=true and -use_fast_math builds.

The profiler settled the instruction-count prediction with no rounding:

kernel             ffma      fmul      fadd
residualFused      1048576   0         0
residualSeparate   0         1048576   1048576

The profiler evidence, and who can capture it

Timing cannot show one instruction removed from a memory-bound kernel, but Nsight Compute counts instructions directly. The three metrics this page will quote, from the report's Instruction Statistics section: smsp__sass_thread_inst_executed_op_ffma_pred_on.sum, smsp__sass_thread_inst_executed_op_fmul_pred_on.sum and smsp__sass_thread_inst_executed_op_fadd_pred_on.sum. Prediction: residualFused executes 1,048,576 FFMAs and no separate FMUL or FADD for its arithmetic; residualSeparate executes zero FFMAs and 1,048,576 of each of the other two.

Running ncu yourself needs root on a stock driver. Plain ncu returns ERR_NVGPUCTRPERM, because the driver default RmProfilingAdminOnly: 1 restricts performance counters to administrators; this project measured that on its own T4 on 2026-08-29. On day 42's node the fix is sudo:

sudo ncu --set full --clock-control base \
     --kernel-name-base demangled \
     --kernel-name regex:'residualFused' \
     --launch-skip 3 --launch-count 1 \
     --export profile/fast-math-residual-fused \
     ./fast_math

That is one kernel's capture; a single unfiltered --launch-count 1 would hold only the first launch in the program, a warm-up. The README carries the full six-report loop, one per kernel, including the regex escaping that ncu's demangler forces on the four softmaxRows instantiations. The existing reports prove that all four patterns matched the intended kernels; their details exports render the same specializations as softmaxRows<0, 0> through <1, 1> after capture.

The permission capture proves ordinary-user NCU was blocked on the measured node with RmProfilingAdminOnly: 1; it does not settle current free-Colab behavior. That hosted-tier question remains open.

This page ships reports regardless: six per-kernel fast-math-<tag>.ncu-rep files from the default build and six fast-math-fast-<tag>.ncu-rep from the -use_fast_math build, tags residual-fused, residual-separate and softmax-v0 to v3, each with a .details.txt export, at content/profiles/fast-math-and-precision/. Opening a shipped report with ncu --import needs no counters and no GPU. You can complete the exercise from the reports and text transcripts alone.

The compiler output needs no profiler or root and gives more detail here. nvcc -ptx and cuobjdump -sass show the substitutions directly, expf's polynomial sequence in the default build against a short MUFU.EX2-based sequence in the fast one, and the FFMA present in one residual kernel and absent from the other. Day 46 taught the reading skill; the README has the exact commands, and both work on a machine with no GPU.

Run it yourself

One source file produces four builds, and the README lists each command. On Compiler Explorer, select the target architecture you want to inspect and use -O3 -std=c++17 plus that -arch value. Toggle -use_fast_math and compare the build-identification block in the two outputs.

The shipped artifact targeted sm_75 with -O3 -arch=sm_75 -std=c++17. Use a target that matches the code you want to inspect when you repeat the exercise.

Exercise

The per-element divide is the wrong place to spend precision. Change the kernel to compute each row's reciprocal once, 1.0f / denom, and multiply by it in the final loop. Then measure: accuracy and time against the shipped IEEE / and __fdividef rows.

Time: 25 to 40 minutes. Submit: the changed kernel, your accuracy and time rows beside the two shipped ones, and one sentence on which shipped row your reciprocal variant sits closer to and why.

Check: the harness checks that every output is finite and in [0, 1], every row summing to 1 within 1e-3, and the no-intrinsic variant within 1e-4 of the double reference. Your reciprocal variant must pass all of them. Compare where its error lands between the two shipped rows.

const float invDenom = 1.0f / denom;

Pitfalls

Setting -use_fast_math globally and forgetting it. The flag applies to every kernel in the translation unit. The tested kernel may pass, while a later numerically sensitive reduction can inherit flush-to-zero and 2 ULP divides. Prefer the source-level spellings, __expf and __fdividef where you measured them, and keep the flag off.

Flush-to-zero is the one piece with no source spelling, so if you need ftz you need the flag, scoped to one file.

Turning off FMA to make the GPU "match the CPU". -fmad=false makes both platforms equally wrong instead of one of them right. If you need cross-platform agreement, state that need and accept the cost; that trade has a name, floating point determinism, and day 68 is built around it. Do not call it a fix.

Trusting the build line instead of the binary. A flag in one CMake target can turn a "precise" baseline into fast math. Run the build-identification kernel and print what the compiler did.

Using __expf on unshifted arguments. The documented bound grows with |x|. After a max-subtraction the arguments are small and the intrinsic is nearly free accuracy-wise; on raw logits in the tens you lose more precision than you measured.

Assuming the error bound is the error. 2 ULP per operation is a ceiling per call. What accumulates across a 1024-term sum, and what cancellation amplifies, depends on the program. This lesson measures a full softmax against a double reference instead of adding up documentation numbers.

Go deeper

• The mathematical functions appendix, with every ULP bound quoted here and the full --use_fast_math substitution table: https://docs.nvidia.com/cuda/cuda-programming-guide/05-appendices/mathematical-functions.html (checked 2026-08-30)
• The nvcc manual's entries for -ftz, -prec-div, -prec-sqrt, -fmad and -use_fast_math: https://docs.nvidia.com/cuda/cuda-compiler-driver-nvcc/index.html (checked 2026-08-30)
• Whitehead and Fit-Florea, "Precision and Performance: Floating Point and IEEE 754 Compliance for NVIDIA GPUs", the source for the FMA section: https://developer.nvidia.com/sites/default/files/akamai/cuda/files/NVIDIA-CUDA-Floating-Point.pdf (checked 2026-08-30)
• Goldberg, "What Every Computer Scientist Should Know About Floating-Point Arithmetic", for rounding, ULP and cancellation from first principles.

Next

Day 48 stays inside one kernel's cost model but changes the question: instead of making each operation cheaper, it asks what a whole extra kernel launch costs, and fuses a four-kernel chain to find out.