//
//  charPoly.cpp
//  fp
//
//  Created by Robert Delaney on 2/4/17.
//
//

/*
Finds the characteristic polynomial of a square matrix using a method described by
Prof. Bumby of Rutgers Univ.

Let M be the matrix of size nxn, and the characteristic polynomial be:

x^n - p[1]*x^n-1 - p[2]*x^n-2 - ... - p[n], and let p[0] = 1

Then define matrices F such that

F_1(M) = M

and for k = 1, 2,...,n

p[k] = (1/k)*Tr(F_k)

F_k+1 = M*(F_k - p[k]*I)

*/

#include "charPoly.hpp"
#include "matc.hpp"
#include "matcConv.hpp"
#include "matfp.hpp"
#include "matcMath.hpp"

// M is input as an n x n matrix; p is created and returned as a n+1 dim. array
// containing the coeffs. of the characterestic polynomial of M with p[n] = 1 being the coeff of x^n
bool charPoly(matc& M, fpComplex*& p)
{
    INT32       i, k, n;
    matc        temp, F;
    fpComplex   temp1;
    
    if(M.nr != M.nc)
        return false;
    
    n = M.nr;
    
    p = (fpComplex*)malloc((n+1)*sizeof(fpComplex));
    
    for(i=0;i<=n;++i)
    {
        init(p[i].re);
        init(p[i].im);
    }
    
    p[0].re = 1.;
    p[0].im = 0.;
    
    F = M;
    for(k=1;k<n;k++)
    {
        temp1 = trace(F);
        p[k] = temp1 / k;
        for(i=0;i<n;i++)
            F.array[i][i] = F.array[i][i] - p[k];
        temp = M * F;
        F  = temp;
    }
    
    temp1 = trace(F);
    p[n] = temp1 / n;
    
    // change signs except for p[0]
    for(i=1;i<=n;i++)
        p[i] = -p[i];
    
    // now reverse p[] to get coefficients in standard form
    for(i=0;i<=n/2;i++)
    {
        temp1 = p[i];
        p[i] = p[n-i];
        p[n-i] = temp1;
    }
    
    return true;
    
}/* charPoly */


// M is input as an n x n matrix; p is created and returned as a n+1 dim. array
// containing the coeffs. of the characterestic polynomial of M with p[n] = 1 being the coeff of x^n
bool charPoly(matfp& M, fpComplex*& p)
{
    matc    z;
    
    conv(z, M);  // has init!
    
    return charPoly(z, p);
    
}/* charPoly */

