polynomial.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_POLYNOMIAL_H
#define THEORETICA_POLYNOMIAL_H
 
#ifndef THEORETICA_NO_PRINT
#include <sstream>
#include <ostream>
#endif
 
#include "../core/real_analysis.h"
#include "../algebra/vec.h"
#include "../complex/complex.h"
#include "../complex/complex_analysis.h"
 
 
namespace theoretica {
 
    template<typename Type = real>
    class polynomial {
 
    private:
 
        vec<Type> coeff {Type(0.0)};
 
    public:
 
        polynomial() = default;
 
        polynomial(Type a) : coeff({a}) {}
 
        template<typename Vector, enable_vector<Vector> = true>
        explicit polynomial(const Vector& c) {
            
            coeff.resize(c.size());
 
            for (size_t i = 0; i < c.size(); ++i)
                coeff[i] = c[i];
        }
        
 
        polynomial(std::initializer_list<Type> l) : coeff(l) {}
 
 
        template<typename ...Args>
        polynomial(Args... args) : coeff({args...}) {}
 
 
        inline const Type& at(unsigned int i) const {
            return coeff.at(i);
        }
 
        inline Type& at(unsigned int i) {
            return coeff.at(i);
        }
 
 
        inline const Type& operator[](unsigned int i) const {
            return coeff[i];
        }
 
 
        inline Type& operator[](unsigned int i) {
            return coeff[i];
        }
 
 
        template<typename InputType = Type>
        inline Type eval(InputType x) const {
 
            Type sum = Type(0.0);
 
            // Evaluate using Horner's method
            for (unsigned int i = 0; i < coeff.size(); ++i)
                sum = coeff[coeff.size() - i - 1] + sum * x;
 
            return sum;
        }
 
 
        template<typename InputType = Type>
        inline Type operator()(InputType x) const {
            return eval(x);
        }
 
 
        inline unsigned int degree(real tolerance = MACH_EPSILON) const {
 
            for (int i = coeff.size() - 1; i >= 0; --i) {
                if(abs(coeff[i]) > MACH_EPSILON)
                    return i;
            }
 
            return 0;
        }
 
 
        inline void trim(real tolerance = MACH_EPSILON) {
 
            size_t new_size = coeff.size();
            for (size_t i = coeff.size() - 1; i > 0; --i) {
 
                if(abs(coeff[i]) > tolerance)
                    break;
 
                new_size--;
            }
 
            if (new_size == 0) {
                coeff.resize(1);
                coeff[0] = Type(0.0);
                return;
            }
 
            if (new_size < coeff.size())
                coeff.resize(new_size);
        }
 
 
        inline const vec<Type>& coeffs() const {
            return coeff;
        }
 
 
        inline vec<Type>& coeffs() {
            return coeff;
        }
 
 
        inline size_t size() const {
            return coeff.size();
        }
 
 
        inline void resize(size_t sz) {
            coeff.resize(sz);
        }
 
 
        inline polynomial operator+(const polynomial& p) const {
 
            polynomial r;
            r.coeff.resize(max(coeff.size(), p.size()));
 
            unsigned int i = 0;
 
            for (; i < min(coeff.size(), p.size()); ++i)
                r[i] = coeff[i] + p[i];
 
            if (coeff.size() > p.size())
                for (; i < coeff.size(); ++i)
                    r[i] = coeff[i];
            else
                for (; i < p.size(); ++i)
                    r[i] = p[i];
 
            return r;
        }
 
 
        inline polynomial operator-(const polynomial& p) const {
            
            polynomial r;
            r.coeff.resize(max(coeff.size(), p.size()));
 
            unsigned int i = 0;
 
            for (; i < min(r.size(), p.size()); ++i)
                r[i] = coeff[i] - p[i];
 
            if (coeff.size() > p.size())
                for (; i < coeff.size(); ++i)
                    r[i] = coeff[i];
            else
                for (; i < p.size(); ++i)
                    r[i] = -p[i];
 
            return r;
        }
 
 
        inline polynomial operator*(const polynomial& p) const {
 
            polynomial r = polynomial();
            r.coeff.resize(this->size() + p.size() - 1);
 
            for (unsigned int i = 0; i < size(); ++i) {
                for (unsigned int j = 0; j < p.size(); ++j) {
                    r[i + j] += coeff[i] * p[j];
                }
            }
 
            return r;
        }
 
 
        inline polynomial operator/(const polynomial& d) const {
 
            const unsigned int d_order = d.degree();
            const unsigned int this_order = degree();
 
            if(d_order == 0 && d[0] == 0) {
                TH_MATH_ERROR("polynomial::operator/", d[0], MathError::DivByZero);
                return polynomial(nan());
            }
 
            // Remainder
            polynomial r = *this;
 
            // Quotient
            polynomial q = polynomial();
            unsigned int i = 0;
 
            while(i <= this_order) {
 
                // Compute only once the degree of the polynomial
                const unsigned int r_order = r.degree();
 
                // Stop execution if the division is complete
                // (when the remainder is 0 or has lower degree)
                if((r_order == 0 && (abs(r[0]) < MACH_EPSILON)) || r_order < d_order)
                    break;
 
                // Simple division between highest degree terms
                const polynomial t = polynomial<Type>::monomial(
                    r[r_order] / d[d_order],
                    r_order - d_order);
 
                // Add monomial to quotient and subtract the
                // monomial times the dividend from the remainder
                q += t;
                r -= t * d;
 
                i++;
            }
 
            // Check if division failed to converge
            if(i > this_order && r.degree() >= d_order) {
                TH_MATH_ERROR("polynomial::operator/", i, MathError::NoConvergence);
                return polynomial(nan());
            }
 
            return q;
        }
 
 
        template<typename ScalarType, disable_vector<ScalarType> = true>
        inline polynomial operator*(ScalarType a) const {
            return polynomial(coeff * a);
        }
 
 
        template<typename ScalarType, disable_vector<ScalarType> = true>
        inline polynomial operator/(ScalarType a) const {
 
            if(abs(a) < MACH_EPSILON) {
                TH_MATH_ERROR("polynomial::operator/", a, MathError::DivByZero);
            }
 
            polynomial r = polynomial(*this);
 
            for (unsigned int i = 0; i < size(); ++i)
                r.coeff[i] /= a;
 
            return r;
        }
 
 
        inline polynomial& operator+=(const polynomial& p) {
 
            // Make room for the new coefficients
            if(coeff.size() < p.size())
                coeff.resize(p.size(), Type(0));
 
            for (unsigned int i = 0; i < min(size(), p.size()); ++i)
                coeff[i] += p[i];
 
            return *this;
        }
 
 
        inline polynomial& operator-=(const polynomial& p) {
 
            // Make room for the new coefficients
            if(coeff.size() < p.size())
                coeff.resize(p.size(), Type(0));
 
            for (unsigned int i = 0; i < min(coeff.size(), p.size()); ++i)
                coeff[i] -= p[i];
 
            return *this;
        }
 
 
        inline polynomial& operator*=(const polynomial& p) {
 
            polynomial r = polynomial();
            r.coeff.resize(this->size() + p.size() - 1, Type(0));
 
            for (unsigned int i = 0; i < coeff.size(); ++i)
                for (unsigned int j = 0; j < p.size(); ++j)
                    r[i + j] += coeff[i] * p[j];
 
            *this = r;
            return *this;
        }
 
 
        template<typename ScalarType, disable_vector<ScalarType> = true>
        inline polynomial& operator*=(ScalarType a) {
 
            for (unsigned int i = 0; i < coeff.size(); ++i)
                coeff[i] *= a;
            
            return *this;
        }
 
 
        inline polynomial& operator/=(const polynomial& a) {
            return (*this = (*this / a));
        }
 
 
        template<typename ScalarType, disable_vector<ScalarType> = true>
        inline polynomial& operator/=(ScalarType a) {
 
            if(abs(a) < MACH_EPSILON) {
                TH_MATH_ERROR("polynomial::operator/=", a, MathError::DivByZero);
                return *this;
            }
 
            for (unsigned int i = 0; i < coeff.size(); ++i)
                coeff[i] /= a;
            
            return *this;
        }
 
 
        inline bool operator==(const polynomial& other) const {
 
            if (other.size() == 0 && size() == 0)
                return true;
 
            const size_t other_degree = other.degree();
            const size_t this_degree = degree();
 
            if (this_degree != other_degree)
                return false;
 
            for (size_t i = 0; i <= this_degree; ++i)
                if (abs(coeff[i] - other[i]) > MACH_EPSILON)
                    return false;
 
            return true;
        }
 
 
        inline auto begin() {
            return coeff.begin();
        }
 
 
        inline auto end() {
            return coeff.end();
        }
 
 
        inline auto cbegin() const {
            return coeff.cbegin();
        }
 
 
        inline auto cend() const {
            return coeff.cend();
        }
 
 
        inline vec<complex<>, 2> quadratic_roots() const {
 
            const int order = degree();
 
            // Check that the polynomial is quadratic
            if(order != 2) {
                TH_MATH_ERROR("quadratic_roots", order, MathError::ImpossibleOperation);
                return vec<complex<>, 2>({nan(), nan()});
            }
 
            complex<> z1, z2;
            const Type a = coeff[2];
            const Type b = coeff[1];
            const Type c = coeff[0];
 
            // Linear case
            if (abs(a) < MACH_EPSILON) {
 
                if (abs(b) < MACH_EPSILON)
                    return {nan(), nan()};
 
                return {-c/b, -c/b};
            }
 
            const Type discriminant = b * b - 4 * a * c;
            const complex<> delta = sqrt(complex<>(discriminant));
 
            const real sign = (b > 0) ? 1 : -1;
            const complex<> q = -0.5 * (b + sign * delta);
 
            if (abs(q) < MACH_EPSILON)
                return {0.0, 0.0};
 
            z1 = q / a;
            z2 = c / q;
 
            return {z1, z2};
        }
 
 
        template<typename Vector, enable_vector<Vector> = true>
        inline static polynomial<Type> from_roots(const Vector& roots) {
 
            polynomial<Type> p;
            p.resize(roots.size() + 1);
            p[0] = Type(1.0);
 
            for (size_t i = 0; i < roots.size(); i++)
                p *= polynomial<Type>({-roots[i], 1});
            
            return p;
        }
 
 
        inline static polynomial<Type> monomial(Type c, unsigned int order) {
            
            polynomial m;
            m.coeff = vec<Type>(order + 1, Type(0.0));
            m.coeff[order] = c;
            return m;
        }
 
 
        // Friend operators to enable equations of the form
        // (Type) op. (polynomial<Type>)
 
