/*
{m:nat, n:nat}
record <'a> matrix(m,n) {
  row: int(m);
  col: int(n);
  data[m][n]: 'a
}
*/

{p:nat, q:nat, r:nat}
<float> matrix(p,r)
matmul (a: <float> matrix(p,q), b: <float> matrix(q,r)) {
  var: nat i, j, k, row, col; float sum; c: <float> matrix;;

  c = matrix { row = a.row;
               col = b.col;
               data = alloc(a.row, alloc(b.col, 0.0)) };

  for (i = 1; i < a.row; i = i + 1) {
    c.data[i] = alloc(b.col, 0.0);
  }

  for (i = 0; i < a.row; i = i + 1) {
    for (j = 0; j < b.col; j = j + 1) {
      sum = 0.0;
      for (k = 0; k < a.col; k = k + 1) {
        sum = sum +. a.data[i][k] *. b.data[k][j];
      }
      c.data[i][j] = sum;
    }
  }
  return c;
}

{n:nat}
unit print_array (size: int(n), a[n]: float) {
  var: nat i;;
  for (i = 0; i < size; i = i + 1) {
      print_float (a[i]); print_string ("\t");
    }
    print_newline ();
}

unit print_matrix (a: <float> matrix) {
  var: nat i;;
  for (i = 0; i < a.row; i = i + 1) {
      print_array (a.col, a.data[i]);
    }
}

unit main () {
  var: nat i, j, n; a: <float> matrix;;

  n = 8;
  a = matrix { row = n; col = n; data = alloc (n, alloc (n, 0.0)) };

  for (i = 1; i < a.row; i = i + 1) {
    a.data[i] = alloc(a.col, 0.0);
  }

  for (i = 0; i < a.row; i = i + 1) {
    for (j = 0; j < a.col; j = j + 1) {
      a.data[i][j] = float_of_int (i - j);
    }
  }
  print_matrix (a);
  a = matmul (a, a);
  print_newline ();
  print_matrix (a);
}