Internals

Double32s.add_dd_(x::Tuple{T,T}, y::Tuple{T,T}) -> Tuple{T,T}

Accurate double-word addition. Branch free; the caller must check the result for non-finiteness.

This is, operation for operation, AccurateDWPlusDW — Algorithm 6 of Joldeş, Muller & Popescu, Tight and rigorous error bounds for basic building blocks of double-word arithmetic, ACM TOMS 44(2), Article 15res, 2017, DOI 10.1145/3121432; the algorithm and its Theorem 3.1 are both on p. 15res:7 — whose relative error bound is

|result − (x + y)| ≤ 3u²/(1 − 4u) · |x + y|,    u = 2⁻²⁴,

valid for every p ≥ 3 (Theorem 3.1), proved in that paper and subsequently machine-checked (Coq/Flocq) and shown asymptotically optimal by Muller & Rideau, Formalization of double-word arithmetic, and comments on "Tight and rigorous error bounds…", ACM TOMS 48(1), Article 9, 2022, DOI 10.1145/3484514, §2.2 Property 2.1 (the algorithm is renumbered Algorithm 4 there). The formalization is generic in the precision, so the bound applies at p = 24, and it covers the two Fast2Sum steps — the algorithm's correctness does not rest on the |a| ≥ |b| precondition.

Validity: the bound assumes no overflow and canonical (nonoverlapping) inputs. It survives into the subnormal range in absolute terms — floating-point addition errors are always exactly representable (Hauser), so every EFT here stays exact — but the relative form of the bound presumes a normal result.

Double32s.add_dd_fast_(x::Tuple{T,T}, y::Tuple{T,T}) -> Tuple{T,T}

Fewer operations and registers than add_dd_, weaker last-bit error bound. For reductions and explicitly requested fast mode.

Double32s.sub_dd_(x::Tuple{T,T}, y::Tuple{T,T}) -> Tuple{T,T}

Accurate double-word subtraction. Its own kernel rather than add_dd_(x, -y), for the signed-zero reason given in design section 7.3.

Error bound: identical to add_dd_ — negation of both words is exact, and two_diff_(a, b) computes the same values as two_sum_(a, -b) apart from the sign of a zero, so the AccurateDWPlusDW bound 3u²/(1-4u) transfers.

Double32s.mul_dd_(x::Tuple{T,T}, y::Tuple{T,T}) -> Tuple{T,T}

Accurate double-word multiplication: exact leading product plus both cross terms plus the low-low term.

fma here is chosen for accuracy, not exactness — each fma retains one rounding a separate multiply-and-add would lose. That is not in tension with the "never muladd in an EFT" rule, which is about muladd's optionality; these are explicit fma calls whose fusion is guaranteed.

Not a proof of correctly rounded Double32 multiplication: residuals of the cross products are not all retained.

Double32s.mul_dd_fast_(x::Tuple{T,T}, y::Tuple{T,T}) -> Tuple{T,T}

Omits the x.lo * y.lo term, whose magnitude is normally of order u² relative to the leading product. Efficient, but it can affect the final double-word bits.

Double32s.div_dd_(x::Tuple{T,T}, y::Tuple{T,T}) -> Tuple{T,T}

One leading quotient and one compensated correction — the division kernel for every policy; see div_dd_refined_ for why the second correction is measured, retained, and not used.

Accuracy: worst observed error 4.00 nominal-48-bit ulps (1.56% of 101,761 lattice-range quotients exceed 1 ulp), same family as Joldeş–Muller–Popescu's DWDivDW2, whose proven bound for their variant is 15u². No formal bound has been derived for this exact kernel; the claim is the measured one, and faithful rounding is not claimed (design section 9's target revised by measurement — see experiment 10's resolution).

The correction is not skipped when y[2] == 0: a one-word divisor such as 3.0f0 still needs it to reach double-word precision. That skip is a tempting "optimization" precisely because it looks like an algebraic identity.

Double32s.div_dd_refined_(x::Tuple{T,T}, y::Tuple{T,T}) -> Tuple{T,T}

div_dd_ plus a second Newton correction.

Measured, and deliberately not selected by any policy. Over 101,761 lattice-range quotients against a 400-bit BigFloat reference (test/core/arithmetic.jl carries the regression):

The second correction roughly doubles the kernel's operation count and register pressure — the costs that set GPU occupancy — and still does not reach faithful rounding. Design section 8.4's contingency ("StrictPolicy may select the refinement if one correction is not faithfully rounded") was written on the assumption that the refinement would close the gap; it does not, so paying 2× in every device kernel buys only a smaller constant. The resolution of design experiment 10 is therefore: one correction for every policy, with the measured bound documented, and this kernel retained under test for any future host-only accuracy policy.

