polynomial.hpp

Ugrás a fájl dokumentációjához.

#ifndef POLYNOMIAL_H
#define POLYNOMIAL_H
 
#include <vector>
#include <iostream>
using namespace std;
 
namespace bme {

    template <typename T>
    class Polynomial {
    private:
        vector<T> coef;
        int dg;
        
        void deg_check() {
            dg = -1;
            for (int k = coef.size()-1; k >= 0 ; k--) {
                if (coef[k] != T(0)) {
                    dg = k;
                    break;
                }
            }
            setCoefVectorSize(dg);
        }
 
        void setCoefVectorSize(int _dg) {
            if (_dg < -1) {
                return;
            }
            else if (_dg == -1) {
                clear();
                return;
            }
            size_t new_deg = static_cast<size_t>(_dg);
            coef.resize(new_deg + 1);
        }
        
        static Polynomial<T> InvalidPoly() {
            Polynomial<T> ans;
            ans.dg = -2;
            return ans;
        }
 
        
    public:

        Polynomial<T>() {
            clear();
        }
        

        template <typename U>
        Polynomial<T>(U a) {
            coef.resize(1);
            coef[0] = T(a);
            deg_check();
        }


        template<typename U, typename V, typename W>
        Polynomial<T>(U a, V b, W c) {
            coef.resize(3);
            coef[2] = T(a);
            coef[1] = T(b);
            coef[0] = T(c);
            deg_check();
        }


        template <typename U>
        Polynomial<T>(int array_size, const U * array) {
            if (array_size <= 0) {
                coef.clear();
                return;
            }
            coef.resize(array_size);
            for (int k = 0; k < array_size; k++) {
                coef[k] = T(array[k]);
            }
            deg_check();
        }


        int deg() const {
            return dg; 
        }


        bool isInvalid()const {
            return dg == -2;
        }
        
        void setInvalidPoly() {
            coef.resize(1);
            coef[0] = T(0);
            dg = -2;
        }


        const T& getCoef(int deg) const {
            if ((deg < 0) || (deg > dg)) { return T(0); }
            return coef[deg];
        }


        const T& operator[](int deg) const {
            return getCoef(deg); 
        }


        void set(int deg, const T& val) {
            if (isInvalid()) {
                return;
            }
            if (deg < 0) { return; }
            const size_t ix = static_cast<size_t>(deg);
            if (ix + 1 > coef.size()) {
                coef.resize(ix + 1);
            }
            coef[ix] = val;
            deg_check();
        }


        void clear() {
            coef.resize(1);
            dg = -1;
            coef[0] = T(0);
        }


        void minimize() {
            if (isInvalid()) {
                return;
            }
            const size_t rsize = static_cast<size_t>(deg()) + 1;
            coef.reserve(rsize);
        }


        bool isZero() const {
            return (deg() == -1); 
        }


        void shift(int n = 1) {
            if (n == 0 || isInvalid()) return;
            if (n < -dg) {
                clear();
                return;
            }
            else if (n < 0) {
                for (int k = 0; k <= dg + n; k++) coef[k] = coef[k - n];
                for (int k = dg + n + 1; k <= dg; k++) coef[k] = T(0);
                deg_check();
                return;
            }
            if (static_cast<size_t>(n) + static_cast<size_t>(deg()) > INT_MAX) {
                setInvalidPoly();
                return;
            }
            const int old_deg = deg();
            setCoefVectorSize(n + old_deg);
            for (int k = old_deg; k >= 0; k--) coef[k + n] = coef[k];
            for (int k = 0; k < n; k++) coef[k] = T(0);
            dg += n;
        }
 
        // FUNCTION APPLICATION //


        T of(const T& a) const {
            if (isInvalid()) {
                return T(0);
            }
            if (deg() <= 0) { return coef[0]; }
            T ans;
            ans = T(0);
            for (int k = deg(); k >= 0; k--) {
                ans *= a;
                ans += coef[k];
            }
            return ans;
        }


        T operator()(const T& a) const { return of(a); }


        Polynomial<T> of(const Polynomial<T>& other) const {
            if (isInvalid() || other.isInvaid()) {
                return InvalidPoly();
            }
            if (deg() <= 0) { return Polynomial<T>(coef[0]); }
            Polynomial<T> ans(T(0));
            for (int k = deg(); k >= 0; k--) {
                ans *= other;
                ans += Polynomial<T>(coef[k]);
            }
            return ans;
        }


        Polynomial<T> operator()(const Polynomial<T>& other) const {
            return of(other);
        }
 
        // COMPARISON //
 