        inline friend polynomial<Type> operator+(
            Type r, const polynomial<Type>& z) {
 
            return z + polynomial(r);
        }
 
        inline friend polynomial<Type> operator-(
            Type r, const polynomial<Type>& z) {
 
            return (z * -1) + polynomial(r);
        }
 
        template<typename ScalarType, disable_vector<ScalarType> = true>
        inline friend polynomial<Type> operator*(ScalarType r, const polynomial<Type>& z) {
            return z * r;
        }
 
 
#ifndef THEORETICA_NO_PRINT
 
        inline std::string to_string(
            const std::string& unknown = "x",
            const std::string& exponentiation = "^") const {
 
            std::stringstream res;
            const int sz = coeff.size();
 
            if (size() == 0) {
                res << Type(0.0);
                return res.str();
            }
 
            res << coeff[sz - 1];
            
            if (size() > 1)
                res << "*" << unknown << exponentiation << (sz - 1) << " ";
 
            for (int i = sz - 2; i >= 0; --i) {
 
                if(abs(coeff[i]) < MACH_EPSILON)
                    continue;
 
                res << (coeff[i] >= 0 ? "+ " : "- ");
                res << abs(coeff[i]);
 
                if(i) {
                    res << "*" << unknown << exponentiation << i;
                    res << " ";
                }
            }
 
            return res.str();
        }
 
 
        inline operator std::string() {
            return to_string();
        }
 
 
        friend std::ostream& operator<<(std::ostream& out, const polynomial& obj) {
            return out << obj.to_string();
        }
 
#endif
 
        
    };
    
}
 
#endif