/*@ axiomatic dirichlet_maths {
  @
  @ logic real c;
  @ logic real p0(real x);
  @ logic real psol(real x, real t);

  @ axiom c_pos: 0 < c;

  @ logic real psol_1(real x, real t);
  @ axiom psol_1_def:
  @   \forall real x; \forall real t;
  @   \forall real eps; \exists real C; 0 < C && \forall real dx;
  @   \abs(dx) < C ==>
  @   \abs((psol(x + dx, t) - psol(x, t)) / dx - psol_1(x, t)) < eps;

  @ logic real psol_11(real x, real t);
  @ axiom psol_11_def:
  @   \forall real x; \forall real t;
  @   \forall real eps; \exists real C; 0 < C && \forall real dx;
  @   \abs(dx) < C ==>
  @   \abs((psol_1(x + dx, t) - psol_1(x, t)) / dx - psol_11(x, t)) < eps;

  @ logic real psol_2(real x, real t);
  @ axiom psol_2_def:
  @   \forall real x; \forall real t;
  @   \forall real eps; \exists real C; 0 < C && \forall real dt;
  @   \abs(dt) < C ==>
  @   \abs((psol(x, t + dt) - psol(x, t)) / dt - psol_2(x, t)) < eps;

  @ logic real psol_22(real x, real t);
  @ axiom psol_22_def:
  @   \forall real x; \forall real t;
  @   \forall real eps; \exists real C; 0 < C && \forall real dt;
  @   \abs(dt) < C ==>
  @   \abs((psol_2(x, t + dt) - psol_2(x, t)) / dt - psol_22(x, t)) < eps;

  @ axiom wave_eq_0: \forall real x; 0 <= x <= 1 ==> psol(x, 0) == p0(x);
  @ axiom wave_eq_1: \forall real x; 0 <= x <= 1 ==> psol_2(x, 0) == 0;
  @ axiom wave_eq_2:
  @   \forall real x; \forall real t;
  @   0 <= x <= 1 ==> psol_22(x, t) - c * c * psol_11(x, t) == 0;
  @ axiom wave_eq_dirichlet_1: \forall real t; psol(0, t) == 0;
  @ axiom wave_eq_dirichlet_2: \forall real t; psol(1, t) == 0;

  @ logic real psol_Taylor_3(real x, real t, real dx, real dt);
  @ logic real psol_Taylor_4(real x, real t, real dx, real dt);

  @ logic real alpha_3; logic real C_3;
  @ logic real alpha_4; logic real C_4;

  @ axiom psol_suff_regular_3:
  @   0 < alpha_3 && 0 < C_3 &&
  @   \forall real x; \forall real t; \forall real dx; \forall real dt;
  @   0 <= x <= 1 ==> \sqrt(dx * dx + dt * dt) <= alpha_3 ==>
  @   \abs(psol(x + dx, t + dt) - psol_Taylor_3(x, t, dx, dt)) <=
  @     C_3 * \abs(\pow(\sqrt(dx * dx + dt * dt), 3));

  @ axiom psol_suff_regular_4:
  @   0 < alpha_4 && 0 < C_4 &&
  @   \forall real x; \forall real t; \forall real dx; \forall real dt;
  @   0 <= x <= 1 ==> \sqrt(dx * dx + dt * dt) <= alpha_4 ==>
  @   \abs(psol(x + dx, t + dt) - psol_Taylor_4(x, t, dx, dt)) <=
  @     C_4 * \abs(\pow(\sqrt(dx * dx + dt * dt), 4));

  @ axiom psol_le:
  @   \forall real x; \forall real t;
  @   0 <= x <= 1 ==> 0 <= t ==> \abs(psol(x, t)) <= 1;

  @ logic real T_max;
  @ axiom T_max_pos: 0 < T_max;

  @ logic real C_conv; logic real alpha_conv;
  @ lemma alpha_conv_pos: 0 < alpha_conv;
  @
  @ } */


/*@ axiomatic dirichlet_prog {
  @
  @  predicate analytic_error{L}
  @    (double **p, integer ni, integer i, integer k, double a, double dt)
  @    reads p[..][..];
  @
  @  lemma analytic_error_le{L}:
  @    \forall double **p; \forall integer ni; \forall integer i;
  @    \forall integer nk; \forall integer k;
  @    \forall double a; \forall double dt;
  @    0 < ni ==> 0 <= i <= ni ==> 0 <= k ==>
  @    0 < \exact(dt) ==>
  @    analytic_error(p, ni, i, k, a, dt) ==>
  @    \sqrt(1. / (ni * ni) + \exact(dt) * \exact(dt)) < alpha_conv ==>
  @    k <= nk ==> nk <= 7598581 ==> nk * \exact(dt) <= T_max ==>
  @    \exact(dt) * ni * c <= 1 - 0x1.p-50 ==>
  @    \forall integer i1; \forall integer k1;
  @    0 <= i1 <= ni ==> 0 <= k1 < k ==>
  @    \abs(p[i1][k1]) <= 2;
  @
  @  predicate separated_matrix{L}(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]);
  @
  @ logic real sqr_norm_dx_conv_err{L}
  @   (double **p, real dx, real dt, integer ni, integer i, integer k)
  @   reads p[..][..];
  @ logic real sqr(real x) = x * x;
  @ lemma sqr_norm_dx_conv_err_0{L}:
  @   \forall double **p; \forall real dx; \forall real dt;
  @   \forall integer ni; \forall integer k;
  @   sqr_norm_dx_conv_err(p, dx, dt, ni, 0, k) == 0;
  @ lemma sqr_norm_dx_conv_err_succ{L}:
  @   \forall double **p; \forall real dx; \forall real dt;
  @   \forall integer ni; \forall integer i; \forall integer k;
  @   0 <= i ==>
  @   sqr_norm_dx_conv_err(p, dx, dt, ni, i + 1, k) ==
  @     sqr_norm_dx_conv_err(p, dx, dt, ni, i, k) +
  @     dx * sqr(psol(1. * i / ni, k * dt) - \exact(p[i][k]));
  @ logic real norm_dx_conv_err{L}
  @   (double **p, real dt, integer ni, integer k) =
  @   \sqrt(sqr_norm_dx_conv_err(p, 1. / ni, dt, ni, ni, k));
  @ } */


