/*
  Author:  Pate Williams (c) 1997

  Reduction of a n by n matrix to Smith normal form.
  See "A Course in Computational Algebraic Number
  Theory" by Henri Cohen 2.4.4 Section pages 75-79
  Algorithm 2.4.14 page 77.
*/

#include <assert.h>
#include <malloc.h>
#include <math.h>
#include <mem.h>
#include <stdio.h>

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

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

void delete_matrix(long m, long **matrix)
{
  long i;

  for (i = 0; i < m; i++)
    free(matrix[i]);
  free(matrix);
}

void extended_euclid(long a, long b, long *x,
                      long *y, long *d)
{
  long q, r, x1 = 0, x2 = 1, y1 = 1, y2 = 0;

  if (b == 0) {
    *d = a;
    *x = 1;
    *y = 0;
    return;
  }
  while (b > 0) {
    q = a / b;
    r = a - q * b;
    *x = x2 - q * x1;
    *y = y2 - q * y1;
    a = b, b = r, x2 = x1, x1 = *x, y2 = y1, y1 = *y;
  }
  *d = a, *x = x2, *y = y2;
}

long gcd(long a, long b)
{
  long r;

  while (b != 0)
    r = a %b, a = b, b = r;
  return a;
}

long invmod(long x, long y)
{
  long a, b, d;

  extended_euclid(x, y, &a, &b, &d);
  assert(d == 1);
  return a;
}

long det(long n, long p, long **m)
{
  int found;
  long ck, d, i, j, k, l, t, x = 1;

  for (j = 0; j < n; j++) {
    found = 0, i = 0;
    while (!found && i < n) {
      found = m[i][j] != 0;
      if (!found) i++;
    }
    assert(found != 0);
    if (i > j) {
      for (l = j; l < n; l++)
        t = m[i][l], m[i][l] = m[j][l], m[j][l] = t;
      x = - x;
    }
    d = invmod(m[j][j], p);
    for (k = j + 1; k < n; k++) {
      ck = d * m[k][j] % p;
      for (l = j + 1; l < n; l++)
        m[k][l] = (m[k][l] - ck * m[j][l]) % p;
    }
    x *= m[j][j];
  }
  return x;
}

void Smith_normal_form(long n, long p, long *D, long **m)
{
  int found;
  long aii, aiid, aij, aijd, aji, ajid;
  long b, c, d, R, i, j, k, l, u, v;
  long *B = calloc(n, sizeof(long));
  long **a = create_matrix(n, n);

  assert(B != 0);
  memcpy(a, m, (size_t) n * (size_t) n * sizeof(long));
  i = n - 1, R = labs(det(n, p, m));
  if (n == 1) D[0] = R;
  else {
    L2:
    c = 0, j = i;
    L3:
    if (j == 0) goto L5;
    j--;
    if (a[i][j] == 0) goto L3;
    aii = a[i][i], aij = a[i][j];
    extended_euclid(aii, aij, &u, &v, &d);
    for (k = 0; k < n; k++)
      B[k] = (u * a[k][i] + v * a[k][j]) % p;
    aiid = aii / d, aijd = aij / d;
    for (k = 0; k < n; k++)
      a[k][j] = (aiid * a[k][j] - aijd * a[k][i]) % R;
    for (k = 0; k < n; k++) a[k][i] = B[k] % R;
    goto L3;
    L5:
    j = i;
    L6:
    if (j == 0) goto L8;
    j--;
    if (a[j][i] == 0) goto L6;
    aii = a[i][i], aji = a[j][i];
    extended_euclid(aii, aji, &u, &v, &d);
    for (k = 0; k < n; k++)
      B[k] = (u * a[i][k] + v * a[j][k]) % p;
    aiid = aii / d, ajid = aji / d;
    for (k = 0; k < n; k++)
      a[j][k] = (aiid * a[j][k] - ajid * a[i][k]) % R;
    for (k = 0; k < n; k++) a[i][k] = B[k] % R;
    c++;
    goto L6;
    L8:
    if (c > 0) goto L2;
    b = a[i][i];
    for (k = 0; k < i; k++) {
      found = 0;
      for (l = 0; l < i && !found; l++)
        found = a[k][l] % b != 0;
      if (found) {
        for (l = 0; l < n; l++)
          a[i][l] = (a[i][l] + a[k][l]) % p;
        break;
      }
    }
    D[i] = gcd(a[i][i], R);
    R /= d;
    if (i == 1) D[0] = gcd(a[0][0], R);
    else {
      i--;
      goto L2;
    }
  }
  delete_matrix(n, a);
}

int main(void)
{
  long d[2], i, n = 2, p = 5;
  long **m = create_matrix(n, n);

  m[0][0] = m[1][1] = 2, m[0][1] = 1;
  Smith_normal_form(n, p, d, m);
  for (i = 0; i < n; i++) printf("%ld\n", d[i]);
  delete_matrix(n, m);
  return 0;
}