/*
  Author:  Pate Williams (c) 1997

  "Algorithm 2.2.3 (Determinant, Using Ordinary
  Elimination). Let M be an n by n matrix. This
  algorithm either outputs the determinant of
  M." -Henri Cohen- See "A Course in Computational
  Algebraic Number Theory" by Henri Cohen pages
  50-51. We specialize to the case where m[i][j] is
  in the field Q (the rationals).
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

double **create_square_matrix(long n)
{
  double **matrix = calloc(n, sizeof(double *));
  long i;

  if (!matrix) {
    fprintf(stderr, "fatal error\ninsufficient memory\n");
    fprintf(stderr, "from create_matrix\n");
    exit(1);
  }
  for (i = 0; i < n; i++) {
    matrix[i] = calloc(n, sizeof(double));
    if (!matrix[i]) {
      fprintf(stderr, "fatal error\ninsufficient memory\n");
      fprintf(stderr, "from create_matrix\n");
      exit(1);
    }
  }
  return matrix;
}

void delete_matrix(long n, double **matrix)
{
  long i;

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

double determinant(long n, double **m)
{
  int found;
  double ck, d, t, x = 1.0;
  long i, j, k, l;

  for (j = 0; j < n; j++) {
    found = 0, i = j;
    while (!found && i < n) {
      found = m[i][j] != 0;
      if (!found) i++;
    }
    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 = 1.0 / m[j][j];
    for (k = j + 1; k < n; k++) {
      ck = d * m[k][j];
      for (l = j + 1; l < n; l++)
        m[k][l] -= ck * m[j][l];
    }
    x *= m[j][j];
  }
  return x;
}

int main(void)
{
  long i, j, n = 6;
  double p1 = 1.0, p2 = 1.0;
  double x[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
  double **m = create_square_matrix(n);

  /* create the Vandermonde's matrix of order n */
  printf("the Vandermonde's matrix of order %ld is:\n", n);
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      m[i][j] = pow(x[j], i + 1);
      printf("%5.0lf ", m[i][j]);
    }
    printf("\n");
  }
  printf("det = %lf\n", determinant(n, m));
  for (j = 0; j < n; j++) p1 *= x[j];
  for (i = 0; i < n; i++)
    for (j = i + 1; j < n; j++)
      p2 *= x[j] - x[i];
  printf("det = %lf\n", p1 * p2);
  delete_matrix(n, m);
  return 0;
}