Public API

Double32s

A GPU-oriented double-word float built from two Float32 words.

The package exports the scalar type Double32, the word accessors HI, LO, and HILO, canonicalize, isnormal, the named arithmetic entry points add, sub, mul, divide, reciprocal, and rsqrt, the compensated reductions, and norm2, matmul, matmul!. The rest of the supported API is marked public and reached as Double32s.name; design section 1.3 records which names carry collision risk and why.

The module is named Double32s and the type Double32 because a Julia module cannot also bind its own name to a type: inside module Double32, the binding Double32 is the module, and using Double32 would import the module rather than the scalar.

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

A double-word number represented as the unevaluated sum hi + lo of two Float32 words, canonicalized so that hi == fl(hi + lo).

The represented value is the exact real sum hi + lo, not a Float32 evaluation of it. A canonical finite pair satisfies |lo| <= ulp(hi)/2, and hi is the nearest Float32 to the pair sum.

Double32 carries roughly 48 significant bits but keeps the Float32 exponent range, so it is not an IEEE format with a uniform significand lattice: the distance to the next representable pair is not eps(Double32) and can be as small as 2f0^-149 beside a high word of 1.0f0. See nominal_eps.

All NaN payloads may be canonicalized to NaN32 by any operation.

Examples

julia> Double32(1.0f0) / 3
0.333333333333333

julia> Float64(Double32(1.0f0) / 3)
0.33333333333333304

HI(x::Double32) -> Float32

High word of x: the nearest Float32 to the represented value, by canonical invariant 3 of design section 5.2.

Exported, and capitalized to match DoubleFloats.jl, where HI/LO/HILO name the same accessors. The lowercase spellings are deliberately not used: hi and lo are too common as local names to inject into a user's namespace.

LO(x::Double32) -> Float32

Low word of x, carrying what the high word cannot represent. See HI.

HILO(x::Double32) -> Tuple{Float32,Float32}

Both words as a tuple, high first. This is the layer-1/layer-2 adapter of design section 8.0 and compiles away entirely.

Double32(HILO(x)...) === x for every canonical x, the sign of a zero included — see canonicalize.

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

Return the canonical Double32 whose value is the exact real sum hi + lo.

This is the normalizing constructor: Double32(hi, lo) calls it. It is the only public route that establishes the invariants of design section 5.2, so a raw pair from anywhere else must be justified by the algorithm that produced it.

Two behaviors are worth stating because they are not what a plain two_sum would give:

• a non-finite high word (including genuine overflow of two finite words) is mapped to the canonical non-finite pair, so lo is never NaN beside an infinite hi;
• when both words are zero the sign of hi is preserved, rather than taking IEEE's (-0.0f0) + (+0.0f0) == +0.0f0. This is what makes Double32(x.hi, x.lo) === x hold for every canonical x, so a value can always be rebuilt from Double32s.HILO without losing the sign of a zero.

This is also the byte-level normal form. A zero low word always comes back +0.0f0, so two pairs that are isequal — which, by design section 5.2 invariant 8, may differ in that sign bit — canonicalize to identical bits. That is the tool for anything that must compare or key on representation rather than value: serialization, byte-wise deduplication, a content hash of a raw buffer. Verified over the signed-zero surface: 0 negative-zero low words emitted, and every isequal pair byte-identical after canonicalization.

The sign is deliberately not normalized at construction. raw is the hot path and negation is raw(-hi, -lo); clearing the sign there would add a select to every negation to fix something that is not a defect, since no value-based operation can observe it (see the show docstring).

Examples

julia> Double32s.canonicalize(1.0f0, 2.0f0^-30)     # already canonical
1.0000000009313226

julia> Double32s.canonicalize(2.0f0^-30, 1.0f0)     # words out of order
1.0000000009313226

julia> Double32s.canonicalize(1.0f0, 0.5f0)         # |lo| > ulp(hi)/2: renormalizes
1.5

canonicalize(x::Double32) -> Double32

Renormalize a possibly weakly-normalized pair. A no-op for canonical input.

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

Whether x is finite, nonzero, and not subnormal in its high word.

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

Whether Float64(x) reproduces the exact pair sum, and therefore whether a Float64 reference may be used for x in validation.

True iff the pair spans at most 53 bits. Design section 4.2 gives the proof; the short version is that |lo| <= ulp(hi)/2 bounds the low word's magnitude and constrains neither its bit position nor the gap between the two words, so both a subnormal lo and a normal lo far below its permitted bound push the pair past 53 bits.

Validated against a 400-bit BigFloat reference over 1,701,602 canonical pairs, 0 mismatches.

Examples

julia> Double32s.float64_reference_is_exact(Double32(1.0f0, 2.0f0^-52))
true

julia> Double32s.float64_reference_is_exact(Double32(1.0f0, 2.0f0^-100))  # normal lo, 101-bit span
false

julia> Double32s.float64_reference_is_exact(Double32(2.0f0^-96, -2.0f0^-149))  # cancellation buys a bit
true

Double32s.lsb_exponent(z::Float32) -> Int

Binary exponent of the lowest set significand bit of a nonzero z, i.e. the k with z == m * 2^k for an odd integer m.

Returns typemin(Int) for zero, so that "no constraint" compares correctly against any bound. Non-finite z is not meaningful here and is not checked.

This is the quantity the Float64-exactness condition of design section 4.2 is stated in: a canonical pair converts exactly iff lo == 0 or lsb_exponent(lo) >= exponent(hi + lo) - 52.

Examples

julia> Double32s.lsb_exponent(1.0f0), Double32s.lsb_exponent(1.5f0)
(0, -1)

