#include <random>
#include <cmath>
#include "host_aberth.h"

using namespace std;


struct AberthState {
	// boundPoly is 's' in the stop condition formulated at p.189-190 of
	// https://link.springer.com/article/10.1007%2FBF02207694

	const Poly &poly;
	Poly deriv;
	Poly boundPoly;
	AApprox approx;
	double radius;

	AberthState(const Poly &poly);
	void regenerate();
	bool step();
	void iterate();
};

static thread_local minstd_rand randgenerator = minstd_rand(random_device()());

AberthState::AberthState(const Poly &poly)
        : poly(poly), deriv(derivative(poly)), radius(maxRootNorm(poly)) {

    regenerate();
    for (int i = 0; i <= N; i++) {
        boundPoly[i] = abs(poly[i]) * (4 * i + 1);
    }
}

void AberthState::regenerate() {
    auto genCoord = [this]() {
        return uniform_real_distribution<double>(-radius, radius)(randgenerator);
    };

    for (int i = 0; i < N; i++) {
        approx[i] = Com(genCoord(), genCoord());
    }
}

// Jacobi-style step where the new elements are computed in parallel from the previous values
bool AberthState::step() {
    array<Com, N * N> pairs;
    for (int i = 0; i < N - 1; i++) {
        for (int j = i + 1; j < N; j++) {
            pairs[N * i + j] = 1.0 / (approx[i] - approx[j]);
        }
    }

    bool allConverged = true;

    AApprox newapprox;
    AApprox offsets;
    for (int i = 0; i < N; i++) {
        Com pval = eval(poly, approx[i]);
        Com derivval = eval(deriv, poly.size() - 1, approx[i]);
        Com quo = pval / derivval;
        Com sum = 0;
        for (int j = 0; j < i; j++) sum -= pairs[N * j + i];
        for (int j = i + 1; j < N; j++) sum += pairs[N * i + j];
        offsets[i] = quo / (1.0 - quo * sum);

        newapprox[i] = approx[i] - offsets[i];

        double sval = eval(boundPoly, abs(newapprox[i]));
        if (abs(pval) > 1e-5 * sval) allConverged = false;
    }

    approx = newapprox;

    return allConverged;
}

void AberthState::iterate() {
    int tries = 1, stepIdx = 1;

    while (!step()) {
        stepIdx++;

        if (stepIdx > tries * 100) {
            regenerate();
            stepIdx = 0;
            tries++;
        }
    }
}

AApprox aberth(const Poly &poly) {
	AberthState state(poly);
	state.iterate();
	return state.approx;
}