using namespace std;

#include <iostream>
#include <iomanip>
#include <cmath>

double f(double x){
    return pow(M_E, (pow(x - 15, 0.25))) * sin(1.2 * x);
}

double *slar(double a, double b, int N){
    double *res = new double[N * (N + 1)] {0};
    double h = (b - a) / (double) N;

    for (int i = 0; i < N; i++){
        res[i * (N + 1) + i] = 4 * h;
        if (i > 0) res[i * (N + 1) + i - 1] = h;
        if (i < N - 1) res[i * (N + 1) + i + 1] = h;
        double x_prev = a + h * (i - 1);
        double x_i = a + h * i;
        double x_next = a + h * (i + 1);
        res[i * (N + 1) + N] = 6 * (f(x_next) - 2 * f(x_i) + f(x_prev)) / h;
    }
    return res;
}


double *gaus_jordan(double *A, int N) {
    double *res = new double[N];
    double coeficient;

    for (int k = 0; k < N; k++) {
        coeficient = A[k * (N + 1) + k];
        for (int j = 0; j < N + 1; j++) {
            A[k * (N + 1) + j] = A[k * (N + 1) + j] / coeficient;
        }

        for (int i = 0; i < N; i++) {
            if (i == k) continue;
            coeficient = A[i * (N + 1) + k];
            for (int j = 0; j < N + 1; j++) {
                A[i * (N + 1) + j] -= A[k * (N + 1) + j] * coeficient;
            }
        }

    }

    for (int i = 0; i < N; i++) {
        res[i] = A[i * (N + 1) + N];
    }

    return res;
}


double* ai(double a, double b, int N){
    double* A = new double[N];
    double h = (b - a) / N;
    double x = a;
    for (int i = 0; i < N; ++i) {
        A[i] = f(x);
        x += h;
    }
    return A;
}

double* di(double a, double b, double* C, int N){
    double* D = new double[N];
    double h = (b - a) / N;

    D[0] = C[0];
    for (int i = 1; i < N; ++i) {
        D[i] = (C[i] - C[i - 1]) / h;
    }
    return D;
}

double* bi(double a, double b, double* C, double* D, int N){
    double* B = new double[N];
    double h = (b - a) / N;
    B[0] = 0;

    for (int i = 1; i < N; ++i) {
        B[i] = (h * C[i] / 2) - (h * h * D[i] / 2) + (f(a+h*i) - f(a+h*(i-1))) / h;
    }

    return B;
}

double spline(double x, double a, double b, double* A, double* B, double* C, double* D, int N){
    double h = (b - a) / N;
    int i = 0;

    while (x >= (a + i * h)) i++;
    i--;

    double xi = a + i * h;

    return A[i] + B[i] * (x - xi) + C[i] * (x - xi) * (x - xi) / 2.0 + D[i] * (x - xi) * (x - xi) * (x - xi) / 6.0;
}

int main(){
    double a = 20, b = 30;
    int N = 50;

    double *M = slar(a, b, N);
    double *A = ai(a, b, N);
    double *C = gaus_jordan(M, N);
    double *D = di(a, b, C, N);
    double *B = bi(a, b, C, D, N);

    int k = (b - a) * N + 1;
    double step = (b - a) / (double)k;
    double x = a, y;

    cout << "   X          Y" << endl;
    for (int i = 0; i <= k; i++){
        cout << setw(6) << x << setw(12) << spline(x, a, b, A, B, C, D, N) << endl;
        x += step;
    }

    return 0;
}