Polynom.hpp

Go to the documentation of this file.

#pragma once
 
// std
#include <cstdint>
#include <ostream>
#include <initializer_list>
#include <iterator>
#include <concepts>
#include <complex>
#include <limits>
 
// eigen
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Eigenvalues> 
 
#include "math.hpp"
 
namespace controlpp
{
 
    template<class T, int Order>
    class Polynom{
        public:
            using value_type = T;
            static constexpr int eigen_vector_size = (Order == Eigen::Dynamic) ? Eigen::Dynamic : Order + 1;
            using vector_type = Eigen::Vector<T, eigen_vector_size>;
 
        private:
            vector_type vector_;
 
        public:
 
            Polynom() = default;
 
            Polynom(const Polynom&) = default;
 
            template<int OtherOrder>
            requires(OtherOrder < Order)
            Polynom(const Polynom<T, OtherOrder>& other){
                size_t i = 0;
                for(; i < other.size(); ++i){
                    this->at(i) = other.at(i);
                }
                for(; i < this->size(); ++i){
                    this->at(i) = static_cast<T>(0);
                }
            }
 
            template<class... Args>
            requires(Order != Eigen::Dynamic && (std::convertible_to<Args, T> && ...))
            explicit Polynom(Args&&... args) : vector_(std::forward<Args>(args)...){}
 
            template<class... Args>
            requires(Order == Eigen::Dynamic && (std::convertible_to<Args, T> && ...))
            explicit Polynom(Args&&... args) : vector_(Eigen::Vector<T, sizeof...(Args)>(std::forward<Args>(args)...)){}
 
 
            Polynom& operator=(const Polynom&) = default;
 
            template<int M>
            requires(M < Order)
            Polynom& operator=(const Polynom<T, M>& other){
                size_t i = 0;
                for(; i < other.size(); ++i){
                    this->at(i) = other.at(i);
                }
                for(; i < this->size(); ++i){
                    this->at(i) = static_cast<T>(0);
                }
                return *this;
            }
 
            template<int M>
            requires((Order != Eigen::Dynamic) && (M == Order+1))
            explicit Polynom(const T (&array)[M]) : vector_(array){}
 
            template<int M>
            requires((Order != Eigen::Dynamic) && (M < Order+1))
            explicit Polynom(const T (&array)[M]){
                size_t i = 0;
                for(; i < M; ++i){
                    this->vector_[i] = array[i];
                } 
                for(; i < this->size(); ++i){
                    this->vector_[i] = static_cast<T>(0);
                } 
            }
            
            explicit Polynom(const T* values, size_t length){
                if constexpr (Order == Eigen::Dynamic){
                    this->vector_.resize(length);
                }
                size_t i = 0;
                for(; i < length && i < this->size(); ++i){
                    this->vector_(i) = values[i];
                }
                for(; i < this->size(); ++i){
                    this->vector_(i) = static_cast<T>(0);
                }
            }
 
            template<int M>
            requires(Order == Eigen::Dynamic)
            explicit Polynom(const T (&array)[M])
                : Polynom(array, M){}
 
            template<int M>
            requires((Order != Eigen::Dynamic) && (M != Eigen::Dynamic) && (M == Order+1))
            explicit Polynom(const Eigen::Vector<T, M>& vector) : vector_(vector){}
 
            template<class U>
            requires(std::constructible_from<Eigen::Vector<T, eigen_vector_size>, U>)
            explicit Polynom(const U& vector_expression) : vector_(vector_expression){}
 
            template<int M>
            requires((Order != Eigen::Dynamic) && (M != Eigen::Dynamic) && (M < Order+1))
            explicit Polynom(const Eigen::Vector<T, M>& vector){
                size_t i = 0;
                for(; i < vector.size(); ++i){
                    this->vector_[i] = vector(i);
                } 
                for(; i < this->size(); ++i){
                    this->vector_[i] = static_cast<T>(0);
                } 
            }
 
            template<std::same_as<T> U>
            requires(Order == Eigen::Dynamic)
            explicit Polynom(const Eigen::Vector<U, Eigen::Dynamic>& vector) : vector_(vector.size()){
                this->vector_ = vector;
            }
 
