#include <iostream>
#include <Cmath>
#include <complex>
#include <vector>

using namespace std;

double getAbs(complex<double> x) {
    return sqrt(x.imag() * x.imag() + x.real() * x.real());
}

vector<complex<double>> derriativeCoefficients(vector<complex<double>> coef) {
    vector<complex<double>> temp;
    int N = coef.size() - 1;
    for (int i = 0; i < N; i++) {
        temp.push_back(coef[i] * (double) (N - i));
    }
    return temp;
}

vector<complex<double>> coefficients = {complex<double>(25, 0), complex<double>(16, 0), complex<double>(18, 0),
                                        complex<double>(1, 0), complex<double>(13, 0), complex<double>(0, 1)};
vector<complex<double>> coefficients1 = derriativeCoefficients(coefficients);

double A(vector<complex<double>> coefficients) {
    double max = getAbs(coefficients[0].real());
    for (int i = 1; i < coefficients.size(); i++) {
        if (getAbs(coefficients[i].real()) > max)
            max = getAbs(coefficients[i].real());
    }
    return max;
}

double an(vector<complex<double>> coefficients) {
    return getAbs(coefficients[0]);
}

vector<complex<double>> horners_method(complex<double> x, vector<complex<double>> coefficients) {
    vector<complex<double>> newcoefficients;
    newcoefficients.push_back(coefficients.operator[](0));
    for (int i = 1; i < coefficients.size(); i++) {
        newcoefficients.push_back(x * newcoefficients[i - 1] + coefficients[i]);
    }
    return newcoefficients;
}

double W(complex<double> z, vector<complex<double>> coefficients) {
    vector<complex<double>> temp = horners_method(z, coefficients);
    return getAbs(temp.operator[](temp.size() - 1));
}

complex<double> grad(complex<double> z, vector<complex<double>> coefficients) {
    complex<double> Pz = horners_method(z, coefficients).back();
    complex<double> Pzconj(Pz.real(), -Pz.imag());
    complex<double> Pz1 = horners_method(z, derriativeCoefficients(coefficients)).back();
    complex<double> gr = Pzconj * Pz1;
    gr = complex<double>(gr.real(), -gr.imag());
    return gr;
}

float failed_attempts = 0;
int N = 0;

complex<double>
coordinate_descent_method(complex<double> z0, double Epsilon, double alpha1, double alpha2,
                          vector<complex<double>> coefficients) {
    complex<double> z;
    double w0, w;
    bool check1 = false;
    bool check2 = false;
    complex<double> delta = complex<double>(Epsilon, Epsilon);
    double deltaW;
    while (getAbs(delta) >= Epsilon || deltaW >= Epsilon) {
        N++;
        if (failed_attempts >= 10) {
            failed_attempts = 0;
            cout << "Root is " << z0 << "  " << "W(z0) =  " << W(z0, coefficients) << '\t' << N << endl;
            return z0;
        }
        complex<double> z1 = z0;
        z = z0 + complex<double>(alpha1, 0);
        w = getAbs(W(z, coefficients));
        w0 = getAbs(W(z0, coefficients));
        if (w <= w0) {
            z0 = z0 + complex<double>(alpha1, 0);
        } else {
            z = z0 - complex<double>(alpha1, 0);
            w = getAbs(W(z, coefficients));
            N++;
            if (w <= w0) {
                z0 = z0 - complex<double>(alpha1, 0);
            } else { check1 = true; }
        }
        z = z0 + complex<double>(0, alpha2);
        w = getAbs(W(z, coefficients));
        if (w <= w0) {
            z0 = z0 + complex<double>(0, alpha2);
        } else {
            z = z0 - complex<double>(0, alpha2);
            w = getAbs(W(z, coefficients));

            if (w <= w0) {
                z0 = z0 - complex<double>(0, alpha2);
            } else { check2 = true; }
        }
        delta = z1 - z0;
        deltaW = getAbs(W(z1, coefficients) - W(z0, coefficients));
    }
    if (check1) {
        alpha1 /= 12.0;
        failed_attempts++;
        return coordinate_descent_method(z0, Epsilon, alpha1, alpha2, coefficients);
    } else if (check2) {
        alpha2 /= 12.0;
        failed_attempts++;
        return coordinate_descent_method(z0, Epsilon, alpha1, alpha2, coefficients);
    } else {
        cout << "Root is " << z0 << "  " << "W(z0) =  " << W(z0, coefficients) << '\t' << N << endl;
        N = 0;
        return z0;
    }
}


complex<double> gradient_method_with_crushing_step(complex<double> z, double Epsilon, double delta, double alpha,
                                                   vector<complex<double>> coefficients) {
    cout << "gradient_method_with_crushing_step" << endl;
    complex<double> zi, zi1;
    complex<double> gr = grad(z, coefficients);
    N++;
    double wi, wi1;
    zi1 = z;
    wi = W(zi1, coefficients);
    while (getAbs(gr) > delta) {
        zi = zi1;
        zi1 = zi - alpha * gr;
        wi1 = W(zi1, coefficients);
        gr = grad(zi1, coefficients);
        N++;
        if (wi1 - wi > -alpha * Epsilon * getAbs(gr)) {
            alpha /= 10.0;
        }
//        cout << wi1 << "  " << zi1 << '\t' << gr << endl;
    }
    cout << "Answer: " << wi1 << "  " << zi1 << '\t' << gr << '\t' << N << endl;
    N = 0;
    return zi1;
}