A faithfully-rounded double-word division needs a different algorithm, not a repeated correction; that remains out of scope (design section 25.4 territory).

Double32s.inv_dd_(y::Tuple{T,T}) -> Tuple{T,T}

Reciprocal as a specialized division, saving the high-word divide's operand setup.

Double32s.sqrt_dd_(x::Tuple{T,T}) -> Tuple{T,T}

One Newton correction on the Float32 square root. Valid for positive finite x whose high word is normal; the layer-2 method handles scaling and the exceptional cases.

Double32s.fma_dd_(x, y, z) -> Tuple{T,T}

Compensated x*y + z computing all major pair terms with one final renormalization. Not correctly rounded; the derived bound is absolute, in units of the operand scale, because under near-total cancellation no relative bound is achievable at any cost in this format.

Derived bound. Let u = 2⁻²⁴ and M = max(|xh·yh|, |zh|). For canonical finite inputs with M ∈ [2⁻⁷⁷, 2¹²⁶] (every intermediate then avoids overflow and keeps its two_prod_ residual representable):

|fma_dd_(x, y, z) − (x·y + z)| ≤ 2u²·(1 + 80u)·M

Derivation. Every two_sum_/two_prod_ is exact, so the chain telescopes: s + a4 + (b1+b2+b3+b4) + e2 + e3 + p4 accounts for every term of the exact result except the xl·yl rounding (≤ u³M). Magnitudes: |t| ≤ 2uM(1+u), |e1| ≤ uM, |p2|,|p3| ≤ uM(1+u), |zl| ≤ uM, whence |b1..b4| ≤ (3,4,5,6)u²M·(1+O(u)) and |e2|,|e3|,|p4| ≤ u²M(1+O(u)), giving |low| ≤ 21u²M(1+O(u)) with its six accumulation roundings ≤ 126u³M. The 2Sum (s, a4) is exact with |l| ≤ u|h| ≤ 2uM(1+5u), so the only O(u²) rounding in the whole kernel is the final l + low: ≤ u(|l| + |low|) ≤ 2u²M(1+16u). Total ≤ 2u²M(1+80u). ∎

Measured worst over structured pools: 2.0000 u²·M — the bound is tight to well under one percent. The ordering question of design section 23.9 is thereby closed: with the exact-telescope structure fixed, ordering moves only the O(u³) terms; the 2u² floor is the final-recombination rounding and is irreducible without a three-word accumulator.

The final normalization uses two_sum_ rather than two_hilo_sum_ deliberately — nothing establishes the FastTwoSum residual condition here, and design section 7.2's rule is to pay the three extra operations rather than rely on an unproven one.

Outside the validity window the result degrades gradually with the underflow taper of design section 5.5 (products below 2⁻⁷⁷ lose cross-term residual bits at subnormal granularity); no rescue is attempted because, unlike *, the three-operand scale structure admits no single exact rescaling.

Double32s.mul_guarded(policy, x, y) -> Double32

Cold path for *, reached in exactly four situations — a non-finite operand, a zero operand, an overflowing product, an underflowing-band product — and the scaled retry in the last two is a correctness requirement (design section 7.5), not tuning. canonical_nonfinite(x.hi * y.hi) — correct for +/-, where no rescue exists — would report Inf for representable products and garbage low words for tiny ones.

Trigger ⟹ the shift is safe. Let es = exponent(x.hi) + exponent(y.hi). With finite nonzero operands the kernel's high word satisfies |h| = |x·y|·(1 ± 3.1u), and 2^es ≤ |xhi·yhi| < 2^(es+2), so:

• h non-finite requires |xhi·yhi| ≥ (2^128 − 2^103)/(1+3.1u), hence es ≥ 126; scaling both operands down by 2^-64 puts their exponents in [-1-64, 127-64] = [-65, 63] — both normal — and the scaled product exponent in [-2, 62], where the kernel is exact-residual clean (a nonzero two_prod_ residual is ≥ 2^(e_p - 47) ≥ 2^-50).
• |h| < MUL_TINY = 2^-98 requires es ≤ -97; each operand exponent is then ≤ -97 + 149 = 52, so scaling both up by 2^64 cannot overflow (≤ 2^116), and the scaled product exponent lands in [-170, 31] — fully accurate wherever the true result is representable at all; a true product below 2^-150 correctly becomes zero on the rescale.

The two bands are disjoint (es ≤ -97 vs es ≥ 126), so one signed shift serves both. Dropped low-word bits from operand scaling cost at most 2^-63 relative — an operand's lo can underflow only if |lo/hi| < 2^(-85-es_min), far below u² = 2^-48.