            template<int M>
            requires(Order == Eigen::Dynamic && M != Eigen::Dynamic)
            explicit Polynom(const Eigen::Vector<T, M>& vector) : vector_(M){
                this->vector_ = vector;
            }
 
            const T& at(size_t i) const {return this->vector_[i];}
 
            T& at(size_t i) {return this->vector_[i];}
 
            const T& operator[](size_t i) const {return this->at(i);}
 
            T& operator[](size_t i) {return this->at(i);}
 
            vector_type& vector(){return this->vector_;}
 
            const vector_type& vector() const {return this->vector_;}
 
            size_t size() const {return this->vector_.size();}
 
            size_t order() const {
                size_t result = this->size()-1;
                for(size_t i = 0; i < this->size(); ++i){
                    const T& elem = this->at(this->size() - 1 - i);
                    const T zero(0);
                    if(elem != zero){
                        break;
                    }
                    --result;
                }
                return result;
            }
 
            void setZero() {this->vector_.setZero();}
            
            T eval(const T& x) const {
                T sum = static_cast<T>(0);
                for(int i = this->size()-1; i >= 0; --i){
                    sum = sum * x + this->at(i);
                }
                return sum;
            }
 
            T operator() (const T& x) const {
                return this->eval(x);
            }
 
            template<int M>
            Eigen::Vector<T, M> eval(const Eigen::Vector<T, M>& x_vec) const {
                Eigen::Vector<T, M> sum_vec;
                if constexpr (M == Eigen::Dynamic) sum_vec.resize(x_vec.size());
                sum_vec.setZero();
 
                for(int i = this->size()-1; i >= 0; --i){
                    sum_vec.array() = sum_vec.array() * x_vec.array() + this->at(i);
                }
 
                return sum_vec;
            }
 
            template<int M>
            Eigen::Vector<T, M> operator() (const Eigen::Vector<T, M>& x_vec) const {
                return this->eval(x_vec);
            }
 
            std::complex<T> eval(const std::complex<T>& x) const {
                std::complex<T> sum(static_cast<T>(0), static_cast<T>(0));
                for(int i = this->size()-1; i >= 0; --i){
                    sum = sum * x + this->at(i);
                }
                return sum;
            }
 
            std::complex<T> operator() (const std::complex<T>& x) const {
                return this->eval(x);
            }
 
            template<int M>
            Eigen::Vector<std::complex<T>, M> eval(const Eigen::Vector<std::complex<T>, M>& x_vec) const {
                Eigen::Vector<std::complex<T>, M> sum_vec;
                if constexpr (M == Eigen::Dynamic) sum_vec.resize(x_vec.size());
                sum_vec.setZero();
 
                for(int i = this->size()-1; i >= 0; --i){
                    sum_vec.array() = sum_vec.array() * x_vec.array() + this->at(i);
                }
 
                return sum_vec;
            }
 
            template<int M>
            Eigen::Vector<std::complex<T>, M> operator() (const Eigen::Vector<std::complex<T>, M>& x_vec) const {
                return this->eval(x_vec);
            }
 
            bool is_zero(T epsilon = std::numeric_limits<T>::min()) const {
                for(size_t i = 0; i < this->size(); ++i){
                    if((this->at(i) >= epsilon)){
                        return false;
                    }
                }
                return true;
            }
 
            void print (std::ostream& stream, std::string_view var="x") const {
                if(this->is_zero()){
                    stream << '0';
                    return;
                }
                for(size_t i = 0; i < this->size(); ++i){
                    if(this->at(i) != static_cast<T>(0)){
                        const auto value = this->at(i);
                        if(value > 0){
                            if(i > 0) stream << " + ";
                            
                            stream << value;
 
                            if(i == 1) stream << ' ' << var;
                            else if(i > 1) stream << ' ' << var << '^' << i;
                        }else if(value < 0){
                            if(i == 0) stream << "- ";
                            else if(i > 0) stream << " - ";
 
                            stream << (-value);
 
                            if(i == 1) stream << ' ' << var;
                            else if(i > 1) stream << ' ' << var << '^' << i;
                        }
                    }
                }
            }
 
