fix(qd_cascade): compress the expansion before renormalizing addition - #1321
Merged
Conversation
…#1317) qd_cascade sqrt/exp/log delivered 48-50 correct decimal digits where the format carries 63.6 and where classic qd delivers all of them. The cause was in none of those functions. expansion_ops::add_cascades() merges the two operands' components by magnitude, accumulates them with two_sum, and returns the running sum plus the collected errors. The VALUE of that expansion is exact, but its components overlap - consecutive terms can be within a factor of two of each other rather than a factor of 2^53: e[0] = 1.743e+00 e[1] = 1.110e-16 <- e[2] is 3x smaller, not 2^53x smaller e[2] = 3.385e-17 compress_8to4() is the QD library's renorm, built entirely from fast_two_sum chains, and fast_two_sum is only error-free when |a| >= |b|. Handed an overlapping expansion it runs out of places to put the tail and returns a result whose fourth component is exactly zero. That is a relative error of 2^-160 in a format carrying 2^-212, on roughly one addition in eight, and everything built on addition inherited it - which is why it surfaced in the 16 squarings of exp() rather than in addition's own test suite. The fix is expansion_ops::compress_expansion<M>(), Shewchuk's COMPRESS (Geometric Predicates, Fig. 22): a downward sweep carrying each rounding remainder into the next component, then an upward sweep leaving the leading components correctly rounded. It is value-preserving - it changes only the SHAPE of the expansion, which is exactly what the precondition is about. Applied to the floatcascade<4> overload only. The <2> and <3> overloads were measured over 5000 random full-precision additions each and stay within 0.49 ulp of their own format with and without it: fewer merged terms compressed into fewer components leaves enough slack that the overlap never reaches a surviving component. Compressing them anyway would have doubled dd_cascade addition for no accuracy gain, so both carry a comment saying so. Verified against an exact oracle (Python decimal at 140 digits, no shared code), 200 random full-precision operand pairs, worst relative error: before after classic qd add 1.7e-49 exact 6.2e-65 sqrt 1.3e-49 6.1e-64 5.5e-65 exp 4.8e-51 4.0e-65 4.0e-65 log 1.3e-50 2.3e-63 2.3e-63 exp and log now match qd to the last bit. Cost: quad-double addition is 46% more expensive (48.4 -> 70.4 nsec/op), and that reverses one of the #1315 benchmark's conclusions - qd_cascade's dot product was 0.83x qd only because the addition was dropping a component; with the correct addition it is 1.15x. The benchmark README is re-measured and its conclusions updated accordingly. Tests: static/highprecision/{qd,td}_cascade/arithmetic/addition_oracle.cpp validate addition against exact dyadic-rational arithmetic (dyadic_exact.hpp, backed by einteger), which shares no code with the expansion arithmetic. The qd_cascade suite fails 825 of 4096 cases at REGRESSION_LEVEL_4 on the old code and passes on the new. The td_cascade suite is a guard for the overload deliberately left uncompressed. Also corrects the 2026-03-16 PRECISION NOTE in qd_cascade/exponent.hpp, which attributed this loss to multiply_cascades(). Multiplication is ~2 digits behind qd's on full-width operands (1.7e-63 vs 1.6e-65) - real, inside the format, and tracked separately. 113 dd/qd/cascade regression tests pass under gcc 13.3 and clang 18.1. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Contributor
|
Caution Review failedThe pull request is closed. ℹ️ Recent review info⚙️ Run configurationConfiguration used: Path: .coderabbit.yaml Review profile: CHILL Plan: Pro Run ID: 📒 Files selected for processing (6)
📝 WalkthroughWalkthroughThe change adds two-pass compression for Changesqd_cascade correctness
Estimated code review effort: 4 (Complex) | ~45 minutes Possibly related issues
Possibly related PRs
Suggested labels: ✨ Finishing Touches📝 Generate docstrings
🧪 Generate unit tests (beta)
Comment |
Ravenwater
marked this pull request as ready for review
August 15, 2026 23:00
Coverage Report for CI Build 31918517357Warning Build has drifted: This PR's base is out of sync with its target branch, so coverage data may include unrelated changes. Coverage decreased (-0.03%) to 85.5%Details
Uncovered ChangesNo uncovered changes found. Coverage Regressions31 previously-covered lines in 3 files lost coverage.
Coverage Stats
💛 - Coveralls |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #1317.
qd_cascadesqrt/exp/logdelivered 48-50 correct decimal digits where the format carries63.6 and classic
qddelivers all of them. The cause was in none of those functions.Root cause
expansion_ops::add_cascades()merges the operands' components by magnitude, accumulates themwith
two_sum, and returns the running sum plus the collected errors. The value of thatexpansion is exact, but its components overlap - consecutive terms can be within a factor of
two of each other rather than a factor of 2^53:
compress_8to4()is the QD library'srenorm, built entirely fromfast_two_sumchains, andfast_two_sum(a, b)is only error-free when|a| >= |b|. Handed an overlapping expansion it runsout of places to put the tail and returns a result whose fourth component is exactly zero - a
relative error of 2^-160 in a format carrying 2^-212, on roughly one addition in eight.
Everything built on addition inherited it, which is why it surfaced in the 16 squarings of
exp()rather than in addition's own test suite. The existing suites all check self-consistency, and a
result that is normalized, round-trips, and matches the other implementation of the same flawed
algorithm passes every one of them.
Fix
expansion_ops::compress_expansion<M>()- Shewchuk's COMPRESS (Geometric Predicates, Fig. 22):a downward sweep carrying each rounding remainder into the next component, then an upward sweep
leaving the leading components correctly rounded. It is value-preserving; it changes only the
shape of the expansion, which is exactly what the precondition is about.
Applied to the
floatcascade<4>overload only. The<2>and<3>overloads were measuredover 5000 random full-precision additions each and stay within 0.49 ulp of their own format with
and without it - fewer merged terms compressed into fewer output components leaves enough slack
that the overlap never reaches a surviving component. Compressing them anyway would have doubled
dd_cascadeaddition for no accuracy gain, so both carry a comment recording that.Verification
Independent oracle (Python
decimalat 140 digits, sharing no code with the library), 200 randomfull-precision operand pairs, worst relative error:
qdexpandlognow matchqdto the last bit.Cost, and a reversed conclusion
Quad-double addition is 46% more expensive (48.4 -> 70.4 nsec/op). That reverses one of the #1315
benchmark's findings:
qd_cascade's dot product measured 0.83xqdonly because the additionwas dropping a component; with the correct addition it is 1.15x. The benchmark README is
re-measured and its tables and conclusions updated accordingly - including the "can we retire
classic qd" answer, whose accuracy blocker is now gone and whose remaining blockers are purely
speed.
Tests
static/highprecision/{qd,td}_cascade/arithmetic/addition_oracle.cppvalidate addition againstexact dyadic-rational arithmetic (
dyadic_exact.hpp, backed byeinteger) - no shared code withthe expansion arithmetic under test, so it cannot be fooled by a result that is well-formed and
wrong. The
qd_cascadesuite fails 825 of 4096 cases atREGRESSION_LEVEL_4on the old codeand passes on the new; the
td_cascadesuite guards the overload deliberately left uncompressed.Also corrects the 2026-03-16
PRECISION NOTEinqd_cascade/exponent.hpp, which attributed thisloss to
multiply_cascades(). Multiplication is ~2 digits behindqd's on full-width operands(1.7e-63 vs 1.6e-65) - real, inside the format, and worth tracking separately.
113 dd/qd/cascade regression tests pass under gcc 13.3.0 and clang 18.1.3.
🤖 Generated with Claude Code
Summary by CodeRabbit
Bug Fixes
Tests
Documentation