gauss.cpp

#include <vector>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <fstream>

using namespace std;

ofstream fout ("out.txt");

void print(const vector<vector<double>> & v) {
    for (auto & row : v) {
        for (auto & it : row)
            fout << setw(12) << it;
        fout << endl;
    }
}

vector<double> solve(vector<vector<double>> & D) {
    int m = D.size(), n = D.front().size(), now = 0;
    for (int j = 0, i; j < n-1; ++j) {
        for (i = now; i < m && !D[i][j]; ++i);
        if ( i == m ) continue;
        swap(D[i], D[now]);
        for (i = now+1; i < m; ++i) {
            if ( !D[i][j] ) continue;
            double c = D[i][j] / D[now][j];
            for (int k = j; k < n; ++k)
                D[i][k] -= c * D[now][k];
        }
        ++now;
    }

    now = m-1;

    for (int j = n-2; j >= 0; --j, --now) {
        if ( !D[now][j] ) continue;

        for (int i = now-1; i >= 0; --i) {
            if ( !D[i][j] ) continue;
            double c = D[i][j] / D[now][j];
            for (int k = n-1; k >= 0; --k)
                D[i][k] -= c * D[now][k];
        }
    }

    vector<double> answer(m, 0);
    for (int i = 0; i < m; ++i)
        if ( D[i][i] ) {
            D[i].back() /= D[i][i];
            D[i][i] = 1;
            answer[i] = D[i].back();
        }
    return answer;
}

int main() {
    vector<vector<double>> f = {
        {-1.0, 0.86603},
        {0.0, 1.0},
        {1.0, 0.86603},
        {2.0, 0.5},
        {3.0, 0.0},
        {4.0, -0.5}
    };

    for (int n = 1; n <= 6; ++n) {
        fout << "--------------------------------------------------------------" << endl;
        fout << "n = " << n << endl;
        vector<vector<double>> s (n+1, vector<double>(n+2, 0.0));
        for (int k = 0; k <= n; ++k) {
            for (int i = 0; i <= n; ++i)
                for (int j = 0; j < f.size(); ++j)
                    s[k][i] += pow(f[j][0], k+i);

            for (int j = 0; j < f.size(); ++j)
                s[k].back() += f[j][1] * pow(f[j][0], k);
        }

        fout << "Before:" << endl;
        print(s);
        fout << endl;

        auto a = solve(s);

        fout << "After:" << endl;
        print(s);
        fout << endl;

        fout << "Answer:" << endl;
        for (auto & it : a)
            fout << setw(12) << it;

        fout << endl;

        double F = 0;
        for (int j = 0; j < f.size(); ++j) {
            double r = 0;
            for (int i = 0; i <= n; ++i)
                r += a[i] * pow(f[j][0], i);
            r -= f[j][1];
            F += r*r;
        }
        fout << endl << "F = " << F << endl;
    }

    return 0;
}