            void print (std::ostream& stream, std::function<std::string(int i)> var) const {
                if(this->is_zero()){
                    stream << '0';
                    return;
                }
                for(size_t i = 0; i < this->size(); ++i){
                    if(this->at(i) != static_cast<T>(0)){
                        const auto value = this->at(i);
                        if(value > 0){
                            if(i > 0) stream << " + ";
                            
                            stream << value;
 
                            if(i == 1) stream << ' ' << var(i);
                            else if(i > 1) stream << ' ' << var(i) << '^' << i;
                        }else if(value < 0){
                            if(i == 0) stream << "- ";
                            else if(i > 0) stream << " - ";
 
                            stream << (-value);
 
                            if(i == 1) stream << ' ' << var(i);
                            else if(i > 1) stream << ' ' << var(i) << '^' << i;
                        }
                    }
                }
            }
 
            friend std::ostream& operator<< (std::ostream& stream, const Polynom& poly){
                poly.print(stream, "x");
                return stream;
            }
 
            template<int M>
            requires(M <= Order)
            Polynom& operator+=(const Polynom<T, M>& other){
                for(size_t i = 0; i < other.size(); ++i){
                    this->at(i) += other.at(i);
                }
                return *this;
            }
 
            template<int M>
            requires(M <= Order)
            Polynom& operator-=(const Polynom<T, M>& other){
                for(size_t i = 0; i < other.size(); ++i){
                    this->at(i) -= other.at(i);
                }
                return *this;
            }
 
            template<std::convertible_to<T> U>
            Polynom& operator*=(const U& num){
                this->vector() *= static_cast<T>(num);
                return *this;
            }
 
            template<std::convertible_to<T> U>
            Polynom& operator/=(const U& num){
                this->vector() /= static_cast<T>(num);
                return *this;
            }
    };
 
    // -----------------------------------------------------------------------------------------------
    //                                        Roots, Zeros
    // -----------------------------------------------------------------------------------------------
    
    template<class T>
    Eigen::Vector<std::complex<T>, 0> zeros([[maybe_unused]] const Polynom<T, 0>& polynom){
        return Eigen::Vector<std::complex<T>, 0>();
    }
 
    template<class T>
    Eigen::Vector<std::complex<T>, 1> zeros(const Polynom<T, 1>& polynom){
        const std::complex x0 = - polynom[0] / polynom[1];
        Eigen::Vector<std::complex<T>, 1> result(x0);
        return result;
    }
 
    template<class T>
    Eigen::Vector<std::complex<T>, 2> zeros(const Polynom<T, 2>& polynom){
        const std::complex c = polynom[0];
        const std::complex b = polynom[1];
        const std::complex a = polynom[2];
        const std::complex x1 = (-b - std::sqrt(b * b - static_cast<T>(4) * a * c)) / (static_cast<T>(2) * a);
        const std::complex x2 = (-b + std::sqrt(b * b - static_cast<T>(4) * a * c)) / (static_cast<T>(2) * a);
        const Eigen::Vector<std::complex<T>, 2> result(x1, x2);
        return result;
    }
 
    template<class T, int N>
    requires(N > 1)
    Eigen::Vector<std::complex<T>, N> zeros(const Polynom<T, N>& polynom){
        const Eigen::Matrix<T, N, N> C = controlpp::companion(polynom.vector());
        const Eigen::Vector<std::complex<T>, N> result = C.eigenvalues();
        return result;
    }
 
    // -----------------------------------------------------------------------------------------------
    //                                      Comparison Operators
    // -----------------------------------------------------------------------------------------------
 
    template<class T, int N>
    bool operator==(const Polynom<T, N>& lhs, const Polynom<T, N>& rhs){
        return lhs.vector() == rhs.vector();
    }
 
    template<class T, int N>
    bool operator!=(const Polynom<T, N>& lhs, const Polynom<T, N>& rhs){
        return lhs.vector() != rhs.vector();
    }
 
    // -----------------------------------------------------------------------------------------------
    //                                      Arithmetic Operators
    // -----------------------------------------------------------------------------------------------
 
