Double32s Accuracy Report

Generated by docs/Reports/accuracies.jl on Julia 1.12.6.

Errors are measured against a 400-bit BigFloat reference in nominal 48-bit ulps of the true value: ulp48(V) = 2^(floor(log2|V|)-47), the design's section 9.1 metric. "Faithful" is the fraction of samples with error strictly below 1 ulp48. Results are charged against the 48-bit lattice only where it exists (|V| ≥ 2^-100); below that the format tapers toward 24 bits by construction (design section 5.5) and accuracy statements belong to the format, not the operations.

Gradual underflow honored on this host: true (Double32s.has_gradual_underflow()). All subnormal-range claims are conditional on it.

Arithmetic

operationdomainsamplesmedian (ulp48)worst (ulp48)faithful
+mixed pools4000.0000.500100.00%
+cancellation4000.0000.000100.00%
-mixed pools4000.0000.750100.00%
*mixed pools4000.0000.687100.00%
/mixed pools4000.0831.50099.25%
invmixed pools4000.0621.000100.00%
sqrtpositive4000.0640.750100.00%
abs2mixed pools4000.0000.957100.00%
hypotmixed pools4000.0261.31199.50%

Documented decisions behind these numbers (measured this development cycle, regressions in test/core/arithmetic.jl):

  • sqrt is empirically faithful — worst 0.75 ulp48, zero violations, including near-perfect-square adversaries — with a single Newton correction.
  • / is not faithful and is documented as such: worst 4.00 ulp48 over 101,761 lattice-range quotients. A second correction (implemented, measured: 2.30 ulp48) still does not reach faithfulness and costs 2× in GPU kernels, so every policy uses one correction.
  • inv: worst 1.16 ulp48, one violation per few hundred structured operands.

Compensated fma

Bounded against the operand scale M = max(|x·y|, |z|), in units of u²·M, u = 2⁻²⁴; the derived bound is 2u²(1+80u)·M and it is tight.

operationsamplesworst (u²·M)bound
fma / muladd16002.00002.0000... (derived: 2u²(1+80u)·M)

Reductions

comparisonnDouble32 error (ulp48)naive error (ulp48)
sum(Double32, A) vs sum(::Vector{Float32})1000001.0703614648.9
compensated_dot(x, y) vs Float32 dot5000015.21116244311.3
norm2 (squares overflow Float32)100003.412Inf (Float32 overflows)

Coverage

Every implemented Double32 surface is measured here: scalar arithmetic, compensated fma, reductions, elementary functions, trigonometry, hyperbolics, and dense matrix multiply. The KernelAbstractions reductions and the tiled device matmul are measured separately in gpu_report.md, which has the GPU packages available.

Elementary functions (roadmap Phase 6)

Worst error in nominal-48-bit ulps against a 400-bit BigFloat reference, restricted to results in [2^-100, 2^100] where the 48-bit lattice exists. "unfaithful" counts results more than half a ulp48 from the true value. "operands" is how many arguments from the test pool actually produced a comparable result — the count after discarding those outside that range.

functionworst ulp48unfaithfuloperands
cbrt0.8947786
rsqrt1.23539786
exp0.982651038
expm10.984301133
exp20.966691200
exp100.936561032
log0.98243784
log1p0.90240814
log20.83816784
log101.16486784
sin1.081881500
cos1.001781500
tan2.7984481500
atan1.09337786
asin1.7052301200
acos1.5981131200
sinh1.454641064
cosh0.907641064
tanh1.28281214
asinh1.33834800
acosh0.94444811
atanh1.031621214

sinh, cosh, asinh, acosh, and atanh are each written in the non-cancelling form; cosh is the exception that proves the rule, since it never cancels and the expm1 form measured worse (3.6 against 1.1).

None of these is faithfully rounded, and none is claimed to be. The floor for anything ending in a division is this package's /, measured at 4.00 ulp48 — which is why tan, asin, and acos sit highest. log1p and log2/log10 were each reworked once to avoid a division that dominated their error; see the Phase 6 entry in docs/Checkpoint.md.

sin, cos, tan, and sincos return NaN beyond Double32s.TRIG_ARG_LIMIT = 2^20, where the Cody–Waite reduction can no longer be trusted. That is a deliberate policy, not a range over which the table above applies.

Reductions in detail

Relative error of the sum against a 400-bit BigFloat reference. sum(Double32, A) accumulates Float32 data in ~48 bits through Base's own mapreduce; compensated_sum pins a sequential order.

The exact zeros are real, and worth understanding. When every element is a Float32 and their exponents span a narrow range, the exact sum needs only about 24 bits plus that span — for 1/k over 1e5 terms that is roughly 44 bits, which fits inside Double32's 48. Those sums are not merely accurate, they are exact. The last row spreads the exponents over 2^-60..2^60 so the exact sum needs far more than 48 bits, and there the format's true error appears.

datasum(Double32, A)compensated_sumsum (Float32)sum(Float64, A)
1e8 then 1e5 ones0.00e+000.00e+004.00e-070.00e+00
1/k, k = 1..1e50.00e+000.00e+001.47e-070.00e+00
randn, 1e52.29e-152.29e-151.61e-070.00e+00
randn x 2^(-60..60), 1e53.07e-143.07e-141.92e-081.20e-16

Dot product and norm (50 000 elements)

quantityDouble32Float32Float64
compensated_dot(x, y)2.74e-131.15e-062.90e-17
norm2(x)1.20e-135.08e-08—

compensated_dot forms every Float32 product exactly as a Double32 before accumulating, so nothing is lost before the summation.

Matrix multiply (128 x 128, Float32 data)

Worst normwise error: max|C - Cref| / max|Cref|. An elementwise relative error would instead measure the conditioning of the near-zero entries randn data produces by cancellation — it reported 2.0e-2 for BLAS Float32, which is not a real result.

productworst relative error
Double32.(A) * Double32.(B) (generic)1.83e-14
A * B (Float32, BLAS)3.94e-07

The device-side tiled kernel Double32s.matmul, which forms every Float32 product exactly with two_prod_, is measured in gpu_report.md.