#include <iostream>
#include <map>
#include <vector>
#include <cmath>
#include <fenv.h>
#include <stdio.h>


using namespace std;

double getf(double x, double t) {
    return -1;//3 * t * t - 6 * x;
}

double phi(double x) {
    return x*x;//x * x * x;
}

double alpha(double t) {
    return t;//t * t * t;
}

double betta(double t) {
    return 1+t;//1 + t * t * t;
}

double geta(double x, double t) {
    return 1;
}

double getb(double x, double t) {
    return 0;
}

double getc(double x, double t) {
    return 0;
}

int N = 100;
int M = 100;
const long double T = 0.1;
double h = 1. / (N + 0.);
double tau = T / (M + 0.);

vector < double > x;
vector < double > t;


double ret[1000][1000];


double L(int i, int k) {
    return geta(x[i], t[k]) * (ret[i + 1][k] - 2 * ret[i][k] + ret[i - 1][k]) / (h * h) +
           getb(x[i], t[k]) * (ret[i + 1][k] - ret[i - 1][k]) / (2 * h) +
           getc(x[i], t[k]) * ret[i][k];
}


struct Val{
    int n, m;
    map < pair < double, double>, double > a;
    Val(int n = 0, int m = 0): n(n), m(m) {
    }

    void print() {
        for (auto tmp : a) {
            cout << tmp.first.first << " " << tmp.first.second << " " << tmp.second << endl;
        }
    }

    double operator()(double x, double t) {
        pair < double, double> mP;
        mP.first = x;
        mP.second = t;
        return a[mP];
    }

    void set_(int i, int j, double val) {
        pair < double, double> mP;
        mP.first = x[i];
        mP.second = t[j];
        a[mP] = val;
    }

};

Val naive_dif_scheme() {
    x.clear();
    t.clear();
    h = 1. / (N + 0.);
    tau = (double)T / (M + 0.);
    if ((tau / (h * h)) > (1.0 / 2.0)) {
        cout << "WARNING --- Система не устойчива\n";
    }
    for (int i = 0; i <= N; i++) {
        x.push_back(i * h);
    }

    for (int k = 0; k <= M; k++) {
        t.push_back(tau * k);
    }

    for (int i = 0; i <= N; i++) {
        ret[i][0] = phi(x[i]);
    }
    for (int k = 1; k <= M; k++) {
        ret[0][k] = alpha(t[k]);
        for (int i = 1; i <= N - 1; i++) {
            ret[i][k] = ret[i][k - 1] + tau * (L(i, k - 1) + getf(x[i], t[k - 1]));
        }
        ret[N][k] = betta(t[k]);
    }
    Val retq(N, M);
    for (int i = 0; i <= N; i++)
        for (int j = 0; j <= M; j++)
            retq.set_(i, j, ret[i][j]);
    return retq;
}

template <typename Q, typename U>
double getNorm(Q a, U b, int n, int m) {
    double ret = 0;
    double h = 1./ (n + 0.);
    double tau = T / (m + 0.);
    for (int i = 0; i < n ; i++) {
        for (int j = 0; j < m ; j++) {
            ret = max(ret, fabs(a(i * h, j * tau) - b(i * h, j * tau)));
       //     if (ret == 0)
       //         cout << i << " " << j << " " << a(i * h, j * tau) << " " << b(i * h, j * tau) << endl;
        }
    }
    return ret;
};

struct neededAnswer {
    double operator()(double x, double t) {
        return x*x+t;//x * x * x + t * t * t;
    }
};



Val with_weight(double w) {
    x.clear();
    t.clear();
    h = 1. / (N + 0.);
    tau = (double)T / (M + 0.);
       if ((tau / (h * h)) > (1.0 / 2.0)) {
        cout << "WARNING --- Система не устойчива\n";
    }
    for (int i = 0; i <= N; i++) {
        x.push_back(i * h);
    }

    for (int k = 0; k <= M; k++) {
        t.push_back(tau * k);
    }

    for (int i = 0; i <= N; i++) {
        ret[i][0] = phi(x[i]);
    }
    double G[N+1];
    for (int k = 1; k <= M; k++) {
        G[0] = alpha(t[k]);
        G[N] = betta(t[k]);
        for (int i = 1; i <= N - 1; i++) {
            G[i] = -1 / tau * ret[i][k - 1] - (1 - w) * L(i, k - 1) - getf(x[i], t[k - 1]);
        }
        ret[0][k] = alpha(t[k]);
        double A[N + 1], B[N + 1], C[N + 1];

        A[0] = 0;
        B[0] = -1;
        C[0] = 0;
        C[N] = 0;
        A[N] = -1. / h;

        B[N] = A[N];
        for (int i = 1; i <= N - 1; i++) {
            A[i] = w * (1 / h / h - 1 / (2 * h));
            B[i] = w * 2 / h / h + 1. / tau;
            C[i] = w * (1 / h / h + 1 / (2 * h));
        }
        double T[N+1], S[N+1];
        S[0] = C[0] / B[0];
        T[0] = -G[0] / B[0];
        for (int i = 1; i <= N; ++i) {
            S[i] = C[i] / (B[i] - A[i] * S[i - 1]);
            T[i] = (A[i] * T[i - 1] - G[i]) / (B[i] - A[i] * S[i - 1]);

        }
        double Y[N+1];
        Y[N] = T[N];
        for (int i = N - 1; i >= 0; --i) {
            Y[i] = S[i] * Y[i + 1] + T[i];

        }
        for (int i = 1; i <= N - 1; i++) {
            ret[i][k] = Y[i];


        }
        ret[N][k] = betta(t[k]);
    }
    Val retq(N, M);
    for (int i = 0; i <= N; i++)
        for (int j = 0; j <= M; j++) {
            retq.set_(i, j, ret[i][j]);


        }
    return retq;
}

int main() {
    neededAnswer answer;
    for (int n = 5; n <= 40; n *= 2) {
         if (n == 40)
                n=30;
        N = n;
        M = n * n;
        Val tmp = naive_dif_scheme();
        cout <<  "N =  " << N << " M = " << M << " ";
        cout << " разница = " << getNorm(answer, tmp, N, M) << endl;
    }

    double sigma = 0.0;
    while (sigma < 1) {
        cout << "вес = " << sigma << endl;
        for (int n = 5  ; n <= 40; n *= 2) {
            if (n == 40)
                n=30;
            N = n;
            M = n * n;
            Val tmp = with_weight(sigma);
             cout <<  "N =  " << N << " M = " << M << " ";
             cout << " разница = " << getNorm(answer, tmp, N, M) << endl;
      }
        sigma += 0.1;
    }
    return 0;
 }