    // operator +
    // ----------
 
    template<class T, int Nl, int Nr>
    Polynom<T, (Nl > Nr) ? Nl : Nr> operator+(const Polynom<T, Nl>& lhs, const Polynom<T, Nr>& rhs){
        if constexpr (Nl > Nr){
            Polynom<T, Nl> result;
            result.vector().head(rhs.size()) = lhs.vector().head(rhs.size()) + rhs.vector();
            result.vector().tail(lhs.size() - rhs.size()) = lhs.vector().tail(lhs.size() - rhs.size());
            return result;
        }else if constexpr (Nl == Nr){
            return Polynom<T, Nl>(lhs.vector() + rhs.vector());
        }else{
            Polynom<T, Nr> result;
            result.vector().head(lhs.size()) = lhs.vector() + rhs.vector().head(lhs.size());
            result.vector().tail(rhs.size() - lhs.size()) = rhs.vector().tail(rhs.size() - lhs.size());
            return result;
        }
    }
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    Polynom<Tpoly, N> operator+(const Tscalar& lhs, const Polynom<Tpoly, N>& rhs){
        Polynom<Tpoly, N> result(rhs);
        result[0] += static_cast<Tpoly>(lhs);
        return result;
    }
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    Polynom<Tpoly, N> operator+(const Polynom<Tpoly, N>& lhs, const Tscalar& rhs){
        return ((static_cast<Tpoly>(rhs) + lhs).eval());
    }
 
    // operator -
    // ----------
 
    template<class T, int N>
    Polynom<T, N> operator-(const Polynom<T, N>& poly){
        return Polynom<T, N>((-poly.vector()).eval());
    }
 
    template<class T, int Nl, int Nr>
    Polynom<T, (Nl > Nr) ? Nl : Nr> operator-(const Polynom<T, Nl>& lhs, const Polynom<T, Nr>& rhs){
        if constexpr (Nl > Nr){
            Polynom<T, Nl> result;
            result.vector().head(rhs.size()) = lhs.vector().head(rhs.size()) - rhs.vector();
            result.vector().tail(lhs.size() - rhs.size()) = lhs.vector().tail(lhs.size() - rhs.size());
            return result;
        }else if constexpr (Nl == Nr){
            return Polynom<T, Nl>(lhs.vector() - rhs.vector());
        }else{
            Polynom<T, Nr> result;
            result.vector().head(lhs.size()) = lhs.vector() + rhs.vector().head(lhs.size());
            result.vector().tail(rhs.size() - lhs.size()) = -rhs.vector().tail(rhs.size() - lhs.size());
            return result;
        }
    }
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    Polynom<Tpoly, N> operator-(const Tscalar& lhs, const Polynom<Tpoly, N>& rhs){
        Polynom<Tpoly, N> result(-rhs);
        result[0] += static_cast<Tpoly>(lhs);
        return result;
    }
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    Polynom<Tpoly, N> operator-(const Polynom<Tpoly, N>& lhs, const Tscalar& rhs){
        Polynom<Tpoly, N> result(lhs);
        result[0] -= static_cast<Tpoly>(rhs);
        return result;
    }
 
    // operator *
    // ----------
 
    template<class T, int lOrder, int rOrder>
    Polynom<T, lOrder + rOrder> operator*(const Polynom<T, lOrder>& lhs, const Polynom<T, rOrder>& rhs){
        Polynom<T, lOrder + rOrder> result;
        result.setZero();
 
        for(size_t i = 0; i < rhs.size(); ++i){
            result.vector().segment(i, lhs.size()) += lhs.vector() * rhs[i];
        }
 
        return result;
    }
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int Order>
    Polynom<Tpoly, Order> operator*(const Polynom<Tpoly, Order>& lhs, const Tscalar& rhs){
        return Polynom<Tpoly, Order>((lhs.vector() * static_cast<Tpoly>(rhs)).eval());
    }
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int Order>
    Polynom<Tpoly, Order> operator*(const Tscalar& lhs, const Polynom<Tpoly, Order>& rhs){
        return Polynom<Tpoly, Order>(rhs * static_cast<Tpoly>(lhs));
    }
 
