/* Discrete derivative function using central differencing, which is
 * second-order accurate, except at the two end-points (where it is
 * first-order).
 * This version is incorrect, because it asserts third-order accuracy.
 * To attempt to verify with CIVL, type:
 *		civl verify derivativeBad.cvl
 */
$input double dx; // delta x
$assume(0<dx && dx<1);
$input int num_elements;
$assume(num_elements >= 1);
$input double in[num_elements];
double out[num_elements];
// the following says rho is a function from R to R which has 3 continuous
// derivatives in the closed interval [-1,1]:
$abstract $differentiable(3, [-1,1]) $real rho($real x);

/* Computes discrete derivative by central differencing.
 * Right end-point is computed by backwards differencing.
 * Left end-point is computed by forward differencing.
 */
void differentiate(int n, double y[], double h, double result[]) {
  $assume(n*h<=1);
  $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]-y[j-1])/(2*h);
    @ loop assigns i, result[1..n-2];
    @*/
  for (int i=1; i<n-1; i++) {
    result[i] = (y[i+1]-y[i-1])/(2*h);
  }
  result[0] = (y[1]-y[0])/h; 
  result[n-1] = (y[n-1] - y[n-2])/h;
  $assert($uniform (int i : 1..n-2) result[i]-$D[rho,{x,1}](i*h) == $O(h*h*h));
}

int main() {
  differentiate(num_elements, in, dx, out);
}