Function error_bound 

pub fn error_bound(point: &PointType, file: &str) -> Result<ErrorBound, Error>

Derives the error bound of the approximate coordinates of point.

The constructor of a point type evaluates its coordinates in f64, one node of its PointGraph at a time, and with each value v an error e, a bound on |v - V| where V is the exact value of the node. This function derives the program that computes e and the factor that makes the result an upper bound in spite of the rounding of that program itself.

The rule

With u = 2^-53, the rule gives each node a real bound R:

• an input (a coordinate, an f64 parameter or an integer that f64 represents exactly) has R = 0;
• a negation keeps the bound of its operand;
• a sum or difference s = x ± y has R_s = R_x + R_y, plus u·|s| unless an operand is exactly zero;
• a product p = x·y has R_p = |x|·R_y + |y|·R_x + R_x·R_y, plus u·|p| + 2^-1074 unless an operand is exactly zero.

Here x, y, s and p are the computed f64 values. R bounds the error: rounding to nearest makes |fl(x ± y) - (x ± y)| at most u·|fl(x ± y)| (a sum whose result is below 2^-1021 is exact), and |fl(x·y) - x·y| at most u·|fl(x·y)| or, when x·y is below 2^-1022, 2^-1075; an operation with an operand that is exactly zero is exact. The errors of the operands then propagate as written, since |x·y - X·Y| <= |x|·|y - Y| + (|y| + |y - Y|)·|x - X|. The condition “exactly zero” is on the operands, never on the result: a product that underflows gives zero without being exact.

The program

The program computes, for each node, the terms of its rule in the order above and adds them from left to right. Terms that are zero for every input (the error of an input) are left out; the others are computed even when they are zero at run time, which costs nothing in the bound because an operation with a zero operand is exact. Every product of the program is a covered product: fl(fl(a·b) + 2^-1074), or 0 when a or b is zero. The rounding terms are the covered product of u and |s| for a sum, and the covered product of u and |p| plus 2^-1074 for a product, both 0 when an operand of the node is zero.

The rounding of the program

All the values of the program are non-negative. For non-negative a and b, fl(a + b) >= (1 - u)·(a + b): above 2^-1022 by the bound on the relative error, and below it because the sum is exact. A covered product satisfies c >= (1 - u)^2·a·b whatever the magnitudes, which is why the program adds 2^-1074 after each product: if fl(a·b) >= 2^-1022, fl(a·b) >= (1 - u)·a·b and the addition loses at most another factor 1 - u; if fl(a·b) < 2^-1022, possibly zero after underflow, then fl(a·b) >= a·b - 2^-1075, and adding 2^-1074 to a value below 2^-1022 is exact, so c >= a·b. That covers the underflow of the program itself, including a product rounded to zero, and the term is added only when both operands are nonzero, so a coordinate whose terms are all zero still gets an error of exactly zero.

Each node is given a count h of roundings: one per addition and per multiplication of the program on the longest chain that leads to its error, so a covered product counts two. An error e(x) with h(x) such roundings satisfies e(x) >= (1 - u)^h(x)·R(x), by induction over the operations: an input has h = 0; the error of an operand has its own count; a covered product of a value and an error e(y) has h(y) + 2, and of two errors h(x) + h(y) + 2; the rounding term of a sum has 2 and that of a product 3; and adding a term to the sum of the previous ones gives one more than the larger of the two counts.

The factor and the final rounding

With H the largest count of the three coordinates, the slack factor S is the smallest f64 with S·(1 - u)^H >= 1, derived here in exact arithmetic. The constructor then stores next_up(fl(e·S)), which is at least e·S >= R, since the value after the nearest f64 is above the exact product; and 0 when e is 0, which happens only when every term is zero, so only when R is zero and the coordinate is exact. If the coordinate or fl(e·S) is not finite, it stores +∞ instead.

Errors

An error at the name of the type if an integer of its expression is not exactly an f64, or if the computation of the error of a coordinate takes more than MAX_ROUNDS roundings.