Sylvie Boldo

Gallery of proved programs

Average of two floating-point numbers, Sterbenz's version

This function computes an approximation of the average of two floating-point numbers. It was written using hints given by Sterbenz in order to prevent overflow. Accuracy is 3/2 ulp.

/*@ ensures \result == \abs(x); */
double abs(double x){
  if (x >= 0) return x;
  else return (-x);
}

/*@ axiomatic Floor {
  @  logic integer floor (real x);
  @  axiom floor_prop: \forall real x; floor(x) <= x < floor(x)+1;
  @ } */

/*@ logic real ulp_d (real x) =
  @ \let e = 1+ floor (\log(\abs(x)) / \log(2)); 
  @    \pow(2,\max(e-53,-1074)); */

/*@ logic real l_average (real x, real y) = 
  @   \let same_sign =
  @     (x >= 0) ? ((y >=0) ? \true : \false) : ((y >=0) ? \false : \true);
  @   (same_sign) ? ((\abs(x) <= \abs(y)) ?
  @     \round_double(\NearestEven, x+\round_double(\NearestEven, 
  @           \round_double(\NearestEven, y-x)/2))  :                   
  @        \round_double(\NearestEven, y+\round_double(\NearestEven, 
  @           \round_double(\NearestEven, x-y)/2))) :
  @     \round_double(\NearestEven,\round_double(\NearestEven, x+y)/2);
  @ */

/*@  lemma average_sym: \forall double x; \forall double y; 
  @        l_average(x,y) == l_average (y,x);
  @  lemma average_sym_opp: \forall double x; \forall double y; 
  @         l_average(-x,-y) == - l_average (x,y);
  @
  @ lemma average_props: \forall double x; \forall double y;
  @        \abs(l_average(x,y) - (x+y)/2) <= 3./2*ulp_d((x+y)/2)
  @   &&  (\min(x,y) <= l_average(x,y) <= \max(x,y))
  @   && (0 <= (x+y)/2 ==> 0 <= l_average(x,y))     
  @   && ((x+y)/2 <= 0 ==> l_average(x,y) <= 0)
  @   && ((x+y)/2==0 ==> l_average(x,y)==0)
  @   && (0x1p-1074 <= \abs((x+y)/2) ==> l_average(x,y) != 0);
  @ */


/*@   ensures \result == l_average(x,y);
  @   ensures \abs((\result - (x+y)/2)) <= 3./2*ulp_d((x+y)/2);
  @   ensures \min(x,y) <= \result <= \max(x,y);
  @   ensures 0 <= (x+y)/2 ==> 0 <= \result;     
  @   ensures (x+y)/2 <= 0 ==> \result <= 0;
  @   ensures (x+y)/2 == 0 ==> \result == 0;
  @   ensures 0x1p-1074 <= \abs((x+y)/2) ==> \result != 0;
  @ */

double average(double x, double y) {
  int same_sign;
  double r;
  if (x >= 0) {
    if (y >=0) same_sign=1;
    else same_sign=0; }
  else {
    if (y >=0) same_sign=0;
    else same_sign=1; }
  if (same_sign ==1) {
    if (abs(x) <= abs(y)) r=x+(y-x)/2;
    else r=y+(x-y)/2; }
  else r=(x+y)/2;
  //@ assert r==l_average(x,y);
  return r;
}