julia> Double32s.lsb_exponent(reinterpret(Float32, 0x00000001))   # 2^-149
-149

Double32s.raw_exponent(z::Float32) -> Int

Binary exponent of z, agreeing with Base.exponent on every finite nonzero argument — and never throwing.

Base.exponent raises DomainError for 0, Inf, and NaN. Every call site in this package excludes those first, so the throw is unreachable; the compiler cannot prove it, and the exception path survives into the generated code. Design section 14 forbids exactly that in device code, and it was measured rather than assumed: before this function existed, @code_llvm for /, inv, norm2, and even sqrt_or_nan — the operation documented as the device-safe one — each contained calls to throw1/throw2 from Base.exponent. After it, +, -, *, /, inv, fma, and sqrt_or_nan contain no throw reference, no alloca, and no double, which test/gates.jl now asserts.

What it bought, measured on both sides. On the host: nothing — an A/B of / and sqrt_or_nan against the Base.exponent versions came out at 0.98× and 0.93×, i.e. inside the noise. The gain is on the device, where the dead exception path still costs registers and spill space:

(RTX 5090, CUDA.registers on an elementwise kernel.) Local memory is spill space, so the 43% cut on sqrt_or_nan is the number that matters. The other reason to do this is categorical rather than quantitative: some backends decline to compile a kernel containing exception machinery at all.

Out-of-domain inputs return a value rather than raising: 0 gives typemin(Int) >> 1, and Inf/NaN give typemax(Int) >> 1. Both are chosen so that arithmetic on the result (the +, -, and max that call sites do with it) cannot overflow, and both are far outside the range any guard band tests. A caller that has not already excluded these gets a meaningless scale factor rather than a fault — which is why every call site still excludes them, and says so.

Examples

julia> Double32s.raw_exponent(1.0f0), Double32s.raw_exponent(0.5f0)
(0, -1)

julia> Double32s.raw_exponent(reinterpret(Float32, 0x00000001))   # 2^-149
-149

julia> Double32s.raw_exponent(0.0f0) == typemin(Int) >> 1
true

Double32s.sqrt_nonneg(z::Float32) -> Float32

sqrt with no domain check: llvm.sqrt.f32 directly, returning NaN for a negative argument instead of throwing.

Base.sqrt(::Float32) tests the sign and raises DomainError, and that branch is unreachable at the one call site here — Double32s.sqrt_dd_ is reached only after the layer-2 method has excluded negatives — but survives into the generated code, dragging throw_complex_domainerror into every kernel that takes a square root. Same reasoning as Double32s.raw_exponent.

Double32s.capability_report(target = nothing) -> NamedTuple

Everything the capability probe can say about an execution target, in one call.

target is nothing for the host thread, an array whose backend should be questioned, or — with KernelAbstractions loaded — a Backend. Loading CUDA or AMDGPU adds vendor fields to the report for their backends.

Fields:

• target — a short description of what was probed;
• gradual_underflow::Bool — see Double32s.has_gradual_underflow;
• fused_multiply_add::Bool — see Double32s.has_fused_multiply_add;
• lanes::NTuple{PROBE_LANES,Bool} — the raw per-lane outcomes;
• functional::Bool — whether there was anything to probe. Always true for the host and for an array prototype; a KernelAbstractions.Backend with no device present reports false, with every other field falsified rather than guessed;
• usable::Bool — fused_multiply_add, without which no error bound in this package holds. Gradual underflow is not required for usability: it narrows the exactness claims near the underflow boundary rather than voiding them.

Vendor extensions append their own fields. The CUDA extension adds device, capability, and math_mode — the last one matters, because CUDA.FAST_MATH enables flush-to-zero and permits contraction, so a report that says gradual_underflow = false on NVIDIA hardware is far more often a task-local math mode than a property of the silicon.

Examples

julia> r = Double32s.capability_report();

julia> r.usable && r.fused_multiply_add
true

Double32s.has_fused_multiply_add() -> Bool
Double32s.has_fused_multiply_add(prototype::AbstractArray) -> Bool
Double32s.has_fused_multiply_add(backend::KernelAbstractions.Backend) -> Bool

Whether fma is a genuine single-rounding fused multiply-add in this execution environment.

This is the package's hardest requirement. Every error-free transformation on a product goes through EFT.two_prod_, whose residual is identically zero when the multiply and the add round separately — so an unfused backend does not produce a slightly worse Double32, it produces a Float32 wearing two words. The check stays in the normal range, so unlike Double32s.has_gradual_underflow it is unaffected by flush-to-zero.

Examples

julia> Double32s.has_fused_multiply_add()
true

Double32s.has_gradual_underflow() -> Bool
Double32s.has_gradual_underflow(prototype::AbstractArray) -> Bool
Double32s.has_gradual_underflow(backend::KernelAbstractions.Backend) -> Bool

Whether an execution environment honors gradual underflow, checked by running the design section 5.9 probe list: subnormal addition, a subnormal multiplication result, a subnormal fma result, an EFT.two_prod_ case whose exact residual is subnormal, and subnormal sign preservation.

Returns false under flush-to-zero — e.g. after set_zero_subnormals(true), or on a device compiled with FTZ — in which case the strict kernels still run but their subnormal-range exactness claims (design sections 7.4, 9) do not hold, and validation must widen its tolerances near the underflow boundary. Reach for Double32s.capability_report when the distinction between FTZ and a missing fused multiply-add matters, which it does when the answer is false.