complex<double>
gradient_method_with_const_step(complex<double> z, double Epsilon, vector<complex<double>> coefficients) {
    cout << "gradient_method_with_const_step:" << endl;
    complex<double> zi, zi1;
    complex<double> gr = grad(z, coefficients);
    N++;
    double wi1;
    double x = 1.0 + A(coefficients) / an(coefficients);
    int L = (int) round(pow(horners_method(complex<double>(x, 0), coefficients1).back().real(), 2) +
                        horners_method(complex<double>(x, 0), coefficients).back().real() *
                        horners_method(complex<double>(x, 0), derriativeCoefficients(coefficients1)).back().real());
    cout << "L = " << L << endl;
    zi1 = z;
    double alpha = (1 - Epsilon) / (double) L;
    while (getAbs(gr) > Epsilon) {
        zi = zi1;
        zi1 = zi - alpha * gr;
        wi1 = W(zi1, coefficients);
        gr = grad(zi1, coefficients);
        N++;
        //cout << wi1 << "  " << zi1 << '\t' << gr << endl;
    }
    complex<double> minimizer, minimum;
    minimizer = zi1;
    minimum = W(minimizer, coefficients);
    cout << "Answer: minimizer " << minimizer << '\t' << "minimum " << minimum << 't' << N << endl;
    N = 0;
    return zi1;
}

complex<double>
gradient_method_with_predictable_step(complex<double> z, double Epsilon, vector<complex<double>> coefficients) {
    cout << "gradient_method_with_predictable_step" << endl;
    complex<double> zi;
    complex<double> gr = grad(z, coefficients);
    N++;
    int k = 10000;
    zi = z;
    while (getAbs(gr) > Epsilon) {
        zi = zi - (1 / (double) (k)) * gr;
        gr = grad(zi, coefficients);
        N++;
//        cout << W(zi, coefficients) << "  " << zi << '\t' << gr << endl;
        k = k + 1;
    }
    complex<double> minimizer, minimum;
    minimizer = zi;
    minimum = W(minimizer, coefficients);
    cout << "Answer: minimizer " << minimizer << '\t' << "minimum " << minimum << '\t' << N << endl;
    N = 0;
    return zi;
}

vector<complex<double>> root_search(vector<complex<double>> coefficients) {
    cout << "root_search by coordinate_descent_method" << endl;
    vector<complex<double>> roots;
    int degree = coefficients.size();
    vector<complex<double>> tempcoeff = coefficients;
    complex<double> root;
    for (int i = 1; i < degree; i++) {
        cout << "W(z0)" << "  " << "z0" << endl;
        double alpha = 10e-5;
        double alpha1 = 0.5;
        double alpha2 = 0.5;
        double Epsilon = 10e-4;
        double delta = 10e-3;
        //root = coordinate_descent_method(complex<double>(0.0, 0.0), Epsilon, alpha1, alpha2, tempcoeff);
        //root = coordinate_descent_method(root, Epsilon, alpha1, alpha2, coefficients);
        root = gradient_method_with_crushing_step(complex<double>(0, 0), Epsilon, delta, alpha, tempcoeff);
        root = gradient_method_with_crushing_step(root, Epsilon, delta, alpha, coefficients);
//        root = gradient_method_with_const_step(complex<double>(0, 0), Epsilon, tempcoeff);
//        root = gradient_method_with_const_step(root, Epsilon, coefficients);
//        root = gradient_method_with_predictable_step(complex<double>(0, 0), Epsilon, tempcoeff);
//        root = gradient_method_with_predictable_step(root, Epsilon, coefficients);
        tempcoeff = horners_method(root, tempcoeff);
        tempcoeff.resize(degree - i);
        cout << endl;
        roots.push_back(root);
    }
    for (int i = 0; i < roots.size(); i++) {
        cout << "Root is " << roots[i] << "  " << "W(root) =  " << W(roots[i], coefficients) << endl;
    }
    return roots;
}

int main() {

    double delta = 10e-10;
    double alpha = 10e-5;
    double alpha1 = 5;
    double alpha2 = 5;
    double Epsilon = 10e-5;
//    coordinate_descent_method(complex<double>(0, 0), Epsilon, alpha1, alpha2, coefficients);
    vector<complex<double>> roots = root_search(coefficients);
//    gradient_method_with_crushing_step(complex<double>(0, 0), Epsilon, delta, alpha, coefficients);
    //gradient_method_with_const_step(complex<double>(0, 0), Epsilon, coefficients);
//    gradient_method_with_predictable_step(complex<double>(0, 0), Epsilon, coefficients);
    return 0;
}