Contact Version française

The ProVal team was stopped at the end of August 2012, and reborn into a new team Toccata
These pages do not evolve anymore, please follow the link above for up-to-date informations about our team.

Veltkamp/Dekker algorithm

Veltkamp/Dekker algorithm to compute the exact error of the multiplication with only floating-point operations.


Authors: Sylvie Boldo

Topics: Floating-Point Computations

Tools: Frama-C / Jessie / Gappa

References: NSV-3 Benchmarks

index


The following C function indeed computes the exact error of the multiplication [4, 3] with only floating-point operations (and no FMA). This requires 16 operations and each of them has a specific role to play to get the exact error. This algorithm was already proved with a generic radix and precise underflow restrictions and overflow restrictions so that no infinity will be created [1]. The question was how to specify it as a program to get the full benefit of the proof.

See also [2]

References

[1]
Sylvie Boldo. Pitfalls of a full floating-point proof: example on the formal proof of the Veltkamp/Dekker algorithms. In Ulrich Furbach and Natarajan Shankar, editors, Third International Joint Conference on Automated Reasoning, volume 4130 of Lecture Notes in Computer Science, pages 52–66, Seattle, USA, August 2006. Springer.
[2]
Sylvie Boldo and Claude Marché. Formal verification of numerical programs: from C annotated programs to mechanical proofs. Mathematics in Computer Science, 2011.
[3]
Theodorus J. Dekker. A floating point technique for extending the available precision. Numerische Mathematik, 18(3):224–242, 1971.
[4]
Gerhard W. Veltkamp. Algolprocedures voor het berekenen van een inwendig product in dubbele precisie. RC-Informatie 22, Technische Hogeschool Eindhoven, 1968.

/*@ requires xy == \round_double(\NearestEven,x*y) && 
  @          \abs(x) <= 0x1.p995 && 
  @          \abs(y) <= 0x1.p995 &&
  @          \abs(x*y) <=  0x1.p1021;
  @ ensures  ((x*y == 0 || 0x1.p-969 <= \abs(x*y)) 
  @                 ==> x*y == xy+\result);
  @*/
double Dekker(double x, double y, double xy) {

  double C,px,qx,hx,py,qy,hy,tx,ty,r2;
  C=0x8000001p0;
  /*@ assert C == 0x1p27+1; */

  px=x*C;
  qx=x-px;
  hx=px+qx;
  tx=x-hx;

  py=y*C;
  qy=y-py;
  hy=py+qy;
  ty=y-hy;

  r2=-xy+hx*hy;
  r2+=hx*ty;
  r2+=hy*tx;
  r2+=tx*ty;
  return r2;
}