Pagini recente » Istoria paginii runda/oji_2006_10/clasament | Cod sursa (job #879164) | Cod sursa (job #2529833) | Monitorul de evaluare | Cod sursa (job #1283954)
#include <fstream>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <cmath>
#define EPS 0.0000001
using namespace std;
int N, M;
double** A;
double* X;
ifstream f("gauss.in");
ofstream g("gauss.out");
void debug()
{
for(int i = 1; i <= N; i++, cout << endl)
for(int j = 1; j <= M + 1; j++)
cout << A[i][j] << " ";
cout << endl;
}
void read_input()
{
cin >> N >> M;
A = new double*[N + 1];
for(int i = 1; i <= N; i++)
{
A[i] = new double[M + 2];
}
for(int i = 1; i <= N; i++)
{
for(int j = 1; j <= M + 1; j++)
{
cin >> A[i][j];
}
}
}
int argmax(const int& col)
{
int row_max = col;
for(int j = col + 1; j <= N; j++)
{
if(abs(A[j][col]) > abs(A[row_max][col]))
{
row_max = col;
}
}
return row_max;
}
void swap_rows(const int& k, const int& pivot_index)
{
for(int i = 1; i <= M + 1; i++)
{
swap(A[k][i], A[pivot_index][i]);
}
}
// Compute in-place row echelon form
bool compute_row_echelon_form()
{
//debug(A, N, M);
for(int k = 1; k < min(N, M); k++)
{
int pivot_index = argmax(k);
//cout << "Pivot index " << pivot_index << endl;
if(A[pivot_index][k] > -EPS && A[pivot_index][k] < EPS)
{
return false;
}
if(k != pivot_index)
{
swap_rows(k, pivot_index);
}
for(int i = k + 1; i <= N; i++)
{
for(int j = k + 1; j <= M + 1; j++)
{
A[i][j] = A[i][j] - A[k][j] * (A[i][k] / A[k][k]);
}
A[i][k] = 0;
}
//debug(A, N, M);
}
return true;
}
bool compute_unknowns()
{
//cout << "compute_unknowns " << endl;
//Calculul necunoscutelor.
for(int i = N; i >= 1; i--)
{
for(int j = 1; j <= M + 1; j++)
{
if(A[i][j] < -EPS || A[i][j] > EPS)
{
// Toti coeficientii sunt 0 iar rezultatul este diferit de 0.
if(j == M + 1) // adica ultima coloana, b-ul
{
return false;
}
// Calculam pe necunoscuta j = rezultatul ecuatiei
// i - necunoscutele j+1...M inmultite cu coeficientii lor din aceasta ecuatie.
X[j] = A[i][M + 1];
for(int k = j + 1; k <= M; k++)
{
X[j] -= X[k] * A[i][k];
}
X[j] /= A[i][i];
break;
}
}
}
return true;
}
void free_allocations()
{
for(int i = 1; i <= N; i++)
delete A[i];
delete A;
delete X;
}
int main()
{
cin.rdbuf(f.rdbuf());
cout.rdbuf(g.rdbuf());
read_input();
X = new double[M + 1];
bool ret = compute_row_echelon_form();
if(!ret)
{
cout << "Imposibil\n";
}
else
{
ret = compute_unknowns();
if(!ret)
{
cout << "Imposibil\n";
}
else
{
for(int i = 1; i <= M; i++)
{
cout << fixed << setprecision(12) << X[i] << " ";
}
}
}
free_allocations();
return 0;
}