The no-argument method allocates nothing, so the identical call is available on a device. The result is deliberately not cached: set_zero_subnormals is a per-thread runtime switch and CUDA's math mode is a per-task setting, so a load-time answer can go stale. The probe is a handful of flops; a backend should cache it per device at launch-configuration time.

Examples

julia> Double32s.has_gradual_underflow()
true

julia> Double32s.has_gradual_underflow(Float32[])   # host array, host answer
true

Double32s.probe_lane(code::Int32, a::Float32, b::Float32) -> Bool

One lane of the backend capability probe. Pure Float32, no allocation, no throw: safe to call from inside a device kernel.

Lane 6 is what separates the two failure modes. Lanes 1–5 all die under flush-to-zero, and lane 4 additionally dies if fma is not a genuine single-rounding fused operation; lane 6 stays in the normal range, so it isolates fusion. A backend that passes 6 and fails 1–5 has FTZ; one that fails both has no fused multiply-add and cannot host this package's error bounds at all.

Operands come from the caller — never from literals inside the lane. With constant operands the compiler folds the whole probe to true at compile time, and a folded probe ignores set_zero_subnormals and device FTZ flags entirely. This is not hypothetical; it is how the first version of this probe failed its own test. Expected values, by contrast, are literals, and must be, as reinterpret bit patterns: under FTZ even ldexp(1.0f0, -145) flushes to zero, so a probe that computes its own expectation with the faulty arithmetic reports the fault as health.

See Double32s.has_gradual_underflow, Double32s.has_fused_multiply_add, and Double32s.capability_report for the host and per-backend entry points.

Double32s.probe_results() -> NTuple{PROBE_LANES,Bool}
Double32s.probe_results(prototype::AbstractArray) -> NTuple{PROBE_LANES,Bool}
Double32s.probe_results(backend::KernelAbstractions.Backend) -> NTuple{PROBE_LANES,Bool}

Run every lane of Double32s.probe_lane and return the outcomes as a host tuple.

The no-argument method answers for the host thread and allocates nothing. The prototype method answers for whatever device owns that array, by copying the probe operands there, broadcasting, and copying six Bools back; it needs no GPU package — only an array type supporting similar, copyto!, and broadcast, which every GPUArrays-based backend does. The KernelAbstractions.Backend method is added by the KernelAbstractions extension and launches a real kernel instead, which is what a backend without an array type still has to offer.

Double32s.probe_operands() -> (codes, a, b)

The three host vectors that drive Double32s.probe_lane, each of length Double32s.PROBE_LANES. Copy them to a device to run the probe there; Double32s.probe_results does exactly that.

Double32s.build_report(target, lanes) -> NamedTuple

Assemble a Double32s.capability_report from lane outcomes that have already been collected. Backend extensions use this to report a device that could not be probed at all, rather than letting the probe throw out of a diagnostic call.

Double32s.PROBE_LANES

Number of lanes in the capability probe. See Double32s.probe_lane.

Double32s.UNDERFLOW_LANES

The lane range that tests gradual underflow; the remaining lane tests fma fusion. See Double32s.probe_lane.

Double32s.NO_LANES

All-false lane outcomes, for reporting a target that could not be probed.

Double32s.nominal_eps(::Type{Double32}) -> Double32

The nominal 48-bit unit in the last place, 2.0^-47.

This is an algorithmic precision constant, not the distance to the next representable pair. Design section 5.6.

Double32s.unit_roundoff(::Type{Double32}) -> Double32

The nominal unit roundoff 2.0^-48, i.e. half nominal_eps.

ArithmeticPolicy

Compile-time selector for the accuracy/performance contract of a Double32 operation. Instances are singletons and are always passed by value as the first argument, so that dispatch — not a runtime flag — chooses the kernel.

Three rules follow from making the policy a type rather than a flag, and they are what buys the performance:

1. Never test a policy at runtime. if policy isa FastPolicy inside a kernel is a code smell even though the compiler usually removes it; write a second method instead.
2. Never pass a policy as a keyword argument into device code. Keyword arguments do not participate in dispatch. Accept keywords at the host API boundary, resolve them to singletons there, pass positionally from then on.
3. Never store a policy in a mutable or abstractly typed field. Use a type parameter if a policy must be carried in a struct.

StrictPolicy()

Preserve the documented operation order, use explicit FMA, guard exceptional values, and require capability tests for claims involving subnormals. The default.

NativePolicy()

Preserve the same public result contract where possible but permit backend-specific instruction selection and launch tuning.

FastPolicy()

Explicitly relax error and exceptional-value guarantees. Never the implicit behavior of +, *, or /.

Double32s.DEFAULT_POLICY

The policy the Base operators use: StrictPolicy().

Double32s.sqrt_or_nan(x::Double32) -> Double32

Square root that returns NaN for a negative argument instead of throwing.

This is the spelling device kernels should use. Base.sqrt throws a DomainError to match sqrt(::Float32), and a throw in a device kernel is not free: it forces the exception machinery into the generated code and can abort the kernel.

The operand is always scaled to an even exponent near zero first, and that is not an optimization to be removed. sqrt_dd_ forms the residual x - s0*s0, whose magnitude is about |x| * 2^-48; below exponent(x.hi) ≈ -101 that residual falls under 2f0^-149 and is silently lost even though the result remains fully representable. Measured before scaling was added: sqrt(Double32(2.0165723f-36, 8.0f-44)) had relative error 2.4e-11, about 35 correct bits instead of 48. An even shift keeps both the scaling and its inverse exact.

Exceptional values follow IEEE and Base:

sqrt(+0)   -> +0        sqrt(-0)   -> -0
sqrt(+Inf) -> +Inf      sqrt(x<0)  -> NaN
sqrt(NaN)  -> NaN

