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
| operation | domain | samples | median (ulp48) | worst (ulp48) | faithful |
|---|---|---|---|---|---|
+ | mixed pools | 400 | 0.000 | 0.500 | 100.00% |
+ | cancellation | 400 | 0.000 | 0.000 | 100.00% |
- | mixed pools | 400 | 0.000 | 0.750 | 100.00% |
* | mixed pools | 400 | 0.000 | 0.687 | 100.00% |
/ | mixed pools | 400 | 0.083 | 1.500 | 99.25% |
inv | mixed pools | 400 | 0.062 | 1.000 | 100.00% |
sqrt | positive | 400 | 0.064 | 0.750 | 100.00% |
abs2 | mixed pools | 400 | 0.000 | 0.957 | 100.00% |
hypot | mixed pools | 400 | 0.026 | 1.311 | 99.50% |
Documented decisions behind these numbers (measured this development cycle, regressions in test/core/arithmetic.jl):
sqrtis 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.
| operation | samples | worst (u²·M) | bound |
|---|---|---|---|
fma / muladd | 1600 | 2.0000 | 2.0000... (derived: 2u²(1+80u)·M) |
Reductions
| comparison | n | Double32 error (ulp48) | naive error (ulp48) |
|---|---|---|---|
sum(Double32, A) vs sum(::Vector{Float32}) | 100000 | 1.070 | 3614648.9 |
compensated_dot(x, y) vs Float32 dot | 50000 | 15.211 | 16244311.3 |
norm2 (squares overflow Float32) | 10000 | 3.412 | Inf (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.
| function | worst ulp48 | unfaithful | operands |
|---|---|---|---|
cbrt | 0.894 | 7 | 786 |
rsqrt | 1.235 | 39 | 786 |
exp | 0.982 | 65 | 1038 |
expm1 | 0.984 | 30 | 1133 |
exp2 | 0.966 | 69 | 1200 |
exp10 | 0.936 | 56 | 1032 |
log | 0.982 | 43 | 784 |
log1p | 0.902 | 40 | 814 |
log2 | 0.838 | 16 | 784 |
log10 | 1.164 | 86 | 784 |
sin | 1.081 | 88 | 1500 |
cos | 1.001 | 78 | 1500 |
tan | 2.798 | 448 | 1500 |
atan | 1.093 | 37 | 786 |
asin | 1.705 | 230 | 1200 |
acos | 1.598 | 113 | 1200 |
sinh | 1.454 | 64 | 1064 |
cosh | 0.907 | 64 | 1064 |
tanh | 1.282 | 8 | 1214 |
asinh | 1.338 | 34 | 800 |
acosh | 0.944 | 44 | 811 |
atanh | 1.031 | 62 | 1214 |
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.
| data | sum(Double32, A) | compensated_sum | sum (Float32) | sum(Float64, A) |
|---|---|---|---|---|
| 1e8 then 1e5 ones | 0.00e+00 | 0.00e+00 | 4.00e-07 | 0.00e+00 |
| 1/k, k = 1..1e5 | 0.00e+00 | 0.00e+00 | 1.47e-07 | 0.00e+00 |
| randn, 1e5 | 2.29e-15 | 2.29e-15 | 1.61e-07 | 0.00e+00 |
| randn x 2^(-60..60), 1e5 | 3.07e-14 | 3.07e-14 | 1.92e-08 | 1.20e-16 |
Dot product and norm (50 000 elements)
| quantity | Double32 | Float32 | Float64 |
|---|---|---|---|
compensated_dot(x, y) | 2.74e-13 | 1.15e-06 | 2.90e-17 |
norm2(x) | 1.20e-13 | 5.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.
| product | worst 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.