The rescaling is exact and maps the overflow boundary exactly. Scaling the result by 2^128 is exact on both words whenever the hi word stays in range, and the boundary case cannot be misclassified: the kernel's h is fl of the scaled pair value, and at every binade the Float32 rounding boundary (pair ≥ (1 - 2^-25)·2^(e+1) rounds hi up to the next power of two) maps under exact scaling to precisely the overflow threshold 2^128 - 2^103, including the round-to-even tie — both 1.0 and 2^128 have even significands. So ldexp of the finite result returns Inf iff the true product rounds to Inf. On the underflow side the rescale rounds into the subnormal taper; a one-quantum double-rounding tie there is inside the documented taper of design section 5.5.

What the underflow rescue can and cannot deliver. The scaled kernel is full-precision, but a result near 2^-126 cannot hold that precision: at exponent(hi) = -126 the canonical bound is |lo| ≤ ½ulp(hi) = 2^-150, below the least subnormal, so the format admits no low word at all — precision tapers from 48 bits at 2^-102 to 24 bits at 2^-126 (design section 5.5). The guard's obligation is the nearest representable pair, i.e. absolute error ≤ few·u²·|V| + 2^-150-class. Measured over the band [2^-126, 2^-98]: worst error 8.2e-8·|V| (2.3e7 u², 23 bits lost) before this guard; after it, within 4u²·|V| + 2^-149 everywhere — the representation optimum. The regression is in test/core/arithmetic.jl.

Double32s.div_guarded(policy, x, y) -> Double32

Overflow rescue for / (design sections 7.5, 8.0): the quotient's leading term can overflow while the true pair quotient is representable, and before this path existed the hot kernel turned that band into NaN, not merely Inf — q0 = Inf poisons the residual x - q0*y with Inf - Inf. Measured: Double32(2f0^127) / Double32(0.5f0, 2f0^-25) has true quotient floatmax + 2^80 (representable) and returned NaN.

Reachability ⟹ the shift is safe. The hot path reaches here only when its h was non-finite with finite nonzero operands, which requires exponent(x.hi) - exponent(y.hi) ≥ 127:

• that forces exponent(x.hi) ≥ 127 + exponent(y.hi) ≥ 127 - 149 = -22, so kx was zero (the tiny-numerator scaling and this rescue are mutually exclusive), and scaling x down by 2^-64 keeps its high word normal (exponent ≥ -86);
• the scaled leading quotient has exponent in [63, 212]: if it still overflows, the true quotient exceeds 2^211 and Inf is correct saturation; the rescuable band exponent ∈ [127, 129] lands near 2^64, where the kernel is fully accurate;
• dropped numerator low-word bits cost at most 2^-63 relative, far below u²;
• the 2^64 rescale of the result is exact on the hi word up to the overflow boundary, and the boundary maps exactly for the same reason as in mul_guarded — h is fl of the scaled pair, and the round-up-to-next-binade boundary at any binade scales exactly onto the 2^128 - 2^103 overflow threshold, ties included.

Double32s.divide(policy, x::Double32, y::Double32) -> Double32

Layer-2 division. A tiny numerator is scaled up first, and this is not an optimization to be removed.

div_dd_ forms the residual x - q0*y, whose magnitude is about |x| * 2^-48. When exponent(x.hi) drops near -101 that residual falls below 2f0^-149 and is lost — silently, since the result is perfectly finite. The quotient itself may be fully representable while the algorithm computing it is not: measured on random operands, Double32(2.0165723f-36, 8.0f-44) / Double32(4.0026174f-18, …) came back with relative error 6.0e-10, roughly 35 correct bits instead of 48.

The scaling only ever goes up, and that direction is load-bearing. Scaling a large numerator down into a unit binade — the obvious symmetric thing — destroys a subnormal low word, and with it x / one(Double32) === x: Double32(3.0f0, 2f0^-149) scaled down by one binade loses its low word outright. Scaling up is always exact.

Only the numerator needs it. q0 cannot overflow as a result of the scaling: after it, exponent(x.hi) >= SCALE_FLOOR, and exponent(y.hi) >= -149, so q0's exponent stays below SCALE_FLOOR + 149 < 127.

Double32s.reciprocal(policy, y::Double32) -> Double32

Layer-2 reciprocal.

inv_dd_ forms 1 - q0*y, whose residuals sit near 2^0 whatever the operand's exponent, so unlike Double32s.divide this needs no scaling for accuracy. It needs it only for range: inv(y.hi) overflows once exponent(y.hi) <= -128, which loses the whole computation. Scaling only ever goes up, and that direction is exact.