Double32s.compensated_dot(x, y) -> Double32

Dot product accumulated in Double32.

For Float32 inputs each product is formed exactly by two_prod_ and the exact pair is accumulated, so no information is lost before accumulation — strictly more accurate than rounding each product to Float32 first. The raw call is justified because two_prod_ returns a magnitude-ordered nonoverlapping pair by construction.

"Exactly" carries two_prod_'s condition: the product residual is representable only for products ≥ 2^-101 (and only under gradual underflow — gate on has_gradual_underflow). Smaller products contribute their correctly rounded value with the residual truncated at 2^-149 — the format's underflow taper, harmless to a sum whose magnitude sits in the normal range.

When LinearAlgebra is loaded, an extension may add LinearAlgebra.dot(::AbstractArray{Double32}, ::AbstractArray{Double32}) delegating here — legitimate, since the element type is ours.

Double32s.compensated_sum(A; policy = DEFAULT_POLICY) -> Double32

Sequential Double32-accumulated sum with a specified summation order.

Prefer sum(Double32, A) unless you need the order pinned: Base's pairwise reduction is more accurate and already uses our +. This exists for reproducibility work, where a fixed tree matters more than the error bound.

Double32s.compensated_mapreduce(f, op, A; init, policy) -> Double32

Sequential Double32 map-reduce with a specified order, the general form of compensated_sum.

Double32s.compensated_norm2(x) -> Double32

Euclidean norm via compensated_dot — sequential, host-ordered, like the rest of the compensated_* family — with an exact power-of-two scaling so the squares neither overflow nor underflow prematurely. The parallel, reproducible-tree counterpart is Double32s.norm2.

SplitDouble32Array(hi::AbstractArray{Float32,N}, lo::AbstractArray{Float32,N})

A Double32 array stored as two Float32 arrays — structure of arrays.

Use it when a kernel reads only high words during a seed phase, when vendor Float32 code produces a high approximation followed by a correction, when a reduction stores high and low partials separately, or when benchmarks show better memory behavior than an array of Double32.

Both backing arrays must have identical axes. Mutating one component alone through the public API is not possible; the fields remain accessible for expert use, with the canonical invariant then the caller's responsibility.

The backing arrays may live on a device. Adapt.adapt recurses into both, so a SplitDouble32Array{2,CuArray{Float32,2},CuArray{Float32,2}} passed to a kernel arrives as a SplitDouble32Array of device arrays; see Double32s.allocate_split to build one for a KernelAbstractions backend. Adaptation deliberately bypasses the axes check — inside a kernel that check would compile to a DimensionMismatch throw and drag the exception machinery in with it.

Examples

julia> A = Double32s.SplitDouble32Array([1.0f0, 2.0f0], [0.0f0, 0.0f0]);

julia> A[2]
2.0

julia> A[1] = Double32(1.0f0) / 3;  A[1]
0.333333333333333

Double32s.split(A::AbstractArray{Double32}) -> SplitDouble32Array

Convert array-of-structures to structure-of-arrays.

Backend generic: the two word arrays are produced by broadcasting, so a CuArray{Double32} splits into two CuArray{Float32} on the device with no host round trip and no scalar indexing.

Double32s.unsplit(A::SplitDouble32Array) -> AbstractArray{Double32}

Convert structure-of-arrays back to array-of-structures, on whichever backend holds the words. Elements are already canonical, so no normalization runs.

Double32s.split!(S::SplitDouble32Array, A::AbstractArray{Double32}) -> S

In-place Double32s.split: write the words of A into the existing backing arrays of S.

For a conversion inside a loop this is the form to use — split allocates two arrays per call, split! allocates nothing. Written as two in-place broadcasts, so it is one fused kernel per word on any backend and needs no GPU package loaded.

Examples

julia> A = Double32[1, 2] ./ 3;

julia> S = Double32s.allocate_split_host(size(A));

julia> Double32s.split!(S, A); S[1] === A[1]
true

Double32s.unsplit!(A::AbstractArray{Double32}, S::SplitDouble32Array) -> A

In-place Double32s.unsplit. Allocation free; see Double32s.split!.

Double32s.allocate_split(backend, dims...) -> SplitDouble32Array

Allocate a Double32s.SplitDouble32Array of the given shape on a KernelAbstractions backend, with both word arrays uninitialized.

Defined by the KernelAbstractions extension; using KernelAbstractions (or any package that loads it, such as CUDA or AMDGPU) brings it into existence.

Double32s.allocate_split_host(dims) -> SplitDouble32Array

Allocate a host-memory SplitDouble32Array with uninitialized words. The backend-generic form is Double32s.allocate_split, which needs KernelAbstractions loaded.

sum(Double32, A::DenseArray; policy = DEFAULT_POLICY, kwargs...) -> Double32
sum(A::AbstractArray{Double32}; policy = DEFAULT_POLICY, kwargs...) -> Double32

Parallel sum in Double32, on whatever backend owns A.

The first form is the accurate-accumulation pattern: A may hold Float32 — an element converts exactly and for free — and the reduction tree is the reproducible AcceleratedKernels one. The second form routes every sum of Double32 data through the same tree by dispatch.

policy selects the addition kernel: StrictPolicy() (the default) gives the 3u² bound of design section 8.1; FastPolicy() gives roughly Float32-level error per step for about 8 flops instead of 20.

Measured on 50 M Float32 on an RTX 5090: 0.439 ms against 0.146 ms for a plain Float32 sum — the ~3× is the arithmetic cost of 48-bit accumulation, not a launch artifact.

