Cod sursa(job #2301365)

Utilizator preda.andreiPreda Andrei preda.andrei Data 12 decembrie 2018 21:09:24
Problema Algoritmul lui Gauss Scor 100
Compilator cpp-64 Status done
Runda Arhiva educationala Marime 3.18 kb
#include <cmath>
#include <fstream>
#include <iomanip>
#include <vector>

using namespace std;

using System = vector<vector<double>>;

void SwapRows(System &system, int row1, int row2)
{
    if (row1 == row2) {
        return;
    }
    for (size_t i = 0; i < system[row1].size(); i += 1) {
        swap(system[row1][i], system[row2][i]);
    }
}

int FindPivotRow(const System &system, int col)
{
    for (size_t i = col; i < system.size(); i += 1) {
        if (system[i][col] != 0) {
            return i;
        }
    }
    return -1;
}

void TransformRow(System &system, int dest, int src, double coeff)
{
    for (size_t i = 0; i < system[src].size(); i += 1) {
        system[dest][i] += system[src][i] * coeff;
    }
}

void ClearDownwards(System &system, int col)
{
    auto coeff = system[col][col];
    for (size_t i = 0; i < system[col].size(); i += 1) {
        system[col][i] /= coeff;
    }

    for (size_t i = col + 1; i < system.size(); i += 1) {
        auto coeff = system[i][col] / system[col][col];
        TransformRow(system, i, col, -coeff);
    }
}

bool MakeSolvable(System &system)
{
    auto equations = system.size();
    auto vars = system[0].size() - 1;

    for (size_t i = 0; i < min(equations, vars); i += 1) {
        auto row = FindPivotRow(system, i);
        if (row < 0) {
            return false;
        }

        SwapRows(system, i, row);
        ClearDownwards(system, i);
    }
    return true;
}

double FindValueAtIndex(const vector<double> &equation, int index,
                        const vector<double> &vals)
{
    double other_sum = 0;
    for (size_t i = index + 1; i + 1 < equation.size(); i += 1) {
        other_sum += equation[i] * vals[i];
    }

    return (equation.back() - other_sum) / equation[index];
}

vector<double> ExtractValues(const System &system)
{
    auto equations = system.size();
    auto vars = system[0].size() - 1;
    vector<double> res(vars, 0);

    for (int i = min(equations, vars) - 1; i >= 0; i -= 1) {
        res[i] = FindValueAtIndex(system[i], i, res);
    }
    return res;
}

bool IsSolution(const vector<double> &equation, const vector<double> &res)
{
    double sum = 0;
    for (size_t i = 0; i < res.size(); i += 1) {
        sum += equation[i] * res[i];
    }

    static constexpr auto kEps = 1e-9;
    return fabs(sum - equation.back()) < kEps;
}

vector<double> Solve(System &system)
{
    if (!MakeSolvable(system)) {
        return vector<double>();
    }

    auto vars = system[0].size() - 1;
    auto res = ExtractValues(system);

    for (auto i = vars; i < system.size(); i += 1) {
        if (!IsSolution(system[i], res)) {
            return vector<double>();
        }
    }
    return res;
}

int main()
{
    ifstream fin("gauss.in");
    ofstream fout("gauss.out");

    int equations, vars;
    fin >> equations >> vars;

    System system(equations, vector<double>(vars + 1));
    for (auto &row : system) {
        for (auto &elem : row) {
            fin >> elem;
        }
    }

    auto res = Solve(system);
    if (res.empty()) {
        fout << "Imposibil\n";
    } else {
        for (const auto &elem : res) {
            fout << setprecision(9) << fixed << elem << " ";
        }
        fout << "\n";
    }

    return 0;
}