Cod sursa(job #2705129)

Utilizator flibiaVisanu Cristian flibia Data 11 februarie 2021 22:51:26
Problema Algoritmul lui Gauss Scor 100
Compilator cpp-64 Status done
Runda Arhiva educationala Marime 2.69 kb
// mizerie pentru n != m
#include <bits/stdc++.h>
#define ll long long
#define EPS 1e-10

using namespace std;

ifstream in("gauss.in");
ofstream out("gauss.out");

int n, m;
vector<vector<double> > a;

vector<double> substitutie_descendenta(vector<vector<double> > a) {
    int n = a.size();
    int m = a[0].size() - 1;
    n = min(n, m);

    vector<double> sol(m, 1);
    for (int k = n - 1; k >= 0; k--) {
        double sum = 0;
        for (int j = k + 1; j < m; j++)
            sum += a[k][j] * sol[j];

        sol[k] = (a[k][m] - sum) / a[k][k];
    }

    return sol;
}

vector<double> gauss_pivotare_totala(vector<vector<double> > a) {
    int n = a.size();
    int m = a[0].size() - 1;
    n = min(n, m);

    vector<int> perm(m, 0);
    for (int i = 0; i < m; i++) {
        perm[i] = i;
    }

    for (int k = 0; k < n; k++) {
        int row = 0, col = 0;
        double mx = 0;

        for (int i = k; i < n; i++) {
            for (int j = k; j < n; j++) {
                if (abs(a[i][j]) > mx) {
                    mx = abs(a[i][j]);
                    row = i;
                    col = j;
                }
            }
        }

        if (abs(mx) < EPS) {
            throw "Imposibil";
        }

        if (k != row) {
            for (int j = 0; j <= m; j++) {
                swap(a[k][j], a[row][j]);
            }
        }
        if (k != col) {
            for (int i = 0; i < n; i++) {
                swap(a[i][k], a[i][col]);
            }
            swap(perm[k], perm[col]);
        }

        for (int i = k + 1; i < n; i++) {
            double ratio = a[i][k] / a[k][k]; 

            for (int j = k; j <= m; j++) {
                a[i][j] -= ratio * a[k][j];
            }
        }
    }

    auto aux = substitutie_descendenta(a);

    auto sol = aux;

    for (int i = 0; i < m; i++) {
        sol[perm[i]] = aux[i];
    }

    return sol;
}

int main() {
    in >> n >> m;
    a = vector<vector<double> >(n, vector<double>(m + 1, 0));
    
    for (int i = 0; i < n; i++) {
        for (int j = 0; j <= m; j++) {
            in >> a[i][j];
        }
    }

    try {
        auto sol = gauss_pivotare_totala(a);

        for (int i = 0; i < n; i++) {
            double sum = 0;
            for (int j = 0; j < m; j++) {
                sum += a[i][j] * sol[j]; 
            }

            if (abs(sum - a[i][m]) > 0.0001) {
                return out << "Imposibil", 0;
            }
        }
        
        for (auto it : sol) {
            out << fixed << setprecision(8) << it << ' ';
        }
    } catch (const char *msg) {
        out << msg;
    }

    return 0;
}