The first method is bounded at DenseArray, not AbstractArray, to stay unambiguous against packages that specialize sum(f, a) on their own array types (StaticArrays was the live collision: typeintersect(StaticArray, DenseArray) is Union{}, so the narrower bound removes the crossing). A view or wrapper falls back to Base's generic sum(Double32, A) — same values, Base's tree.

Examples

julia> sum(Double32, Float32[1, 2, 3, 4])
10.0

LinearAlgebra.dot(Double32, x, y; policy = DEFAULT_POLICY, kwargs...) -> Double32

Dot product accumulated in Double32, for inputs of any element type, on whatever backend owns x. The first argument is the accumulator type, exactly as in sum(Double32, A) — a dot(::Vector{Float32}, ::Vector{Float32}) method returning Double32 would change the meaning of every existing Float32 dot, so the widening must be asked for.

When x and y hold Float32, each product is formed exactly as a Double32 by EFT.two_prod_ before it is accumulated, so nothing is lost before the summation itself — strictly better than rounding each product to Float32 first, and better than any Float32 dot. This is the parallel counterpart of Double32s.compensated_dot, which is sequential and host only.

It is also the faster spelling on device: on 50 M Float32 pairs on an RTX 5090 it beats sum(x .* y) by 1.5× (0.418 ms against 0.630 on 50 M pairs), because the fused reduction never materializes the product array that the broadcast form writes and reads back.

Products below 2⁻¹⁰¹ lose their residual to subnormal underflow; see EFT.two_prod_ and gate on Double32s.has_gradual_underflow.

dot(x::StridedVector{Double32}, y::StridedVector{Double32}) -> Double32

Dot product of Double32 vectors.

On the host: the channel-compensated accumulation the matrix multiply uses, one renormalization at the end. Sequential in k (the TwoSum chain is a loop-carried dependency), which still measures far ahead of the generic pairwise reduction over pair arithmetic.

On any other backend: the AcceleratedKernels reduction — see dot(Double32, x, y). The signature and the branch are both load-bearing: CuVector is a StridedVector, so array-type dispatch cannot separate host from device — the backend test does, and it is a compile-time constant. And the bound must stay StridedVector, not AbstractVector: sparse vectors are non-strided, and an AbstractVector method is ambiguous against SparseArrays' own dot specializations.

Double32 is real, so dot needs no conjugation.

Double32s.norm2(x; kwargs...) -> Double32

Euclidean norm on whatever backend owns x, with an exact power-of-two prescaling so the squares neither overflow nor underflow. Accepts Float32 or Double32 elements. (LinearAlgebra has no norm2 generic to extend, so the name is package-local.)

The scale is 2^-k with k taken from the largest |x[i]|, which costs one extra reduction pass to find. Scaling by a power of two is exact, so the only error is the summation and the final sqrt.

mapreduce(f, op, A::AbstractArray, neutral::Double32; kwargs...) -> Double32

Parallel reduction of op(f(A[1]), f(A[2]), …) on whatever backend owns A, returning a host scalar. Wrapper over AcceleratedKernels.mapreduce.

neutral must be the identity of op and fixes the accumulator type. It is positional rather than a keyword because AK's own default, neutral_element(op, eltype(src)), is the wrong type for this package's headline pattern: reducing a Vector{Float32} into a Double32 accumulator needs zero(Double32), not zero(Float32). Requiring it makes the accumulator type impossible to get wrong — the same rule design section 14.4 states as advice, expressed as a signature. It is also not Base's init: AK seeds every block with it, so a non-identity value would be added once per block.

Determinism. For a fixed block_size and length the tree shape is fixed, so the result is reproducible across runs and across backends — a CPU and a CUDA reduction of the same array return the same bits, which the test suite asserts. It is not the same answer as Base's pairwise tree; Double32 addition is not associative, so those differ in the last places. Neither is wrong; only one is reproducible.

Extra keywords pass through to AcceleratedKernels.mapreduce — block_size, temp, switch_below, dims.

The backend is always KernelAbstractions.get_backend(A) — the backend that owns the memory is the only one that can run the kernel, so it is derived, never a parameter.

reduce(op, A::AbstractArray, neutral::Double32; kwargs...) -> Double32

mapreduce with f = identity. eltype(A) must equal typeof(neutral).

map!(f, dest::AbstractArray{Double32}, a)              -> dest
map!(f, dest::AbstractArray{Double32}, a, b)           -> dest
map!(f, dest::AbstractArray{Double32}, a, b, c)        -> dest

Elementwise dest[i] = f(a[i], …) through AcceleratedKernels.foreachindex, on whatever backend owns dest.

Prefer broadcast. dest .= f.(a, b) compiles to one fused device kernel through GPUArrays, handles mixed shapes and scalars, and is the spelling this package's documentation uses everywhere; measured on an RTX 5090 over 50 M elements the two are a tie. These methods exist for backends that provide KernelAbstractions but not GPUArrays broadcast, and for callers who want to pass AK's tuning keywords (block_size, max_tasks, min_elems) explicitly.

Arities are fixed at one, two, and three rather than variadic: a variadic kernel would splat a tuple of device arrays inside the body, which is the construct most likely to defeat inference on a backend compiler, and these three cover every use in this package.

Stricter than Base.map! in one respect, on purpose: all arrays must share axes (Base silently stops at the shortest), and dest must be distinct from the inputs.

Double32s.DEFAULT_GROUPSIZE

Workgroup size used by the ka_* kernels when none is given: 256.

Chosen, not inherited. It must be a power of two for the tree reduction, and 256 is the smallest power of two that satisfies all four constraints at once:

