Pagini recente » Cod sursa (job #2048930) | Cod sursa (job #889285) | Cod sursa (job #265214) | Cod sursa (job #2437871) | Cod sursa (job #1246171)
#include <fstream>
#include <iomanip>
#define Nmax 310
#define eps 0.000000001
#define nonZero(x) ((x) < -eps || eps < (x))
using namespace std;
int N,M;
double A[Nmax][Nmax], X[Nmax];
bool Solution;
void findSolution() {
int i, j, k;
Solution = true;
// Find a solution by substitution
// Where X[i] is variable, X[i] = 0
for(i = N; i >= 1; i--)
for(j = 1; j <= M+1; j++)
if(nonZero(A[i][j])) {
if(j == M+1) {
// No solution ( 0 = non zero number )
Solution = false;
return;
}
else {
// We calculate X[j]
X[j] = A[i][M+1];
for(k = j+1; k <= M; k++)
X[j] -= X[k] * A[i][k];
break;
}
}
}
void reduceEquations(int i, int j) {
int l, c;
// Divide equation i by the coefficient of X[j] ( A[i][j] )
for(c = j+1; c <= M+1; c++)
A[i][c] /= A[i][j];
A[i][j] = 1;
// Eliminate X[j] from the next equations
for(l = i+1; l <= N; l++) {
for(c = j+1; c <= M+1; c++)
A[l][c] -= A[i][c] * A[l][j];
A[l][j] = 0;
}
}
void swapEquations(int a, int b) {
if(a == b)
return;
for(int i = 1; i <= M+1; i++)
swap(A[a][i], A[b][i]);
}
int findXj(int start, int j) {
for(int i = start; i <= N; i++)
if(nonZero(A[i][j]))
return i;
return -1;
}
void Read() {
int i, j;
ifstream in("gauss.in");
in >> N >> M;
for(i = 1; i <= N; i++)
for(j = 1; j <= M+1; j++)
in >> A[i][j];
in.close();
}
void writeSolution(ofstream & out) {
if(!Solution)
out << "Imposibil";
else
for(int i = 1; i <= M; i++)
out << fixed << setprecision(11) << X[i] << ' ';
out << '\n';
}
void Solve() {
int i, j, k;
// Output file
ofstream out("gauss.out");
// Gauss elimination
for(i = 1, j = 1; (i <= N && j <= M); i++, j++) {
// Try to find an equation which contains X[j]
if((k = findXj(i, j)) == -1)
continue; // Not found
// Swap equations (if necessary) to preserve the row echelon form
swapEquations(i, k);
// reduced row echelon
reduceEquations(i, j);
}
// Try to find a solution
findSolution();
// Write solution
writeSolution(out);
out.close();
}
int main() {
// Read A and b (A * x = b)
Read();
// Solve the linear system
Solve();
return 0;
}