Рассчитать характеристический многочлен матрицы 3-го порядка

//----------------------------------------------------------------------------
#include <iostream>
#include <iomanip>
#include <ctime>
//----------------------------------------------------------------------------
double **Create3x3()
{
    double **arr = new double *[3];
    for (int i = 0; i<3; i++)
        arr[i] = new double [3];
 
    for (int i = 0; i<3; i++)
        for (int j = 0; j<3; j++)
            arr[i][j] = 1.0 *(rand()%100) / (rand()%7 + 1);
 
    return arr;
}
 
void Kill3x3 (double **arr)
{
    for (int i = 0; i<3; i++)
        delete [] arr[i];
    delete []arr;
}
 
void Print3x3 (double **arr)
{
    for (int i = 0; i<3; i++)
    {   for (int j = 0; j<3; j++)
            std::cout << std::fixed << std::setw(10) << std::setprecision(3) << arr[i][j];
        std::cout << "\n";
    }
}
//----------------------------------------------------------------------------
double Det2 (double a, double b, double c, double d)
{
    return a*c-b*d;
}
 
double Det3 (double **arr)
{
    return  arr[0][0]*Det2(arr[1][1], arr[1][2], arr[2][1], arr[2][2]) - 
            arr[0][1]*Det2(arr[1][0], arr[1][2], arr[2][0], arr[2][2]) +
            arr[0][2]*Det2(arr[1][0], arr[1][1], arr[2][0], arr[2][1]);
}
 
//----------------------------------------------------------------------------
int main()
{   
    srand (unsigned int (time));
    
    double **arr = Create3x3();
    std::cout << "\nArray:\n";
    Print3x3 (arr);
    
    double P[4];
    P[0] = -1.0;
    P[1] = arr[0][0] + arr[1][1] + arr[2][2];
    P[2] = -(   Det2(arr[0][0], arr[0][1], arr[1][0], arr[1][1]) + 
                Det2(arr[1][1], arr[1][2], arr[2][1], arr[2][2]) +
                Det2(arr[0][0], arr[0][2], arr[2][0], arr[2][2]));
 
    P[3] = Det3 (arr);
 
    std::cout << "\nPolynomial:\n";
    for (int i = 0; i<4 ; i++)
    {   std::cout << std::fixed << std::setw(10) << std::setprecision(3) << P[i] << " L^" << 3-i;
        if (i<3 && P[i+1] >0)
            std::cout << " +";
    }
 
    Kill3x3 (arr);
    std::cin.get();
    return 0;
}
//-----------------------------------------------------------------------------