• it is a whole number of 32-lane NVIDIA warps and of 64-lane AMD wavefronts, so neither vendor runs a partially-populated wave;
• 256 × sizeof(Double32) is 2 KiB of local memory per workgroup, against 48–228 KiB of shared memory per SM on CUDA and 64 KiB of LDS per CU on AMD — small enough that local memory never limits occupancy;
• the reduction tree is log2(256) = 8 barriers deep, and each stage consumes 2 × 256 elements, so a 100-million-element array finishes in 4 passes;
• it is large enough to amortize the launch, and small enough that a short array still fills several workgroups.

Measured on an RTX 5090, sum(Double32, A) over 50 M Float32:

256 wins, and the whole spread is 4.5% — the reduction is dominated by the arithmetic of Double32 addition, not by launch shape, so this is a default worth having and not worth agonizing over. CUDA's launch_configuration suggests 768 instead; see Double32s.launch_config for why the measurement is preferred.

Override per call with the groupsize keyword when a measurement disagrees.

Double32s.PolicyOp(f, policy) -> callable

Bind a policy to a policy-taking operation, giving a plain binary callable (a, b) -> f(policy, a, b).

This is rule 2 of ArithmeticPolicy made concrete for kernels. A reduction needs an op of two arguments, but every accuracy-selecting operation in this package takes the policy first. Wrapping in a closure would work on the host and is the wrong thing to hand a kernel: a closure over a local is an ordinary struct with a field, and its type is not the same across two call sites. PolicyOp is zero-size — both fields are singletons — so it costs no register, no argument-buffer slot, and no launch-configuration difference, and the reduction specializes on it completely.

Examples

julia> op = Double32s.PolicyOp(Double32s.add, Double32s.FastPolicy());

julia> sizeof(op)
0

julia> op(Double32(1.0f0), Double32(2.0f0))
3.0

Double32s.matmul(A, B; kwargs...) -> AbstractMatrix{Double32}

Allocating form of Double32s.matmul!.

Examples

julia> A = Float32[1 2; 3 4]; B = Float32[5 6; 7 8];

julia> Float64.(Double32s.matmul(A, B))
2×2 Matrix{Float64}:
 19.0  22.0
 43.0  50.0

Double32s.matmul!(C, A, B; kwargs...) -> C

Tiled matrix multiply with Double32 accumulation, on any KernelAbstractions backend. Design section 19.5.

Both Float32 × Float32 → Double32 and Double32 × Double32 → Double32 are supported, and they differ in more than the element type:

• Float32 operands: each product is formed exactly by EFT.two_prod_ before it enters the accumulator, so a dot product of length K loses nothing until the additions. Shared tiles hold Float32, half the memory of the Double32 case.
• Double32 operands: the ordinary pair multiply and its 2u² bound. Nothing exact is available — the operands have already spent their bits.

Measured on an RTX 5090 (normwise error against a Float64 reference):

The split-channel inner op is what made those times: replacing the pair multiply + pair add (~30 flops, long dependency chain) with the fused product-accumulate (~11–14 flops, short chain) took 1024³ from 5.62 ms to 0.33 ms for Float32 operands and from 22.1 ms to 0.42 ms for Double32 operands — 17× and 53× — while the error stays at the accuracy model (~28 nominal ulps normwise at K = 1024, against ~4.4e-07 for Float32 BLAS).

A square problem whose size is an exact multiple of the output block (tile × block, 16 by default) takes a bounds-free instantiation: every range test is statically true there and compiles away. Everything else takes the checked path, which pads out-of-range loads with zeros.

At these times the kernel is within 7–9× of cuBLAS SGEMM on the same device — for ~7 more decimal digits. No vendor library accumulates Double32, so there is nothing faster to call at this accuracy.

Results are bit-identical across backends for a given tile and block — the tile loop fixes the summation order. Since the reductions moved to AcceleratedKernels, this is the one remaining cross-backend-reproducible routine in the package.

Keyword arguments

• tile — workgroup edge, default MATMUL_TILE; a power of two.
• block — register-blocking factor, default MATMUL_BLOCK.

The backend is always KernelAbstractions.get_backend(C): the backend that owns the output is the only one that can run the kernel, so it is derived, never a parameter.

There is no policy keyword: the split-channel accumulation is the algorithm, and it is both faster than the StrictPolicy pair chain and holds its measured error. A knob that selected a slower, equally-accurate path would be surface area without a decision behind it.

Double32s.MATMUL_TILE

Tile edge for Double32s.matmul!: 16.

16 × 16 is 256 work items — the same DEFAULT_GROUPSIZE the reductions use, and for the same reasons: whole warps on NVIDIA, whole wavefronts on AMD. Local memory is two Float32 tiles plus one Double32 accumulator tile, 16² × (4 + 4 + 8) = 4 KiB per workgroup, which leaves occupancy unconstrained on both vendors.

Double32s.MATMUL_BLOCK

Blocking factor for Double32s.matmul!: 1.

Each work item computes a MATMUL_BLOCK × MATMUL_BLOCK sub-block of C, so MATMUL_BLOCK > 1 reuses each staged value that many times — the textbook fix for a shared-memory-bound GEMM.

Measured twice, with opposite-looking results, and the second measurement is the one that stands. Under the original pair-arithmetic inner op, block = 2 was 3× slower (16.08 vs 5.62 ms at 1024³): KA's @private maps to per-thread local memory — global-backed spill space, not registers — so blocking multiplied local-memory traffic instead of removing shared-memory traffic. Under the split-channel inner op that replaced it, the accumulator traffic is a fraction of the old cost and block = 2 measures within ±5% of block = 1 (0.327 vs 0.340 ms at 1024³, 2.48 vs 2.62 ms at 2048³).