A large operand is left alone deliberately. The correction q1 then underflows to zero, but the result's low word would sit below 2f0^-149 anyway — that is the documented precision taper of design section 5.5, not a defect.

Double32s.compensated_fma(policy, x, y, z) -> Double32

High-accuracy compensated x*y + z. See fma_dd_: no single-rounding or correct-rounding claim is made.

Double32s.square_root(policy, x::Double32) -> Double32

Layer-2 square root. Non-throwing; see sqrt_or_nan, which this is.

Raw()

Internal tag selecting the non-normalizing Double32 constructor. Passing Raw() asserts that the caller has already established the canonical invariant of design section 5.2. Not exported, not public API.

raw(hi::Float32, lo::Float32) -> Double32

Construct without canonicalizing. Every call site must be able to justify the canonical invariant from the algorithm that produced hi and lo. Internal; use canonicalize or the Double32 constructor instead.

Double32s.to_odd_hi(x::Double32) -> Float32

The high word rounded to odd against the low word: hi itself when the pair is exact at Float32 precision or hi's significand is already odd, otherwise the Float32 neighbor of hi on lo's side.

This is the required first step of every conversion narrower than Float32. hi alone is RN₂₄(x), and rounding it again at p < 24 bits misrounds exactly when hi lies on a tie of the p-bit lattice and lo breaks that tie — RN_p(RN₂₄(v)) ≠ RN_p(v) on those operands. Rounding the intermediate to odd instead of to nearest makes the composition exact for every p ≤ 22 (Boldo–Melquiond), which covers Float16 (11) and BFloat16 (8).

Derivation that the neighbor selection is round-to-odd: lo ≠ 0 means the value is not a Float32, so RO₂₄(v) is whichever of the two enclosing Float32s has an odd significand. hi = RN₂₄(v) is one of the two, and the other lies on lo's side of hi because v = hi + lo. If hi is odd it is the answer; otherwise the neighbor nextfloat(hi, sign(lo)) is adjacent and therefore odd.

At hi = ±floatmax(Float32) with lo pointing outward the neighbor is ±Inf32; the narrower formats overflow at those magnitudes anyway, so the final rounding is unaffected.

Double32s.iscanonical(x::Double32) -> Bool

Whether x satisfies the canonical invariants of design section 5.2. Host and test code only — this is the assertion that raw call sites are audited against, not something a kernel should call.

Note that invariant 8 makes the sign of a zero low word unspecified, so this compares low words with == rather than ===.

Double32s.canonical_nonfinite(x::Float32) -> Double32

Map a non-finite Float32 to its canonical Double32 pair. Branch free: NaN maps to the package quiet NaN, ±Inf pass through, both get +0.0f0 low.

Double32s.reference(x::Double32) -> BigFloat

The exact represented value. Host and test code only — always correct, unlike a Float64 conversion. Design section 22.1.

The zero low word is dropped rather than added, because of design section 5.2 invariant 8: the sign of a zero lo carries no information. Adding it would leak that sign into the reference — big(-0.0f0) + big(-0.0f0) is -0.0 while big(-0.0f0) + big(0.0f0) is +0.0, so two isequal pairs would get BigFloat references that are not isequal, and every test that treats this function as ground truth would inherit the discrepancy. Adding a nonzero lo is exact by canonicality, so no other case needs care.

Double32s.zero_sign(h::Float32, reference::Float32) -> Float32

Restore the IEEE sign of a zero result.

A pair kernel that ends in s = p + t loses the sign whenever p and t are zeros of opposite sign: (-0.0f0) + (+0.0f0) is +0.0f0. IEEE fixes the sign of a zero result from the operands — (-0) + (-0) is -0, 0 * -0 is -0 — so the layer-2 methods pass the corresponding Float32 high-word operation as the reference. The reference is already computed inside the kernel, so this costs a select rather than an arithmetic op.

Complex{Double32} branch cuts and every copysign-based algorithm depend on this, which is why it is not treated as a rounding detail.

Both operands must be zero for the substitution to fire, and that condition is load-bearing rather than defensive. A pair-level cancellation can be exact while the high-word operation is not: fma on the high words of an operand triple whose pair product-sum is zero returns a nonzero value, and substituting it would replace a correct zero with garbage — measured at 7.9e-8 relative before the iszero(reference) conjunct was added.

Double32s.scale2(x::Double32, n::Int) -> Double32

Unguarded exact scaling by 2^n, for call sites that have already established that neither word leaves range. Use ldexp anywhere else.