Scalar product of vectors using floating-point numbers
Bounding the rounding error made in the computation of the scalar product of the vectors given as arrays of double precision floating-point numbers.
Authors: Claude Marché / Thi Minh Tuyen Nguyen
Topics: Floating-Point Computations
Tools: Frama-C / Jessie / Gappa
see also the index (by topic, by tool, by reference, by year)
Description
We consider the following C code which computes the scalar product of two vectors x and y of a given size n.
double scalar_product(double x[], double y[], int n) { double p = 0.0; for (int i=0; i < n; i++) { p = p + x[i]*y[i]; } return p; }
In order to specify some properties of this code with respect to the exact scalar product, in particular to bound the rounding errors, we introduce an axiomatic ScalarProduct which models of the mathematical scalar product. The specification of the C function then provides a bound on the rounding errors made in the floating-point computations. The complete annotated code is show below.
Proofs
Verification Conditions are generated using the Frama-C/Jessie tool chain and its support for floating-point computations [1]. They are proved automatically, most of them by SMT solvers like Alt-Ergo, Z3 or CVC3. The code below indicated the only 2 assertions which are proved by Gappa and only Gappa: the first is the bound on the rounding errors made when computing p + x_{i} × y_{i}. The second is the show that |x_{i}|×|y_{i}| ≤ 1 whenever |x_{i}|≤ 1 and |y_{i}|≤ 1.
The table below summaries the bound found for different values of the size NMAX of the vectors.
NMAX | 10 | 100 | 1000 |
B | 0x1.1p-50 | 0x1.02p-47 | 0x1.004p-44 |
See also [2] for analogous examples.
References
- [1]
- Ali Ayad and Claude Marché. Multi-prover verification of floating-point programs. In Jürgen Giesl and Reiner Hähnle, editors, Fifth International Joint Conference on Automated Reasoning, Lecture Notes in Artificial Intelligence, Edinburgh, Scotland, July 2010. Springer.
- [2]
- Sylvie Boldo and Claude Marché. Formal verification of numerical programs: from C annotated programs to mechanical proofs. Mathematics in Computer Science, 2011.
Complete Annotated Code
#define NMAX 1000 #define NMAXR 1000.0 #define B 0x1.004p-44 /*@ axiomatic ScalarProduct { @ // exact_scalar_product(x,y,n) = x[0]*y[0] + ... + x[n-1] * y[n-1] @ logic real exact_scalar_product{L}(double *x, double *y, integer n) @ reads x[..], y[..]; @ axiom A1{L}: \forall double *x,*y; @ exact_scalar_product(x,y,0) == 0.0; @ axiom A2{L}: \forall double *x,*y; \forall integer n ; @ n >= 0 ==> @ exact_scalar_product(x,y,n+1) == @ exact_scalar_product(x,y,n) + x[n]*y[n]; @ } @*/ /*@ lemma bound_int_to_real: @ \forall integer i; i <= NMAX ==> i <= NMAXR; @*/ /*@ lemma triangular_inequality: @ \forall real x,y,z; @ \abs(x-z) <= \abs(x-y) + \abs(y-z); @*/ /*@ requires 0 <= n <= NMAX; @ requires \valid_range(x,0,n-1) && \valid_range(y,0,n-1) ; @ requires \forall integer i; 0 <= i < n ==> @ \abs(x[i]) <= 1.0 && \abs(y[i]) <= 1.0 ; @ ensures @ \abs(\result - exact_scalar_product(x,y,n)) <= n * B; @*/ double scalar_product(double x[], double y[], int n) { double p = 0.0; /*@ loop invariant 0 <= i <= n ; @ loop invariant \abs(exact_scalar_product(x,y,i)) <= i; @ loop invariant \abs(p - exact_scalar_product(x,y,i)) <= i * B; @ loop variant n-i; @*/ for (int i=0; i < n; i++) { // bounds, needed by Gappa //@ assert \abs(x[i]) <= 1.0; //@ assert \abs(y[i]) <= 1.0; //@ assert \abs(p) <= NMAXR*(1+B) ; L: p = p + x[i]*y[i]; // bound on the rounding errors in the statement above, proved by gappa /*@ assert \abs(p - (\at(p,L) + x[i]*y[i])) <= B; */ // the proper instance of triangular inequality to show the main invariant /*@ assert \abs(p - exact_scalar_product(x,y,i+1)) <= \abs(p - (\at(p,L) + x[i]*y[i])) + \abs((\at(p,L) + x[i]*y[i]) - (exact_scalar_product(x,y,i) + x[i]*y[i])) ; */ // a lemma to show the invariant \abs(exact_scalar_product(x,y,i)) <= i /*@ assert \abs(exact_scalar_product(x,y,i+1)) <= \abs(exact_scalar_product(x,y,i)) + \abs(x[i]) * \abs(y[i]); */ // a necessary lemma, proved only by gappa //@ assert \abs(x[i]) * \abs(y[i]) <= 1.0; } return p; }