#include <bits/stdc++.h>
///#include <boost/rational.hpp>
///#include <boost/multiprecision/cpp_int.hpp>
using namespace std;
///using namespace boost;
#include <iostream>
#include <memory>
#include <algorithm>
#include <vector>
#include <cassert>
#include <ostream>
#include <sstream>
template<typename T>
class Matrix
{
public:
Matrix();
Matrix(size_t, size_t);
void resize(size_t, size_t);
typename std::vector<T>::iterator begin();
typename std::vector<T>::iterator end();
typename std::vector<T>::const_iterator begin() const;
typename std::vector<T>::const_iterator end() const;
/*Retrieves a single element in the Matrix. Bounds-checking is performed.*/
T& at(size_t, size_t);
const T& at(size_t, size_t) const;
T& operator()(size_t, size_t);
const T& operator()(size_t, size_t) const;
/* Compound addition and subtraction. */
Matrix& operator += (const Matrix&);
Matrix& operator -= (const Matrix&);
/* Compound scalar multiplication and division. */
Matrix& operator *= (const T&);
Matrix& operator /= (const T&);
size_t rows() const;
size_t columns() const;
size_t size() const;
void addToLine(size_t, const T&);
void addToColumn(size_t, const T&);
void multiplyLine(size_t, const T&);
void multiplyColumn(size_t, const T&);
void swapLines(size_t, size_t);
void swapColumns(size_t, size_t);
void linearCombinationLine(size_t, size_t, size_t, const T&, const T&);
void linearCombinationColumn(size_t, size_t, size_t, const T&, const T&);
private:
std::vector<T> data;
size_t M, N;
};
template<typename T>
Matrix<T>::Matrix() : data{}, M{0}, N{0} {
}
template<typename T>
Matrix<T>::Matrix(size_t i, size_t j) : data{}, M{i}, N{j} {
this->data.resize(i * j);
}
template<typename T>
void Matrix<T>::resize(size_t i, size_t j)
{
this->M = i;
this->N = j;
this->data.resize(i * j);
}
template <typename T>
typename std::vector<T>::iterator Matrix<T>::begin()
{
return this->data.begin();
}
template <typename T>
typename std::vector<T>::iterator Matrix<T>::end()
{
return this->data.end();
}
template <typename T>
typename std::vector<T>::const_iterator Matrix<T>::begin() const
{
return this->data.begin();
}
template <typename T>
typename std::vector<T>::const_iterator Matrix<T>::end() const
{
return this->data.end();
}
template <typename T>
T& Matrix<T>::at(size_t i, size_t j)
{
if (i >= this->M || j >= this->N)
throw std::out_of_range{"Matrix::at"};
return this->data[i * N + j];
}
template <typename T>
const T& Matrix<T>::at(size_t i, size_t j) const
{
if (i >= this->M || j >= this->N)
throw std::out_of_range{"Matrix::at"};
return this->data[i * N + j];
}
template <typename T>
T& Matrix<T>::operator ()(size_t i, size_t j)
{
return this->at(i, j);
}
template <typename T>
const T& Matrix<T>::operator ()(size_t i, size_t j) const
{
return this->at(i, j);
}
template <typename T>
Matrix<T>& Matrix<T>::operator += (const Matrix<T>& rhs)
{
assert(this->M == rhs.rows());
assert(this->N == rhs.columns());
std::transform(this->begin(), this->end(), rhs.begin(), this->begin(), std::plus<T>());
return *this;
}
template <typename T>
Matrix<T>& Matrix<T>::operator -= (const Matrix<T>& rhs)
{
assert(this->M == rhs.rows());
assert(this->N == rhs.columns());
std::transform(this->begin(), this->end(), rhs.begin(), this->begin(), std::minus<T>());
return *this;
}
template <typename T>
Matrix<T>& Matrix<T>::operator *= (const T& scalar)
{
std::transform(this->begin(), this->end(), this->begin(),
std::bind2nd(std::multiplies<T>(), scalar));
return *this;
}
template <typename T>
Matrix<T>& Matrix<T>::operator /= (const T& scalar)
{
std::transform(this->begin(), this->end(), this->begin(),
std::bind2nd(std::divides<T>(), scalar));
return *this;
}
template <typename T>
size_t Matrix<T>::rows() const
{
return this->M;
}
template <typename T>
size_t Matrix<T>::columns() const
{
return this->N;
}
template <typename T>
size_t Matrix<T>::size() const
{
return this->M * this->N;
}
template<typename T>
std::istream& operator >> (std::istream& stream, Matrix<T>& lhs)
{
size_t n, m;
stream >> n >> m;
lhs.resize(n, m);
for (size_t i = 0; i < n; ++i)
for (size_t j = 0; j < m; ++j)
stream >> lhs(i, j);
return stream;
}
template<typename T>
std::ostream& operator << (std::ostream& stream, const Matrix<T>& lhs)
{
std::stringstream sstream;
for (size_t i = 0; i < lhs.rows(); ++i)
{
for (size_t j = 0; j < lhs.columns(); ++j)
{
sstream << lhs(i, j);
if (j + 1 != lhs.columns())
sstream << " ";
}
if (i + 1 != lhs.rows())
sstream << "\n";
}
stream << sstream.str();
return stream;
}
/* Binary addition and subtraction. */
template <typename T>
const Matrix<T> operator + (const Matrix<T>& lhs, const Matrix<T>& rhs)
{
assert(lhs.rows() == rhs.rows());
assert(lhs.columns() == rhs.columns());
return Matrix<T>(lhs) += rhs;
}
template <typename T>
const Matrix<T> operator - (const Matrix<T>& lhs, const Matrix<T>& rhs)
{
assert(lhs.rows() == rhs.rows());
assert(lhs.columns() == rhs.columns());
return Matrix<T>(lhs) -= rhs;
}
/* Binary operators implemented by compound assignment operators. */
template <typename T>
const Matrix<T> operator * (const Matrix<T>& lhs, const T& scalar)
{
return Matrix<T>(lhs) *= scalar;
}
template <typename T>
const Matrix<T> operator / (const Matrix<T>& lhs, const T& scalar)
{
return Matrix<T>(lhs) /= scalar;
}
template <typename T>
const Matrix<T> operator * (const T& scalar, const Matrix<T>& lhs)
{
return Matrix<T>(lhs) *= scalar;
}
template <typename T>
const Matrix<T> operator / (const T& scalar, const Matrix<T>& lhs)
{
return Matrix<T>(lhs) /= scalar;
}
/* Unary + just returns a copy of its argument. */
template <typename T>
const Matrix<T> operator + (const Matrix<T>& lhs)
{
return lhs;
}
/* Unary minus negates its argument. */
template <typename T>
const Matrix<T> operator - (const Matrix<T>& lhs)
{
return Matrix<T>(lhs) *= T(-1);
}
template<typename T>
const Matrix<T> operator * (const Matrix<T>& lhs, const Matrix<T>& rhs)
{
size_t N = lhs.rows();
size_t M = lhs.columns();
size_t P = rhs.columns();
assert(lhs.columns() == rhs.rows());
Matrix<T> result(N, P);
for (size_t i = 0; i < M; ++i)
for (size_t k = 0; k < N; ++k)
for (size_t j = 0; j < P; ++j)
result(i, j) += lhs(i, k) * rhs(k, j);
return result;
}
/* Special-case: Multiplication with assignment. */
template <typename T>
Matrix<T>& operator *= (Matrix<T>& lhs, const Matrix<T>& rhs)
{
return lhs = lhs * rhs;
}
template<typename T>
bool operator == (const Matrix<T>& lhs, const Matrix<T>& rhs)
{
return std::equal(lhs.begin(), lhs.end(), rhs.begin());
}
template<typename T>
bool operator != (const Matrix<T>& lhs, const Matrix<T>& rhs)
{
return !(lhs == rhs);
}
template<typename T>
const Matrix<T> transpose(const Matrix<T>& lhs)
{
size_t N = lhs.rows();
size_t M = lhs.columns;
Matrix<T> result(N, M);
for (size_t i = 0; i < M; ++i)
for (size_t j = 0; j < N; ++j)
result(j, i) = lhs(i, j);
return result;
}
template<typename T>
const Matrix<T> identity(size_t N)
{
Matrix<T> result(N, N);
for (size_t i = 0; i < N; ++i)
result(i, i) = T(1);
return result;
}
template<typename T>
const Matrix<T> matrix_exponential(const Matrix<T>& lhs, size_t expo)
{
assert(lhs.rows() == lhs.columns());
size_t M = lhs.rows();
Matrix<T> result = identity<T>(M);
Matrix<T> A = lhs;
while (expo > 0)
{
if (expo & 1)
result *= A;
A *= A;
expo >>= 1;
}
return result;
}
template<typename T>
void Matrix<T>::addToLine(size_t x, const T& scalar)
{
for (size_t i = 0; i < this->columns(); ++i)
this->at(x, i) += scalar;
}
template<typename T>
void Matrix<T>::addToColumn(size_t x, const T& scalar)
{
for (size_t i = 0; i < this->rows(); ++i)
this->at(i, x) += scalar;
}
template<typename T>
void Matrix<T>::multiplyLine(size_t x, const T& scalar)
{
for (size_t i = 0; i < this->columns(); ++i)
this->at(x, i) *= scalar;
}
template<typename T>
void Matrix<T>::multiplyColumn(size_t x, const T& scalar)
{
for (size_t i = 0; i < this->rows(); ++i)
this->at(i, x) *= scalar;
}
template<typename T>
void Matrix<T>::swapLines(size_t x, size_t y)
{
for (size_t i = 0; i < this->columns(); ++i)
std::swap(this->at(x, i), this->at(y, i));
}
template<typename T>
void Matrix<T>::swapColumns(size_t x, size_t y)
{
for (size_t i = 0; i < this->rows(); ++i)
std::swap(this->at(i, x), this->at(i, y));
}
template<typename T>
void Matrix<T>::linearCombinationLine(size_t dest, size_t src1, size_t src2, const T& scalar1, const T& scalar2)
{
for (size_t i = 0; i < this->columns(); ++i)
this->at(dest, i) = scalar1 * this->at(src1, i) + scalar2 * this->at(src2, i);
}
template<typename T>
void Matrix<T>::linearCombinationColumn(size_t dest, size_t src1, size_t src2, const T& scalar1, const T& scalar2)
{
for (size_t i = 0; i < this->rows(); ++i)
this->at(i, dest) = scalar1 * this->at(i, src1) + scalar2 * this->at(i, src2);
}
namespace LinearAlgebra
{
///------------------------------------------------------------
/// using T = boost::multiprecision::cpp_int;
///using TYPE = rational<T>;
using TYPE = double;
const TYPE EPSILON = 1e-8;///T(0);
const TYPE ONE = 1;///T(1);
///------------------------------------------------------------
///------------------------------------------------------------
bool checkIfInvertible(const Matrix<TYPE> &A)
{
assert(A.rows() == A.columns());
for (size_t i = 0; i < A.rows(); ++i)
if (abs(A(i, i)) <= EPSILON)
return false;
return true;
}
bool checkIfIdentity(const Matrix<TYPE> &A)
{
return A == identity<TYPE>(A.rows());
}
template<typename S>
Matrix<TYPE> createMatrix(const vector<vector<S>> &coefs)
{
int n = coefs.size();
assert(n >= 1);
int m = coefs[0].size();
assert(m >= 1);
Matrix<TYPE> A(n, m);
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j)
A(i, j) = coefs[i][j];
return A;
}
///------------------------------------------------------------
void rowEchelonForm(Matrix<TYPE> &A, Matrix<TYPE> &B, TYPE &determinant, vector<pair<int,int>> &inversions)
{
int N = A.rows();
int M = A.columns();
int i = 0, j = 0;
while (i < N && j < M)
{
int k = i;
while (k < N && abs(A(k, j)) <= EPSILON)
k++;
if (k == N)
{
j++;
continue;
}
if (k != i)
{
if (B.size() > 0)
B.swapLines(k, i);
A.swapLines(k, i);
determinant *= -ONE;
inversions.push_back({k, i});
}
determinant *= A(i, j);
if (abs(A(i, j) - ONE) > EPSILON)
{
TYPE coef = A(i, j);
if (B.size() > 0)
B.multiplyLine(i, ONE / coef);
A.multiplyLine(i, ONE / coef);
}
for (int l = i + 1; l < N; ++l)
if (abs(A(l, j)) > EPSILON)
{
TYPE coef = A(l, j);
if (B.size() > 0)
B.linearCombinationLine(l, l, i, ONE, -coef);
A.linearCombinationLine(l, l, i, ONE, -coef);
}
i++;
j++;
}
}
void reducedRowEchelonForm(Matrix<TYPE> &A, Matrix<TYPE> &B)
{
int N = A.rows();
int M = A.columns();
for (int i = N - 1; i >= 0; i--)
{
int j = 0;
while (j <= M - 1 && abs(A(i, j)) <= EPSILON)
j++;
if (j < M)
{
for (int k = i - 1; k >= 0; k--)
if (abs(A(k, j)) > EPSILON)
{
TYPE coef = A(k, j);
B.linearCombinationLine(k, k, i, ONE, -coef);
A.linearCombinationLine(k, k, i, ONE, -coef);
}
}
}
}
template<typename S>
double computeDeterminant(const vector<vector<S>> &coef)
{
assert(coef.size() >= 1);
assert(coef[0].size() >= 1);
int N = coef.size();
int M = coef[0].size();
if (N != M)
{
cerr << "Matrix not invertible : N != M (non-square)\n";
return {};
}
Matrix<TYPE> A = createMatrix(coef);
Matrix<TYPE> B;
TYPE determinant = ONE;
vector<pair<int,int>> invs{};
rowEchelonForm(A, B, determinant, invs);
///return boost::rational_cast<double>(determinant);
return ONE;
}
template<typename S>
vector<vector<double>> computeInverse(const vector<vector<S>> &coef)
{
assert(coef.size() >= 1);
assert(coef[0].size() >= 1);
int N = coef.size();
int M = coef[0].size();
if (N != M)
{
cerr << "Matrix not invertible : N != M (non-square)\n";
return {};
}
Matrix<TYPE> A = createMatrix(coef);
Matrix<TYPE> B = identity<TYPE>(N);
TYPE determinant = ONE;
vector<pair<int,int>> invs{};
rowEchelonForm(A, B, determinant, invs);
if (!checkIfInvertible(A))
{
cerr << "Matrix not invertible\n";
return {};
}
reducedRowEchelonForm(A, B);
vector<vector<double>> inverse(N);
for (int i = 0; i < N; ++i)
{
inverse[i].resize(M);
///for (int j = 0; j < M; ++j)
///inverse[i][j] = boost::rational_cast<double>(B(i, j));
}
cout << B << endl;
///CHECKING PART-----------------------------
Matrix<TYPE> tmp = createMatrix(coef);
Matrix<TYPE> D = tmp * B;
assert(checkIfIdentity(D));
Matrix<TYPE> E = B * tmp;
assert(checkIfIdentity(E));
///------------------------------------------
return inverse;
}
template<typename S>
vector<double> solveSystemOfEquations(const vector<vector<S>> &coef, const vector<S> &bs)
{
assert(coef.size() >= 1);
assert(coef[0].size() >= 1);
int N = coef.size();
int M = coef[0].size();
assert(static_cast<int>(bs.size()) == N);
Matrix<TYPE> A = createMatrix(coef);
Matrix<TYPE> B(N, 1);
for (int i = 0; i < N; ++i)
B(i, 0) = bs[i];
TYPE determinant = ONE;
vector<pair<int,int>> invs{};
rowEchelonForm(A, B, determinant, invs);
reducedRowEchelonForm(A, B);
Matrix<TYPE> solutions(1, M);
for (int i = N - 1; i >= 0; i--)
{
int j = 0;
while (j < M && abs(A(i, j)) <= EPSILON)
j++;
if (j == M)
{
if (abs(B(i, 0)) > EPSILON)
{
cerr << "Impossible!" << endl;
return {};
}
}
else
solutions(0, i) = B(i, 0);
}
vector<double> solution(M);
for (int j = 0; j < M; ++j)
///solution[j] = boost::rational_cast<double>(solutions(0, j));
solution[j] = solutions(0, j);
return solution;
}
}
int main()
{
ifstream in("gauss.in");
ofstream out("gauss.out");
int n, m;
in >> n >> m;
vector<vector<int>> A(n);
vector<int> bs(n);
for (int i = 0; i < n; ++i)
{
A[i].resize(m);
for (int j = 0; j < m; ++j)
in >> A[i][j];
in >> bs[i];
}
auto it = LinearAlgebra::solveSystemOfEquations(A, bs);
if (it.size() == 0)
out << "Imposibil\n";
else
{
for (auto x : it)
out << fixed << setprecision(10) << x << " ";
out << endl;
}
return 0;
}