//#####################################################################
// Copyright 2002-2005, Ronald Fedkiw, Geoffrey Irving.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Core/Arrays/ARRAY_VIEW.h>
#include <Core/Matrices/LAPACK.h>
#include <cassert>
namespace PhysBAM{
template<class T>
void Polynomial_Roots(ARRAY_VIEW<const T> coeffs,ARRAY_VIEW<std::complex<T> > roots)
{
    assert(roots.m==coeffs.m-1);
    int n=roots.m;
    MATRIX_MXN<T> M(n,n);
    for(int i=0;i<n-1;i++) M(i+1,i)=1;
    for(int i=0;i<n;i++) M(i,n-1)=-coeffs(i)/coeffs(n);
    Compute_Eigenvalues_Raw(M,roots);
}

template<class T>
int Polynomial_Roots(ARRAY_VIEW<const T> coeffs,ARRAY_VIEW<T> roots)
{
    assert(roots.m==coeffs.m-1);
    ARRAY<std::complex<T> > r(coeffs.m-1);
    Polynomial_Roots(coeffs,ARRAY_VIEW<std::complex<T> >(r));
    int k=0;
    for(int i=0;i<r.m;i++)
        if(!imag(r(i)))
            roots(k++)=real(r(i));
    return k;
}
template void Polynomial_Roots<double>(ARRAY_VIEW<double const>,
    ARRAY_VIEW<std::complex<double>>);
template void Polynomial_Roots<float>(ARRAY_VIEW<float const>,
    ARRAY_VIEW<std::complex<float>>);
template int Polynomial_Roots<double>(ARRAY_VIEW<double const>,
    ARRAY_VIEW<double>);
template int Polynomial_Roots<float>(ARRAY_VIEW<float const>,
    ARRAY_VIEW<float>);
}
