Sylvie Boldo

Gallery of proved programs

All those programs were proved using the Bore version of Frama-C and the 2.26 version of Why. The Coq proofs use both the PFF and the gappalib libraries for Coq results about floating-point arithmetic.

dirichlet: discretization of the 1D acoustic wave equation

This example was developed among the FOST project. It is a finite difference numerical scheme for the resolution of the one-dimensional acoustic wave equation. A reasonable bound on the rounding error requires a complex property that expresses each rounding error (see here for more details).

The bound on the method error of this scheme has been formally certified in Coq (see here and here for more details). As the bound on the method error (and its properties) has not yet been linked to the bound rounding error, this proof is incomplete (2 admits are left, corresponding to the bounds on the exact values taken by the scheme, and another one corresponding to a complex mathematical property).

#pragma JessieIntegerModel(math)

/*@ axiomatic dirichlet {
  @  predicate separated_matrix(double **p,integer leni) =
  @     \forall integer i; \forall integer j; 
  @          0 <= i < leni ==> 0 <= j < leni ==> i != j ==>
  @          \base_addr(p[i]) != \base_addr(p[j]);
  @
  @  predicate analytic_error{L}
  @          (double **p, integer ni, integer i, integer k, double a)
  @          reads p[..][..];
  @ } */




/*@ requires leni >= 1 && lenj >= 1;
  @  ensures  \valid_range(\result,0,leni-1) &&
  @           (\forall integer i; 0 <= i < leni ==> 
  @                \valid_range(\result[i],0,lenj-1)) &&
  @            separated_matrix(\result,leni);
  @ */
double **array2d_alloc(int leni, int lenj);


/*@ requires (l != 0);
  @ ensures \round_error(\result) <= 14*0x1.p-52 && \abs(\result) <= 2;
  @ */
double p_zero(double xs, double l, double x);


/*@ requires ni >= 2 && nk >= 2
  @          && l != 0
  @          && dx > 0. && dt > 0.  && v > 0.
  @          && \exact(dx) > 0. && \exact(dt) > 0.
  @          && \exact(v)==v
  @          && \abs(\exact(dx) - dx) / dx <= 0x1.p-53
  @          && \abs(\exact(dt) - dt) / dt <= 0x1.p-51
  @          && 3./5. <= \exact(dt)/\exact(dx) * \exact(v) <= 1-0x1.p-50
  @          && 0x1.p-1000 <= v <= 0x1.p1000
  @          && ni <= 0x1.p64 && nk <= 4194304
  @          && \exact(dx) <= 1;
  @
  @
  @ ensures \forall int i; \forall int k;
  @     0 <= i <= ni ==> 0 <= k <= nk ==> 
  @     \round_error(\result[i][k]) <= 78./2*0x1.p-52*(k+1)*(k+2); 
  @ */

double **forward_prop(int ni, int nk, double dx, double dt, double v,
                      double xs, double l) {

  /* output variable */
  double **p;
  
  /* local variables */
  int i, k;
  double a1, a, dp;

  
  /* assemble diagonal of stiffness matrix M[v]^{-1}A */
  a1 =  dt/dx*v;
  a  = a1*a1;
  
  /*@ assert 1./4 <= a <= 1 && 
    @         0 < \exact(a) <= 1 &&
    @         \round_error(a) <= 0x1.p-49;
    @  */

  p = array2d_alloc(ni+1, nk+1);
  /* initial conditions (includes boundary conditions) */

  p[0][0]=0.;
  /*@ loop invariant 1 <= i <= ni && analytic_error(p,ni,i-1,0,a);
    @ loop variant ni-i; */
  for (i=1; i<ni; i++) {
    p[i][0] = p_zero(xs, l, i*dx);
  }
  p[ni][0] =0.;

  /*@ assert analytic_error(p,ni,ni,0,a); */


  /* initial time derivative is supposed null + boundary conditions */

  p[0][1] = 0.;

  /*@ loop invariant 1 <= i <= ni && analytic_error(p,ni,i-1,1,a);
    @ loop variant ni-i; */
  for (i=1; i<ni; i++) {
    dp = p[i+1][0] - 2.*p[i][0] + p[i-1][0];
    p[i][1] = p[i][0] + 0.5*a*dp;
  }
  p[ni][1] = 0.;
  
  /*@ assert analytic_error(p,ni,ni,1,a); */


  /* propagation = time loop */
  /*@ loop invariant 1 <= k <= nk && analytic_error(p,ni,ni,k,a);
    @ loop variant nk-k; */
  for (k=1; k<nk; k++) {
    /* boundary condition 1 */
    p[0][k+1] = 0.;

    /* time iteration = space loop */
    /*@ loop invariant 1 <= i <= ni && analytic_error(p,ni,i-1,k+1,a);
      @ loop variant ni-i; */
 
    for (i=1; i<ni; i++) {
      dp = p[i+1][k] - 2.*p[i][k] + p[i-1][k];
      p[i][k+1] = 2.*p[i][k] - p[i][k-1] + a*dp;
    }
    
    /* boundary condition 2 */
    p[ni][k+1] = 0.;
    /*@ assert analytic_error(p,ni,ni,k+1,a); */
    
  }

  return p;

}