// ConsoleApplication1.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "iostream"
#include "cmath"


using namespace std;
void xI(double nach, double kon, double *a, double n)
{
    a[0] = nach;
    for (int i = 1; i < n; ++i)
        a[i] = a[i - 1] + ((kon - nach) / (n - 1));

}

void pxI(double *p, double *a, double n)
{

    for (int i = 0; i < n; ++i)
        p[i] = -a[i]*a[i] - 1;

}

void qxI(double *q, double *a, double n)
{

    for (int i = 0; i < n; ++i)
        q[i] = a[i] *4;

}

void fxi(double *a, double *f, double m, double n)
{
    for (int i = 0; i < n; ++i)
        f[i] = 4 - 4 * a[i] * a[i] * a[i] + 4 * m*a[i];
}


void matrix(double *y0, double *y1, double *y2, double *y3, double *y4, double *q, double *p, double *f, double *svch, double n, double h, double o, double l, double *ksi, double *eta)
{
    y0[1] = 0;
    for (int i = 2; i < n - 1; ++i)
        y0[i] = 1 - (p[i] * h) / 2;
    //////////////////
    for (int i = 1; i < n - 1; ++i)
        y1[i] = q[i] * h*h - 2;
    /////////////////
    for (int i = 1; i < n - 2; ++i)
        y2[i] = 1 + (p[i] * h) / 2;
    y2[3] = 0;
    /////////////////
    svch[1] = h*h*f[1] - (1 - (h*p[1]) / 2)*o;
    for (int i = 2; i < n - 2; ++i)
        svch[i] = h*h*f[i];
    svch[3] = svch[3] = h*h*f[3] - (1 + (h*p[3]) / 2)*l;

        ////////////////
    ksi[1] = -y2[1] / y1[1];
    for (int i = 2; i < n - 1; ++i)
        ksi[i] = -y2[i] / (y0[i] * ksi[i - 1] + y1[i]);

    //////////////////////////////
    eta[1] = svch[1] / y1[1];
    for (int i = 2; i < n - 1; ++i)
        eta[i] = (svch[i] - y0[i] * eta[i - 1]) / (y0[i] * ksi[i - 1] + y1[i]);
}


//дописать нахождение обратной матрицы
void obrmatr(double *yk, double *ytoch, double *dy, double *a, double n, double m, double o, double l, double *eta, double *ksi)
{

    yk[0] = o;
    yk[3] = eta[3];
    yk[4] = l;
    for (int i = 2; i > 0; --i)
        yk[i] = ksi[i] * yk[i + 1] + eta[i];


    for (int i = 0; i < n; ++i)
        ytoch[i] = (a[i])*(a[i])*(a[i])*(a[i]) - 4 * a[i] + m;


    for(int i=0;i<n;++i)
    dy[i] = abs(ytoch[i] - yk[i]);

}



void matrix2(double *y0, double *y1, double *y2, double *q, double *p, double *f, double *svch, double n, double h, double o, double l, double *ksi, double *eta )
{

    y0[1] = 0;
    for (int i = 2; i < n-1; ++i)
        y0[i] = 1 - (p[i] * h) / 2;
    //////////////////
    for (int i = 1; i < n-1; ++i)
        y1[i] = q[i] * h*h - 2;
    /////////////////
    for (int i = 1; i < n-2; ++i)
    y2[i] = 1 + (p[i] * h) / 2;
    y2[7] = 0;
    /////////////////
    svch[1] = h*h*f[1] - (1 - (h*p[1]) / 2)*o;
    for (int i = 2; i < n - 2; ++i)
        svch[i] = h*h*f[i];
    svch[7]= svch[7] = h*h*f[7] - (1 + (h*p[7]) / 2)*l;
    /////////////////
    ksi[1] = -y2[1] / y1[1];
    for (int i = 2; i < n-1; ++i)
        ksi[i] = -y2[i]/(y0[i]*ksi[i-1]+y1[i]);

    //////////////////////////////
    eta[1] = svch[1]/y1[1];
    for (int i = 2; i < n-1; ++i)
        eta[i] = (svch[i]-y0[i]*eta[i-1])/(y0[i]*ksi[i-1]+y1[i]);
}