float *Eig(float **a,long Size) //релизация функции Eig(X) для квадратной матрицы (Size*Size)
//реализуется расчет коэффициентов при степенях полинома методом Данилевского
{
 float *Vec = NULL;
 long S = Size-1;
 float **Buff = NULL; //буферная матрица
 float **Buff1 = NULL; //еще одна буферная матрица
 float **A = NULL; //буферная матрица, в которую складываю промежуточные результаты. Матрица заменяет матрицу а на второй и последующих итерациях
 A = (float**)malloc(Size*sizeof(float*));
 for (long i = 0; i < Size; i++)
  {
   A[i] = (float*)malloc(Size*sizeof(float));
   for (long j = 0; j < Size; j++) { A[i][j] = 0.0; }
  }
 float **M = NULL; //матрица, на которую постоянно домножать придется
 long Sa,Sa1;
 for (long i = S; i > 1; i--) //операция итеративная
  {
   //Первый блок приведения к матрице Фронебиуса
   if (i == S) M = Diag(Vec,Size,1.0,1); //диагональная матрица с 1 в качестве значений по диагонали
   if (i == S) //матрица M^(-1)n-1 заполняется
    {
     for (long j = 0; j < Size; j++) { M[i-1][j] = a[i][j]; }
    }
   else
    {
     for (long j = 0; j < Size; j++) { M[i-1][j] = A[i][j]; }
    }
   if (i == S) { Buff = MultiplyMatrix(M,a,Size,Size,Size,Size,Sa,Sa1); }
    else
      {
       Buff = MultiplyMatrix(M,Buff1,Size,Size,Size,Size,Sa,Sa1);
       if (Buff1 != NULL) { for (long i = 0; i < Size; i++) { free(Buff1[i]); Buff1[i] = 0; } free(Buff1); Buff1 = 0; }
      }
   for (long j = 0; j < Size; j++) //опять матрицу к диагональному виду привожу
    {
     if (j != i-1) { M[i-1][j] = 0.0; }
      else { M[i-1][j] = 1.0; }
    }
   if (i == S) //заполнение матрицы M(n-1)
    {
     for (long j = 0; j < Size; j++)
      {
       if (j == (i-1))
        {
         if (a[i][i-1] == 0.0) M[i-1][j] = 0.0;
          else M[i-1][j] = 1/a[i][i-1];
        }
       else
        {
         if (a[i][i-1] == 0.0) M[i-1][j] = 0.0;
          else M[i-1][j] = -1*(a[i][j]/a[i][i-1]);
        }
      }
    }
   else
    {
     for (long j = 0; j < Size; j++)
      {
       if (j == (i-1))
        {
         if (A[i][i-1] == 0.0) M[i-1][j] = 0.0;
          else M[i-1][j] = 1/A[i][i-1];
        }
       else
        {
         if (A[i][i-1] == 0.0) M[i-1][j] = 0.0;
          else M[i-1][j] = -1*(A[i][j]/A[i][i-1]);
        }
      }
    }
   Buff1 = MultiplyMatrix(Buff,M,Size,Size,Size,Size,Sa,Sa1);
   for (long j = 0; j < Size; j++) { free(Buff[j]); Buff[j] = 0; } free(Buff); Buff = 0; //чищу матрицу Buff вместе с памятью
   //Второй блок приведения к матрице Фронебиуса
   for (long j = 0; j < Size; j++) //обратное приведение матрицы к диагональному виду
    {
     if (j != i-1) { M[i-1][j] = 0.0; }
      else { M[i-1][j] = 1.0; }
    }
   if (i == S) //заполнение матрицы M^(-1)(n-2)
    {
     for (long j = 0; j < Size; j++) { M[i-2][j] = a[i-1][j]; }
    }
   else
    {
     for (long j = 0; j < Size; j++) { M[i-2][j] = A[i-1][j]; }
    }
   Buff = MultiplyMatrix(M,Buff1,Size,Size,Size,Size,Sa,Sa1);
   if (Buff1 != NULL) { for (long i = 0; i < Size; i++) { free(Buff1[i]); Buff1[i] = 0; } free(Buff1); Buff1 = 0; } //чищу память от буферной матрицы Buff1
   for (long j = 0; j < Size; j++) //обратное приведение матрицы к диагональному виду
    {
     if (j != (i-2)) { M[i-2][j] = 0.0; }
      else { M[i-2][j] = 1.0; }
    }
   if (i == S) //заполнение матрицы M(n-2)
    {
     for (long j = 0; j < Size; j++)
      {
       if (j == (i-2))
        {
         if (a[i-1][i-2] == 0.0) { M[i-2][j] = 0.0; }
          else { M[i-2][j] = 1/a[i-1][i-2]; }
        }
       else
        {
         if (a[i-1][i-2] == 0.0) { M[i-2][j] = 0.0; }
          else { M[i-2][j] = -1*(a[i-1][j])/a[i-1][i-2]; }
        }
      }
    }
   else
    {
     for (long j = 0; j < Size; j++)
      {
       if (j == (i-2))
        {
         if (A[i-1][i-2] == 0.0) { M[i-2][j] = 0.0; }
          else { M[i-2][j] = 1/A[i-1][i-2]; }
        }
       else
        {
         if (A[i-1][i-2] == 0.0) { M[i-2][j] = 0.0; }
          else { M[i-2][j] = -1*(A[i-1][j])/A[i-1][i-2]; }
        }
      }
    }
   Buff1 = MultiplyMatrix(Buff,M,Size,Size,Size,Size,Sa,Sa1);
   if (Buff != NULL) { for (long i = 0; i < Size; i++) { free(Buff[i]); Buff[i] = 0; } free(Buff); Buff = 0; }
   for (long j = 0; j < Size; j++)  //снова матрицу к диагональному виду привожу
    {
     if (j != (i-2)) { M[i-2][j] = 0.0; }
      else { M[i-2][j] = 1.0; }
    }
   //промежуточную матрицу А заполняю
   if (i != 2) //лишние телодвижения пресекаю
    {
     for (long j = 0; j < Size; j++)
      {
       for (long k = 0; k < Size; k++)
        {
         A[j][k] = Buff1[j][k];
        }
      }
    }
  }
 //чищу память
 if (M != NULL) { for (long i = 0; i < Size; i++) { free(M[i]); M[i] = 0; } free(M); M = 0; }
 if (A != NULL) { for (long i = 0; i < Size; i++) { free(A[i]); A[i] = 0; } free(A); A = 0; }
 //результат работы функции
 Vec = (float*)malloc(sizeof(float)*Size);
 for (long i = 0; i < Size; i++)
  {
   Vec[i] = Buff1[0][i];   
  } 
 if (Buff1 != NULL) { for (long i = 0; i < Size; i++) { free(Buff1[i]); Buff1[i] = 0; } free(Buff1); Buff1 = 0; }
 return Vec;
}