exp_key Function

public function exp_key(u) result(e)

-log(u) for u in [2**-53, 1], using only IEEE + - * /. Not a general-purpose log.

The construction. Split u = m * 2**k with m in [1/sqrt(2), sqrt(2)), so log(u) = log(m) + k*log(2). Then with f = (m-1)/(m+1), at most 0.1716 in magnitude, log(m) = 2*atanh(f) = 2f*(1 + f**2/3 + f**4/5 + ...). Twelve terms leave a truncation error near 6e-21, four orders below one ulp of the result -- and the series stays relatively accurate as m approaches 1, where f -> 0 and the answer is essentially 2f, which a naive log(1+x) would reach only through cancellation. Measured worst error against a real128 reference over every exponent a key can carry: 2.0 ulp (tools/check_exp_key.f90 --accuracy), not the 1 ulp originally aimed at; the floor is f itself, since m+1 rounds and the division rounds again. Two ulp is 2e-16 relative, twelve orders below anything a distributional test could see, so closing it would buy nothing and would cost the two-part f that fdlibm needs for it.

volatile is what makes this reproducible, and it is the whole point of the procedure. A Horner step is a*b + c, the exact shape a compiler fuses into an FMA -- which rounds once where IEEE rounds twice, and so answers differently. Storing each product through a volatile forces it to be rounded to real64 before the add, which no compiler may skip. This is not hypothetical and it is not exotic: gfortran -O3 -march=native changes the fingerprint of the unbarriered form, and -march=native is an ordinary thing to build with. An earlier design note concluded FMA was harmless on the strength of four builds that happened not to enable it on gfortran; tools/check_exp_key.sh now sweeps for it.

The cost is about 2x on the polynomial and it buys a contract that holds under every IEEE-conforming flag set. Measured over 5M calls at -O3 -funroll-loops -march=native on machine C, best of 3: gfortran 10.58 ns unbarriered -> 21.86 ns (2.07x), flang 12.46 -> 23.80 ns (1.91x). (An earlier figure here, 6.58 -> 16.68 ns / 2.53x, was the noinline barrier on machine B; the shape of the conclusion is unchanged.) If the key loop ever dominates a real workload, the way to recover it is a table-driven reduction -- a 16-entry table of log(m_i) brings |f| under 2**-5 and the series down to six terms, halving the barriers -- not deleting the barriers.

-ffast-math / -Ofast / ifx's default -fp-model=fast ARE in scope, and the barrier on the final sum is what brings them in. An earlier version of this note said the opposite -- that no library could survive them -- on the strength of a build that broke. The diagnosis behind that was wrong: the failure was blamed on reciprocal approximation of (m-1)/(m+1), and it is not division at all. Measured on ifx 2026.1.1 at -O2 with the default fast model: -prec-div does NOT fix it, -no-fma does NOT fix it, and -assume protect_parens DOES. The transformation was reassociation, and it had exactly one place to bite -- e = big + (small - logm), where the parentheses were the only thing holding the grouping and fast-math is licensed to ignore them. big is about k * 0.693 and small about k * 1.9e-10, so regrouping to (big + small) - logm annihilates the low-order correction. Wrapping that subtraction in ek_rnd makes the grouping a memory fact rather than a syntactic one, and the fingerprint then holds at every level on both compilers, -Ofast included. It also repaired gfortran -Ofast, which was silently wrong before.

The lesson generalises past this file: a barrier on every PRODUCT is not a barrier on the expression. Each ek_rnd here stops an FMA from spanning an add; none of them stops the adds being re-grouped among themselves. When a value is built as a large term plus a small correction, the grouping IS the algorithm, and it needs its own barrier.

A SUM feeding a PRODUCT is the second shape, and ek_c(0) being exactly 1.0 is what makes the optimiser want it. The last Horner step is poly = t + ek_c(0), which logm then multiplies by f + f. Under reassoc, g * (t + 1.0) distributes to g*t + g -- free, because multiplying by one vanishes -- and that form is a single FMA. The barriers on the two products either side did not help, because the unbarriered add BETWEEN them was the target; ek_rnd(poly) is what closes it. Measured on machine A (arm64): flang -O3 -ffast-math and -Ofast moved the fingerprint to 7108262605301466839, differing on 148 of 13824 inputs by 1 ulp each. Changing ek_c(0) from 1.0 to 0.9 in a scratch copy removes the FMA, which is what identifies the mechanism rather than merely correlating with it.

