Clever Geek Handbook

Exact solutions for lib (ic)

Using ECLiPSe Prolog lib(ic)

I came across the following problem from David H. Bailey, “Solving Numerical Anomalies in Scientific Computing”. which I mentioned in the book Unum. Actually, this is only part of it. First, let me put the equation in terms (is)/2

. Also, note that all of these decimal numbers are exactly represented in floix 2 fields (including 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 that's true 0.0 ( no rounding at all ). But now it's the same with $=

instead 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 know that interval arithmetic is often a little approximate, as in:

[eclipse 4]: 1 $= sqrt(1).

Delayed goals:
    0 $= -1.1102230246251565e-16__2.2204460492503131e-16
Yes (0.00s cpu)

But at least there is equality! However, in the first case, zero is no longer included. Obviously, I misunderstood something. I have also tried 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 unexpected answer)

Here is a quote from the book on page 187, which I interpreted to mean that the numbers are represented accurately (now stroked through ).

Use the {3,5} environment, an environment that can mimic IEEE single precision. The input values ​​are accurately represented ....
{-1, 2}


This did the job, calculating the exact answer with less than half the bits used ...

Otherwise, the request page 184 is executed:

...

     

0.80143857 x + 1.65707065 y = 2.51270273

The equations certainly look innocent. Assuming exact decimal inputs, this system is solved exactly with x = -1 and y = 2.

Here it is 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 exact solutions.

Exact wording

Here is a reformulation that has no rounding problems in the input coefficients, yet the solution is correct - & infin; ... + & infin ;. So trivially correct but not used.

[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)
+3


source to share


1 answer


Several problems arise here to create confusion:

  • apart from the declared ones, the three constants in the example do not have exact representations in the form of double floats.

  • it is not true that the initial example does not require rounding.

  • The seemingly correct result in the first example is actually related to a lucky round-off error. Other settlement orders produce different results.

  • the exact result given the closest double float representation of the constant is indeed not zero, but 2.2204460492503131e-16.

  • interval arithmetic can only give accurate results when the inputs are accurate, which is not the case here. Constants must be expanded into ranges that include the desired decimal fraction.

  • relational arithmetic like the one offered by lib (ic) does not, by its very nature, guarantee a specific order of evaluation. For this reason, rounding errors may differ from those encountered during functional evaluation. However, the results will be accurate for the given constants.

See below for more details. As I'll show some using ECLiPSe queries, a quick word in advance about the syntax:

  • two floats separated by double underlining, for example, 0.99__1.01

    denote an interval constant with a lower and upper limit, in this case a number in the vicinity of 1.

  • two integers separated by one underscore, for example, 3_4

    denote a rational constant with a numerator and a denominator, in this case three quarters.

To demonstrate point (1), convert the float representation 0.80143857 to rational. This gives the exact fraction 3609358445212343/4503599627370496, which is close, but not identical, to the implied decimal 80143857/100000000. Floating point so the representation is not exact:

?- F is rational(0.80143857), F =\= 80143857_100000000.
F = 3609358445212343_4503599627370496
Yes (0.00s cpu)

Below is how the result depends on the order of evaluation (point 3 above, note that I simplified the original example to get rid of non-local multiplications):

?- Null is -0.80143857 + 3.3141413 - 2.51270273.
Null = 0.0
Yes (0.00s cpu)

?- Null is -2.51270273 + 3.3141413 - 0.80143857.
Null = 2.2204460492503131e-16
Yes (0.00s cpu)


The order dependency proves round-off errors occur (point 2). For those familiar with floating point operations, it is actually easy to see that when added, -0.80143857 + 3.3141413

two bits of precision from 0.80143857

are lost when adjusting the exponents of the operands. It is actually this lucky rounding error that gives the OP his seemingly correct result!

In fact, the second result is more accurate in terms of floating point representation of constants. We can prove this by repeating the calculations using precise rational arithmetic:

?- Null is rational(-0.80143857) + rational(3.3141413) - rational(2.51270273).
Null = 1_4503599627370496
Yes (0.00s cpu)

?- Null is rational(-2.51270273) + rational(3.3141413) - rational(0.80143857).
Null = 1_4503599627370496
Yes (0.00s cpu)

Since the additions are performed with exact rationalities, the result is now order independent, but because, 1_4503599627370496 =:= 2.2204460492503131e-16

it confirms the nonzero floating point result obtained above (point 4).

How can interval arithmetic help? It works by calculating using intervals that enclose the true value, so that the results will always be accurate with respect to the inputs. Therefore, it is imperative to have input intervals (limited verifications in ECLiPSe terminology) that enclose the desired true meaning. You can get them by writing them explicitly, for example 0.80143856__0.80143858

; by converting from an exact number, such as rational use breal(80143857_100000000)

; or by automatically instructing the parser to expand all floating point numbers to limited real intervals as follows:

?- set_flag(syntax_option, read_floats_as_breals).
Yes (0.00s cpu)

?- Null is -0.80143857 + 3.3141413 - 2.51270273.
Null = -8.8817841970012523e-16__1.3322676295501878e-15
Yes (0.00s cpu)

?- Null is -2.51270273 + 3.3141413 - 0.80143857.
Null = -7.7715611723760958e-16__1.2212453270876722e-15
Yes (0.00s cpu)

Both results now contain zero, and it becomes obvious how the accuracy of the result depends on the order of the estimate.

+5


source






More articles:


All Articles