lab6

using namespace std;

#include <cmath>
#include <vector>
#include <cstdio>

template <int N>
struct Splain {
public:
    Splain(double a, double b, double (*f)(double)) : a(a), b(b), h((b - a) / (double) N), f(f) {
        coef_a();
        coef_b();
        coef_c();
        coef_d();
    }

    double operator()(double x) {
        int i = 0;
        double h = (b - a) / N;
        double xi = a + i * h;

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

        return A[i] + B[i] * (x - xi) + C[i] * (x - xi) * (x - xi) / 2 + D[i] * (x - xi) * (x - xi) * (x - xi) / 6;
    }
private:
    double A[N]{ 0 }, B[N]{ 0 }, C[N]{ 0 }, D[N]{ 0 }, a, b, h;
    double(*f)(double);

    void coef_a() {
        double x = a;
        for (int i = 0; i < N; i++) {
            A[i] = f(x);
            x += h;
        }
    }

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

    void coef_b() {
        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;
        }
    }

    void coef_c() {
        double *matrix = new double[N * (N + 1)];

        double h = (b - a) / N;
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N + 1; j++) {
                matrix[i * (N + 1) + j] = 0;
            }
        }

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

        for (int j = 0; j < N; j++) {
            double m = 1.0 / matrix[j * (N + 1) + j];
            for (int i = j; i < N + 1; i++) {
                matrix[j * (N + 1) + i] *= m;
            }
            for (int k = j + 1; k < N; k++) {
                double a = matrix[k * (N + 1) + j];
                for (int l = j; l < N + 1; l++) {
                    matrix[k * (N + 1) + l] -= matrix[j * (N + 1) + l] * a;
                }
            }
        }
        for (int j = N - 1; j >= 0; j--) {
            double s = 0;
            for (int i = N - 1; i > j; i--) {
                s += matrix[j * (N + 1) + i] * C[i];
            }
            C[j] = matrix[j * (N + 1) + N] - s;
        }
    }
};

int main() {
    auto f = [](double x) -> double {
        return 1 / x - 0.1 * x * x * sin(2 * x);
    };

    auto splain = Splain<1000>(2, 11, f);

    printf("|   x   |      y      |     error     |\n");
    for (double x = 2; x <= 11; x += 0.2) {
        double y = splain(x);

        printf("| %5.2f | %11.7f | %13.6e |\n", x, y, f(x) - y);

    }
}