#include <stdio.h>
#ifdef _CIVL
#include <civlc.cvh>
#pragma CIVL ACSL
#endif

#define UP_ONE_REM(up, step) ((up) + (1 - (up) % (step)))
#define UP_MUL_REM(up, step) ((up) - ((up) % (step) - 1) + (step))

#ifdef _CIVL
$input int ni,nj;
$input int bi,bj;  //TODO: bounds of block sizes ?
$assume(ni > 1);
$assume(nj > 1);
$assume(bi > 0);
$assume(bj > 0);
$input double uin[ni][nj];
#endif
double u1[ni][nj];
double u2[ni][nj];

int main(int argc, char** argv) {
    int i,j,ir,jr,ib,jb;
    
    /*
    for(i=0;i<ni;i++) {
      for(j=0;j<nj;j++) {
        u1[i][j] = 0;
        u2[i][j] = 0;
      }
    }

    // straightforward loop nest
    for(i=1;i<ni-1;i++) {
      for(j=1;j<nj-1;j++) {
        u1[i][j] = uin[i][j]*uin[i][j];
      }
    }
    */
    $havoc(&u1);
    $havoc(&u2); 
    $assume($forall (int i : 0 .. ni-1) u1[i][0] == 0 &&  u1[i][nj-1] == 0 && u2[i][0] == 0 &&  u2[i][nj-1] == 0);
    $assume($forall (int j : 0 .. nj-1) u1[0][j] == 0 &&  u1[ni-1][j] == 0 && u2[0][j] == 0 &&  u2[ni-1][j] == 0);
    $assume($forall (int i : 1 .. ni-2) ($forall (int j : 1 .. nj-2) u1[i][j] == uin[i][j] * uin[i][j]));

    if (ni > 2) {
      i = ni - 1;	
      if (nj > 2)
	j = nj - 1;
      else
	j = 1;
    } else {
      i = 1;
      j = nj;
    }

    /*@ predicate noRemain(int upper, int step) = upper % step == 1;
      @
      @ predicate oneRemain(int upper, int step) = upper % step == 0;
      @*/
    
    $choose {
      $when (1) $assume(noRemain(ni-bi, bi));
      $when (1) $assume(oneRemain(ni-bi, bi));
      $when (1) $assume(!noRemain(ni-bi, bi) && !oneRemain(ni-bi, bi));
    }

    $choose {    
      $when (1) $assume(noRemain(nj-bj, bj));
      $when (1) $assume(oneRemain(nj-bj, bj));
      $when (1) $assume(!noRemain(nj-bj, bj) && !oneRemain(nj-bj, bj));
    }
    
    // blocked loop nest

    /*@ loop invariant noRemain(ni-bi, bi) ==> (1 <= i <= ni-bi);
      @ loop invariant oneRemain(ni-bi, bi) ==> 1 <= i <= UP_ONE_REM(ni-bi, bi);
      @ loop invariant (!noRemain(ni-bi, bi) && !oneRemain(ni-bi, bi)) ==> 1 <= i <= UP_MUL_REM(ni-bi,bi);
      @ loop invariant i % bi == 1;
      @ loop invariant noRemain(nj-bj, bj) ==> 
      @                  (\forall int t; 1 <= t < i ==> 
      @                    (\forall int k; 1 <= k < nj-bj ==> 
      @                      u2[t][k] == uin[t][k]*uin[t][k]));
      @ loop invariant oneRemain(nj-bj, bj) ==>
      @                  (\forall int t; 1 <= t < i ==> 
      @                    (\forall int k; 1 <= k < UP_ONE_REM(nj-bj, bj) ==> 
      @                      u2[t][k] == uin[t][k]*uin[t][k]));
      @ loop invariant !oneRemain(nj-bj, bj) && !noRemain(nj-bj, bj) ==>
      @                  (\forall int t; 1 <= t < i ==> 
      @                    (\forall int k; 1 <= k < UP_MUL_REM(nj-bj, bj) ==> 
      @                      u2[t][k] == uin[t][k]*uin[t][k]));
      @ loop assigns i, j, ib, jb, u2[1 .. ni-2][1 .. nj-2];
      @*/
    for(i=1;i<ni-bi;i=i+bi) {
      /*@ loop invariant noRemain(nj-bj, bj) ==> 1 <= j <= nj-bj;
	@ loop invariant oneRemain(nj-bj, bj) ==> 1 <= j <= UP_ONE_REM(nj-bj,bj);
	@ loop invariant !noRemain(nj-bj, bj) && !oneRemain(nj-bj, bj) ==> 1 <= j <= UP_MUL_REM(nj-bj,bj);
	@ loop invariant j % bj == 1;
	@ loop invariant noRemain(nj-bj, bj) ==> 
	@                  (\forall int t; 0 <= t < bi ==> 
	@                    (\forall int k; 1 <= k < nj-bj ==> 
	@                      u2[i+t][k] == uin[i+t][k]*uin[i+t][k]));
	@ loop invariant oneRemain(nj-bj, bj) ==> 
	@                  (\forall int t; 0 <= t < bi ==> 
	@                    (\forall int k; 1 <= k < UP_ONE_REM(nj-bj, bj) ==> 
	@                      u2[i+t][k] == uin[i+t][k]*uin[i+t][k]));
	@ loop invariant !noRemain(nj-bj, bj) && !oneRemain(nj-bj, bj) ==>
	@                  (\forall int t; 0 <= t < bi ==> 
	@                    (\forall int k; 1 <= k < UP_MUL_REM(nj-bj, bj) ==> 
	@                      u2[i+t][k] == uin[i+t][k]*uin[i+t][k]));
	@ loop assigns j, ib, jb, u2[i .. i + bi - 1][1 .. nj - 2];
	@*/
      for(j=1;j<nj-bj;j=j+bj) {
	/*@ loop invariant 0 <= ib <= bi;
	  @ loop invariant \forall int t; 0 <= t < ib ==> 
	  @                   (\forall int k; 0 <= k < bj ==> u2[i+t][j+k] == uin[i+t][j+k]*uin[i+t][j+k]);
	  @ loop assigns   ib, jb, u2[i .. (i + bi - 1)][j .. (j + bj - 1)];
	  @*/
        for(ib=0;ib<bi;ib++) {
	  /*@ loop invariant 0 <= jb <= bj;
	    @ loop invariant \forall int k; 0 <= k < jb ==> u2[i+ib][j+k] == uin[i+ib][j+k]*uin[i+ib][j+k];
	    @ loop assigns   jb, u2[i + ib][j .. (j + bj - 1)];
	    @*/
          for(jb=0;jb<bj;jb++) {
            u2[i+ib][j+jb] = uin[i+ib][j+jb]*uin[i+ib][j+jb];
          }
        }
      }
    }
}