#ifndef POLYNOMIAL_H
#define POLYNOMIAL_H 1

#include <vector>
#include <string>
#include <iostream>
#include <utility>
#include <cmath>
#include <cassert>

using std::size_t;

constexpr double POLYNOMIAL_EPS = 1e-8;

class Polynomial: public std::vector<double> {
private:
    // The vector contains coefficients of the polynomial
    // in ascending order.
    // For example, the polynomial
    // x^3 - x - 1 has coefficients
    // {-1, -1, 0, 1}
public:
    Polynomial(int deg = 0):
        std::vector<double>(deg + 1)
    {
    }

    Polynomial(double freeTerm = 0.):
        std::vector<double>(1, freeTerm)
    {
    }

    Polynomial(const std::vector<double>& coeffs):
        std::vector<double>(coeffs)
    {
    }

    Polynomial(std::vector<double>&& coeffs):
        std::vector<double>(coeffs)
    {
    }

    // Coefficients in this constructor are given in ascending order
    Polynomial(int deg, const double* coeffs):
        std::vector<double>(coeffs, coeffs + deg + 1)
    {
        trim();
    }

    int degree() const {
        return int(size() - 1);
    }

    // Remove upper zero coefficients
    void trim() {
        int i = degree();
        while (i > 0) {
            if (fabs((*this)[i]) > POLYNOMIAL_EPS)
                break;
            --i;
        }
        assert(i >= 0);
        resize(i + 1);
    }

    Polynomial& operator+=(const Polynomial& f) {
        int d0 = degree();
        int d1 = f.degree();
        int d = d0;
        if (d1 > d)
            d = d1;
        resize(d + 1);
        for (int i = 0; i <= d; ++i) {
            if (i <= d0) {
                if (i <= d1)
                    (*this)[i] += f[i];
            } else {
                assert(i <= d1);
                (*this)[i] = f[i];
            }
        }
        trim();
        return *this;
    }

    Polynomial& operator-=(const Polynomial& f) {
        int d0 = degree();
        int d1 = f.degree();
        int d = d0;
        if (d1 > d)
            d = d1;
        resize(d + 1);
        for (int i = 0; i <= d; ++i) {
            if (i <= d0) {
                if (i <= d1)
                    (*this)[i] -= f[i];
            } else {
                assert(i <= d1);
                (*this)[i] = (-f[i]);
            }
        }
        trim();
        return *this;
    }

    Polynomial& operator*=(const Polynomial& f) {
        int d0 = degree();
        int d1 = f.degree();
        int d = d0 + d1;
        std::vector<double> product(d + 1, 0.);
        for (int i = 0; i <= d0; ++i) {
            for (int j = 0; j <= d1; ++j) {
                product[i + j] += (*this)[i]*f[j];
            }
        }
        (*this) = std::move(product);
        return *this;
    }

    std::string toString() const;
};

inline Polynomial operator+(const Polynomial& f, const Polynomial& g) {
    Polynomial res = f;
    res += g;
    return res;
}

inline Polynomial operator-(const Polynomial& f, const Polynomial& g) {
    Polynomial res = f;
    res -= g;
    return res;
}

inline Polynomial operator*(const Polynomial& f, const Polynomial& g) {
    Polynomial res = f;
    res *= g;
    return res;
}

inline std::ostream& operator<<(std::ostream& s, const Polynomial& f) {
    s << f.toString();
    return s;
}

std::istream& operator>>(std::istream& s, Polynomial& f);

#endif /* POLYNOMIAL_H */