Using ECLiPSe Prolog's lib(ic) I stumbled upon the following problem from David H. Bailey, "Resolving numerical anomalies in scientific computation." which I was referred to by the Unum book. Actually, it is only part of it. First, let me formulate the equation in terms of (is)/2. Also, please note that all these decimal numbers have an exact representation in radix 2 floats (which comprises IEEE):
ECLiPSe Constraint Logic Programming System [kernel]
...
Version 6.2development #21 (x86_64_linux), Wed May 27 20:58 2015
[eclipse 1]: lib(ic).
...
Yes (0.36s cpu)
[eclipse 2]: X= -1, Y = 2, Null is 0.80143857*X+1.65707065*Y-2.51270273.
X = -1
Y = 2
Null = 0.0
Yes (0.00s cpu)
So this is truly 0.0 (no rounding at all). But now the same with $= in place of is:
[eclipse 3]: X= -1, Y = 2, Null $= 0.80143857*X+1.65707065*Y-2.51270273.
X = -1
Y = 2
Null = 2.2204460492503131e-16__2.2204460492503131e-16
Yes (0.00s cpu)
This interval does not contain 0.0. I am aware that interval arithmetic often is a bit too approximate as in:
[eclipse 4]: 1 $= sqrt(1).
Delayed goals:
0 $= -1.1102230246251565e-16__2.2204460492503131e-16
Yes (0.00s cpu)
But at least the equation holds! However, in the first case zero is no longer included. Evidently I have not understood something. I tried also eval/1 but to no avail.
[eclipse 5]: X= -1, Y = 2, Null $= eval(0.80143857*X+1.65707065*Y-2.51270273).
X = -1
Y = 2
Null = 2.2204460492503131e-16__2.2204460492503131e-16
Yes (0.00s cpu)
What is the reason for Null not including 0.0?
(Edit after @jschimpf's surprising answer)
Here is the quotation from the book page 187 that I interpreted as meaning that the numbers are represented exactly (now stroked through).
Use a {3,5}, environment, the one that can simulate IEEE single precision. The input values are exactly representable. ...
{-1, 2}
...
That did the job, computing the exact answer with fewer than half the bits used by ...
Otherwise the statement page 184 holds:
...
0.80143857 x + 1.65707065 y = 2.51270273
The equations certainly look innocent enough. Assuming exact decimal inputs, this
system is solved exactly by x = -1 and y = 2.
Here is it rechecked with SICStus' library(clpq):
| ?- {X= -1,Y=2,
A = 80143857/100000000,
B = 165707065/100000000,
C = 251270273/100000000,
Null = A*X+B*Y-C}.
X = -1,
Y = 2,
A = 80143857/100000000,
B = 33141413/20000000,
C = 251270273/100000000,
Null = 0 ?
yes
So -1, 2 are the exact solutions.
A precise formulation
Here is a reformulation that does not have rounding problems in the input coefficients, still the solution is just -∞...+∞. Thus trivially correct, but not useable.
[eclipse 2]: A = 25510582, B = 52746197, U = 79981812,
C = 80143857, D = 165707065, V = 251270273,
A*X+B*Y$=U,C*X+D*Y$=V.
A = 25510582
B = 52746197
U = 79981812
C = 80143857
D = 165707065
V = 251270273
X = X{-1.0Inf .. 1.0Inf}
Y = Y{-1.0Inf .. 1.0Inf}
Delayed goals:
52746197 * Y{-1.0Inf .. 1.0Inf} + 25510582 * X{-1.0Inf .. 1.0Inf} $= 79981812
80143857 * X{-1.0Inf .. 1.0Inf} + 165707065 * Y{-1.0Inf .. 1.0Inf} $= 251270273
Yes (0.00s cpu)
Several issues conspire here to create confusion:
other than claimed, the three constants in the example do not have exact representations as double floats.
it is not true that the initial example involves no rounding.
the seemingly correct result in the first example is actually due to a fortunate rounding error. Other computation orders give different results.
the exact result, given the nearest double float representation of the constants, is indeed not zero but 2.2204460492503131e-16.
interval arithmetic can only give accurate results when the inputs are accurate, which isn't the case here. The constants have to be widened into intervals that include the desired decimal fraction.
relational arithmetic like the one offered by lib(ic) does by nature not guarantee a particular evaluation order. For that reason the rounding errors can differ from the ones encountered during functional evaluation. The results will however be accurate with respect to the given constants.
The following goes into a bit more detail. As I will demonstrate some points using ECLiPSe queries, a quick word beforehand about the syntax:
two floats separated by a double underscore, such as
0.99__1.01denote an interval constant with lower and upper bound, in this case a number in the vicinity of 1.two integers separated by a single underscore, such as
3_4denote a rational constant with numerator and denominator, in this case three quarters.To demonstrate point (1), convert the float representation of 0.80143857 into a rational. This gives the precise fraction 3609358445212343/4503599627370496, which is close, but not identical, to the intended decimal fraction 80143857/100000000. The floating point representation is therefore not exact:
The following shows how the result depends on the evaluation order (point 3 above; note that I have simplified the original example by getting rid of the irrelevant multiplications):
The order dependence proves that rounding errors occur (point 2). For those familiar with floating point operations, it is in fact easy to see that when adding
-0.80143857 + 3.3141413, two bits of precision from0.80143857get lost while adjusting the exponents of the operands. In fact it is this lucky rounding error that gives the OP his seemingly correct result!In reality, the second result is more accurate with respect to the floating point representations of the constants. We can prove this by repeating the computations using precise rational arithmetic:
As the additions are done with precise rationals, the result is now order-independent, and because
1_4503599627370496 =:= 2.2204460492503131e-16, this confirms the non-zero floating point result obtained above (point 4).How can interval arithmetic help here? It works by computing with intervals that enclose the true value, such that results will always be accurate with respect to the inputs. So it is crucial to have input intervals (bounded reals in ECLiPSe terminology) that enclose the desired true value. These can be obtained by either writing them down explicitly, such as
0.80143856__0.80143858; by converting from a precise number, for instance a rational usingbreal(80143857_100000000); or by instructing the parser to automatically widen all floating point numbers into bounded real intervals, as follows:Both results now contains zero, and it becomes apparent how the precision of the result depends on the evaluation order.