Cod sursa(job #2497086)

Utilizator Kln1000Ciobanu Bogdan Kln1000 Data 22 noiembrie 2019 00:39:12
Problema Algoritmul lui Gauss Scor 100
Compilator cpp-64 Status done
Runda Arhiva educationala Marime 3.41 kb
#include <fstream>
#include <iomanip>

#include <sys/mman.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <fcntl.h>
#include <unistd.h>

#define NMAX 310
#define EPSILON 1e-8

template <typename T>
using Matrix = T[NMAX][NMAX];
template <typename T>
using Vector = T[NMAX];

class mapping_parser {
public:
    inline mapping_parser() {
        /// default c-tor
    }

    inline mapping_parser(const char* file_name) {
        int fd = open(file_name, O_RDONLY);
        index &= 0;
        fstat(fd, &sb);
        buffer = (char*)mmap(0, sb.st_size, PROT_READ, MAP_PRIVATE | MAP_NORESERVE, fd, 0);
        close(fd);
    }

    inline mapping_parser& operator >> (double& n) {
        aux = 0;
        sign &= 0;
        for (; buffer[index] < '0' or buffer[index] > '9'; ++index);
        sign |= (buffer[index - 1] == '-');
        for (; '0' <= buffer[index] and buffer[index] <= '9'; ++index)
            aux = (aux << 3) + (aux << 1) + buffer[index] - '0';
        aux ^= ((aux ^ -aux ) & -sign);
        n = aux;
        return *this;
    }

    inline mapping_parser& operator >> (int& n) {
        n = 0;
        for (; buffer[index] < '0' or buffer[index] > '9'; ++index);
        for (; '0' <= buffer[index] and buffer[index] <= '9'; ++index)
            n = (n << 3) + (n << 1) + buffer[index] - '0';
        return *this;
    }

    ~mapping_parser() {
        munmap(buffer, sb.st_size);
    }

private:
    int aux;
    struct stat sb;
    int index, sign;
    char* buffer;
};

std::ofstream t("gauss.out");
mapping_parser f("gauss.in");

Matrix<double> A;
Vector<double> b;
Vector<double> x;

char buffer[0x40000];

int main() {
    t.sync_with_stdio(false);
    t.rdbuf()->pubsetbuf(buffer, 0x40000);
    int n, m;

    f >> n >> m;

    auto abs = [](const double &target) {
        return target < 0 ? -target : target;
    };

    auto subtract_line_with_scalar = [&A, &b, m](int dst, int src, double scalar, int start = 0) mutable {
#pragma GCC ivdep
        for (int i = start; i < m; ++i)
            A[dst][i] -= A[src][i] * scalar;

        b[dst] -= b[src] * scalar;
    };

    auto divide_line_with_scalar = [&A, &b, m](int dst, double scalar, int start = 0) mutable {
#pragma GCC ivdep
        for (int i = start; i < m; ++i)
            A[dst][i] /= scalar;

        b[dst] /= scalar;
    };

#pragma GCC ivdep
    for (int i = 0; i < n; ++i) {
#pragma GCC ivdep
        for (int j = 0; j < m; ++j)
            f >> A[i][j];
        f >> b[i];
    }

#pragma GCC ivdep
    for (int i = 0; i < m; ++i) {
        divide_line_with_scalar(i, A[i][i], i);
#pragma GCC ivdep
        for (int j = i + 1; j < n; ++j) {
            subtract_line_with_scalar(j, i, A[j][i], i);
        }
    }

    double sigma;
    int target = n < m ? n - 1 : m - 1;

    for (int i = n - 1, aux = target; i >= 0; --i) {
        sigma = .0;
        for (int j = target; j > aux; --j)
            sigma += x[j] * A[i][j];
        if (abs(sigma + A[i][aux]) < EPSILON && abs(b[i]) > EPSILON) {
            t << "Imposibil\n";
            return 0;
        }
        else if (A[i][aux] > EPSILON) {
            x[aux] = (b[i] - sigma) / A[i][aux];
            --aux;
        }
    }

    t << std::setprecision(9) << std::fixed;
#pragma GCC ivdep
    for (int i = 0; i < m; ++i)
        t << x[i] << " ";
    t.close();
    return 0;
}