#include "polynomial.h"
#include <limits>

Term::Term(double a, int n) : a(a), n(n) {}

Complex::Complex(double a, double b) : a(a), b(b) {}
Complex::Complex(double a) : a(a), b(0) {}
Complex::Complex() : a(0), b(0) {}

Complex operator + (Complex A, Complex B){
    return Complex(A.a+B.a, A.b+B.b);
}
Complex operator - (Complex A, Complex B){
    return Complex(A.a-B.a, A.b-B.b);
}
Complex operator * (Complex A, Complex B){
    return Complex(A.a*B.a - A.b*B.b, A.a*B.b + A.b*B.a);
}
Complex operator / (Complex A, Complex B){
    return Complex( (A.a*B.a+A.b*B.b) / (B.a*B.a+B.b*B.b), (A.b*B.a-A.a*B.b) / (B.a*B.a+B.b*B.b) );
}
bool operator == (Complex A, Complex B){
    return abs(A.a-B.a) < std::numeric_limits<double>::epsilon()*256
        && abs(A.b-B.b) < std::numeric_limits<double>::epsilon()*256;
}

std::ostream& operator << (std::ostream& out, Complex c){
    out << "(" << c.a << ", " << c.b << "i)";
    return out;
}

Complex pow(Complex c, int n){
    if(n == 0) return 1;
    if(n == 1) return c;
    return pow(c, n/2) * pow(c, n/2 + n%2);
}

//f(x)
Complex Polynomial::eval(Complex x){
    Complex sum = 0;
    for(unsigned i = 0; i < p.size(); i++)
        sum = sum + p[i].a * pow(x, p[i].n);
    return sum;
}

//highest power of x
int Polynomial::degree(){
    int max = 0;
    for(unsigned i = 0; i < p.size(); i++)
        max = max > p[i].n ? max : p[i].n;
    return max;
}

//Durand–Kerner method
std::vector<Complex> Polynomial::solve(){
    std::vector<Complex> approx;

    Term max(0, 0);
    for(unsigned i = 0; i < p.size(); i++)
        if(i == 0 || max.n < p[i].n) max = p[i];

    for(unsigned i = 0; i < p.size(); i++)
        p[i].a /= max.a;

    for(unsigned i = 0; i < max.n; i++)
        approx.push_back( Complex( rand()%32768/65536.0, rand()%32768/65536.0 ) );

    while(1){
        std::vector<Complex> root(approx.size());

        bool end = true;

        for(unsigned i = 0; i < approx.size(); i++){
            root[i] = 1;
            for(unsigned j = 0; j < approx.size(); j++)
                if(i != j) root[i] = root[i] * (approx[i]-approx[j]);
            root[i] = approx[i] - eval(approx[i])/root[i];

            end &= root[i] == approx[i];
        }

        if(end) return root;
        else approx = root;
    }
}

std::ostream& operator << (std::ostream& out, Polynomial p){
    for(unsigned i = 0; i < p.p.size(); i++){
        if(i != 0 && p.p[i].a > 0) out << "+";
        out << p.p[i].a << "*x^" << p.p[i].n << " ";
    }
    return out;
}