    // operator /
    // ----------
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int Order>
    Polynom<Tpoly, Order> operator/(const Polynom<Tpoly, Order>& lhs, const Tscalar& rhs){
        return Polynom<Tpoly, Order>(lhs.vector() / static_cast<Tscalar>(rhs));
    }
 
    
    namespace polynom{
 
        template<class T>
        inline const Polynom<T, 1> x(0, 1);
    }
 
// -----------------------------------------------------------------------------------------------------------------
// -----------------------------------------------------------------------------------------------------------------
// -----------------------------------------------------------------------------------------------------------------
 
    template<class T, int Order>
    class FixedPolynom{
        public:
            using value_type = T;
            using vector_type = Eigen::Vector<T, Order+1>;
 
        private:
            vector_type vector_;
 
        public:
 
            FixedPolynom() = default;
 
            FixedPolynom(const FixedPolynom&) = default;
 
            template<std::convertible_to<T> U>
            FixedPolynom(const U& value) : vector_(Eigen::Vector<T, Order+1>::Zero()){
                vector_(0) = static_cast<T>(value);
            }
 
            template<int M>
            requires(M < Order)
            FixedPolynom(const FixedPolynom<T, M>& other){
                size_t i = 0;
                for(; i < other.size(); ++i){
                    this->at(i) = other.at(i);
                }
                for(; i < this->size(); ++i){
                    this->at(i) = static_cast<T>(0);
                }
            };
 
            FixedPolynom& operator=(const FixedPolynom&) = default;
 
            template<int M>
            requires(M < Order)
            FixedPolynom& operator=(const FixedPolynom<T, M>& other){
                size_t i = 0;
                for(; i < other.size(); ++i){
                    this->at(i) = other.at(i);
                }
                for(; i < this->size(); ++i){
                    this->at(i) = static_cast<T>(0);
                }
                return *this;
            };
 
 
            explicit FixedPolynom(const T (&values)[Order+1]){
                for(size_t i = 0; i < this->size(); ++i){
                    this->vector_[i] = values[i];
                }
            }
 
            template<int M>
            requires(M < Order+1)
            explicit FixedPolynom(const T (&values)[M]){
                size_t i = 0;
                for(; i < M; ++i){
                    this->vector_[i] = values[i];
                }
                for(; i < this->size(); ++i){
                    this->vector_[i] = static_cast<T>(0);
                }
            }
 
            explicit FixedPolynom(const Eigen::Vector<T, Order+1>& vector) : vector_(vector){}
 
            template<int M>
            requires(M < Order+1)
            explicit FixedPolynom(const Eigen::Vector<T, M>& vector){
                size_t i = 0;
                for(; i < vector.size(); ++i){
                    this->vector_[i] = vector[i];
                }
                for(; i < this->size(); ++i){
                    this->vector_[i] = static_cast<T>(0);
                }
            }
            
            FixedPolynom& operator=(const T (&values)[Order+1]) {
                this->vector_ = values;
                return *this;
            }
 
            template<int M>
            requires(M < Order+1)
            FixedPolynom& operator=(const T (&values)[M]) {
                size_t i = 0;
                for(; i < M; ++i){
                    this->vector_[i] = values[i];
                }
                for(; i < this->size(); ++i){
                    this->vector_[i] = static_cast<T>(0);
                }
                return *this;
            }
            
            explicit FixedPolynom(const Polynom<T, Order>& poly) : FixedPolynom(poly.vector()){}
 
            template<int M>
            requires(M < Order)
            explicit FixedPolynom(const Polynom<T, M>& poly){
                size_t i = 0;
                for(; i < M; ++i){
                    this->vector_[i] = poly[i];
                }
                for(; i < this->size(); ++i){
                    this->vector_[i] = static_cast<T>(0);
                }
            }
 
            operator Polynom<T, Order>() const {return Polynom<T, Order>(this->vector_);}
 
            FixedPolynom& operator=(const Eigen::Vector<T, Order+1>& vector){
                this->vector_ = vector;
                return *this;
            }
 
