/* A discrete approximation to the second derivative of a function from R to R.
 * It is second-order accurate, except at the two endpoints. 
 * To verify with CIVL, type:
 *		civl verify secondDerivative.cvl
 *
 * Note: based on Quarteroni, Sacco, Saleri. "Numerical Mathematics" 2nd ed. sec 10.10.1
 */

$input double dx;
$assume(0<dx && dx<1);
$input int num_elements;
$assume (num_elements > 0);
$input double in[num_elements];
double out[num_elements];
// assume rho:R->R has 4 continuous derivatives on [-1,1]:
$abstract $differentiable(4, [-1,1]) $real rho($real x);

void secondDerivative(int n, double y[], double h, double result[]) {
  $assume($forall (int i:0..n-1) y[i] == rho(i*h));
  /*@ loop invariant 1<=i && i<=n-1;
    @ loop invariant $forall (int j:1..i-1)
    @   result[j] == (y[j+1] - 2*y[j] + y[j-1])/(h*h);
    @ loop assigns i, result[1..n-2];
    @*/
  for (int i = 1; i < n-1; i++) {
    result[i] = (y[i+1] - 2*y[i] + y[i-1])/(h*h);
  }
  result[0] = (y[2] - 2*y[1] + y[0])/h; 
  result[n-1] = (y[n-3] - 2*y[n-2] - y[n-1])/h;
  $assert($uniform (int i:1..n-2) result[i] - $D[rho,{x,2}](i*h) == $O(h*h));
}

void main() {
  secondDerivative(num_elements, in, dx, out);
}