[C++] 다항식 클래스와 라그랑주 보간법(Lagrange Interpolation) 구현

다음 문제를 풀다가 다이나믹 프로그래밍으로 푸는 법이 생각이 안 나서 Lagrange Interpolation으로 밀어붙이려고 했다.

사실 생긴 것만 봐도 다이나믹 어쩌고저쩌고보다는 수학 문제를 닮아 있는 것 같지 않은가? 사실 이 문제에 대한 일반화는 이미 존재하고, 베르누이 수열인가 뭔가를 쓴다고 알고 있지만 그걸 굳이 알고 싶지는 않았다. 그리고 도전 중인 문제가 하나 더 있는데 그 문제에서도 Lagrange Interpolation을 쓰기 때문에 어차피 한 번은 이걸 구현해야 해서, 그냥 오늘 구현하기로 했다.

특이사항으로, class Polynomial을 구현할 때 std::vector를 쓰지 않았다. C++의 이동 연산에 대해 공부하고 있었기 때문에 이동 연산을 구현해야 했기 때문이다.

시간복잡도는 O(n^2).

참고로 Default Constructor는 덧셈의 항등원 0을 만드는 것이 아니라 곱셈의 항등원 1을 만든다.

using namespace std;

template<class Type>
class Polynomial {
    Type* p;
    int size;
public:
    explicit Polynomial(Type* t, int deg) : size(deg + 1) {
        p = new Type[deg + 1];
        for (int i = 0; i <= deg; ++i) {
            p[i] = t[i];
        }
    }
    Polynomial() : size(1) {
        p = new Type[1];
        p[0] = 1;
    }
private:
    explicit Polynomial(int _size) : size(_size) {
        p = new Type[_size];
    }
public:
    Polynomial(const Polynomial& rhs) : size(rhs.size) {
        p = new Type[rhs.size];
        for (int i = 0; i < rhs.size; ++i) {
            p[i] = rhs.p[i];
        }
    }
    Polynomial(Polynomial&& rhs) noexcept : size(rhs.size) {
        p = rhs.p;
        rhs.p = NULL;
    }
    Polynomial& operator=(const Polynomial& rhs) {
        size = rhs.size;
        p = new Type[rhs.size];
        for (int i = 0; i < rhs.size; ++i) {
            p[i] = rhs.p[i];
        }
        
        return *this;
    }
    Polynomial& operator=(Polynomial&& rhs) {
        size = rhs.size;
        p = rhs.p;
        rhs.p = NULL;
        
        return *this;
    }
    ~Polynomial() {
        if (p) {
            delete[] p;
        }
    }
    
    Type operator()(const Type& x) const {
        Type a = 1, ret = 0;
        for (int i = 0; i < size; ++i) {
            ret += p[i] * a;
            a *= x;
        }

        return ret;
    }
    
    Polynomial operator+(const Polynomial& rhs) const {
        if (size < rhs.size) {
            Polynomial temp(rhs.size);
            for (int i = 0; i < size; ++i) {
                temp.p[i] = p[i] + rhs.p[i];
            }
            for (int i = size; i < rhs.size; ++i) {
                temp.p[i] = rhs.p[i];
            }
            temp.resize();
            return temp;
        } else {
            Polynomial temp(size);
            for (int i = 0; i < rhs.size; ++i) {
                temp.p[i] = p[i] + rhs.p[i];
            }
            for (int i = rhs.size; i < size; ++i) {
                temp.p[i] = p[i];
            }
            temp.resize();
            return temp;
        }
    }
    Polynomial& operator+=(const Polynomial& rhs) {
        *this = *this + rhs;
        return *this;
    }
    
    Polynomial operator*(const Type& k) const {
        Polynomial temp(*this);
        temp *= k; return temp;
    }
    Polynomial& operator*=(const Type& k) {
        for (int i = 0; i < size; ++i) {
            p[i] *= k;
        }
        
        return *this;
    }
    Polynomial operator*(const Polynomial& rhs) const {
        Polynomial temp(size + rhs.size - 1);
        for (int i = 0; i < temp.size; ++i) {
            temp.p[i] = 0;
        }
        for (int i = 0; i < size; ++i) {
            for (int j = 0; j < rhs.size; ++j) {
                temp.p[i + j] += p[i] * rhs.p[j];
            }
        }
        
        return temp;
    }
    Polynomial& operator*=(const Polynomial& rhs) {
        *this = *this * rhs;
        return *this;
    }
    
    Polynomial operator/(const Type& k) const {
        Polynomial temp(*this);
        temp /= k; return temp;
    }
    Polynomial& operator/=(const Type& k) {
        for (int i = 0; i < size; ++i) {
            p[i] /= k;
        }
        
        return *this;
    }

    
    Polynomial synthetic_div(const Type& rhs) const {
        if (size == 1) {
            Polynomial ret(1);
            ret.p[0] = 0;
            return ret;
        }
        Polynomial temp(size - 1);
        temp.p[size - 2] = p[size - 1];
        for (int i = size - 2; i > 0; --i) {
            temp.p[i - 1] = p[i] + temp.p[i] * rhs;
        }
        
        return temp;
    }
    Polynomial& apply_synthetic_div(const Type& rhs) {
        *this = this->synthetic_div(rhs);
        return *this;
    }

    int degree() const {
        return size - 1;
    }
    
    const Polynomial& print() const {
        if (size == 0) {
            cout << "null polynomial" << endl;
        } else {
            cout << p[0];
            for (int i = 1; i < size; ++i) {
                cout << " + " << p[i] << 'x' << i;
            } cout << endl;
        }
        
        return *this;
    }
    
    static Polynomial get_linear(const Type& t) {
        Type arr[2];
        arr[0] = -t; arr[1] = 1;
        return Polynomial<Type>(arr, 1);
    }
    
private:
    Polynomial& resize() {
        while (p[size - 1] == 0 && size > 1) {
            --size;
        }
        
        return *this;
    }
};

template<class Type>
Polynomial<Type> products(Type* arr, int size) {
    Polynomial<Type> ret;
    for_each(arr, arr + size, [&ret](const Type& r){ret *= Polynomial<Type>::get_linear(r);});
    
    return ret;
}

template<class Type>
Type coeffs(Type* x, int size, int i) {
    Type ret = 1, missed = x[i];
    for_each(x, x + size, [&ret, missed](const Type& r){ret *= (r == missed ? 1 : missed - r);});
    
    return ret;
}

template<class Type>
Polynomial<Type> lag_interpol(Type* x, Type* y, int size) {
    const Polynomial<Type> product = products(x, size);
    
    Polynomial<Type> ret = product.synthetic_div(x[0]) * (y[0] / coeffs(x, size, 0));
    for (int i = 1; i < size; ++i) {
        ret += product.synthetic_div(x[i]) * (y[i] / coeffs(x, size, i));
    }
    
    return ret;
}