#include <iostream>
#include <cmath>
#include <list>
using namespace std;
float eps = 0.0001;

float fun(float x, float y, float dy, float d2y);
float* runge_kutt(float x, float* y, float h);
float* double_recalculation(float a, float* y, float h);
float mt_shoot(float* y, float h);
float bisect(float x, float* y);
float* get_interpol(float* x, float* y);
float Aitken(float* X, float* y, float x);
float Simpson(float a, float b, float n, float* interpol);

int main()
{
    float* interpol = new float[6];
    float* x_n = new float[6];
    float** answer = new float* [6];
    float* Koshi = new float[2]{ 1, 6 };
    float x = 0.0, h = 0.2, dy = mt_shoot(Koshi, h);
    float* y = new float[2]{ Koshi[0], dy };
    for (int i = 0; i < 6; i++) answer[i] = new float[3];

    for (int i = 0; i < 6; i++)
        x_n[i] = i * h;

    int i = 0;
    while (x <= 1)
    {
        answer[i][0] = x;
        answer[i][1] = y[0];
        answer[i][2] = y[1];
        interpol[i] = y[0];
        i++; x += h;
        y = runge_kutt(x, y, h);
    }
    answer[i][0] = x; answer[i][1] = Koshi[1]; answer[i][2] = y[1];
    interpol[i] = Koshi[1];
    float* interpol_res = get_interpol(x_n, interpol);
    cout << "Integral: " << Simpson(0, 1, 1000, interpol_res) << endl;
    for (i = 0; i < 6; i++)
    {
        for (int j = 0; j < 2; j++)
            cout << answer[i][j];
        cout << answer[i][2];
        cout << endl;
    }

    return 0;
}

float fun(float x, float y, float dy, float d2y)
{
    return pow(d2y, 5) - 7 * cos(x) * d2y - 10 * x - 4 * exp(x) * dy + (5 * y) / (x + 6);
}

float f(float x, float* y)
{
    float z = 0; int i;
    for (i = 0; z < x; i++) z += 0.01;

    return y[i];
}

float bisect(float x, float* y)
{
    float a = 0, b = 1, fun_a = 1, fun_b = 1, fun_differ, differ;

    while (fun_a * fun_b > 0)
    {
        a--; b++;
        fun_a = fun(x, y[0], y[1], a);
        fun_b = fun(x, y[0], y[1], b);
    }

    while (true)
    {
        fun_a = fun(x, y[0], y[1], a);
        differ = (a + b) / 2;
        fun_differ = fun(x, y[0], y[1], differ);

        if (fabs(fun_differ) < eps)
            break;

        if (fun_a * fun_differ < 0)
            b = differ;
        else a = differ;

    }

    return differ;
}

float* get_interpol(float* x, float* y)
{
    float* res = new float[102];
    int k = 0;
    for (float i = 0; i <= 1.01; i += 0.01)
    {
        res[k] = (Aitken(x, y, i));
        k++;
    }

    return res;
}

float Aitken(float* X, float* y, float x)
{
    float* P = new float[sizeof(y)];
    for (int i = 0; i < sizeof(y); i++) P[i] = y[i];

    for (int j = 1; j < sizeof(y); j++)
        for (int i = 0; i < sizeof(y) - j; i++)
            P[i] = (P[i] * (X[i + j] - x) - P[i + 1] * (X[i] - x)) / (X[i + j] - X[i]);

    return P[0];
}

float Simpson(float a, float b, float n, float* interpol)
{
    int k;
    float h = (b - a) / n;
    float answer = f(a, interpol) + f(b, interpol);

    for (int i = 1; i < n; i++)
    {
        if (i % 2 == 0) k = 2;
        else k = 4;
        answer += k * f(a + i * h, interpol);
    }
    answer *= h / 3;

    return answer;
}

float* runge_kutt(float x, float* y, float h)
{
    float* y_current = new float[2];
    y_current[0] = y[0]; y_current[1] = y[1];

    for (int i = 0; i < 2; i++)
    {
        y_current[0] = y_current[0] + (h / 2) * y_current[1];
        y_current[1] = y_current[1] + (h / 2) * bisect(x + h, y_current);
    }

    return y_current;
}

float* double_recalculation(float a, float* y, float h)
{
    float* y_1 = new float[2];
    float* y_2 = new float[2];
    float x_i = 0;
    for (int i = 0; i < 2; i++)
    {
        y_1[0] = y_2[0] = 1;
        y_1[1] = y_2[1] = a;
    }

    while (x_i <= 1)
    {
        y_1 = runge_kutt(x_i, y_1, h);
        y_2 = runge_kutt(x_i, y_2, h / 2);
        y_2 = runge_kutt(x_i + h / 2, y_2, h / 2);

        if (2 * eps <= fabs(y_1[0] - y_2[0]))
        {
            x_i = 0;
            h /= 2;
            y_1 = new float[2]{ 1, a };
            y_2 = new float[2]{ 1, a };
        }
        else x_i += h;
    }

    y[0] = y_2[0]; y[1] = y_2[1];

    return y;
}

float mt_shoot(float* y, float h)
{
    float a = 0, b = 1;
    float* y_a = new float[2];
    float* y_b = new float[2];
    float* y_c = new float[2];
    y_a = double_recalculation(a, y_a, h);

    bool f = y_a[1] < y[1] ? true : false;
    while (f == (y_a[1] < y[1]))
    {
        if (y_a[1] < y[1])
        {
            f = true;
            a = b;
            b += h;
        }
        else
        {
            f = false;
            b = a;
            a -= h;
        }
        y_a = double_recalculation(a, y_a, h);
    }

    b = a; a -= h;
    y_a = double_recalculation(a, y_a, h);
    y_b = double_recalculation(b, y_b, h);

    float c = 0;
    while (fabs(y[1] - y_c[1]) > 0.2)
    {
        c = (a + b) / 2;
        y_c = double_recalculation(c, y_c, h);
        if ((y_a[1] - y[1]) * (y_b[1] - y[1]) < 0)
        {
            b = c;
            y_b = y_c;
        }
        else
        {
            a = c;
            y_a = y_c;
        }
    }

    return c;
}