1 stays the default: simpler indexing, fastest at 512³, and never worse than 5% behind. The keyword stays tunable, and both measurements stay recorded so a future reader knows the 3× penalty was a property of the old inner op, not of blocking.

Double32s.launch_config(backend; kind = :reduction) -> NamedTuple

The workgroup size this package uses on backend, and what the vendor's occupancy heuristic suggests instead.

Fields: groupsize (what is used), suggested (the occupancy API's answer, or nothing), kind, and any vendor detail. Design section 15.3 asked for occupancy-informed launch tuning and called its candidate numbers "hypotheses, not fixed constants" — this is where the hypothesis is tested.

The measured answer and the heuristic disagree, and the measurement wins. On an RTX 5090, CUDA's launch_configuration suggests 768 threads for the reduction stage; measured times over 50 M elements are 128 → 0.413 ms, 256 → 0.395 ms, 512 → 0.407 ms, 1024 → 0.406 ms. The heuristic maximizes resident warps to hide memory latency, but this kernel is arithmetic bound on Double32 addition — extra resident warps buy nothing, while a larger block deepens the intra-block tree by log2(W) barriers. Both numbers are reported so the disagreement stays visible instead of being buried in a constant.

Requires KernelAbstractions; vendor detail requires CUDA or AMDGPU.

Double32s.kernel_report(backend) -> Vector{NamedTuple}

Register, local-memory, and occupancy figures for representative Double32 kernels on backend — design section 15.1's "PTX and register-count test hooks", and the Phase 3 occupancy report.

One row per operation (+, *, /, fma, sqrt_or_nan, and the reduction stage), each with registers, local_bytes, shared_bytes, max_threads, and occupancy at the package's default workgroup size. The Phase 3 gate item "no spills or occupancy collapse" is checked against exactly this: local_bytes is spill space, so a nonzero figure that grows with a code change is the warning.

Requires CUDA (implemented in the CUDA extension).

Double32s.wavefront_size(backend) -> Int

Lanes per hardware execution unit on backend: a CUDA warp, an AMD wavefront, or 1 for CPU().

Design section 16.2 is explicit that 64 must not be hardcoded for AMD — RDNA parts run wave32, CDNA parts wave64, and the same device can differ by compilation mode. The number is therefore always queried, never assumed, and this is the one place to query it.

The package's own kernels do not branch on it: DEFAULT_GROUPSIZE is 256, a whole multiple of 32 and 64, so no partially populated wave occurs on either vendor at the default. It is exposed because a caller choosing a non-default groupsize needs it to keep that property.

Requires KernelAbstractions; the vendor answers require CUDA or AMDGPU.

Double32s.decimal(x::Double32) -> String

The shortest decimal string s with isequal(parse(Double32, s), x) — verified by parsing back before returning, not promised by digit-count analysis. This is what show, repr, string, and interpolation all print.

A short decimal is not always available, and the length is not slack. Wide-gap pairs (design section 23.11) need up to ~90 significant digits, because a high word near floatmax(Float32) above a subnormal low word pins ~277 bits. Even a pair whose low word is zero can be long: Double32(1.0f-7) needs 39 digits, since the pair lattice is fine enough that no shorter decimal selects that exact value over its neighbors.

The 97-digit fallback is unreachable — 97 digits over-determine 384 bits — and is asserted rather than silently trusted.

The sign of a zero low word never appears, because the rendering is computed from the represented value. Design section 5.2 invariant 8 makes that sign non-information: - flips it unconditionally, so -one(Double32) and Double32(-1.0f0) are bit-distinct pairs holding one value. A rendering that exposed the sign would make repr a function of the bit pattern rather than of the value, giving isequal(a, b) while repr(a) != repr(b). Measured against the signed-zero surface before the value-based rendering existed: 18 disagreements over its 36 isequal pairs; now 0. The sign of a zero high word does appear, because there it is the value's own sign.

Use Double32s.HILO to inspect the representation.

Examples

julia> Double32s.decimal(sqrt(Double32(5)))
"2.23606797749979"

julia> Double32s.HILO(sqrt(Double32(5)))         # the representation
(2.236068f0, -3.283041f-8)

parse(Double32, s::AbstractString) -> Double32

Throwing form of tryparse.

tryparse(Double32, s::AbstractString) -> Union{Double32, Nothing}

Decimal (or MPFR-accepted) string to the nearest canonical pair, via a 384-bit BigFloat intermediate — host-only, per design section 11.5. Returns nothing when s is not a number. "NaN", "Inf", "-Inf", and signed zeros parse to their canonical pairs.

The intermediate is exact for every string this package prints. For arbitrary input the composite rounding is correct unless the written value lies within 2^-384 relative of a tie of the pair lattice without being that tie — below any decimal a user can plausibly produce by accident.

Double32s.rsqrt(x::Double32) -> Double32

Reciprocal square root, 1/sqrt(x).

It is sqrt(inv(x)), and the order is the whole point. Design section 20.2 lists reciprocal square root among the seed-and-refine candidates, so a dedicated Newton kernel was written and measured against both compositions:

The composition wins, and the order of the composition wins again, for a reason worth keeping: sqrt halves the relative error of its argument, while inv passes it through. Doing the reciprocal first puts inv's 1.16 ulp through that halving and adds sqrt's own 0.75 — about 1.33, which is what is observed. Doing it the other way adds the two errors instead.

So there is no dedicated kernel and no new code path: fewer operations, better accuracy, and nothing extra to maintain.

rsqrt(0) is +Inf, rsqrt(-0.0) is -Inf, and a negative argument gives NaN rather than throwing — this is a kernel-safe function.

Examples

julia> Double32s.rsqrt(Double32(4.0f0))
0.5

julia> Double32s.rsqrt(Double32(0.0f0))
Inf

julia> isnan(Double32s.rsqrt(Double32(-1.0f0)))
true

Double32s.log_or_nan(x::Double32) -> Double32

Natural logarithm returning NaN for a negative argument instead of throwing. The kernel-safe spelling, exactly as Double32s.sqrt_or_nan is for sqrt.

log(+0)   -> -Inf       log(-0)   -> -Inf
log(+Inf) -> +Inf       log(x<0)  -> NaN
log(1)    -> +0         log(NaN)  -> NaN

Double32s.TRIG_ARG_LIMIT

Largest |x| for which the trigonometric functions reduce their argument exactly: 2^20.

This is the large-argument policy design section 20.4 asks for, and it is a limit rather than a promise. Reduction subtracts k·π/2 using the three-word Cody–Waite split Double32s.PIO2_C1 / PIO2_C2 / PIO2_C3, which together carry about 72 bits of π/2. A Double32 argument near 2^n needs about n + 48 correct bits of π/2 for the reduced value to keep full precision, so 72 bits covers n <= 24 with margin; the limit is set at 2^20 to leave room for the reduction's own rounding.

Beyond it the functions return NaN. An earlier draft of this docstring claimed they would "still return a value in [-1, 1]". They did not: the Cody–Waite subtraction cancels catastrophically once q outgrows the split's precision, and sin(2^30) came back as 9.88e8 — a silently wrong answer with no indication anything had gone wrong, which is the worst possible failure for a numerical routine.

NaN is the policy because it is loud. A caller can test for it; a caller cannot test for "this number in [-1, 1] is fictional". Correct reduction past this point needs Payne–Hanek with a multi-thousand-bit 2/π — a table and a loop that do not belong in a GPU kernel — so it is recorded as open work rather than approximated.

Callers who need large arguments should reduce in a wider type first.

Double32s.POS_ZERO

Canonical +0.0, the pair (+0.0f0, +0.0f0).

Double32s.NEG_ZERO

Canonical -0.0, the pair (-0.0f0, +0.0f0).

Double32s.POS_INF

Canonical +Inf, the pair (+Inf32, +0.0f0).

Double32s.NEG_INF

Canonical -Inf, the pair (-Inf32, +0.0f0).

Double32s.D32_NAN

The package's canonical quiet NaN, the pair (NaN32, +0.0f0). NaN32 === reinterpret(Float32, 0x7fc00000).

Double32s.EFT

Error-free transformations on IEEE floats. Reached as Double32s.EFT.two_sum and friends; nothing here enters the namespace of using Double32s.

Trailing-underscore names are the unguarded kernels: branch free, no finiteness check, and the caller owns the preconditions. The names without the underscore add the guard.

two_sum(a::T, b::T) -> (s::T, e::T)

Guarded two_sum_: returns (s, zero(T)) when s is not finite.

two_sum_(a::T, b::T) -> (s::T, e::T)

Knuth's TwoSum. s + e == a + b exactly whenever s is finite, for any a and b. Six add/sub, no branch, no precondition.

two_hilo_sum(a::T, b::T) -> (s::T, e::T)

Guarded two_hilo_sum_.

two_hilo_sum_(a::T, b::T) -> (s::T, e::T)

Dekker's FastTwoSum. Three add/sub instead of six.

Exactness condition. s + e == a + b exactly, for finite s, iff fl(a + b) - a is itself representable. abs(a) >= abs(b) is the familiar sufficient condition, not a necessary one; a == 0, and more generally "a + b is exact", are further sufficient cases that the pair kernels of design section 8 actually rely on.

Stating the sufficient condition as though it were the requirement is what let the section 8.1 justification stand unchallenged for so long. The uses inside add_dd_/sub_dd_ are covered end-to-end by the machine-checked AccurateDWPlusDW bound — see Double32s.add_dd_ — and every other call site pairs this with a magnitude argument (|correction| ≤ u|leading|) in a comment. The condition is not checked here. Use two_hilo_sum_checked in host tests.

two_hilo_sum_checked(a::T, b::T) -> (s::T, e::T)

Differentially-checked FastTwoSum for host code and tests. Not for kernels.

Checks the property that actually matters — agreement with two_sum_ — rather than the merely sufficient abs(a) >= abs(b), which the pair kernels violate legitimately under cancellation. @boundscheck makes the check elidable under @inbounds at the call site.

two_diff(a::T, b::T) -> (s::T, e::T)

Guarded two_diff_.

two_diff_(a::T, b::T) -> (s::T, e::T)

Knuth's TwoDiff. s + e == a - b exactly whenever s is finite.

Written out symmetrically rather than as two_sum_(a, -b): the two are algebraically identical and disagree on the sign of a zero result.

two_prod(a::T, b::T) -> (p::T, e::T)

Guarded two_prod_.

two_prod_(a::T, b::T) -> (p::T, e::T)

FMA-based TwoProduct. p + e == a * b exactly, provided p is finite, the exact residual is representable (no underflow below the least subnormal), and fma is a genuine single-rounding fused multiply-add.

fma, never muladd. Base.muladd lowers to llvm.fmuladd, which permits but does not require fusion. An EFT built on muladd is silently wrong on whichever backend declines to fuse.

nonfinite_pair(s::T) -> Tuple{T,T}

Canonical non-finite EFT result: the non-finite word and a positive zero.