Interval arithmetic with outward rounding, and how to check it exactly

An interval [lo, hi] is a pair of doubles with lo <= hi that is known to contain a real number. Every operation here returns an interval that contains every result the exact operation can give on members of its inputs. The integral enclosures of lib/hf/iv_ints are built from these operations, so their soundness rests on this module.

Representation. An interval is a 2-entry array [lo, hi]. A table of intervals is an array whose last axis has length 2, so flattening it puts entry k's ends at 2k and 2k + 1.

Outward rounding without a rounding mode. Viv has no control over the rounding mode, and both engines round to nearest. So each operation computes the rounded result r and then finds out on which side of r the exact result lies, from its error term:

• a + b: Knuth's TwoSum gives e with a + b = s + e exactly.
• a * b: Dekker's TwoProduct (Veltkamp split by 2^27 + 1) gives e with a * b = p + e exactly.
• a / b: with q = fl(a / b), the residual r = a - b q has the sign of (a / b - q) b. b q = P + E by TwoProduct, a - P is exact because P is within a few ulps of a (Sterbenz), and a rounded (a - P) - E has the sign of the exact one.
• sqrt(a): with s = fl(sqrt(a)), a - s^2 has the sign of sqrt(a) - s, computed the same way.

When the error is zero the result is exact and both ends are r. Otherwise one end is r and the other is r moved one step outward. The step is Rump, Zimmermann, Boldo and Melquiond's successor in round-to-nearest, c + (phi |c| + eta) with phi = 2^-53 (1 + 2^-52) and eta = 2^-1074, which gives the next double or the one after it, never c itself. A result whose inputs are exact doubles is therefore at most one or two ulps wide.

Dekker's product and the residuals are exact only when nothing underflows, so below IV_TINY = 1e-290 an operation skips the error term and moves both ends out. Overflow and NaN are not handled; nothing in lib/hf comes near them.

Why the error terms are exact in both engines. The generated C starts with #pragma STDC FP_CONTRACT OFF, and a probe where a fused multiply-add would have given -2^-60 printed 0 in both engines. Without that, the C compiler could fuse ah * bh - p and the product error would be wrong.

The libm assumption. exp is the one operation whose rounding is not under our control. iv_exp assumes the library exp is within 2 ulps of the exact value and widens the result by 4 units of 2^-53 relative on each side, plus one outward step. Every Viv engine takes exp from the same code, fdlibm's e_exp (lang/viv/src/vmath), whose error analysis gives under 1 ulp. This is an assumption, not a proof: that bound is fdlibm's own analysis, not something this library has checked.

How the tests check containment exactly. An interval can be checked against the exact result without trusting interval code. The exact result of a + b or a * b is the unevaluated sum s + e of two doubles, and lo <= s + e is decided by the sign of the three-term sum s + e - lo. That sign is found exactly by Shewchuk's expansion growth (iv_esign): TwoSum turns any list of doubles into a nonoverlapping expansion with the same exact sum, and the sign of a nonoverlapping expansion is the sign of its largest component. For / and sqrt the tests check the residual at each end: a - b lo and a - b hi must have the signs that put a / b between them, and a - lo^2 >= 0 >= a - hi^2.

Used with import hf.interval.