            template<int M>
            requires(M < Order+1)
            FixedPolynom& operator=(const Eigen::Vector<T, M>& vector){
                size_t i = 0;
                for(; i < M; ++i){
                    this->vector_[i] = vector[i];
                }
                for(; i < this->size(); ++i){
                    this->vector_[i] = static_cast<T>(0);
                }
                return *this;
            }
 
            const T& operator[](size_t i) const {return this->vector_[i];}
 
            T& operator[](size_t i) {return this->vector_[i];}
 
            const T& at(size_t i) const {return this->vector_[i];}
 
            T& at(size_t i) {return this->vector_[i];}
 
            vector_type& vector(){return this->vector_;}
 
            const vector_type& vector() const {return this->vector_;}
 
            size_t size() const {return Order+1;}
 
            size_t order() const {
                size_t result = Order;
                const T zero(0);
                for(size_t i = 0; i < this->size(); ++i){
                    if(this->at(this->size() - i - 1) != zero){
                        break;
                    }
                    --result;
                }
                return result;
            }
 
            void setZero() {this->vector_.setZero();}
            
            void print (std::ostream& stream, std::string_view var="x") const {
                std::string_view plus = "";
                for(size_t i = 0; i < this->size(); ++i){
                    if(this->at(i) != static_cast<T>(0)){
                        stream << plus << this->at(i) << ' ' << var << '^' << i;
                        plus = " + ";
                    }
                }
            }
 
            friend std::ostream& operator<< (std::ostream& stream, const FixedPolynom& poly){
                poly.print(stream, "x");
                return stream;
            }
 
            template<int M>
            requires(M <= Order)
            FixedPolynom& operator+=(const FixedPolynom<T, M>& other){
                for(size_t i = 0; i < other->size(); ++i){
                    this->at(i) += other.at(i);
                }
                return *this;
            }
 
            template<int M>
            requires(M <= Order)
            FixedPolynom& operator-=(const FixedPolynom<T, M>& other){
                for(size_t i = 0; i < other->size(); ++i){
                    this->at(i) -= other.at(i);
                }
                return *this;
            }
    };
 
    // -----------------------------------------------------------------------------------------------
    //                                      Comparison Operators
    // -----------------------------------------------------------------------------------------------
 
    template<class T, int N>
    bool operator==(const FixedPolynom<T, N>& lhs, const FixedPolynom<T, N>& rhs){
        return lhs.vector() == rhs.vector();
    }
 
    template<class T, int N>
    bool operator!=(const FixedPolynom<T, N>& lhs, const FixedPolynom<T, N>& rhs){
        return lhs.vector() != rhs.vector();
    }
 
    // -----------------------------------------------------------------------------------------------
    //                                      Arithmetic Operators
    // -----------------------------------------------------------------------------------------------
 
    // operator +
    // ----------
 
    template<class T, int N>
    FixedPolynom<T, N> operator+(const FixedPolynom<T, N>& lhs, const FixedPolynom<T, N>& rhs){
        return FixedPolynom<T, N>(lhs.vector() + rhs.vector());
    }
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    requires (N >= 1)
    FixedPolynom<Tpoly, N> operator+(const Tscalar& lhs, const FixedPolynom<Tpoly, N>& rhs){
        FixedPolynom<Tpoly, N> result(rhs);
        result[0] += static_cast<Tpoly>(lhs);
        return result;
    }
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    requires (N >= 1)
    FixedPolynom<Tpoly, N> operator+(const FixedPolynom<Tpoly, N>& lhs, const Tscalar& rhs){
        return (static_cast<Tpoly>(rhs) + lhs);
    }
 
    // operator -
    // ----------
 
    template<class T, int N>
    FixedPolynom<T, N> operator-(const FixedPolynom<T, N>& lhs, const FixedPolynom<T, N>& rhs){
        return FixedPolynom<T, N>(lhs.vector() - rhs.vector());
    }
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    requires (N >= 1)
    FixedPolynom<Tpoly, N> operator-(const Tscalar& lhs, const FixedPolynom<Tpoly, N>& rhs){
        FixedPolynom<Tpoly, N> result(-rhs);
        result[0] += static_cast<Tpoly>(lhs);
        return result;
    }
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    requires (N >= 1)
    FixedPolynom<Tpoly, N> operator-(const FixedPolynom<Tpoly, N>& lhs, const Tscalar& rhs){
        FixedPolynom<Tpoly, N> result(lhs);
        result[0] -= static_cast<Tpoly>(rhs);
        return result;
    }
 