void obrmatr2(double *yk, double *ytoch, double *dy, double *a, double n, double m, double o, double l, double *eta, double *ksi)
{

    yk[0] = o;
    yk[7] = eta[7];
    yk[8] = l;
    for (int i = n-3; i > 0; --i)
    {
        yk[i] = ksi[i] * yk[i + 1] + eta[i];
    }

    ytoch[0] = o;
    ytoch[8] = l;
    for (int i = 1; i < n-1; ++i)
        ytoch[i] = (a[i])*(a[i])*(a[i])*(a[i]) - 4 * a[i] + m;


    for (int i = 0; i<n; ++i)
        dy[i] = abs(ytoch[i] - yk[i]);

}

void Dyi(double *e, double *d, double *b, double n, double *ytoch)
{
    d[0] = b[0];
    for (int h = 1, i = 2; i < n; i += 2)
    {
        d[i] = b[i - h];
        h++;
    }

    for (int i = 1; i < n; i += 2)
        d[i] = (d[i - 1] + d[i + 1]) / 2;


    for (int i = 0; i < n; ++i)

        e[i] = fabs(d[i] - ytoch[i]);

}


void Dyk(double *g, double *yk, double *ytoch, double n)
{
    for (int i = 0; i < n; ++i)
        g[i] = fabs(yk[i] - ytoch[i]);
}

void Rk(double *r, double *d, double *yk, double n)
{

    for (int i = 1; i < n; i+=2)
        r[i] = fabs(yk[i] - d[i]);


}

double max;

