#define _CRT_SECURE_NO_WARNINGS
#define _USE_MATH_DEFINES

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include <chrono>
#include <ctime>

#include "graphic.h"

#define A 0.0
#define B 1.0

#define STARTPA 1.0
#define STARTPB 1.0

#define P 1.0

#define NSTART 5

struct config {
    int n;
    int nseg;
    double h;
    double* x;
    double* u;
};

/*double c2(double x) {
    double denum = P * (exp(P) + exp(-P)) - STARTPB * (exp(-P) - exp(P));
    double mid = (2 * exp(P) / P + 1 + 6.0 / pow(P, 2)) / pow(P, 2);

    double bpart = STARTPB * (exp(P) * (STARTPA + 2.0 / pow(P, 4)) + 4.0 / pow(P, 4));
    double apart = STARTPA * P * exp(P);

    return (mid + apart + bpart) / denum;
}
double c1(double x) {
    return STARTPA + 2.0 / pow(P, 4) - c2(x);
}*/
double u(double x) {
    double mul1 = STARTPA * pow(P, 4) + 2;
    double mul2 = P - STARTPB;
    double mul3 = P + STARTPB;
    double mul4 = pow(P, 2) * pow(x, 3) - pow(P, 2) * pow(x, 2) + 6 * x - 2;
    double mul5 = 4 * STARTPB + pow(P, 2) + 6;

    double add1 = mul1 * mul2 * exp(2*P*x);
    double add2 = mul1 * mul3 * exp(2*P);
    double add3 = mul2 * mul4 * exp(P*x);
    double add4 = mul3 * mul4 * exp(P*(x + 2));
    double add5 = mul5 * exp(2*P*x + P);
    double add6 = mul5 * exp(P);

    double denum = pow(P, 4) * (STARTPB * (exp(2*P) - 1) + P * (exp(2*P) + 1));

    double sol = exp(-P*x) * (add1 + add2 + add3 + add4 - add5 + add6) / denum;
    return sol;
    //return c1(x) * exp(P * x) + c2(x) * exp(-P * x) + (pow(x, 3) - pow(x, 2)) / pow(P, 2) + (6 * x - 2) / pow(P, 4);
}

double f(double x) {
    return pow(x, 2) - pow(x, 3);
}
void solve(int n, double* a, double* c, double* b, double* f, double* x)
{
    double m;
    for (int i = 1; i < n; i++)
    {
        m = a[i] / c[i - 1];
        c[i] = c[i] - m * b[i - 1];
        f[i] = f[i] - m * f[i - 1];
    }

    x[n - 1] = f[n - 1] / c[n - 1];

    for (int i = n - 2; i >= 0; i--)
    {
        x[i] = (f[i] - b[i] * x[i + 1]) / c[i];
    }
}

config* create(int n) {
    config* con = new config();

    con->n = n;
    con->nseg = n - 1;
    con->h = (B - A) / con->nseg;

    con->x = new double[n];
    con->x[0] = A;

    for (int i = 1; i < n; i++) {
        con->x[i] = con->x[i - 1] + con->h;
    }

    return con;
}

void methodTridig(config* con) {
    int n = con->n;

    double* a = new double[n];
    double* b = new double[n];
    double* c = new double[n];
    double* d = new double[n];

    for (int i = 0; i < n; i++) {
        a[i] = 1;
        b[i] = 1;
        c[i] = -2.0 - pow(con->h * P, 2);
        d[i] = f(con->x[i]) * pow(con->h, 2);
    }

    double k = 1.0 / con->h;

    a[0] = 0;
    a[n - 1] = -k;

    c[0] = 1;
    c[n - 1] = STARTPB + k + con->h * pow(P, 2) / 2;

    b[0] = 0;
    b[n - 1] = 0;

    d[0] = STARTPA;
    d[n - 1] = -con->h/2 * f(con->x[n - 1]);

    solve(n, a, c, b, d, con->u);
}
double calc(config* con) {
    auto start = std::chrono::system_clock::now();

    con->u = new double[con->n];

    methodTridig(con);

    std::chrono::duration<double> elapsed = std::chrono::system_clock::now() - start;
    return elapsed.count();
}

void axis(int size, Graphic* g) {
    g->start();

    g->moveTo(-size, 0);
    g->lineTo(size, 0);

    g->moveTo(0, -size);
    g->lineTo(0, size);
}
void presets(Graphic* g) {
    g->setOfsx(250);
    g->setOfsy(250);
    g->setZoom(200);

    g->setColor(RGB(255, 255, 255));
    axis(100, g);
}
void chart(config* con, Graphic* g) {
    g->start();

    g->moveTo(con->x[0], con->u[0]);
    for (int i = 1; i < con->n; i++) {
        g->lineTo(con->x[i], con->u[i]);
    }
}
void chart(double (*func)(double x), COLORREF color, Graphic* g) {
    g->setColor(color);
    g->start();

    g->moveTo(A, func(A));
    for (double i = A; i <= B; i += 0.01) {
        g->lineTo(i, func(i));
    }
}

void chartEval(int nsteps, Graphic* g) {
    int n = NSTART;
    int colorSat = 10;
    config* con;

    for (int i = 0; i < nsteps; i++) {
        con = create(n + 1);
        calc(con);

        g->setColor(RGB(50, 50, 100 + colorSat));
        chart(con, g);
        delete con;

        n *= 2;
        colorSat *= 1.5;
    }
}

double getErr(config* con) {
    double locErr, maxErr = 0;

    for (int i = 0; i < con->n; i++) {
        locErr = abs(con->u[i] - u(con->x[i]));

        if (locErr > maxErr) {
            maxErr = locErr;
        }
    }

    return maxErr;
}
void printEval(int nsteps) {
    int n = NSTART;
    config* con;

    for (int i = 0; i < nsteps; i++) {
        con = create(n + 1);
        printf("[N: %d] ", con->nseg);
        printf("duration: %.4lf ", calc(con));
        printf("err: %.16lf h: %lf\n\n", getErr(con), con->h);

        delete con;
        n *= 2;
    }
}

void main() {
    Graphic* g = new Graphic();
    int power = 10;

    presets(g);
    axis(100, g);

    chartEval(power, g);

    _getch();

    chart(u, RGB(230, 230, 255), g);

    _getch();

    printEval(power);

    _getch();
}