    template<class T, int N>
    FixedPolynom<T, N> operator-(const FixedPolynom<T, N>& values){
        return FixedPolynom<T, N>(-values.vector());
    }
 
    // operator *
    // ----------
 
    template<class T, int N>
    FixedPolynom<T, N> operator*(const FixedPolynom<T, N>& lhs, const FixedPolynom<T, N>& rhs){
        FixedPolynom<T, N> result;
        result.setZero();
 
        for(size_t i = 0; i < lhs.size(); ++i){
            result.vector().tail(lhs.size()-i) += lhs.vector().head(lhs.size()-i) * rhs[i];
        }
 
        return result;
    }
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    FixedPolynom<Tpoly, N> operator*(const FixedPolynom<Tpoly, N>& lhs, const Tscalar& rhs){
        return FixedPolynom<Tpoly, N>(lhs.vector() * static_cast<Tpoly>(rhs));
    }
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    FixedPolynom<Tpoly, N> operator*(const Tscalar& lhs, const FixedPolynom<Tpoly, N>& rhs){
        return (rhs * static_cast<Tpoly>(lhs));
    }
 
    // operator /
    // ----------
 
    template<class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    FixedPolynom<Tpoly, N> operator/(const FixedPolynom<Tpoly, N>& lhs, const Tscalar& rhs){
        return FixedPolynom<Tpoly, N>(lhs.vector() / static_cast<Tscalar>(rhs));
    }
} // namespace controlpp
 
namespace Eigen {
    template<typename T, int N>
    struct NumTraits<controlpp::FixedPolynom<T, N>> {
        using Self = controlpp::FixedPolynom<T, N>;
        using Real = Self; // or T if you want to treat real parts differently
        using NonInteger = Self;
        using Literal = Self;
        using Nested = Self;
 
        enum {
            IsComplex = 0,
            IsInteger = 0,
            IsSigned = 1,
            RequireInitialization = 1,
            ReadCost = 1,
            AddCost = 5,
            MulCost = 10
        };
 
        static inline Self epsilon() { return Self(); }
        static inline Self dummy_precision() { return Self(); }
        static inline Self highest() { return Self(); }
        static inline Self lowest() { return Self(); }
    };
 
    // ResultType for operator +
    // -------------------------
 
    template <class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    struct ScalarBinaryOpTraits<
        Tscalar,
        controlpp::FixedPolynom<Tpoly, N>,
        internal::scalar_sum_op<Tscalar, controlpp::FixedPolynom<Tpoly, N>>> 
    {
        using ReturnType = decltype(std::declval<Tscalar>() + std::declval<controlpp::FixedPolynom<Tpoly, N>>());
    };
 
    template <class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    struct ScalarBinaryOpTraits<
        controlpp::FixedPolynom<Tpoly, N>,
        Tscalar,
        internal::scalar_sum_op<controlpp::FixedPolynom<Tpoly, N>, Tscalar>> 
    {
        using ReturnType = decltype(std::declval<controlpp::FixedPolynom<Tpoly, N>>() + std::declval<Tscalar>());
    };
 
    template <class Tpoly, int N>
    struct ScalarBinaryOpTraits<
        controlpp::FixedPolynom<Tpoly, N>,
        controlpp::FixedPolynom<Tpoly, N>,
        internal::scalar_sum_op<controlpp::FixedPolynom<Tpoly, N>, controlpp::FixedPolynom<Tpoly, N>>> 
    {
        using ReturnType = decltype(std::declval<controlpp::FixedPolynom<Tpoly, N>>() + std::declval<controlpp::FixedPolynom<Tpoly, N>>());
    };
 
    // ResultType for operator -
    // -------------------------
 
