/*
  Author:  Pate Williams (c) 1997

  Reduction of m by n matrix to Hermite normal
  form. See "A Course in Computational Algebraic
  Number Theory" by Henri Cohen 2.4.2 Section
  Algorithm 2.4.5 page 69.
*/

#include <assert.h>
#include <malloc.h>
#include <stdio.h>
#include "lip.h"

verylong  **create_matrix(long m, long n)
{
  long i;
  verylong **zmatrix = calloc(m, sizeof(verylong *));

  assert(zmatrix != 0);
  for (i = 0; i < m; i++) {
    zmatrix[i] = calloc(n, sizeof(verylong));
    assert(zmatrix[i] != 0);
  }
  return zmatrix;
}

void delete_matrix(long m, long n, verylong **zmatrix)
{
  long i, j;

  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++)
      zfree(&zmatrix[i][j]);
    free(zmatrix[i]);
  }
  free(zmatrix);
}

void Hermite_normal_form(long m, long n, verylong zn,
                         verylong **za, verylong **zw,
                         long *nk)
{
  long i, j, k, l, p;
  verylong zaij = 0, zaijd = 0, zaik = 0, zaikd = 0;
  verylong zb = 0, zd = 0, zq = 0, zr = 0, zs = 0;
  verylong zu = 0, zv = 0;
  verylong *zB = calloc(m, sizeof(verylong));

  assert(zB != 0);
  i = m - 1, j = n - 1, k = n - 1;
  if (m <= n) l = 0; else l = m - n;
  L2:
  if (j == 0) goto L4;
  j--;
  if (zscompare(za[i][j], 0l) == 0) goto L2;
  zcopy(za[i][j], &zaij);
  zcopy(za[i][k], &zaik);
  zexteucl(zaik, &zu, zaij, &zv, &zd);
  zdiv(zaik, zd, &zaikd, &zr);
  zdiv(zaij, zd, &zaijd, &zr);
  for (p = 0; p < m; p++) {
    zmulmod(zu, za[p][k], zn, &zr);
    zmulmod(zv, za[p][j], zn, &zs);
    zaddmod(zr, zs, zn, &zB[p]);
    zmulmod(zaikd, za[p][j], zn, &zr);
    zmulmod(zaijd, za[p][k], zn, &zs);
    zsubmod(zr, zs, zn, &za[p][j]);
  }
  for (p = 0; p < m; p++) zcopy(zB[p], &za[p][k]);
  goto L2;
  L4:
  zcopy(za[i][k], &zb);
  if (zscompare(zb, 0l) < 0) {
    for (p = 0; p < m; p++) znegate(&za[p][k]);
    znegate(&zb);
  }
  if (zscompare(zb, 0l) == 0) {
    k++;
    goto L5;
  }
  for (j = k + 1; j < n; j++) {
    zdiv(za[i][j], zb, &zq, &zr);
    for (p = 0; p < m; p++) {
      zmul(zq, za[p][k], &zr);
      zsubmod(za[p][j], zr, zn, &zs);
      zcopy(zs, &za[p][j]);
    }
  }
  L5:
  if (i == l) {
    *nk = n - k;
    for (j = 0; j < *nk; j++)
      for (p = 0; p < m; p++)
        zcopy(za[p][j + k], &zw[p][j]);
  }
  else {
    i--;
    k--;
    j = k;
    goto L2;
  }
  zfree(&zaij);
  zfree(&zaijd);
  zfree(&zaik);
  zfree(&zaikd);
  zfree(&zb);
  zfree(&zd);
  zfree(&zq);
  zfree(&zr);
  zfree(&zs);
  zfree(&zu);
  zfree(&zv);
  free(zB);
}

int main(void)
{
  long i, j, m = 6, n = 11, nk;
  verylong zn = 0;
  verylong **za = create_matrix(m, n);
  verylong **zw = create_matrix(m, n);

  zintoz(2l, &za[0][0]);
  zintoz(2l, &za[0][1]);
  zone(&za[0][2]);
  zintoz(4l, &za[1][0]);
  zone(&za[1][4]);
  zone(&za[2][1]);
  zone(&za[2][2]);
  zone(&za[2][4]);
  zone(&za[3][0]);
  zone(&za[3][3]);
  zone(&za[3][4]);
  zone(&za[4][0]);
  zintoz(2l, &za[4][1]);
  zone(&za[4][4]);
  zone(&za[5][0]);
  zone(&za[5][1]);
  zone(&za[5][2]);
  zone(&za[5][3]);
  zintoz(228l, &za[0][5]);
  zintoz(228l, &za[1][6]);
  zintoz(228l, &za[2][7]);
  zintoz(228l, &za[3][8]);
  zintoz(228l, &za[4][9]);
  zintoz(228l, &za[5][10]);
  zintoz(228l, &zn);
  printf("the original adjoined matrix is as follows:\n\n");
  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++)
      printf("%3ld ", ztoint(za[i][j]));
    printf("\n");
  }
  printf("\n");
  Hermite_normal_form(m, n, zn, za, zw, &nk);
  printf("the Hermite normal form matrix is as follows:\n\n");
  for (i = 0; i < m; i++) {
    for (j = 0; j < nk; j++)
      printf("%3ld ", ztoint(zw[i][j]));
    printf("\n");
  }
  delete_matrix(m, n, za);
  delete_matrix(m, n, zw);
  return 0;
}