#include <iostream>
#include <conio.h>
#include <cmath>

using namespace std;

class Polynom
{
    int power;
    double* gaps = nullptr;
    int rootCount = 0;
public:
    double* coefs = nullptr;
    double* roots = nullptr;

    ~Polynom()
    {
        if (coefs != nullptr)
            delete[] coefs;
        if (gaps != nullptr)
            delete[] gaps;
        if (roots != nullptr)
            delete[] roots;
    }
    double getVal(double X)
    {
        double val = coefs[power];
        for (int i = 0; i < power; i++)
            val += coefs[i] * pow(X, power - i);

        return val;
    }
    void setCoefs(int power, double* coefs)
    {
        this->power = power;
        this->coefs = new double[power + 1];
        for (int i = 0; i <= power; i++)
            this->coefs[i] = coefs[i];
    }
    void setGaps(double* extrems)
    {
        gaps = new double[power + 1];
        double maxCoef = fabs(coefs[1]);
        for (int i = 1; i <= power; i++)
            if (fabs(coefs[i]) > maxCoef)
                maxCoef = fabs(coefs[i]);

        double R = 1 + maxCoef / fabs(coefs[0]);
        gaps[0] = -R;
        gaps[power] = R;
        for (int i = 1; i < power; i++)
            gaps[i] = extrems[i - 1];

        /*DEBUG
        cout << "GAPS" << endl;
        for (int i = 0; i < power + 1; i++)
            cout << gaps[i] << "   ";
        cout << endl;*/
    }
    void solve(double eps)
    {
        if (power == 1)
        {
            rootCount = 1;
            roots = new double[1];
            roots[0] = -coefs[1] / coefs[0];
        }
        else
        {

            for (int i = 0; i < power; i++)
            {
                if (getVal(gaps[i])*getVal(gaps[i + 1]) < 0)
                    rootCount++;
                else if (getVal(gaps[i]) == 0)
                    rootCount++;
            }
            roots = new double[rootCount];
            int curRootC = 0;
            for (int i = 0; i < power; i++)
            {
                if (getVal(gaps[i])*getVal(gaps[i + 1]) < 0)
                {
                    double a = gaps[i];
                    double b = gaps[i + 1];
                    double curX = a;
                    double prevX;

                    do
                    {
                        prevX = curX;
                        curX = a - getVal(a) * (b - a) / (getVal(b) - getVal(a));
                        if (getVal(a)*getVal(curX) < 0)
                        {
                            b = curX;
                        }
                        else
                        {
                            a = curX;
                        }
                    } while (fabs(curX - prevX) > eps);
                    roots[curRootC] = curX;
                    curRootC++;
                }
                else if (getVal(gaps[i]) == 0)
                {
                    roots[curRootC] = gaps[i];
                    curRootC++;
                }
            }
        }
        /*DEBUG
        cout << "ROOTS" << endl;
        for (int i = 0; i < rootCount; i++)
        {
            cout << roots[i] << "   ";
        }
        cout << endl; */
    }

    int getRootQt()
    {
        return rootCount;
    }

    double * getRoots()
    {
        return roots;
    }
    double * getCoefs()
    {
        return coefs;
    }
};

int main()
{
    int power;
    double* eqCoefs;
    cout << "Solving polynomical equtsion" << endl;
    cout << "Set power of polynom: ";
    cin >> power;
    eqCoefs = new double[power + 1];
    cout << "Input coefficients" << endl;
    for (int i = 0; i <= power; i++)
    {
        cin >> eqCoefs[i];
    }

    Polynom* dervs = new Polynom[power];
    dervs[0].setCoefs(power, eqCoefs);
    for (int i = 0; i < power-1; i++)
    {
        double *dervCoefs = new double[power - i];
        for (int j = 0; j < power - i; j++)
            dervCoefs[j] = dervs[i].coefs[j] * (power - i - j);

        dervs[i + 1].setCoefs(power - 1 - i, dervCoefs);
        delete[] dervCoefs;
    }

    cout << "Polynomical equation: ";
    for (int i = 0; i < power; i++)
        cout << dervs[0].coefs[i] << "*X^" << power - i << " + ";
    cout << dervs[0].coefs[power] << " = 0" << endl;

    double accur;
    cout << "Set accuracy: ";
    cin >> accur;

    for (int i = power - 1; i > 0; i--)
    {
        dervs[i].solve(accur);
        dervs[i - 1].setGaps(dervs[i].roots);
    }
    dervs[0].solve(accur);

    double* roots = dervs[0].getRoots();
    cout << "Roots: ";
    for (int i = 0; i < dervs[0].getRootQt(); i++)
        cout << roots[i] << "       ";
    cout << endl;

    system("pause");
    delete[] eqCoefs;
    delete[] dervs;
    return 0;
}