Cod sursa(job #1704769)

Utilizator AlexandruValeanuAlexandru Valeanu AlexandruValeanu Data 19 mai 2016 12:02:13
Problema Algoritmul lui Gauss Scor 100
Compilator cpp Status done
Runda Arhiva educationala Marime 12.8 kb
#include <bits/stdc++.h>
///#include <boost/rational.hpp>
///#include <boost/multiprecision/cpp_int.hpp>

using namespace std;
///using namespace boost;

namespace LinearAlgebra
{
    ///------------------------------------------------------------
    ///using T = boost::multiprecision::cpp_int;
    using TYPE = double;///rational<T>;
    using Vector = vector<TYPE>;
    using Matrix = vector<Vector>;

    const TYPE EPSILON = 1e-8;///T(0);
    ///------------------------------------------------------------

    ///------------------------------------------------------------
    bool checkIfInvertible(Matrix &A)
    {
        assert(A.size() >= 1);
        assert(A.size() == A[0].size());

        for (size_t i = 0; i < A.size(); ++i)
            if (abs(A[i][i]) <= EPSILON)
                return false;

        return true;
    }

    bool checkIfIdentity(const Matrix &A)
    {
        assert(A.size() >= 1);
        assert(A.size() == A[0].size());

        for (size_t i = 0; i < A.size(); ++i)
            for (size_t j = 0; j < A.size(); ++j)
            {
                if (i == j && abs(A[i][j] - 1) > EPSILON)
                    return false;

                if (i != j && abs(A[i][j]) > EPSILON)
                    return false;
            }

        return true;
    }

    template<typename S>
    Vector createVector(const vector<S> &a)
    {
        int n = a.size();
        assert(n >= 1);

        Vector sol;
        sol.resize(n);

        for (size_t i = 0; i < a.size(); ++i)
            sol[i] = a[i];

        return sol;
    }

    Vector createVector(int n)
    {
        assert(n >= 1);

        Vector sol;
        sol.resize(n);

        return sol;
    }

    Matrix createIdentity(int n)
    {
        assert(n >= 1);

        Matrix sol;
        sol.resize(n);

        for (int i = 0; i < n; ++i)
        {
            sol[i].resize(n);
            sol[i][i] = static_cast<TYPE>(1);
        }

        return sol;
    }

    Matrix createTranspose(const Matrix &A)
    {
        int n = A.size();
        assert(n >= 1);

        int m = A[0].size();
        assert(m >= 1);

        Matrix sol;
        sol.resize(m);

        for (int i = 0; i < m; ++i)
        {
            sol[i].resize(n);

            for (int j = 0; j < n; ++j)
                sol[i][j] = A[j][i];
        }

        return sol;
    }

    Matrix createMatrix(int n, int m)
    {
        assert(n >= 1);
        assert(m >= 1);

        Matrix sol;
        sol.resize(n);

        for (int i = 0; i < n; ++i)
            sol[i].resize(m);

        return sol;
    }

    template<typename S>
    Matrix createMatrix(const vector<vector<S>> &coefs)
    {
        int n = coefs.size();
        assert(n >= 1);

        int m = coefs[0].size();
        assert(m >= 1);

        Matrix sol;
        sol.resize(n);

        for (int i = 0; i < n; ++i)
        {
            sol[i].resize(m);

            for (int j = 0; j < m; ++j)
                sol[i][j] = coefs[i][j];
        }

        return sol;
    }

    Matrix matrix_multiplication(const Matrix &A, const Matrix &B)
    {
        int N = A.size();
        assert(N >= 1);

        int M = A[0].size();
        assert(M >= 1);

        int P = B.size();
        assert(P >= 1);

        int Q = B[0].size();
        assert(Q >= 1);

        assert(M == P);

        Matrix D = createMatrix(N, Q);

        for (int i = 0; i < N; ++i)
            for (int k = 0; k < M; ++k)
                for (int j = 0; j < Q; ++j)
                    D[i][j] += A[i][k] * B[k][j];

        return D;
    }

    void printMaxtrix(const Matrix &A)
    {
        int N = A.size();
        assert(N >= 1);

        int M = A[0].size();
        assert(M >= 1);

        for (auto v : A)
        {
            for (auto x : v)
                cerr << x << " ";

            cerr << endl;
        }
    }
    ///------------------------------------------------------------

    ///------------------------------------------------------------
    void swapLines(Matrix &A, int x, int y) /// swap(X, Y)
    {
        assert(x != y);
        swap(A[x], A[y]);

        ///cerr << "SWAP LINES : " << x << " " << y << endl;
    }

    void divideLine(Matrix &A, int line, TYPE rat) /// X = X / rat
    {
        int N = A.size();
        assert(N >= 1);

        int M = A[0].size();
        assert(M >= 1);

        assert(abs(rat) > EPSILON);

        for (int j = 0; j < M; ++j)
            A[line][j] /= rat;

        ///cerr << "MULTIPLY :: " << "L" << line << " => " << "L" << line << " x " << 1 / rat << endl;
    }

    void changeLines(Matrix &A, int x, int y, TYPE rat) /// X = X - rat * Y
    {
        if (A.size() == 0)
            return;

        int N = A.size();
        assert(N >= 1);

        int M = A[0].size();
        assert(M >= 1);

        assert(abs(rat) > EPSILON);

        for (int i = 0; i < M; ++i)
            A[x][i] -= rat * A[y][i];

        ///cerr << "ADD :: " << "L" << x << " => " << "L" << x << " + (" << "L" << y << " x " << -rat << ")" << endl;
    }
    ///------------------------------------------------------------

    void rowEchelonForm(Matrix &A, Matrix &B, TYPE &determinant, vector<pair<int,int>> &inversions)
    {
        int N = A.size();
        assert(N >= 1);

        int M = A[0].size();
        assert(M >= 1);

        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())
                    swapLines(B, k, i);

                swapLines(A, k, i);

                determinant *= -1;
                inversions.push_back({k, i});
            }

            determinant *= A[i][j];

            if (abs(A[i][j] - 1) > EPSILON)
            {
                if (B.size())
                    divideLine(B, i, A[i][j]);

                divideLine(A, i, A[i][j]);
            }

            for (int l = i + 1; l < N; ++l)
                if (abs(A[l][j]) > EPSILON)
                {
                    if (B.size())
                        changeLines(B, l, i, A[l][j]);

                    changeLines(A, l, i, A[l][j]);
                }

            i++;
            j++;
        }
    }

    void reducedRowEchelonForm(Matrix &A, Matrix &B)
    {
        int N = A.size();
        assert(N >= 1);

        int M = A[0].size();
        assert(M >= 1);

        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)
                    {
                        changeLines(B, k, i, A[k][j]);
                        changeLines(A, k, i, A[k][j]);
                    }
            }
        }
    }

    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 A = createMatrix(coef);
        Matrix B;

        TYPE determinant = static_cast<TYPE>(1);
        vector<pair<int,int>> invs{};

        rowEchelonForm(A, B, determinant, invs);

        ///return boost::rational_cast<double>(determinant);
        return determinant;
    }

    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 A = createMatrix(coef);
        Matrix B = createIdentity(N);

        TYPE determinant = static_cast<TYPE>(1);
        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]);
                inverse[i][j] = B[i][j];
        }

        ///CHECKING PART-----------------------------
        Matrix tmp = createMatrix(coef);
        Matrix D = matrix_multiplication(tmp, B);
        assert(checkIfIdentity(D));
        Matrix E = matrix_multiplication(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 A = createMatrix(coef);
        Matrix B = createMatrix(N, 1);

        for (int i = 0; i < N; ++i)
            B[i][0] = bs[i];



        TYPE determinant = static_cast<TYPE>(1);
        vector<pair<int,int>> invs{};

        rowEchelonForm(A, B, determinant, invs);
        reducedRowEchelonForm(A, B);

        Vector solutions = createVector(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[i] = B[i][0];
        }

        vector<double> solution(M);

        for (int j = 0; j < M; ++j)
///            solution[j] = boost::rational_cast<double>(solutions[j]);
            solution[j] = solutions[j];

        return solution;
    }

    /*

    template<typename S>
    void LU_decomposition(const vector<vector<S>> &coef)
    {
        assert(coef.size() >= 1);
        assert(coef[0].size() >= 1);

        GaussianElimination::N = coef.size();
        GaussianElimination::M = coef[0].size();

        TYPE C[MAX_VAR + 1][MAX_VAR + 1];

        assert(N == M);

        clear();

        for (int i = 1; i <= N; ++i)
            for (int j = 1; j <= M; ++j)
            {
                A[i][j] = coef[i - 1][j - 1];
                C[i][j] = 0;
            }

        for (int i = 1; i <= N; ++i)
        {
            for (int j = 1; j <= M; ++j)
                cout << A[i][j] << " ";

            cout << endl;
        }

        int i = 1, j = 1;

        while (i <= N && j <= M)
        {
            C[i][j] = 1;
            i++;
            j++;
        }

        i = 1; j = 1;

        while (i <= N && j <= M)
        {
            int k = i;

            while (k <= N && A[k][j] == 0)
                k++;

            if (k == N + 1)
            {
                j++;
                continue;
            }

            if (k != i)
            {
                assert(false);
                swapLines(k, i);
            }

            if (A[i][j] != 1)
            {
                C[i][j] = A[i][j];
                divideLine(i, A[i][j]);
            }

            for (int l = i + 1; l <= N; ++l)
                if (A[l][j] != 0)
                {
                    C[l][j] = A[l][j];
                    changeLines(l, i, A[l][j]);
                }

            i++;
            j++;
        }

        cout << endl;
        cout << "Lower matrix :\n";

        for (int i = 1; i <= N; ++i)
        {
            for (int j = 1; j <= M; ++j)
                cout << C[i][j] << " ";

            cout << endl;
        }

        cout << endl;

        cout << "Upper matrix :\n";

        for (int i = 1; i <= N; ++i)
        {
            for (int j = 1; j <= M; ++j)
                cout << A[i][j] << " ";

            cout << endl;
        }

        assert(N == M);
        matrix_multiplication(C, A, N);
    }
    */
}


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;
}