    template <class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    struct ScalarBinaryOpTraits<
        Tscalar,
        controlpp::FixedPolynom<Tpoly, N>,
        internal::scalar_difference_op<Tscalar, controlpp::FixedPolynom<Tpoly, N>>> 
    {
        using ReturnType = decltype(std::declval<Tscalar>() - std::declval<controlpp::FixedPolynom<Tpoly, N>>());
    };
 
    template <class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    struct ScalarBinaryOpTraits<
        controlpp::FixedPolynom<Tpoly, N>,
        Tscalar,
        internal::scalar_difference_op<controlpp::FixedPolynom<Tpoly, N>, Tscalar>> 
    {
        using ReturnType = decltype(std::declval<controlpp::FixedPolynom<Tpoly, N>>() - std::declval<Tscalar>());
    };
 
    template <class Tpoly, int N>
    struct ScalarBinaryOpTraits<
        controlpp::FixedPolynom<Tpoly, N>,
        controlpp::FixedPolynom<Tpoly, N>,
        internal::scalar_difference_op<controlpp::FixedPolynom<Tpoly, N>, controlpp::FixedPolynom<Tpoly, N>>> 
    {
        using ReturnType = decltype(std::declval<controlpp::FixedPolynom<Tpoly, N>>() - std::declval<controlpp::FixedPolynom<Tpoly, N>>());
    };
 
    // ResultType for operator *
    // -------------------------
 
    template <class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    struct ScalarBinaryOpTraits<
        Tscalar,
        controlpp::FixedPolynom<Tpoly, N>,
        internal::scalar_product_op<Tscalar, controlpp::FixedPolynom<Tpoly, N>>> 
    {
        using ReturnType = decltype(std::declval<Tscalar>() * std::declval<controlpp::FixedPolynom<Tpoly, N>>());
    };
 
    template <class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    struct ScalarBinaryOpTraits<
        controlpp::FixedPolynom<Tpoly, N>,
        Tscalar,
        internal::scalar_product_op<controlpp::FixedPolynom<Tpoly, N>, Tscalar>> 
    {
        using ReturnType = decltype(std::declval<controlpp::FixedPolynom<Tpoly, N>>() * std::declval<Tscalar>());
    };
 
    template <class Tpoly, int N>
    struct ScalarBinaryOpTraits<
        controlpp::FixedPolynom<Tpoly, N>,
        controlpp::FixedPolynom<Tpoly, N>,
        internal::scalar_product_op<controlpp::FixedPolynom<Tpoly, N>, controlpp::FixedPolynom<Tpoly, N>>> 
    {
        using ReturnType = decltype(std::declval<controlpp::FixedPolynom<Tpoly, N>>() * std::declval<controlpp::FixedPolynom<Tpoly, N>>());
    };
 
    // ResultType for operator /
    // -------------------------
 
    template <class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    struct ScalarBinaryOpTraits<
        Tscalar,
        controlpp::FixedPolynom<Tpoly, N>,
        internal::scalar_quotient_op<Tscalar, controlpp::FixedPolynom<Tpoly, N>>> 
    {
        using ReturnType = decltype(std::declval<Tscalar>() / std::declval<controlpp::FixedPolynom<Tpoly, N>>());
    };
 
    template <class Tpoly, std::convertible_to<Tpoly> Tscalar, int N>
    struct ScalarBinaryOpTraits<
        controlpp::FixedPolynom<Tpoly, N>,
        Tscalar,
        internal::scalar_quotient_op<controlpp::FixedPolynom<Tpoly, N>, Tscalar>> 
    {
        using ReturnType = decltype(std::declval<controlpp::FixedPolynom<Tpoly, N>>() / std::declval<Tscalar>());
    };
 
    template <class Tpoly, int N>
    struct ScalarBinaryOpTraits<
        controlpp::FixedPolynom<Tpoly, N>,
        controlpp::FixedPolynom<Tpoly, N>,
        internal::scalar_quotient_op<controlpp::FixedPolynom<Tpoly, N>, controlpp::FixedPolynom<Tpoly, N>>> 
    {
        using ReturnType = decltype(std::declval<controlpp::FixedPolynom<Tpoly, N>>() / std::declval<controlpp::FixedPolynom<Tpoly, N>>());
    };
}// Eigen