        bool operator==(const Polynomial<T>& other) const {
            if (deg() != other.deg()) return false;
            for (int k = 0; k <= dg; k++) {
                if (coef[k] != other.coef[k]) return false;
            }
            return true;
        }
 
        bool operator!=(const Polynomial<T>& other) const {
            return !(*this == other);
        }


        bool operator<(const Polynomial<T>& other) const {
            if (isInvalid()) {
                return false;
            }
            else if (deg() < other.deg()) {
                return true;
            }
            else if (deg() > other.deg()) {
                return false;
            }
            for (int i = 0; i <= deg(); ++i) {
                const T& c1 = getCoef(i);
                const T& c2 = other.getCoef(i);
                if (c1 < c2) {
                    return true;
                }
                else if (c1 > c2) {
                    return false;
                }
            }
            return false;
        }
 
        bool operator<=(const Polynomial<T>& other) const {
            return (*this) < other || (*this) == other;
        }
 
        bool operator>(const Polynomial<T>& other) const {
            return !((*this) <=  other);
        }
 
        bool operator>=(const Polynomial<T>& other) const {
            return !((*this) < other);
        }
 
        // ARITHMETIC //
 
        Polynomial<T> operator+(const Polynomial<T>& other) const {
            if (isInvalid() || other.isInvalid()) {
                return InvalidPoly();
            }
            Polynomial<T> ans;
            const int dmax = (deg() > other.deg()) ? deg() : other.deg();
            ans.setCoefVectorSize(dmax);
            for (int k = 0; k <= dmax; k++) {
                ans.coef[k] = getCoef(k) + other.getCoef(k);
            }
            ans.deg_check();
            return ans;
        }
 
        Polynomial<T>& operator+=(const Polynomial<T>& other) {
            (*this) = (*this) + other;
            return *this;
        }
 
        Polynomial<T> operator-() const {
            if (isInvalid()) {
                return InvalidPoly();
            }
            Polynomial<T> ans;
            ans.setCoefVectorSize(deg());
            for (int k = 0; k <= dg; k++) { ans.coef[k] = -coef[k]; }
            ans.deg_check();
            return ans;
        }
 
        Polynomial<T> operator-(const Polynomial<T>& other) const {
            if (isInvalid() || other.isInvalid()) {
                return InvalidPoly();
            }
            Polynomial<T> ans;
            const int dmax = (deg() > other.deg()) ? deg() : other.deg();
            ans.setCoefVectorSize(dmax);
            for (int k = 0; k <= dmax; k++) {
                ans.coef[k] = getCoef(k) - other.getCoef(k);
            }
            ans.deg_check();
            return ans;
        }
 
        Polynomial<T>& operator-=(const Polynomial<T>& other) {
            (*this) = (*this) - other;
            return *this;
        }


        Polynomial<T> operator*(const Polynomial<T>& other) const {
            if (isInvalid() || other.isInvalid()) {
                return InvalidPoly();
            }
            Polynomial<T> ans;
            if (isZero() || other.isZero()) { return ans; }
            if (static_cast<size_t>(deg()) + static_cast<size_t>(other.deg()) > INT_MAX) {
                return InvalidPoly();
            }
            const int dans = deg() + other.deg();
            ans.setCoefVectorSize(dans);
            for (int k = 0; k <= dans; k++) {
                T c(0);
                for (int j = 0; j <= k; j++) {
                    if ((j <= dg) && (k - j <= other.dg)) {
                        c += coef[j] * other.coef[k - j];
                    }
                }
                ans.coef[k] = c;
            }
            ans.deg_check();
            return ans;
        }
 
        Polynomial<T>& operator*=(const Polynomial<T>& other) {
            *this = (*this) * other;
            return *this;
        }


        Polynomial<T> pow(int k) const {
            if (k < 0 || isInvalid()) {
                Polynomial<T> ans;
                return InvalidPoly();
            }
            if (k == 0) { return Polynomial<T>(1); }
            if (k == 1) { return *this; }
 
            if (k % 2 == 0) {
                int half_k = k / 2;
                Polynomial<T> ans = (*this).pow(half_k);
                ans *= ans;
                return ans;
            }
            int half_k = (k - 1) / 2;
            Polynomial<T> ans = (*this).pow(half_k);
            ans *= ans;
            ans *= *this;
            return ans;
        }


        Polynomial<T> operator/(const Polynomial<T>& other) const {
            if (other.isZero() || isInvalid() || other.isInvalid()) {
                return InvalidPoly();
            }
            Polynomial<T> Q, R;
            quot_rem(*this, other, Q, R);
            return Q;
        }
 
