Fixed-point Newton division

The algorithm and its verification

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)}]

Adding hints

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.

Full listing

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”.