LINP V2.0

#include <algorithm>
#include <iostream>
#include <vector>
#include <complex>
#include <functional>
#include <stdexcept>

using namespace std;
class LINP
{
public:
    // Index => Coeff
    typedef vector<double> cont_type;

public:
    LINP(): cont() {}
    LINP(double rhs): cont(1) {
        cont[0] = rhs;
    }

    LINP operator [](size_t rhs) const {
        if (rhs == 0)
            return LINP(1);
        LINP tmp(*this);
        while (--rhs != 0)
            tmp *= *this;
        return tmp;
    }

    LINP& operator *=(double rhs) {
        for (auto& v : cont)
            v *= rhs;
        return *this;
    }
    LINP operator *(double rhs) const {
        LINP tmp(*this);
        tmp *= rhs;
        return tmp;
    }

    LINP& operator *=(const LINP& rhs) {
        if (rhs.cont.empty()) {
            cont.clear();
            return *this;
        } else if (this == &rhs) {
            return *this *= LINP(*this);
        }
        shrinkToFit();
        cont.reserve(rhs.maxIndex() + maxIndex());

        LINP tmp(*this);
        cont.clear();
        for (auto& v : rhs.cont) {
            if (v != 0)
                *this += tmp * v;
            tmp.cont.insert(tmp.cont.begin(), 0);
        }
        return *this;
    }
    LINP operator *(const LINP& rhs) const {
        LINP tmp(*this);
        tmp *= rhs;
        return tmp;
    }

    LINP& operator +=(const LINP& rhs) {
        // Seems safe when this == &rhs
        if (cont.size() < rhs.cont.size())
            cont.resize(rhs.cont.size());
        auto it = cont.begin();
        for (auto& v : rhs.cont)
            *it++ += v;
        return *this;
    }
    LINP operator +(const LINP& rhs) const {
        LINP tmp(*this);
        tmp += rhs;
        return tmp;
    }

    LINP& operator -=(const LINP& rhs) {
        if (this == &rhs) {
            cont.clear();
        } else {
            if (cont.size() < rhs.cont.size())
                cont.resize(rhs.cont.size());
            auto it = cont.begin();
            for (auto& v : rhs.cont)
                *it++ -= v;
        }
        return *this;
    }
    LINP operator -(const LINP& rhs) const {
        LINP tmp(*this);
        tmp -= rhs;
        return tmp;
    }
    LINP operator -() const {
        LINP tmp(*this);
        for (auto& v : tmp.cont)
            v = -v;
        return tmp;
    }

    friend ostream& operator <<(ostream& os, const LINP& rhs) {
        if (rhs.cont.size() > 1) {
            for (cont_type::size_type idx = rhs.cont.size() - 1; idx != 0; --idx) {
                double coeff = rhs.cont[idx];
                if (coeff == 0)
                    continue;
                if (coeff < 0) {
                    os << "-";
                    if (coeff != -1)
                        os << -coeff;
                } else {
                    os << "+";
                    if (coeff != 1)
                        os << coeff;
                }
                os << 'X';
                if (idx != 1)
                    os << '[' << idx << ']';
            }
        }
        if (!rhs.cont.empty()) {
            double coeff = rhs.cont.front();
            if (coeff != 0) {
                if (coeff > 0)
                    os << '+';
                os << coeff;
            }
        } else {
            os << '0';
        }
        return os;
    }

    void shrinkToFit() {
        cont.resize(maxIndex() + 1);
    }

    size_t maxIndex() const {
        return distance(find_if(cont.rbegin(), cont.rend(), [] (double coeff) {
            return coeff != 0;
        }), cont.rend()) - 1;
    }

    double getCoeffOfIndex(size_t Index) const {
        return Index > maxIndex() ? cont[Index] : 0;
    }

    vector<complex<double> > solve() const throw(invalid_argument)  {
        size_t maxIdx = maxIndex();
        if (maxIdx > 2)
            throw invalid_argument("只支持一次和二次多项式");
        vector<complex<double> > result;
        result.reserve(maxIdx);
        if (maxIdx == 1) {
            result.push_back(-cont[0] / cont[1]);
        } else {
            complex<double> delta=pow(cont[1], 2) - 4 * cont[2] * cont[0];
            result.push_back((+sqrt(delta) - cont[1]) / (2 * cont[2]));
            result.push_back((-sqrt(delta) - cont[1]) / (2 * cont[2]));
        }
        return result;
    }

    static LINP getX() {
        LINP x;
        x.cont.resize(2);
        x.cont.back() = 1;
        return x;
    }

private:
    cont_type cont;
};

int main()
{
    const static LINP X = LINP::getX();
    LINP expression = (X + 1)[3];
    cout << expression;

    LINP linp=X-(X*0.5+1)-1;
    cout << linp << " = 0" << endl;
    cout << "X = " << linp.solve().front().real() << endl << endl;

    LINP linp2=X[2] + X*2 + 1;
    cout << linp2 << " = 0" << endl;
    cout << "X1 = X2 = " << linp2.solve().front().real() << endl << endl;

    LINP linp3=X[2] + 1;
    cout << linp3 << " = 0" << endl;
    vector<complex<double> > result=linp3.solve();
    cout << "x1 = " << result[0].real() << " + " << result[0].imag() << 'i' << endl
         << "x2 = " << result[1].real() << " + " << result[1].imag() << 'i' << endl
         << endl;
}