That one is ARCHITECTURE-gated, so the same compiler and flags can pass on one machine and fail on another. The rewrite fires only where the result becomes one instruction: the same -O3 -ffast-math emits the FMA on aarch64, where FMA is baseline, and on x86-64 under -march=haswell -- and does NOT emit it on baseline x86-64, which has no FMA to contract into. That is why machine C reported clean while machine A failed, and it is the same trap as the original -mfma discovery: a hazard hidden by builds that happen not to enable FMA. Note -ffp-contract=off does not undo it, because flang stamps contract on the IR from -ffast-math and the later flag does not clear it.

NOT pure, and the impurity is the price of the barrier being GUARANTEED. A volatile local is standard Fortran and every compiler must honour it; !GCC$/!DIR$ ATTRIBUTES NOINLINE are directives a compiler is free to ignore. This transform was built on the directives for a while, precisely so that exp_key could stay pure -- and flang 22.1.8 then proved that unsound: it warns on the !DIR$ spelling, ignores the !GCC$ one in silence, inlines the helper and fuses. Measured on machine C (x86-64, i7-10700K), flang -O3 -march=native moved the full fingerprint to -8585622607960331921, differing on 4 of 13824 swept inputs by 1 ulp each. So the directive form did not merely risk losing the barrier; it had already lost it on a compiler this project builds under.

Purity cost nothing to give up and the barrier got FASTER. Nothing needs exp_key to be pure: its only caller is wperm_impl's ordinary serial do loop, and the two other entry points here (exp_key_contract_ok, parquet_debug_exp_key) were never pure either. A stack store/reload also beats an out-of-line call -- measured over 5M calls at -O3 -funroll-loops -march=native on machine C, best of 3:

The transfer(transfer(x, 0_int64), 0.0_real64) row is recorded because it looks like the obvious way to keep pure: it does block the fusion on both compilers, and it is unaffordable on flang, which lowers it through memory per call. Do not re-adopt it.

Inlining is now harmless -- an inlined volatile store and reload is still a store and a reload -- so the directives are gone, and with them flang's -Wignored-directive warning. tools/check_exp_key.sh remains the thing that would notice a lost barrier: it compares the fingerprint across every configuration a compiler supports. Verified on gfortran 15.2.1 (7/7, machines B and C), ifx 2026.1.1 (4/4, machine B) and flang 22.1.8 (6/6, machine C).

exp_key_contract_ok is NOT a second line of defence against this -- see its own doc-comment. It samples 32 inputs and reproduced the frozen value under the diverging flang build, so it did not see the divergence at all. The script is the check that works.

The domain is not general and the guard is the caller's. The race feeds this 1 - u for a uniform u in [0, 1), so the argument lies in [2**-53, 1] and is never zero, never denormal and never above 1. exponent/fraction are exact bit operations on a normal number, so nothing rounds before the polynomial does. u = 1 gives exactly 0.

Arguments

Return Value real(kind=real64)

-log(u), in [0, 36.74]

    function exp_key(u) result(e)
        real(real64), intent(in) :: u   !! a uniform in `[2**-53, 1]`; nothing validates this
        real(real64) :: e               !! `-log(u)`, in `[0, 36.74]`
        real(real64) :: m, f, s, poly, en, big, small, logm
        integer :: k, i

        k = exponent(u)
        m = fraction(u)                 ! in [0.5, 1)
        if (m < ek_sqrt_half) then
            m = m + m
            k = k - 1
        end if
        f = (m - 1.0_real64) / (m + 1.0_real64)
        s = f * f
        poly = ek_c(11)
        do i = 10, 0, -1
            poly = ek_rnd(poly * s) + ek_c(i)
        end do
        logm = ek_rnd((f + f) * ek_rnd(poly))
        en = real(-k, real64)
        big = ek_rnd(en * ek_log2_hi)
        small = ek_rnd(en * ek_log2_lo)
        e = big + ek_rnd(small - logm)
    end function exp_key