int main()
{
    setlocale(LC_ALL, "Russian");

    //Объявление массивов под значения
    double* a;
    double* b;
    double* d;
    double* e;
    double* g;
    double* r;
    double* p;
    double* q;
    double* f;
    double* y0;
    double* y1;
    double* y2;
    double* y3;
    double* y4;
    double* svch;
    double* ytoch;
    double* dy;
    double* ksi;
    double* eta;
    double* yk;



    //Объявляем количество на которое делим отрезок от а до b
    int n;

    //Объявление под переменные а b h
    double nach, kon, o, l, m;
    double h;
    double Eps;



    bool flag = false;

    cout << "Введите n в формате (n+1)" << endl;
    cout << "n= ";
    cin >> n;

    cout << "Введите а" << endl;
    cout << "а= ";
    cin >> nach;

    cout << "Введите b" << endl;
    cout << "b= ";
    cin >> kon;

    cout << "Введите M" << endl;
    cout << "M= ";
    cin >> m;

    cout << "Введите Y_-2" << endl;
    cout << "Y_-2= ";
    cin >> o;

    cout << "Введите Y_-1" << endl;
    cout << "Y_-1= ";
    cin >> l;

    cout << "Введите точность" << endl;
    cout << "Eps= ";
    cin >> Eps;



    do
    {

    //Выделение памяти под массивы
    a = new double[n];
    b = new double[n];
    p= new double[n];
    q = new double[n];
    f = new double[n];
    y0 = new double[n];
    y1 = new double[n];
    y2 = new double[n];
    y3 = new double[n];
    y4 = new double[n];
    svch = new double[n];
    yk = new double[n];
    ytoch = new double[n];
    dy = new double[n];
    eta = new double[n];
    ksi  = new double[n];



    h = (kon - nach) / (n-1);



    xI(nach, kon, a, n);
    pxI(p, a, n);
    qxI(q, a, n);
    fxi(a, f, m, n);
    matrix(y0, y1, y2, y3, y4, q, p, f, svch, n, h, o, l, ksi, eta);
    obrmatr(yk, ytoch, dy, a, n, m, o, l, eta, ksi);
    /*

    for (int i = 0; i < n; ++i)
        cout << "xi 0.25[" << i << "]= " << a[i] << endl;


    for (int i = 0; i < n; ++i)
        cout << "pxi 0.25[" << i << "]= " << p[i] << endl;


    for (int i = 0; i < n; ++i)
        cout << "qxi 0.25[" << i << "]= " << q[i] << endl;


for (int i = 0; i < n; ++i)
        cout << "fxi 0.25[" << i << "]= " << f[i] << endl;

    for (int i = 0; i < n; ++i)
        cout << "y0 0.25[" << i << "]= " << y0[i] << endl;

for (int i = 1; i < n-1; ++i)
        cout << "y1 0.25[" << i << "]= " << y1[i] << endl;

    for (int i = 1; i < n-1; ++i)
        cout << "y2 0.25[" << i << "]= " << y2[i] << endl;


    for (int i = 0; i < n; ++i)
        cout << "свободный член 0.25[" << i << "]= " << svch[i] << endl;

    for (int i = 1; i < n-1; ++i)
        cout << "ksi 0.25[" << i << "]= " << ksi[i] << endl;

    for (int i = 1; i < n-1; ++i)
        cout << "eta 0.25[" << i << "]= " << eta[i] << endl;

    for (int i = 0; i < n; ++i)
        cout << "yi 0.25[" << i << "]= " << yk[i] << endl;

    for (int i = 0; i < n; ++i)
        cout << "y точное 0.25[" << i << "]= " << ytoch[i] << endl;

    for (int i = 0; i < n; ++i)
        cout << "погрешность 0.25[" << i << "]= " << dy[i] << endl;
        */

    for (int i = 0; i < n; ++i)
        b[i] = yk[i];


    n = n * 2 - 1;
    h = (kon - nach) / (n-1);



    a = new double[n];
    p = new double[n];
    q = new double[n];
    f = new double[n];
    yk = new double[n];
    d = new double[n];
    e = new double[n];
    g = new double[n];
    r = new double[n];
    f = new double[n];
    y0 = new double[n];
    y1 = new double[n];
    y2 = new double[n];
    svch = new double[n];
    eta = new double[n];
    ksi = new double[n];
    ytoch = new double[n];
    dy = new double[n];


    xI(nach, kon, a, n);
    pxI(p, a, n);
    qxI(q, a, n);
    fxi(a, f, m, n);
    matrix2(y0, y1, y2, q, p, f, svch, n, h, o, l, ksi, eta);
    obrmatr2(yk, ytoch, dy, a, n, m, o, l, eta, ksi);

    Dyi(e, d, b, n, ytoch);
    Dyk(g, yk, ytoch, n);
    Rk(r, d, yk, n);



    max = 0;

    for (int i = 1; i < n; i+=2)
    {
        if (r[i] > max)
            max = r[i];
    }

    for (int i = 1; i < n; i+=2)
    {
        if (max < Eps)
        {
            flag = true;

            cout << "otvet=" << "r[" << i << "]=" << r[i] << endl;

            //cout << setiosflags(ios::fixed) << setprecision(3) << r[i] << endl;

        }
    }
    if (!flag)
        n = n * 2 - 1;

    } while (!flag);

    /*

    for (int i = 0; i < n; ++i)
        cout << "xi 0.125[" << i << "]= " << a[i] << endl;


    for (int i = 0; i < n; ++i)
        cout << "pxi 0.125[" << i << "]= " << p[i] << endl;


    for (int i = 0; i < n; ++i)
        cout << "qxi 0.125[" << i << "]= " << q[i] << endl;


    for (int i = 0; i < n; ++i)
        cout << "fxi 0.125[" << i << "]= " << f[i] << endl;

    for (int i = 1; i < n-1; ++i)
        cout << "y0 0.25[" << i << "]= " << y0[i] << endl;

    for (int i = 1; i < n - 1; ++i)
        cout << "y1 0.25[" << i << "]= " << y1[i] << endl;

    for (int i = 1; i < n - 1; ++i)
        cout << "y2 0.25[" << i << "]= " << y2[i] << endl;

    for (int i = 1; i < n - 1; ++i)
        cout << "свободный член 0.125[" << i << "]= " << svch[i] << endl;

    for (int i = 1; i < n - 1; ++i)
        cout << "ksi 0.125[" << i << "]= " << ksi[i] << endl;

    for (int i = 1; i < n - 1; ++i)
        cout << "eta 0.125[" << i << "]= " << eta[i] << endl;

    for (int i = 0; i < n; ++i)
        cout << "yk 0.125[" << i << "]= " << yk[i] << endl;

    for (int i = 0; i < n; ++i)
        cout << "y точное 0.125[" << i << "]= " << ytoch[i] << endl;

    for (int i = 0; i < n; ++i)
        cout << "погрешность 0.125[" << i << "]= " << dy[i] << endl;




    for (int i = 0; i < n; ++i)
        cout << "Dyi[" << i << "]= " << e[i] << endl;

    for (int i = 0; i < n; ++i)
        cout << "Dyk[" << i << "]= " << g[i] << endl;
*/
    for (int i = 1; i < n; i += 2)
        cout << "Rung[" << i << "]= " << r[i] << endl;

    delete[] a;
    delete[] b;
    delete[] d;
    delete[] e;
    delete[] g;
    delete[] r;
    delete[] p;
    delete[] q;
    delete[] f;
    delete[] y0;
    delete[] y1;
    delete[] y2;
    delete[] y3;
    delete[] y4;
    delete[] svch;
    delete[] ytoch;
    delete[] dy;
    delete[] ksi;
    delete[] eta;
    delete[] yk;

    system("pause");
}