polynomial_root.h

Go to the documentation of this file.

/* === S Y N F I G ========================================================= */
/* ========================================================================= */

/* === S T A R T =========================================================== */

#ifndef __SYNFIG_POLYNOMIAL_ROOT_H
#define __SYNFIG_POLYNOMIAL_ROOT_H

/* === H E A D E R S ======================================================= */

#include <complex>
#include <vector>

/* === M A C R O S ========================================================= */

/* === T Y P E D E F S ===================================================== */

/* === C L A S S E S & S T R U C T S ======================================= */
template < typename T = float, typename F = float >
class Polynomial : public std::vector<T> //a0 + a1x + a2x^2 + ... + anx^n
{
public:

    //Will maintain all lower constants
    void degree(unsigned int d, const T & def = (T)0) { resize(d+1,def); }
    unsigned int degree()const { return this->size() - 1; }

    const Polynomial & operator+=(const Polynomial &p)
    {
        if(p.size() > this->size())
            resize(p.size(), (T)0);

        for(int i = 0; i < p.size(); ++i)
        {
            (*this)[i] += p[i];
        }
        return *this;
    }

    const Polynomial & operator-=(const Polynomial &p)
    {
        if(p.size() > this->size())
            resize(p.size(), (T)0);

        for(int i = 0; i < p.size(); ++i)
        {
            (*this)[i] -= p[i];
        }
        return *this;
    }

    const Polynomial & operator*=(const Polynomial &p)
    {
        if(p.size() < 1)
        {
            this->resize(0);
            return *this;
        }

        unsigned int i,j;
        std::vector<T> nc(*this);

        //in place for constant stuff
        for(i = 0; i < nc.size(); ++i)
        {
            (*this)[i] *= p[0];
        }

        if(p.size() < 2) return *this;

        this->resize(this->size() + p.degree());
        for(i = 0; i < nc.size(); ++i)
        {
            for(j = 1; j < p.size(); ++j)
            {
                nc[i+j] += nc[i]*p[j];
            }
        }

        return *this;
    }
};

class RootFinder
{
    std::vector< std::complex<float> >  workcoefs;
    //int   its;

public:
    std::vector< std::complex<float> >  coefs; //the number of coefficients determines the degree of polynomial

    std::vector< std::complex<float> >  roots;

    void find_all_roots(bool polish);
};



/* === E N D =============================================================== */

#endif