Surfaced while fixing #1317 (PR #1321).
With full-width operands, expansion_ops::multiply_cascades<4>() is about two decimal digits less
accurate than classic qd's hand-tuned multiplication, and roughly 8x slower.
Accuracy
Independent oracle (Python decimal at 160 digits, sharing no code with the library), 200 random
operand pairs where both operands carry all four components (the shape that arises inside every
Taylor series, Newton iteration and Horner evaluation):
operation type worst rel err median worst digits (2^-212 = 1.5e-64)
mul4 qd 1.604e-65 7.823e-67 64.8
mul4 qdc 1.654e-63 4.687e-65 62.8
62.8 digits is inside the format's 63.6, so this is a quality gap rather than a defect - unlike the
addition bug in #1317, which cost 15 digits. Worth noting that operands built from a single
double do not show it: their product is exactly representable and both implementations return
it bit-for-bit. Any test that only multiplies simple values will see nothing.
The likely mechanism is in the accumulation, not the partial products. After sorting the 32-term
expansion, multiply_cascades folds each term into a 4-component result and absorbs whatever falls
below the last component:
// Sub-ULP carry below the last component is precision loss
// beyond what 4 doubles can represent; absorb into result[3].
if (carry != 0.0) {
double s = 0.0, e = 0.0;
two_sum(result[3], carry, s, e);
result[3] = s; // e is discarded
}
Those discarded e terms are individually sub-ulp but there are many of them, and they are
discarded rather than accumulated. Classic qd instead carries a running three-term accumulation
and renormalizes once at the end.
Performance
From benchmark_hp_scalar (i7-12700K, gcc 13.3, -O3, pinned), nsec/op:
operation qd qd_cascade ratio
multiply 73.01 581.26 7.96
multiply_cascades<4> computes 16 two_prods into a 32-term expansion and then bubble sorts
it - up to 496 compare-and-swap steps - before accumulating. The sort is very likely the
dominant cost; classic qd does no sorting at all.
Why it matters
qd_cascade multiply is the single worst row in the #1315 benchmark suite, and #1317 removed the
accuracy blocker to qd_cascade replacing classic qd, leaving this as the main remaining
obstacle. exp (16 squarings) and log (Newton on exp) are 6x qd, essentially all of it here.
Suggested approach
Port classic qd's accurate_multiplication (qd_impl.hpp) to the cascade framework, the way
#1317 ported Shewchuk's COMPRESS: keep the generic template for other widths, add a specialization
for floatcascade<4> that mirrors the proven algorithm. Validate against the same exact-oracle
harness added in static/highprecision/qd_cascade/arithmetic/addition_oracle.cpp - extending it to
multiplication is a few lines, since dyadic multiplication is exact.
Surfaced while fixing #1317 (PR #1321).
With full-width operands,
expansion_ops::multiply_cascades<4>()is about two decimal digits lessaccurate than classic
qd's hand-tuned multiplication, and roughly 8x slower.Accuracy
Independent oracle (Python
decimalat 160 digits, sharing no code with the library), 200 randomoperand pairs where both operands carry all four components (the shape that arises inside every
Taylor series, Newton iteration and Horner evaluation):
62.8 digits is inside the format's 63.6, so this is a quality gap rather than a defect - unlike the
addition bug in #1317, which cost 15 digits. Worth noting that operands built from a single
doubledo not show it: their product is exactly representable and both implementations returnit bit-for-bit. Any test that only multiplies simple values will see nothing.
The likely mechanism is in the accumulation, not the partial products. After sorting the 32-term
expansion,
multiply_cascadesfolds each term into a 4-component result and absorbs whatever fallsbelow the last component:
Those discarded
eterms are individually sub-ulp but there are many of them, and they arediscarded rather than accumulated. Classic
qdinstead carries a running three-term accumulationand renormalizes once at the end.
Performance
From
benchmark_hp_scalar(i7-12700K, gcc 13.3, -O3, pinned), nsec/op:multiply_cascades<4>computes 16two_prods into a 32-term expansion and then bubble sortsit - up to 496 compare-and-swap steps - before accumulating. The sort is very likely the
dominant cost; classic
qddoes no sorting at all.Why it matters
qd_cascademultiply is the single worst row in the #1315 benchmark suite, and #1317 removed theaccuracy blocker to
qd_cascadereplacing classicqd, leaving this as the main remainingobstacle.
exp(16 squarings) andlog(Newton onexp) are 6xqd, essentially all of it here.Suggested approach
Port classic
qd'saccurate_multiplication(qd_impl.hpp) to the cascade framework, the way#1317 ported Shewchuk's COMPRESS: keep the generic template for other widths, add a specialization
for
floatcascade<4>that mirrors the proven algorithm. Validate against the same exact-oracleharness added in
static/highprecision/qd_cascade/arithmetic/addition_oracle.cpp- extending it tomultiplication is a few lines, since dyadic multiplication is exact.