Cod sursa(job #1705865)

Utilizator AlexandruValeanuAlexandru Valeanu AlexandruValeanu Data 21 mai 2016 00:48:42
Problema Algoritmul lui Gauss Scor 60
Compilator cpp Status done
Runda Arhiva educationala Marime 16.33 kb
#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;
}