anonymous

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

inline void skip_ws(const char*& str){
    while(*str == ' ') str++;
}

//== COMPLEX ==//

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;
}


//== ELEMENT ==//

Element::Element() : b(0), ref(0) {}

Element::Element(Base* b) : b(b), ref(new int(1)) {}

Element::Element(const Element& e){
    b = e.b;
    ref = e.ref;
    if(ref) (*ref)++;
}

Element::~Element(){
    rem();
}

Element& Element::operator = (const Element& e){
    if(&e == this) return *this;
    rem();
    b = e.b;
    ref = e.ref;
    if(ref) (*ref)++;
    return *this;
}

void Element::rem(){
    if(!ref) return;
    (*ref)--;
    if(*ref == 0){
        delete b;
        delete ref;
    }
}

std::ostream& operator << (std::ostream& out, Element e){
    if(e.b) e.b->print(out);
    else out << " [null] ";
    return out;
}

//== NUMBER ==//
Number::Number(int i) : num(i), den(1) {}

Number::Number() : num(0), den(1) {}

Number::Number(int n, int d){
    int g = gcd(n, d);
    num = n/g;
    den = d/g;
}

Number::Number(const char*& str){
    num = 0;
    den = 1;

    bool dot = 0;
    while(*str >= '0' && *str <= '9' || *str == '.'){
        if(*str == '.') dot = 1;
        else{
            num = num*10 + *str - '0';
            if(dot) den *= 10;
        }
        str++;
    }
}

int Number::gcd(int a, int b){
    if(a != 0 && b != 0) while(a % b != 0){
        a = a % b;
        std::swap(a, b);
    }
    return b;
}

Element Number::simplify(){
    return new Number(num, den);
}

Base::Type Number::type() { return Base::number; }

void Number::print(std::ostream& out){
    out<< " [" << num << "/" << den << "] ";
}

Complex Number::eval(Complex x){
    return double(num)/den;
}
std::vector<Complex> Number::solve(){
    return std::vector<Complex>();
}
Element Number::coefficient(){
    return new Number(num, den);
}
unsigned Number::degree(){
    return 0;
}
//== X ==//
X::X() {}

X::X(const char*& str){
    str++;
}

Element X::simplify(){
    return new X();
}

Base::Type X::type() { return Base::x; }

void X::print(std::ostream& out){
    out << " X ";
}

Complex X::eval(Complex x){
    return x;
}
std::vector<Complex> X::solve(){
    std::vector<Complex> root;
    root.push_back(0);
    return root;
}
Element X::coefficient(){
    return new Number(1);
}
unsigned X::degree(){
    return 1;
}
//== MONO ==//
Mono::Mono() {}
Mono::Mono(const char*& str){
    bool mul = 1;
    if(*str == '-') num.push_back(new Number(-1));
    else if(*str != '+') str--;
    str++;

    skip_ws(str);

    while(1){
        if(*str >= '0' && *str <= '9'){
            if(mul) num.push_back(new Number(str));
            else den.push_back(new Number(str));
        }
        else if(*str == '('){
            str++;
            if(mul) num.push_back(new Poly(str));
            else den.push_back(new Poly(str));
        }
        else if(*str == 'x'){
            if(mul) num.push_back(new X(str));
            else den.push_back(new X(str));
        }
        else if(*str == '/'){
            mul = 0;
            str++;
        }
        else if(*str == '*'){
            mul = 1;
            str++;
        }
        else break;
        skip_ws(str);
    }
}

Element Mono::simplify(){
    Element e;
    for(unsigned i = 0; i < num.size(); i++)
        if(e.b == 0) e = num[i].b->simplify();
        else e = e * num[i];

    if(e.b == 0) e = new Number(1);
    for(unsigned i = 0; i < den.size(); i++)
        e = e / den[i];

    return e;
}
Base::Type Mono::type() { return Base::mono; }

void Mono::print(std::ostream& out){
    out << " {";
    for(unsigned i = 0; i < num.size(); i++){
        if(i != 0) out << "*";
        num[i].b->print(out);
    }
    for(unsigned i = 0; i < den.size(); i++){
        out << "/";
        den[i].b->print(out);
    }
    out << "} ";
}

Complex Mono::eval(Complex x){
    Complex n = 1, d = 1;
    for(unsigned i = 0; i < num.size(); i++) n = n * num[i].b->eval(x);
    for(unsigned i = 0; i < den.size(); i++) d = d * den[i].b->eval(x);
    return n/d;
}

Element Mono::coefficient(){
    for(unsigned i = 0; i < num.size(); i++) if(num[i].b->type() == Base::number) return num[i];
    return new Number(1);
}
unsigned Mono::degree(){
    unsigned deg = 0;
    for(unsigned i = 0; i < num.size(); i++) if(num[i].b->type() == Base::x) deg++;
    return deg;
}

std::vector<Complex> Mono::solve(){
    std::vector<Complex> root;

    if(den.size()){

        Mono mn, md;
        mn.num = num;
        md.num = den;

        std::vector<Complex> nroot = mn.simplify().b->solve();
        for(unsigned i = 0; i < nroot.size(); i++)
            if(! (md.eval(nroot[i]) == 0) ) root.push_back(nroot[i]);

    }else root.push_back(0);

    return root;
}

//== POLY ==//
Poly::Poly() {}
Poly::Poly(const char*& str){
    while(*str == ' ') str++;

    bool sgn = 1;
    while(*str != 0 && *str != ')' && *str != '='){
        if(!sgn && *str != '+' && *str != '-'); //throw MADNESS;
        sgn = 0;
        num.push_back(new Mono(str));
        while(*str == ' ') str++;
    }
    if(*str == ')') str++;
}
Element Poly::simplify() {
    Element e;
    for(unsigned i = 0; i < num.size(); i++)
        if(e.b == 0) e = num[i].b->simplify();
        else e = e + num[i];

    if(e.b == 0) e = new Number(0);
    return e;
}
Base::Type Poly::type() { return Base::poly; }

void Poly::print(std::ostream& out){
    out << " (";
    for(unsigned i = 0; i < num.size(); i++){
        if(i != 0) out << "+";
        num[i].b->print(out);
    }
    out << ") ";
}

Complex Poly::eval(Complex x){
    Complex c = 0;
    for(unsigned i = 0; i < num.size(); i++) c = c + num[i].b->eval(x);
    return c;
}
std::vector<Complex> Poly::solve(){
    unsigned deg=0;
    Element coef;

    for(unsigned i = 0; i < num.size(); i++){
        unsigned tdeg = num[i].b->degree();
        if(tdeg > deg){
            deg = tdeg;
            coef = num[i].b->coefficient();
        }
    }
    for(unsigned i = 0; i < num.size(); i++) num[i] = num[i]/coef;

    std::vector<Complex> approx(deg);

    for(unsigned i = 0; i < deg; i++)
        approx[i] = Complex( rand()%32768/65536.0, rand()%32768/65536.0 );

    while(1){
        std::vector<Complex> root(deg);

        bool end = true;

        for(unsigned i = 0; i < deg; 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 &= eval(root[i]) == 0;
        }

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

//== EQUATION ==//

Equation::Equation(const char *str){
    e = new Poly(str);

    while(*str == ' ') str++;
    if(*str == '='){
        str++;
        e = e + Element(new Number(-1))*Element(new Poly(str));
    }
}

std::vector<Complex> Equation::solve(){
    e = e.b->simplify();
    std::cout << e << '\n';
    return e.b->solve();
}

//these two will never be called
Element Poly::coefficient(){
    return Element();
}
unsigned Poly::degree(){
    return 0;
}