/*
  Author:  Pate Williams (c) 1997

  The following program determines the characteristic
  polynomial of a matrix and computes the adjoint
  matrix. See "A Course in Computational Algebraic
  Number Theory" by Henri Cohen page 54.
*/

#include <iostream.h>
#include <math.h>

double **CreateSquareDoubleMatrix(int n)
{
  double **matrix = new double*[n];

  for (int i = 0; i < n; i++) {
    matrix[i] = new double[n];
    for (int j = 0; j < n; j++)
      matrix[i][j] = 0.0;
  }
  return matrix;
}

void DeleteSquareMatrix(int n, double **matrix)
{
  for (int i = 0; i < n; i++) delete[] matrix[i];
  delete[] matrix;
}

void MatMul(int m, int n, int p, double **a, double **b, double **c)
{
  double sum;
  int i, j, k;

  for (i = 0; i < m; i++) {
    for (j = 0; j < p; j++) {
      sum = 0.0;
      for (k = 0; k < n; k++)
        sum += a[i][k] * b[k][j];
      c[i][j] = sum;
    }
  }
}

double Trace(int n, double **matrix)
{
  double sum = 0.0;

  for (int i = 0; i < n; i++) sum += matrix[i][i];
  return sum;
}

void CharPoly(int n, double *a, double **m, double **madj)
//returns the characteristic polynomial and adjoint matrix
{
  double **c = CreateSquareDoubleMatrix(n), **p = CreateSquareDoubleMatrix(n),
         b;
  int i = 0, j, k;

  for (j = 0; j < n; j++) c[j][j] = 1.0;
  a[0] = 1.0;
  for (;;) {
    i++;
    if (i == n) break;
    MatMul(n, n, n, m, c, p);
    for (j = 0; j < n; j++)
      for (k = 0; k < n; k++)
        c[j][k] = p[j][k];
    b = a[i] = - Trace(n, c) / i;
    for (j = 0; j < n; j++) c[j][j] += b;
  }
  MatMul(n, n, n, m, c, p);
  a[n] = - Trace(n, p) / n;
  b = pow(- 1.0, n - 1);
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      madj[i][j] = b * c[i][j];
  DeleteSquareMatrix(n, c);
}

void main(void)
{
  double **m = CreateSquareDoubleMatrix(2),
         **p = CreateSquareDoubleMatrix(2), a[100];
  int i, j, n = 2;

  m[0][0] = 1.0, m[0][1] = 1.0, m[1][0] = 2.0, m[1][1] = 2.0;
  cout << "the matrix is as follows:" << endl << endl;
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++)
      cout << m[i][j] << ' ';
    cout << endl;
  }
  CharPoly(n, a, m, p);
  cout << endl;
  cout << "the characteristic polynomial is as follows:" << endl << endl;
  for (i = 0; i <= n; i++) cout << a[i] << ' ';
  cout << endl << endl;
  cout << "the adjoint matrix is as follows:" << endl << endl;
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++)
      cout << p[i][j] << ' ';
    cout << endl;
  }
  DeleteSquareMatrix(n, m);
  DeleteSquareMatrix(n, p);
}