/*@ 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;
  int i;
  C=1;
  /*@ loop invariant C== \pow(2., i) && 0 <= i <= 27;
    @ loop variant (27-i); */
  for (i=0; i<27; i++) 
    C*=2;
  C++;
  /*@ assert C == \pow(2.,i) + 1. && i==27; */

  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;
}