Pagini recente » Cod sursa (job #17609) | Cod sursa (job #3143133) | Cod sursa (job #1355108)
#include <fstream>
#include <algorithm>
#include <iomanip>
using namespace std;
ifstream in("gauss.in");
ofstream out("gauss.out");
const double err = 0.0000001;
int n, m;
double **mat;
double *solutions;
void read() {
in >> n >> m;
solutions = new double[m + 1];
mat = new double*[n + 1];
for (int i = 1; i <= n; ++i)
mat[i] = new double[m + 2];
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= m + 1; ++j)
in >> mat[i][j];
}
void solve_system() {
int i = 1;
int j = 1;
int x = 1;
// modify the system
while (i <= n && j <= m) {
// search for line with coeficient non zero
x = i;
while (x <= n && (mat[x][j] >= -err && mat[x][j] <= err)) ++x;
// variable is free
if (x == n + 1) {
++j;
continue;
}
// swap lines
if (x != i)
swap(mat[i], mat[x]);
// divide the line by mat[i][j], making him 1 (easier calculations)
// coefients from 1 to j - 1 are 0, so no need for division
for (int k = j + 1; k <= m + 1; ++k)
mat[i][k] /= mat[i][j];
mat[i][j] = 1;
// modify lower lines
// the coef is mat[k][j] because coef = mat[k][j] / mat[i][j] (=1)
for (int k = i + 1; k <= n; ++k) {
double coef = mat[k][j];
for (int t = j; t <= m + 1; ++t)
mat[k][t] -= (coef * mat[i][t]);
}
++i; ++j;
}
// extract the solutions
for (int i = n; i > 0; --i) {
int x = 1;
while (x <= m + 1 && !(mat[i][x] < -err || mat[i][x] > err)) ++x;
// found a line (0, 0, 0 ..., 0, a)
if (x == m + 1) {
out << "Imposibil\n";
return;
}
else if (x == m + 2)
continue;
// calculate solutions[x], no need to divide by mat[i][x] (=1)
solutions[x] = mat[i][m + 1];
for (int j = x + 1; j <= m; ++j)
solutions[x] -= (mat[i][j] * solutions[j]);
}
// print solutions
for (int i = 1; i <= m; ++i)
out << fixed << setprecision(11) << solutions[i] << ' ';
}
void cleanup() {
delete solutions;
for (int i = 1; i <= n; ++i)
delete mat[i];
delete mat;
}
int main() {
read();
solve_system();
cleanup();
return 0;
}