Cod sursa(job #2301267)

Utilizator preda.andreiPreda Andrei preda.andrei Data 12 decembrie 2018 19:55:23
Problema Algoritmul lui Gauss Scor 100
Compilator cpp-64 Status done
Runda Arhiva educationala Marime 3.09 kb
#include <algorithm>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <numeric>
#include <vector>

using namespace std;

template <class T>
using Matrix = vector<vector<T>>;

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

void SwapRows(Matrix<double> &system, size_t row1, size_t row2)
{
    swap_ranges(system[row1].begin(), system[row1].end(), system[row2].begin());
}

void MultiplyAll(vector<double> &vec, double val)
{
    auto mult = [val] (double elem) -> double {
        return elem * val;
    };
    transform(vec.begin(), vec.end(), vec.begin(), mult);
}

void ClearDownwards(Matrix<double> &system, size_t col)
{
    for (auto i = col + 1; i < system.size(); i += 1) {
        auto coeff =  system[i][col] / system[col][col];
        for (auto j = col; j < system[i].size(); j += 1) {
            system[i][j] -= system[col][j] * coeff;
        }
    }
}

bool MakeSolvable(Matrix<double> &system)
{
    for (size_t i = 0; i + 1 < system.size(); i += 1) {
        if (system[i][i] == 0) {
            auto row = FindPivotRow(system, i);
            if (row < 0) {
                return false;
            }
            SwapRows(system, i, row);
        }
        MultiplyAll(system[i], 1.0 / system[i][i]);
        ClearDownwards(system, i);
    }
    return true;
}

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

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

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(Matrix<double> &system)
{
    if (!MakeSolvable(system)) {
        return vector<double>();
    }

    vector<double> res(system[0].size() - 1, 0);
    auto equations = min(system.size(), system[0].size() - 1);

    for (int i = equations - 1; i >= 0; i -= 1) {
        res[i] = FindValue(system[i], i, res);
    }

    if (equations < system.size()) {
        for (auto i = equations; 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;

    Matrix<double> 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;
}