Let us suppose we want to invert a floating-point number on a processor without a floating-point unit. The 24-bit mantissa has to be inverted from a value between 0.5 and 1 to a value between 1 and 2. For the sake of this example, the transformation is performed by Newton's iteration with fixed-point arithmetic.

The mantissa is noted `d`

and its exact reciprocal is
`R`

. Newton's iteration is started with a first approximation
`r0`

taken from a table containing reciprocals at precision
[math]. Two iterations are then performed. The result
`r1`

of the first iteration is computed on 16-bit words in order to
speed up computations. The result `r2`

of the second iteration is
computed on full 32-bit words. We want to prove that this second result is close
enough to the infinitely precise reciprocal `R = 1/d`

.

First, we define `R`

as the reciprocal, and
`d`

and `r0`

as two fixed-point numbers that are
integer multiples of [math] and
[math] respectively. Moreover, `r0`

is
an approximation of `R`

and `d`

is between 0.5 and
1.

R = 1 / d; { @FIX(d,-24) /\ d in [0.5,1] /\ @FIX(r0,-8) /\ r0 - R in [-1b-8,1b-8] -> ... }

Next we have the two iterations. Gappa's representation of fixed-point arithmetic is high-level: the tool is only interested in the weight of the least significant bit. The shifts that occur in an implementation only have an impact on the internal representation of the values, not on the values themselves.

r1 fixed<-14,dn>= r0 * (2 - fixed<-16,dn>(d) * r0); r2 fixed<-30,dn>= r1 * (2 - d * r1);

The property we are looking for is a bound on the absolute error
between `r2`

and `R`

.

{ ... -> r2 - R in ? }

We expect Gappa to prove that `r2`

is [math]. Unfortunately, this is not the case.

Results: r2 - R in [-1320985b-18 {-5.03916, -2^(2.33318)}, 42305669b-23 {5.04323, 2^(2.33435)}]

With the previous script, Gappa computes a range so wide for ```
r2 -
R
```

that it is useless. This is not surprising: The tool does not know what
Newton's iteration is. In particular, Gappa cannot guess that such an iteration
has a quadratic convergence. Testing for `r1 - R`

instead does not
give results any better.

Gappa does not find any useful relation between `r1`

and
`R`

, as the first one is a rounded multiplication while the second one is
an exact division. So we have to split the absolute error into two terms: a
round-off error we expect Gappa to compute, and the convergence due to Newton's
iteration.

{ ... -> r1 - r0 * (2 - d * r0) in ? /\ r0 * (2 - d * r0) - R in ? }

Gappa now gives the answer below. Notice that the range of the round-off error almost matches the precision of the computations.

Results: r1 - r0 * (2 - d * r0) in [-1b-14 {-6.10352e-05, -2^(-14)}, 788481b-32 {0.000183583, 2^(-12.4113)}] r0 * (2 - d * r0) - R in [-131585b-16 {-2.00783, -2^(1.00564)}, 131969b-16 {2.01369, 2^(1.00984)}]

So Gappa computes correct bounds for the round-off error, but not for the algorithmic one. We can help Gappa by providing an expression of the latter one. So we add a rule describing the quadratic convergence of Newton's iteration:

r0 * (2 - d * r0) - R -> (r0 - R) * (r0 - R) * -d; r1 * (2 - d * r1) - R -> (r1 - R) * (r1 - R) * -d;

Gappa answers that [math].

Warning: although present in a quotient, the expression (d) may have not been tested for non-zeroness. Results: r2 - R in [-638882156545b-64 {-3.46339e-08, -2^(-24.7832)}, 32771b-44 {1.86282e-09, 2^(-28.9999)}]

While the answer is the expected one, there is this warning message
about `d`

possibly
being zero. Indeed, `R`

is the reciprocal of `d`

and we
are using the fact that `R * d = 1`

. So the rewriting rules cannot
be proved on their own. (But they can be proved in the context of the
problem, so there is no correctness issue.) In order to eliminate this
warning, we can give the precise hypotheses such that the left hand sides
of the rewriting rules are equal to their right hand sides without any
other assumption. This is indicated at the end of the rule.

r0 * (2 - d * r0) - R -> (r0 - R) * (r0 - R) * -d { d <> 0 };

When generating a script for an external proof checker, Gappa will add
this rewriting rule as a global hypothesis. For example, when selecting the Coq
back-end with the option `-Bcoq`

, the output contains the line
below.

Hypothesis a1 : (_d <> 0)%R -> r9 = r2.

In this hypothesis, `_d`

is the `d`

variable of the
example, while `r9`

and `r2`

are short notations for
`r0 * (2 - d * r0) - R`

and `(r0 - R) * (r0 - R) * -d`

respectively. In order to access the generated proof, the user has to prove
this hypothesis, which can be trivially done with Coq's `field`

tactic.

R = 1 / d; r1 fixed<-14,dn>= r0 * (2 - fixed<-16,dn>(d) * r0); r2 fixed<-30,dn>= r1 * (2 - d * r1); { @FIX(d,-24) /\ d in [0.5,1] /\ @FIX(r0,-8) /\ r0 - R in [-1b-8,1b-8] -> r2 - R in ? } r0 * (2 - d * r0) - R -> (r0 - R) * (r0 - R) * -d { d <> 0 }; r1 * (2 - d * r1) - R -> (r1 - R) * (r1 - R) * -d { d <> 0 };

The answer is the same as before, since Gappa easily proves that
`d`

is not zero.

Results: r2 - R in [-638882156545b-64 {-3.46339e-08, -2^(-24.7832)}, 32771b-44 {1.86282e-09, 2^(-28.9999)}]

Another example of a Newton iteration is given in the section called “Why and Gappa”.