stillwater-sc / universal / 32093743701

85%

master: 84%

push

github

fix(qd): bound the sqrt intermediates across the range (#1332) (#1339)

sqrt(maxpos) returned inf or NaN in dd and qd. Both algorithms square a value of
the argument's own magnitude -- dd's Karp trick forms (a*x)^2, qd's iteration
forms r^2 with r ~ 1/sqrt(a) and closes with r*a -- so within a rounding of
maxpos the square leaves the representable range.

Both now scale the argument into [0.5, 2) before iterating and scale the result
back by half the exponent, which is what the three cascade types have done since
#1331. The scaling is by a power of two, so it is exact and costs no accuracy,
and it bounds every intermediate wherever the argument sits.

Scoring sqrt end to end over the exponent range against an exact oracle showed
the defect was wider than the overflow the issue reported. The same squares that
overflow at the top go subnormal at the bottom, and the trailing limbs of the
result go with them:

              worst over the range      non-finite results
    dd        52 -> 105 bits of 106     0 -> 0
    qd        60 -> 214 bits of 212     28 -> 0

sqrt(1e-300) in a dd held 80 bits of its 106 before this, and sqrt(1e300) in a
qd held 81 of its 212.

qd::sqrt(0) also returned NaN: the seed is 1/sqrt(a[0]), which is inf at zero,
and qd never guarded the case. dd always had. Both now guard zero, NaN and
infinity.

The new suite walks the exponent range end to end for all five multi-component
types -- the cascades included, since #1331 is what keeps them correct. It
checks the residual rather than a reference value: sqrt is irrational and has no
exact reference, but squaring the result is exact, because a and r are sums of
doubles and dyadic rationals are closed under multiplication. |r*r - a| <=
2^-budget * a is decided in exact integer arithmetic. The budgets sit four bits
under each significand, calibrated from the measured worst residual over ~2,800
arguments per type. Against the implementation before this change it fails 14
cases in dd and... (continued)

18 of 18 new or added lines in 2 files covered. (100.0%)

43011 of 50366 relevant lines covered (85.4%)

6901397.81 hits per line