Wiki Agenda Contact Version française

Floating-point square root using Newton iteration

Computing 43 bits of a square root by using three floating-point iterations of Newton's algorithm.


Authors: Guillaume Melquiond

Topics: Floating-Point Computations / Arithmetic / Non-linear Arithmetic

Tools: Gappa / Frama-C

see also the index (by topic, by tool, by reference, by year)


The code below computes the square root of a floating-point number between 0.5 and 2. (Other inputs can be handled by manipulating the exponent.) It starts from a rough approximation of the inverse square root (6 bits of accuracy), performs three steps of floating-point Newton iteration to refine it, and finally multiplies it by the input and returns it. The Gappa tool formally proves that each iteration almost doubled the accuracy (6, 11, 21, 41) and that the result has 43 bits of accuracy.

The three iterations and their annotations are identical, since they are in fact meant to be the body of a loop. Unfortunately, Gappa requires the code to be unrolled and Frama-C fails to do it correctly, so it was done by hand in this example.

/*@
 requires 0.5 <= x <= 2;
 ensures \abs(\result - 1/\sqrt(x)) <= 0x1p-6 * \abs(1/\sqrt(x));
*/
double sqrt_init(double x);

/*@
 requires 0.5 <= x <= 2;
 ensures \abs(\result - \sqrt(x)) <= 0x1p-43 * \abs(\sqrt(x));
*/
double sqrt(double x)
{
  double t, u;
  t = sqrt_init(x);

  u = 0.5 * t * (3 - t * t * x);
  //@ assert \abs(u - 0.5 * t * (3 - t * t * x)) <= 1;
  /*@ assert (0.5 * t * (3 - t * t * x) - 1/\sqrt(x)) / (1/\sqrt(x)) ==
        - (1.5 + 0.5 * ((t - 1/\sqrt(x)) / (1/\sqrt(x)))) *
        (((t - 1/\sqrt(x)) / (1/\sqrt(x))) * ((t - 1/\sqrt(x)) / (1/\sqrt(x)))); */
  //@ assert \abs(u - 1/\sqrt(x)) <= 0x1p-10 * \abs(1/\sqrt(x));
  t = u;

  u = 0.5 * t * (3 - t * t * x);
  //@ assert \abs(u - 0.5 * t * (3 - t * t * x)) <= 1;
  /*@ assert (0.5 * t * (3 - t * t * x) - 1/\sqrt(x)) / (1/\sqrt(x)) ==
        - (1.5 + 0.5 * ((t - 1/\sqrt(x)) / (1/\sqrt(x)))) *
        (((t - 1/\sqrt(x)) / (1/\sqrt(x))) * ((t - 1/\sqrt(x)) / (1/\sqrt(x)))); */
  //@ assert \abs(u - 1/\sqrt(x)) <= 0x1p-10 * \abs(1/\sqrt(x));
  t = u;

  u = 0.5 * t * (3 - t * t * x);
  //@ assert \abs(u - 0.5 * t * (3 - t * t * x)) <= 1;
  /*@ assert (0.5 * t * (3 - t * t * x) - 1/\sqrt(x)) / (1/\sqrt(x)) ==
        - (1.5 + 0.5 * ((t - 1/\sqrt(x)) / (1/\sqrt(x)))) *
        (((t - 1/\sqrt(x)) / (1/\sqrt(x))) * ((t - 1/\sqrt(x)) / (1/\sqrt(x)))); */
  //@ assert \abs(u - 1/\sqrt(x)) <= 0x1p-10 * \abs(1/\sqrt(x));
  t = u;

  //@ assert x * (1/\sqrt(x)) == \sqrt(x);
  return x * t;
}