using namespace std;

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

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

double L(int n, double x){
    double Ln = x, Ln_first = 1;

    if (n == 0) return Ln_first;

    int i = 1;
    for (int i = 1; i < n; i++) {
        double Ln_next = (1.0 * (2 * i + 1) / (double)(i + 1)) * x * Ln - (1.0 * i / (double)(i + 1)) * Ln_first;
        Ln_first = Ln;
        Ln = Ln_next;
    }
    return Ln;
}

double func(int N, double t, int k1, int k2, double a, double b){
    double x = (2 * t - a - b) / (b - a);
    return L(k1, x) * (k2 == N ? f(t) : L(k2, x));
}

double simpson(int k1, int k2, double a, double b, int n, int N) {
    double sig1 = 0, sig2 = 0, h = (b - a) / n, x = a;
    for (int i = 1; i < n; i++) {
        x += h;
        if (i % 2 == 0) sig2 += func(N, x, k1, k2, a, b);
        else sig1 += func(N, x, k1, k2, a, b);
    }

    return h / 3 * (2 * sig2 + 4 * sig1 + func(N, a, k1, k2, a, b) + func(N, b, k1, k2, a, b));
}

double Integral(int k1, int k2, double a, double b, int N) {
    double eps = 1e-8;

    int n = 2.0 * ceil((b - a) / sqrt(sqrt(eps)) / 2.0);
    double In, I2n = simpson(k1, k2, a, b, n, N);
    do {
        In = I2n;
        I2n = simpson(k1, k2, a, b, n *= 2, N);
    } while (fabs(In - I2n) / 15.0 > eps);

    return I2n;
}

double* complete_pivoting(double *matrix, int N){
    int p, q, current_line = 0;
    double t[(N + 1) * (N + 1)], M[N], *x = new double[N];

    for (int k = 0; k < N - 1; k++){
        double max = 0;
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                float a = fabs(matrix[i * (N + 1) + j]);
                if (a > max) {
                    max = a;
                    p = i;
                    q = j;
                }
            }
        }


        for (int j = 0; j < N + 1; j++) {
            t[current_line * (N + 1) + j] = matrix[p * (N + 1) + j];
        }
        current_line++;

        for (int i = 0; i < N; i++) {
            M[i] = -matrix[i * (N + 1) + q] / max;
        }

        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N + 1; j++) {
                if (i != p) {
                    matrix[i * (N + 1) + j] += matrix[p * (N + 1) + j] * M[i];
                }
            }
        }

        for (int j = 0; j < N + 1; j++) {
            matrix[p * (N + 1) + j] = 0;
        }

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

    for (int i = 0; i < N; i++) {
        if (matrix[i * (N + 1) + N] != 0) {
            p = i;
        }
    }

    for (int j = 0; j < N + 1; j++) {
        t[current_line * (N + 1) + j] = matrix[p * (N + 1) + j];
    }

    for (int j = 0; j < N; j++) {
        x[j] = 0;
    }

    for (int i = N - 1; i >= 0; i--){
        double R = t[i * (N + 1) + N];
        for (int j = 0; j < N; j++) {
            if ((t[i * (N + 1) + j] != 0) && (x[j] != 0)) {
                R -= t[i * (N + 1) + j] * x[j];
            }
        }
        for (int j = 0; j < N; j++) {
            if ((t[i * (N + 1) + j] != 0) && (x[j] == 0)) {
                x[j] = R / t[i * (N + 1) + j];
            }
        }
    }

    return x;
}

double *Ak(double a, double b, int N){
    double A[(N + 1) * (N + 1)];

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N + 1; j++) {
            A[i * (N + 1) + j] = A[j * (N + 1) + i] = Integral(i, j, a, b, N);
        }
        A[i * (N + 1) + N] = Integral(i, N, a, b, N);
    }

    return complete_pivoting(A, N);
}

double Pm(double x, double a, double b, double *ak, int n){
    double t = (2 * x - b - a) / (b - a);
    double y = 0;
    for (int j = 0; j < n; ++j) {
        y += ak[j] * L(j, t);
    }
    return y;
}

int main(){
    double a = 20, b = 30, eps = 0.01;

    int k = (b - a) * 2, n = 2;
    double step = (b - a) / k, err;
    do {
        n *= 2;
        err = 0;
        double x = a;
        double *ak = Ak(a, b, n);
        for (int i = 0; i <= k; ++i){
            double delta = f(x) - Pm(x, a, b, ak, n);
            err += delta * delta;
            x += step;
        }
    } while (sqrt(err / (k + 1)) > eps);

    cout << "N = " << n << endl << endl;

    double *ak = Ak(a, b, n);
    cout << "   X          P(x)" << endl;
    for (double x = a; x <= b; x += step){
        printf("%7.2f%15.8f\n", x, Pm(x, a, b, ak, n));
    }

    return 0;
}