lab5

using namespace std;

#include <math.h>
#include <iomanip>
#include <iostream>

const double a = 0, b = 7;
const int k = (b - a) * 50;
const double step = (b - a) / (double)k;
const double sqrEpsilon = 1.0e-6 * (k + 1);

double Function(double x) {
    return (x * x - sin(x)) * cos(2.5 * x);
}

double Legendre(double x, int N) {
    double Tn = x, Tnm1 = 1;
    if (N == 0) return Tnm1;

    for (int i = 1; i < N; i++) {
        double T = 2 * x * Tn - Tnm1;

        Tnm1 = Tn;
        Tn = T;
    }

    return Tn;
}

double Square(double t, int k1, int k2, int N) {
    double x = (2 * t - a - b) / (b - a);
    if (k2 == N) return Function(t) * Legendre(x, k1);
    else return Legendre(x, k1) * Legendre(x, k2);
}

double Simpson(int k1, int k2, int n, int N) {
    double h = (b - a) / (double) n;
    double sig1 = 0, sig2 = 0;
    double x = a + h;

    for (int i = 1; i < n; i++, x += h) {
        if (i % 2 == 0) sig2 += Square(x, k1, k2, N);
        else sig1 += Square(x, k1, k2, N);
    }

    return h / 3.0 * (2.0 * sig2 + 4.0 * sig1 + Square(a, k1, k2, N) + Square(b, k1, k2, N));
}

void Ak(double *X, int N) {
    double matrix[N][N + 1], tmatrix[N][N + 1];

    for (int i = 0; i < N; i++) {
        for (int j = 0; j <= N; j++) {
            int n = (int((b - a) * 100) >> 1) << 1;
            double In, I2n = Simpson(i, j, n, N);
            do {
                In = I2n;
                I2n = Simpson(i, j, n *= 2, N);
            } while (fabs(In - I2n) / 15 > 1e-8);

            matrix[i][j] = I2n;
        }
    }

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

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

        ln++;

        for (int i = 0; i < N; i++) {
            if (i == p) continue;
            double m = -matrix[i][q] / a;
            for (int j = 0; j < N + 1; j++) {
                if (j == q) continue;
                matrix[i][j] += matrix[p][j] * m;
            }
        }

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

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

    }

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

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

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

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

double Polinom(double x, double *A, int N) {
    double t = (2 * x - b - a) / (b - a), s = 0;
    for (int j = 0; j < N; j++) {
        s += A[j] * Legendre(t, j);
    }
    return s;
}

int main() {
    int N = 2;

    double sqrErr = 0;
    do {
        N *= 2;
        sqrErr = 0;
        double *ak = new double[N];
        Ak(ak, N);

        double x = a;
        for (int i = 0; i <= k; i++, x += step) {
            sqrErr += pow(Function(x) - Polinom(x, ak, N), 2);
        }
        delete[] ak;
    } while (sqrErr >= sqrEpsilon);

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

    double *ak = new double[N];
    Ak(ak, N);
    cout << "X     P(x)" << endl;
    double x = a;
    for (int i = 0; i <= k; i++, x += step) {
        printf("%-5g %g\n", x, Polinom(x, ak, N));
    }
    delete[] ak;
}