        Polynomial<T>& operator/=(const Polynomial<T>&other) {
            *this = (*this) / other;
            return *this;
        }


        Polynomial<T> operator%(const Polynomial<T>&other) const {
            if (isInvalid() || other.isInvalid() || other.isZero()) {
                return InvalidPoly();
            }
            Polynomial<T> Q, R;
            quot_rem(*this, other, Q, R);
            return R;
        }
 
        Polynomial<T>& operator%=(const Polynomial<T>& other) {
            (*this) = (*this) % other;
            return *this;
        }


        void make_monic() {
            if (deg() < 0) { return; }
            T lead = coef[dg];
            for (int j = 0; j <= dg; j++) { coef[j] /= lead; }
        }

        bool is_monic() const {
            if (deg() < 0) { return false; } 
            return coef[dg] == T(1);
        }
 
    }; // end of Polynomial<T> class template

    template<typename T>
    ostream& operator<<(ostream & os, const Polynomial<T>& P) {
        if (P.isInvalid()) {
            os << "INVALID" << endl;
            return os;
        }
        if (P.deg() <= 0) {
            os << "(" << P[0] << ")";
            return os;
        }
        for (int k = P.deg(); k >= 0; k--) {
            if (P[k] != T(0)) {
                if (k < P.deg()) os << " + ";
                os << "(" << P[k] << ")";
                if (k > 1) {
                    os << "X^" << k;
                    continue;
                }
                if (k == 1) {
                    os << "X";
                }
            }
        }
        return os;
    }
 
    template<typename U, typename T>
    Polynomial<T> operator+(U x, const Polynomial<T>& P){ 
        return P + T(x); 
    }
 
    template<typename U, typename T>
    Polynomial<T> operator-(U x, const Polynomial<T>& P){
        return (-P) + T(x); 
    }
 
    template<typename U, typename T>
    Polynomial<T> operator*(U x, const Polynomial<T>& P){
        return P * T(x); 
    }
    

    template<typename T>
    void quot_rem(const Polynomial<T>& A, const Polynomial<T>& B,
                    Polynomial<T>& Q, Polynomial<T>& R) {
        if (B.isZero() || B.isInvalid() || A.isInvalid()) {
            Q.setInvalidPoly();
            R.setInvalidPoly();
            return;
        }
        Q.clear();
        R.clear();
        Polynomial<T> AA(A); // az A másolata
        while (AA.deg() >= B.deg()) {
            int k = AA.deg() - B.deg();
            Polynomial<T> BB = B;
            BB.shift(k);
            T a_lead = AA[AA.deg()];
            T b_lead = BB[BB.deg()];
            for (int j = 0; j <= BB.deg(); j++) {
                BB.set(j, BB[j] * a_lead / b_lead);
            }
            AA -= BB;
            Q.set(k, a_lead / b_lead);
        }
        R = A - Q * B;
    }


    template<typename T>
    Polynomial<T> gcd(const Polynomial<T>& A, const Polynomial<T>& B) {
        if (A.isInvalid() || B.isInvalid() || (A.isZero() && B.isZero())) {
            Polynomial<T> ans;
            ans.setInvalidPoly();
            return ans;
        }
        if (B.isZero()) {
            if (A.is_monic()) { return A; }
            Polynomial<T> AA(A);
            AA.make_monic();
            return AA;
        }
        Polynomial<T> C;
        C = A % B;
        return gcd(B, C);
    }


    template<typename T>
    Polynomial<T> gcd(const Polynomial<T>& A, const Polynomial<T>& B,
                        Polynomial<T>& X, Polynomial<T>& Y) {
        Polynomial<T> D; // ebbe kerül a visszatérési érték
 
        if ((A.isZero() && B.isZero()) || A.isInvalid || B.isInvalid()) {
            X.setInvalidPoly();
            Y.setInvalidPoly();
            D.setInvalidPoly();
            return D;
        }
 
        // Ha A nem 0, de B igen, D = A/a_lead, X = a_lead, Y = 0
        if (B.isZero()) {
            D = A;
            T a_lead = A[A.deg()];
            D.make_monic();
            X = Polynomial<T>(T(1) / a_lead);
            Y.clear();
            return D;
        }
 
        // ha sem A, sem B nem 0, akkor rekurzívan számolunk
        Polynomial<T> Q;
        Polynomial<T> R;
        quot_rem(A, B, Q, R);
        Polynomial<T> XX;
        Polynomial<T> YY;
        D = gcd(B, R, XX, YY);
        
        X = YY;
        Y = XX - Q * YY;
        
        return D;
    }
 
}//namespace bme 
 
#endif