stillwater-sc / universal / 30049709164

85%

master: 84%

Build Type

push

github

Committed by web-flow

Commit Message

feat(qd): implement the declared-but-missing fma(a,b,c) (#1190) (#1200)

qd's fma was declared (qd_fwd.hpp) but never defined -- calling it was a link
error. Implement it as a genuine single-rounding fused multiply-add:

  - the product a*b is formed EXACTLY as the expansion of all sixteen a[i]*b[j]
    partials (each split by two_prod -> 32 doubles), c's four limbs are appended,
    and the exact 36-double bag is renormalized to a quad-double in ONE pass
    (Shewchuk grow-then-compress + the qd 5->4 renorm). Because the product is
    never rounded before c is added, the low-order product bits that a plain
    `a*b + c` discards survive into the result -- which is what matters under
    cancellation (c ~ -a*b).
  - IEEE-754 special-value handling done explicitly (nan; inf*0 -> nan;
    inf + (-inf) -> nan; else the signed product-infinity), because qd's own
    operator+ does not always propagate infinities.

Validated against an independent exact dyadic-rational oracle
(static/highprecision/qd/math/fma.cpp, target qd_math_fma):
  - correctly-rounded: full-precision random operands agree with the exact value
    to >= 60 decimal digits (qd carries ~63); worst observed 62 over 200k.
  - fused: under cancellation fma beats the two-rounding a*b+c in 20000/20000
    cases (naive rounds the residual to ~0 digits; fma keeps it).
  - identities fma(a,1,c)==a+c, fma(a,b,0)==a*b, fma(a,0,c)==c, + non-finite.
Green on gcc (CMake) and clang (-Wall -Wpedantic, no warnings).

Resolves #1190. Part of epic #1189.

Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Coverage Stats

38 of 38 new or added lines in 1 file covered. (100.0%)

18 existing lines in 2 files now uncovered.

50135 of 59136 relevant lines covered (84.78%)

7017491.34 hits per line