Pagini recente » Cod sursa (job #814902) | Cod sursa (job #2530596) | Cod sursa (job #2224715) | Cod sursa (job #2111507) | Cod sursa (job #2301263)
#include <algorithm>
#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 = -1e9;
return abs(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;
}