/*@ 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 &&
  @   \exact(\result) == p0(\exact(x));
  @ */
double p_zero(double xs, double l, double x);


/*@ requires
  @   ni >= 2 && nk >= 2 && l != 0 &&
  @   dt > 0. && \exact(dt) > 0. &&
  @   \exact(v) == c && \exact(v) == v &&
  @   0x1.p-1000 <= \exact(dt) &&
  @   ni <= 2147483646 && nk <= 7598581 &&
  @   nk * \exact(dt) <= T_max &&
  @   \abs(\exact(dt) - dt) / dt <= 0x1.p-51 &&
  @   0x1.p-500 <= \exact(dt) * ni * c <= 1 - 0x1.p-50 &&
  @   \sqrt(1. / (ni * ni) + \exact(dt) * \exact(dt)) < alpha_conv;
  @
  @ ensures
  @   \forall integer i; \forall integer k;
  @   0 <= i <= ni ==> 0 <= k <= nk ==>
  @   \round_error(\result[i][k]) <= 78. / 2 * 0x1.p-52 * (k + 1) * (k + 2);
  @
  @ ensures
  @   \forall integer k; 0 <= k <= nk ==>
  @   norm_dx_conv_err(\result, \exact(dt), ni, k) <=
  @     C_conv * (1. / (ni * ni) + \exact(dt) * \exact(dt));
  @ */
double **forward_prop(int ni, int nk, double dt, double v,
                      double xs, double l) {

  /* Output variable. */
  double **p;

  /* Local variables. */
  int i, k;
  double a1, a, dp, dx;

  dx = 1./ni;
  /*@ assert
    @   dx > 0. && dx <= 0.5 &&
    @   \abs(\exact(dx) - dx) / dx <= 0x1.p-53;
    @ */

  /* Compute the constant coefficient of the stiffness matrix. */
  a1 = dt/dx*v;
  a = a1*a1;
  /*@ assert
    @   0 <= a <= 1 &&
    @   0 < \exact(a) <= 1 &&
    @   \round_error(a) <= 0x1.p-49;
    @ */

  /* Allocate space-time variable for the discrete solution. */
  p = array2d_alloc(ni+1, nk+1);

  /* First initial condition and boundary conditions. */
  /* Left boundary. */
  p[0][0] = 0.;
  /* Time iteration -1 = space loop. */
  /*@ loop invariant
    @   1 <= i <= ni &&
    @   analytic_error(p, ni, i - 1, 0, a, dt);
    @ loop variant ni - i; */
  for (i=1; i<ni; i++) {
    p[i][0] = p_zero(xs, l, i*dx);
  }
  /* Right boundary. */
  p[ni][0] = 0.;
  /*@ assert analytic_error(p, ni, ni, 0, a, dt); */

  /* Second initial condition (with p_one=0) and boundary conditions. */
  /* Left boundary. */
  p[0][1] = 0.;
  /* Time iteration 0 = space loop. */
  /*@ loop invariant
    @   1 <= i <= ni &&
    @   analytic_error(p, ni, i - 1, 1, a, dt);
    @ loop variant ni - i; */
  for (i=1; i<ni; i++) {
    /*@ assert \abs(p[i-1][0]) <= 2; */
    /*@ assert \abs(p[i][0]) <= 2; */
    /*@ assert \abs(p[i+1][0]) <= 2; */
    dp = p[i+1][0] - 2.*p[i][0] + p[i-1][0];
    p[i][1] = p[i][0] + 0.5*a*dp;
  }
  /* Right boundary. */
  p[ni][1] = 0.;
  /*@ assert analytic_error(p, ni, ni, 1, a, dt); */

  /* Evolution and boundary conditions. */
  /* Propagation = time loop. */
  /*@ loop invariant
    @   1 <= k <= nk &&
    @   analytic_error(p, ni, ni, k, a, dt);
    @ loop variant nk - k; */
  for (k=1; k<nk; k++) {
    /* Left boundary. */
    p[0][k+1] = 0.;
    /* Time iteration k = space loop. */
    /*@ loop invariant
      @   1 <= i <= ni &&
      @   analytic_error(p, ni, i - 1, k + 1, a, dt);
      @ loop variant ni - i; */
    for (i=1; i<ni; i++) {
      /*@ assert \abs(p[i-1][k]) <= 2; */
      /*@ assert \abs(p[i][k]) <= 2; */
      /*@ assert \abs(p[i+1][k]) <= 2; */
      /*@ assert \abs(p[i][k-1]) <= 2; */
      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;
    }
    /* Right boundary. */
    p[ni][k+1] = 0.;
    /*@ assert analytic_error(p, ni, ni, k + 1, a, dt); */
  }

  return p;

}