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

#define EPSILON 0.0000001

double *eigen(double **matrix, int order) {
  double mxerr = 0;
  double *prev = calloc(order, sizeof(double));
  *prev = 1;

  do {
    double *eigen = calloc(order, sizeof(double));

    for (int i = 0; i < order; i++)
      for (int j = 0; j < order; j++)
        eigen[i] += matrix[i][j] * prev[j];

    double mx = fabs(eigen[0]);
    for (int i = 1; i < order; i++)
      mx = fabs(eigen[i]) > mx ? fabs(eigen[i]) : mx;

    for (int i = 0; i < order; i++)
      eigen[i] /= mx;

    double error[order];

    for (int i = 0; i < order; i++)
      error[i] = fabs(eigen[i] - prev[i]);

    mxerr = error[0];
    for (int i = 1; i < order; i++)
      mxerr = error[i] > mxerr ? error[i] : mxerr;

    free(prev);
    prev = eigen;

  } while (mxerr > EPSILON);

  return prev;
}

int main() {
  const int ORDER = 3;

  double **matrix = malloc(ORDER * sizeof(double *));
  for (int i = 0; i < ORDER; i++)
    matrix[i] = malloc(ORDER * sizeof(double));

  matrix[0][0] = 21, matrix[0][1] = -99, matrix[0][2] = 0;
  matrix[1][0] = -12, matrix[1][1] = 64, matrix[1][2] = 0;
  matrix[2][0] = 9, matrix[2][1] = -0.5, matrix[2][2] = 4;

  double *egn = eigen(matrix, ORDER);

  for (int i = 0; i < ORDER; i++) {
    printf("%lf ", egn[i]);
    free(matrix[i]);
  }
  free(egn);
  free(matrix);
}