cpp/matrix/matrix2.h

#include <iostream>
#include <vector>
#include <cassert>
#include <string>
#include <type_traits>

/*
 * TODO: assert and invert
 * порядок методов
 * проверка на N != M
 */

class BigInteger{
private:
    static const int BASE = 10;
    bool isPositive; //Zero is positive
    std::vector<short> digits;
    void normilize();
    void popBackZeros();
    void fixMinusZero();
    void reverseDigits();
    void addDigit(short);
    friend std::istream& operator>>(std::istream&, BigInteger&);
public:
    BigInteger(int);
    BigInteger();
    BigInteger& operator+=(const BigInteger&);
    BigInteger operator+() const;
    BigInteger& operator++();
    BigInteger operator++(int);
    BigInteger& operator-=(const BigInteger&);
    BigInteger operator-() const;
    BigInteger& operator--();
    BigInteger operator--(int);
    BigInteger& operator*=(const BigInteger&);
    BigInteger& operator/=(const BigInteger&);
    BigInteger& operator%=(const BigInteger&);
    bool operator<(const BigInteger&) const;
    bool operator>(const BigInteger&) const;
    bool operator==(const BigInteger&) const;
    bool operator!=(const BigInteger&) const;
    bool operator<=(const BigInteger&) const;
    bool operator>=(const BigInteger&) const;
    explicit operator int() const;
    explicit operator bool() const;
    std::string toString() const;
};


// if something weird in istream => return interpreted, istream is on the first weird symbol
// also interpret number = sign as 0
std::istream& operator>>(std::istream& in, BigInteger& num) {
    char c;
    bool signDefined = false;
    num.isPositive = true;
    num.digits.clear();
    while(in.get(c)) {
        if (c == '-' && !signDefined) {
            num.isPositive = false;
            signDefined = true;
        } else if (c == '+' && !signDefined) {
            num.isPositive = true;
            signDefined = true;
        } else if (c >= '0' && c <= '9') {
            num.digits.push_back(c - '0');
            signDefined = true;
        } else {
            if (c != ' ' && c != '\n')
                in.putback(c);
            break;
        }
    }
    if (num.digits.size() == 0)
        num.digits.push_back(0);

    num.reverseDigits();
    num.normilize();
    return in;
}

std::ostream& operator<<(std::ostream& out, const BigInteger& num) {
    out << num.toString();
    return out;
}

void BigInteger::normilize() {
    popBackZeros();
    fixMinusZero();
}

void BigInteger::popBackZeros() {
    while(digits.size() > 1 && digits.back() == 0) {
        digits.pop_back();
    }
}

void BigInteger::fixMinusZero() {
    isPositive = !isPositive;
    if (*this == BigInteger(0))
        isPositive = !isPositive;
    isPositive = !isPositive;
}

void BigInteger::reverseDigits() {
    int size = digits.size();
    for (int i = 0; i < size / 2; ++i) {
        std::swap(digits[i], digits[size-1-i]);
    }
}

void BigInteger::addDigit(short digit) {
    reverseDigits();
    digits.push_back(digit);
    reverseDigits();
}

BigInteger::BigInteger(int value) {
    isPositive = (value >= 0);
    if (!isPositive)
        value = -value;
    if (!value) digits.push_back(0);
    while (value) {
        digits.push_back(value % BASE);
        value /= BASE;
    }
}

BigInteger::BigInteger() : BigInteger(0) {
}

BigInteger& BigInteger::operator+=(const BigInteger& num) {
    if (this == &num) {
        return *this += BigInteger(num);
    }
    if (num.isPositive == isPositive) {
        for (unsigned int i = 0; i < std::max(num.digits.size(), digits.size()); ++i) {
            short add = (i < num.digits.size() ? num.digits[i] : 0);
            if (digits.size() == i)
                digits.push_back(0);
            if (digits[i] + add >= BASE) {
                if (digits.size() == i + 1)
                    digits.push_back(0);
                ++digits[i + 1];
            } else if (i >= num.digits.size()) {
                digits[i] = (digits[i] + add) % BASE;
                break;
            }
            digits[i] = (digits[i] + add) % BASE;
        }
    } else {
        isPositive = !isPositive;
        *this -= num;
        isPositive = !isPositive;
    }
    normilize();
    return *this;
}

BigInteger operator+(const BigInteger& num1, const BigInteger& num2) {
    BigInteger copy = num1;
    return copy += num2;
}

BigInteger BigInteger::operator+() const {
    return *this;
}

BigInteger& BigInteger::operator++() {
    return (*this += 1);
}

BigInteger BigInteger::operator++(int) {
    BigInteger copy = *this;
    ++(*this);
    return copy;
}

BigInteger& BigInteger::operator-=(const BigInteger& num) {
    if (this == &num) {
        return *this -= BigInteger(num);
    }
    if (num.isPositive == isPositive) {
        bool needChange = isPositive ^ (num <= *this);
        const BigInteger& a = needChange ? num : *this;
        const BigInteger& b = needChange ? *this : num;
        bool need = false;
        for (unsigned int i = 0; i < a.digits.size(); ++i) {
            bool needNext = (a.digits[i] - need - (i < b.digits.size() ? b.digits[i] : 0) < 0);
            if(i == digits.size()) digits.push_back(0);
            digits[i] = (a.digits[i] - need - (i < b.digits.size() ? b.digits[i] : 0) + BASE) % BASE;
            need = needNext;
        }
        isPositive ^= needChange;
    } else {
        isPositive = !isPositive;
        *this += num;
        isPositive = !isPositive;
    }
    normilize();
    return *this;
}

BigInteger operator-(const BigInteger& num1, const BigInteger& num2) {
    BigInteger copy = num1;
    return copy -= num2;
}

BigInteger BigInteger::operator-() const {
    return BigInteger(0) -= *this;
}

BigInteger& BigInteger::operator--() {
    return (*this -= 1);
}

BigInteger BigInteger::operator--(int) {
    BigInteger copy = *this;
    --(*this);
    return copy;
}

BigInteger& BigInteger::operator*=(const BigInteger& num) {
    BigInteger ans;
    ans.digits.resize(digits.size() + num.digits.size());
    for(unsigned int i = 0; i < digits.size(); ++i) {
        int add = 0;
        for (unsigned int j = 0; j <= num.digits.size(); ++j) {
            short numDigit = (j == num.digits.size() ? 0 : num.digits[j]);
            short addNext = (ans.digits[i + j] + digits[i] * numDigit + add) / BASE;
            ans.digits[i + j] = (ans.digits[i + j] + digits[i] * numDigit + add) % BASE;
            add = addNext;
        }
    }
    ans.isPositive = !(isPositive ^ num.isPositive);
    ans.normilize();
    return *this = ans;
}

BigInteger operator*(const BigInteger& num1, const BigInteger& num2) {
    BigInteger copy = num1;
    return copy *= num2;
}

BigInteger& BigInteger::operator/=(const BigInteger& num) {
    if(num == 0)
        return *this;
    BigInteger ans, val(digits.back());
    ans.isPositive = !(isPositive ^ num.isPositive);
    int ptr = digits.size();
    ptr -= 2;
    val.isPositive = num.isPositive;
    while (ptr >= 0 && (num.isPositive ? val < num : num < val)) {
        val.addDigit(digits[ptr]);
        --ptr;
    }
    val.isPositive = true;
    for (; ptr >= -1; --ptr) {
        int number = 0;
        while (val >= 0) {
            num.isPositive ? val -= num : val += num;
            ++number;
        }
        num.isPositive ? val += num : val -= num;
        --number;

        ans.digits.push_back(number);
        if(ptr != -1)
            val.addDigit(digits[ptr]);
    }
    ans.reverseDigits();
    ans.normilize();
    return *this = ans;
}

BigInteger operator/(const BigInteger& num1, const BigInteger& num2) {
    BigInteger copy = num1;
    return copy /= num2;
}

BigInteger& BigInteger::operator%=(const BigInteger& num) {
    if(num == 0)
        return *this;
    BigInteger copy = (*this) / num;
    return *this -= (copy *= num);
}

BigInteger operator%(const BigInteger& num1, const BigInteger& num2) {
    BigInteger copy = num1;
    return copy %= num2;
}

bool BigInteger::operator<(const BigInteger& num) const {
    if (isPositive != num.isPositive)
        return num.isPositive;
    for (int j = std::max(static_cast<int>(digits.size()), static_cast<int>(num.digits.size())) - 1; j >= 0; --j) {
        unsigned int i = static_cast<unsigned int>(j);
        if ((i < digits.size() ? digits[i] : 0) < (i < num.digits.size() ? num.digits[i] : 0)) return isPositive;
        if ((i < digits.size() ? digits[i] : 0) > (i < num.digits.size() ? num.digits[i] : 0)) return !isPositive;
    }
    return false;
}

bool BigInteger::operator>(const BigInteger& num) const {
    return num < *this;
}

bool BigInteger::operator==(const BigInteger& num) const {
    return !(num < *this || *this < num);
}

bool BigInteger::operator!=(const BigInteger& num) const {
    return num < *this || *this < num;
}

bool BigInteger::operator<=(const BigInteger& num) const {
    return !(num < *this);
}

bool BigInteger::operator>=(const BigInteger& num) const {
    return !(*this < num);
}

std::string BigInteger::toString() const {
    std::string num = "";
    if (!isPositive) num += '-';
    for (int i = digits.size() - 1; i >= 0; --i) {
        num += ('0' + digits[i]);
    }
    return num;
}

BigInteger::operator int() const {
    int value = 0;
    for (int i = digits.size() - 1; i >= 0; --i)
        value = value * 10 + digits[i];
    if(!isPositive)
        value = -value;
    return value;
}

BigInteger::operator bool() const {
    return int(*this);
}


class Rational{
private:
    BigInteger numerator, denominator;
    void normilize();
    static BigInteger gcd(const BigInteger&, const BigInteger&);
    static void reverse(std::string&);
public:
    Rational();
    Rational(const BigInteger&);
    Rational(const BigInteger& numerator, const BigInteger& denominator): numerator(numerator), denominator(denominator) {}
    Rational(int);
    Rational& operator+=(const Rational&);
    Rational operator+() const;
    Rational& operator-=(const Rational&);
    Rational operator-() const;
    Rational& operator*=(const Rational&);
    Rational& operator/=(const Rational&);
    bool operator<(const Rational&) const;
    bool operator>(const Rational&) const;
    bool operator==(const Rational&) const;
    bool operator!=(const Rational&) const;
    bool operator<=(const Rational&) const;
    bool operator>=(const Rational&) const;
    explicit operator double() const;
    std::string asDecimal(size_t) const;
    std::string toString() const;
};

BigInteger Rational::gcd(const BigInteger& a, const BigInteger& b) {
    if (b == 0)
        return a;
    return gcd(b, a % b);
}

void Rational::normilize() {
    BigInteger _gcd = gcd(numerator, denominator);
    numerator /= _gcd;
    denominator /= _gcd;
    if (denominator < 0) {
        numerator *= -1;
        denominator *= -1;
    }
}

void Rational::reverse(std::string& s) {
    for(size_t i = 0; i < s.size()/2; ++i) {
        std::swap(s[i], s[s.size() - 1 - i]);
    }
}

Rational::Rational() : Rational(0) {
}

Rational::Rational(const BigInteger& val) {
    numerator = val;
    denominator = 1;
}

Rational::Rational(int val) {
    numerator = val;
    denominator = 1;
}

Rational& Rational::operator+=(const Rational& frac) {
    if (this == &frac)
        return *this += Rational(frac);
    numerator *= frac.denominator;
    numerator += denominator * frac.numerator;
    denominator *= frac.denominator;
    normilize();
    return *this;
}

Rational operator+(const Rational& frac1, const Rational& frac2) {
    Rational copy = frac1;
    return copy += frac2;
}

Rational Rational::operator+() const {
    return *this;
}

Rational& Rational::operator-=(const Rational& frac) {
    if (this == &frac)
        return *this -= Rational(frac);
    numerator *= frac.denominator;
    numerator -= denominator * frac.numerator;
    denominator *= frac.denominator;
    normilize();
    return *this;
}

Rational operator-(const Rational& frac1, const Rational& frac2) {
    Rational copy = frac1;
    return copy -= frac2;
}

Rational Rational::operator-() const {
    return Rational(0) -= *this;
}

Rational& Rational::operator*=(const Rational& frac) {
    if (this == &frac)
        return *this *= Rational(frac);
    numerator *= frac.numerator;
    denominator *= frac.denominator;
    normilize();
    return *this;
}

Rational operator*(const Rational& frac1, const Rational& frac2) {
    Rational copy = frac1;
    return copy *= frac2;
}

Rational& Rational::operator/=(const Rational& frac) {
    if (frac.numerator == 0)
        return *this;
    if (this == &frac)
        return *this /= Rational(frac);
    numerator *= frac.denominator;
    denominator *= frac.numerator;
    normilize();
    return *this;
}

Rational operator/(const Rational& frac1, const Rational& frac2) {
    Rational copy = frac1;
    return copy /= frac2;
}

bool Rational::operator<(const Rational& frac) const {
    return numerator * frac.denominator < frac.numerator * denominator;
}

bool Rational::operator>(const Rational& frac) const {
    return frac < *this;
}

bool Rational::operator==(const Rational& frac) const {
    return !(frac < *this || *this < frac);
}

bool Rational::operator!=(const Rational& frac) const {
    return frac < *this || *this < frac;
}

bool Rational::operator<=(const Rational& frac) const {
    return !(frac < *this);
}

bool Rational::operator>=(const Rational& frac) const {
    return !(*this < frac);
}

Rational::operator double() const {
  return std::stod(asDecimal(200));
}

std::string Rational::asDecimal(size_t precision = 0u) const {
    BigInteger ten = 1;
    for (size_t i = 0; i < precision; ++i) {
        ten *= 10;
    }
    BigInteger val = (numerator < 0 ? -1 : 1) * numerator * ten / denominator;
    std::string integer = (val / ten).toString();
    std::string frac = (val % ten).toString();
    reverse(frac);
    while(frac.size() < precision) frac += '0';
    frac += '.';
    reverse(frac);
    if(precision == 0u) frac = "";
    integer = (numerator < 0 ? "-" : "") + integer;
    return integer + frac;
}

std::string Rational::toString() const {
    std::string num = numerator.toString();
    if (denominator != 1) {
        num += "/" + denominator.toString();
    }
    return num;
}

std::ostream& operator<<(std::ostream& out, const Rational& value) {
    out << value.asDecimal(5);
    return out;
}

std::istream& operator>>(std::istream& in, Rational& value) {
    BigInteger numerator(0), denominator(0);
    bool was = false;
    char c;
    while (in.get(c)) {
        if (c == ' ' || c == '\n')
            break;
        if (c == '/') {
            was = true;
        } else {
            if (!was) {
                numerator *= 10;
                numerator += (c - '0');
            } else {
                denominator *= 10;
                denominator += (c - '0');
            }
        }
    }
    if (!was) {
        denominator = BigInteger(1);
    }
    value = Rational(numerator, denominator);
    return in;
}

bool isPrime(int N) {
    for (long long i = 2; i * i <= N; ++i) {
        if (N % i == 0)
            return false;
    }
    return true;
}

int fastPow(int a, int b, int module) {
    if (b == 0)
        return 1;
    if (b % 2 == 0) {
        long long val = fastPow(a, b / 2, module);
        return (val * val) % module;
    } else {
        return (1ll * fastPow(a, b - 1, module) * a) % module;
    }
}

int inverseValue(int value, int module) {
    return fastPow(value, module - 2, module);
}

template<int N>
class Finite {
public:
    Finite<N>(int value = 0): value((value % N + N) % N), isPrimeN(isPrime(N)) {}
    Finite<N>& operator+=(const Finite<N>& number) {
        value = (value + number.value) % N;
        return *this;
    }
    Finite<N> operator+(const Finite<N>& number) const {
        Finite<N> copy = *this;
        return copy += number;
    }
    Finite<N>& operator-=(const Finite<N>& number) {
        value = (value - number.value + N) % N;
        return *this;
    }
    Finite<N> operator-(const Finite<N>& number) const {
        Finite<N> copy = *this;
        return copy -= number;
    }
    Finite<N> operator-() const {
        return Finite<N>(0 - value);
    }
    Finite<N>& operator*=(const Finite<N>& number) {
        value = (1ll * value * number.value) % N;
        return *this;
    }
    Finite<N> operator*(const Finite<N>& number) const {
        Finite<N> copy = *this;
        return copy *= number;
    }
    Finite<N>& operator/=(const Finite<N>& number) {
        assert(isPrimeN);
        value = (1ll * value * inverseValue(number.value, N)) % N;
        return *this;
    }
    Finite<N> operator/(const Finite<N>& number) const {
        Finite<N> copy = *this;
        return copy /= number;
    }
    bool operator==(const Finite<N>& number) const {
        return value == number.value;
    }
    bool operator!=(const Finite<N>& number) const {
        return value != number.value;
    }
    Finite<N>& operator++() {
        return (*this += 1);
    }
    Finite<N> operator++(int) {
        Finite<N> copy;
        ++(*this);
        return copy;
    }
    template<int U>
    friend std::ostream& operator<<(std::ostream&, const Finite<U>&);
private:
    int value;
    bool isPrimeN;
};

template<int N>
std::ostream& operator<<(std::ostream& out, const Finite<N>& value){
    out << value.value;
    return out;
}

template<unsigned M, unsigned N, typename Field = Rational>
class Matrix {
public:
    Matrix<M, N, Field>();

    template<typename T>
    Matrix(std::initializer_list<std::initializer_list<T>> a) {
      for (auto i: a) {
        table.push_back({});
        for (auto j: i) {
          table.back().push_back(Field(j));
        }
      }
    }
    Matrix<M, N, Field>(std::vector<std::vector<Field>>&);
    Matrix<M, N, Field>(std::vector<std::vector<int>>&);
    template<unsigned M1, unsigned N1, typename Field1>
    bool operator==(const Matrix<M1, N1, Field1>&) const;
    template<unsigned M1, unsigned N1, typename Field1>
    bool operator!=(const Matrix<M1, N1, Field1>&) const;
    Matrix<M, N, Field>& operator+=(const Matrix<M, N, Field>&);
    Matrix<M, N, Field> operator+(const Matrix<M, N, Field>&) const;
    Matrix<M, N, Field>& operator-=(const Matrix<M, N, Field>&);
    Matrix<M, N, Field> operator-(const Matrix<M, N, Field>&) const;
    Matrix<M, N, Field> operator*(const Field&) const;
    Matrix<M, N, Field>& operator*=(const Matrix<M, N, Field>&);
    Field det() const;
    Matrix<N, M, Field> transposed() const;
    unsigned rank() const;
    Field trace() const;
    void invert();
    Matrix<N, M, Field> inverted() const;
    std::vector<Field> getRow(unsigned) const;
    std::vector<Field> getColumn(unsigned) const;
    std::vector<Field>& operator[](size_t index) {
        return table[index];
    }
    const std::vector<Field>& operator[](size_t index) const {
        return table[index];
    }
    template<unsigned N1, unsigned M1, typename Field1>
    friend class Matrix;
private:
    std::vector<std::vector<Field>> table;
    void swapStrings(unsigned, unsigned);
    void divideString(unsigned, const Field);
    void subtractStrings(unsigned, unsigned, const Field);
    Field makeSimple();
};

template<unsigned int M, unsigned int N, typename Field>
Matrix<M, N, Field>::Matrix() {
    if (M != N)
        assert(0);
    table.resize(M, std::vector<Field>(N));
    for (unsigned i = 0; i < N; ++i) {
        table[i][i] = Field(1);
    }
}

template<unsigned int M, unsigned int N, typename Field>
Matrix<M, N, Field>::Matrix(std::vector<std::vector<Field>>& matrix): table(matrix) {}

template<unsigned int M, unsigned int N, typename Field>
Matrix<M, N, Field>::Matrix(std::vector<std::vector<int>>& matrix) {
    table.resize(M, std::vector<Field>(N));
    for (unsigned i = 0; i < M; ++i) {
        for (unsigned j = 0; j < N; ++j){
            table[i][j] = Field(matrix[i][j]);
        }
    }
}

template<unsigned int M, unsigned int N, typename Field>
template<unsigned int M1, unsigned int N1, typename Field1>
bool Matrix<M, N, Field>::operator==(const Matrix<M1, N1, Field1>& matrix) const {
    if (!std::is_same<Matrix<M, N, Field>, Matrix<M1, N1, Field1>>::value)
        return false;
    for (unsigned i = 0; i < M; ++i) {
        for (unsigned j = 0; j < N; ++j) {
            if (table[i][j] != matrix[i][j])
                return false;
        }
    }
    return true;
}

template<unsigned int M, unsigned int N, typename Field>
template<unsigned int M1, unsigned int N1, typename Field1>
bool Matrix<M, N, Field>::operator!=(const Matrix<M1, N1, Field1>& matrix) const {
    return !(*this == matrix);
}

template<unsigned int M, unsigned int N, typename Field>
Matrix<M, N, Field>& Matrix<M, N, Field>::operator+=(const Matrix<M, N, Field>& matrix) {
    for (unsigned i = 0; i < M; ++i) {
        for (unsigned j = 0; j < N; ++j) {
            table[i][j] += matrix[i][j];
        }
    }
    return *this;
}

template<unsigned int M, unsigned int N, typename Field>
Matrix<M, N, Field> Matrix<M, N, Field>::operator+(const Matrix<M, N, Field>& matrix) const {
    Matrix<M, N, Field> copy = *this;
    return copy += matrix;
}

template<unsigned int M, unsigned int N, typename Field>
Matrix<M, N, Field>& Matrix<M, N, Field>::operator-=(const Matrix<M, N, Field>& matrix) {
    for (unsigned i = 0; i < M; ++i) {
        for (unsigned j = 0; j < N; ++j) {
            table[i][j] -= matrix[i][j];
        }
    }
    return *this;
}

template<unsigned int M, unsigned int N, typename Field>
Matrix<M, N, Field> Matrix<M, N, Field>::operator-(const Matrix<M, N, Field>& matrix) const {
    Matrix<M, N, Field> copy = *this;
    return copy -= matrix;
}

template<unsigned int M, unsigned int N, typename Field>
Matrix<M, N, Field> Matrix<M, N, Field>::operator*(const Field& multiplier) const {
    Matrix<M, N, Field> copy = *this;
    for (unsigned i = 0; i < M; ++i) {
        for (unsigned j = 0; j < N; ++j) {
            copy[i][j] *= multiplier;
        }
    }
    return copy;
}

template<unsigned int M, unsigned int N, typename Field>
Matrix<M, N, Field> operator*(const Field& multiplier, const Matrix<M, N, Field> matrix) {
    return matrix * multiplier;
}

template<unsigned int M, unsigned int N, unsigned int K, typename Field>
Matrix<M, K, Field> operator*(const Matrix<M, N, Field>& matrix1, const Matrix<N, K, Field>& matrix2) {
    std::vector<std::vector<Field>> resultTable(M, std::vector<Field>(K));
    for (unsigned i = 0; i < M; ++i) {
        for (unsigned j = 0; j < K; ++j) {
            for (unsigned k = 0; k < N; ++k) {
                resultTable[i][j] += matrix1[i][k] * matrix2[k][j];
            }
        }
    }
    return Matrix<M, K, Field>(resultTable);
}

template<unsigned int M, unsigned int N, typename Field>
Matrix<M, N, Field>& Matrix<M, N, Field>::operator*=(const Matrix<M, N, Field>& matrix) {
    if (M != N)
        assert(0);
    return (*this = (*this) * matrix);
}

template<unsigned int M, unsigned int N, typename Field>
Field Matrix<M, N, Field>::det() const {
    if (N != M)
        assert(0);
    Matrix<M, N, Field> copy = *this;
    return copy.makeSimple();
}

template<unsigned int M, unsigned int N, typename Field>
Matrix<N, M, Field> Matrix<M, N, Field>::transposed() const {
    std::vector<std::vector<Field>> result(N, std::vector<Field>(M));
    for (unsigned i = 0; i < N; ++i) {
        for (unsigned j = 0; j < M; ++j) {
            result[i][j] = table[j][i];
        }
    }
    return Matrix<N, M, Field>(result);
}

template<unsigned int M, unsigned int N, typename Field>
unsigned Matrix<M, N, Field>::rank() const {
    Matrix<M, N, Field> copy = *this;
    copy.makeSimple();
    unsigned rank = 0;
    for (unsigned i = 0; i < M; ++i) {
        bool empty = true;
        for (unsigned j = 0; j < N; ++j) {
            if (copy[i][j] != 0)
                empty = false;
        }
        rank += !empty;
    }
    return rank;
}

template<unsigned int M, unsigned int N, typename Field>
Field Matrix<M, N, Field>::trace() const {
    if (N != M)
        assert(0);
    Field result(0);
    for (unsigned i = 0; i < N; ++i) {
        result += table[i][i];
    }
    return result;
}

template<unsigned int M, unsigned int N, typename Field>
void Matrix<M, N, Field>::invert() {
    if (M != N)
        assert(0);
    if (N > 5){
        std::cerr << "was:\n";
        for(unsigned i = 0; i < 10; ++i) {
            for(unsigned j = 0; j < 20; ++j){
                std::cerr << table[i][j] << " ";
            }
            std::cerr << std::endl;
        }
    }
    std::vector<std::vector<Field>> complexTable(N, std::vector<Field>(2 * N));
    for (unsigned i = 0; i < N; ++i) {
        for (unsigned j = 0; j < N; ++j) {
            complexTable[i][j] = table[i][j];
        }
    }
    for (unsigned i = 0; i < N; ++i)
        complexTable[i][i + N] = Field(1);
    Matrix<N, 2 * N, Field> complexMatrix(complexTable);
    complexMatrix.makeSimple();
    std::vector<std::vector<Field>> result(N, std::vector<Field>(N));
    for (unsigned i = 0; i < M; ++i) {
        for (unsigned j = 0; j < N; ++j) {
            result[i][j] = complexMatrix[i][j + N];
        }
    }
    if (N > 100){
        std::cerr << "now:\n";
        for(unsigned i = 0; i < M; ++i) {
            for(unsigned j = 0; j < N; ++j){
                std::cerr << result[i][j] << " ";
            }
            std::cerr << std::endl;
        }
    }
    *this = Matrix<N, M, Field>(result);
}

template<unsigned int M, unsigned int N, typename Field>
Matrix<N, M, Field> Matrix<M, N, Field>::inverted() const {
    Matrix<N, M, Field> matrix = *this;
    matrix.invert();
    return matrix;
}

template<unsigned int M, unsigned int N, typename Field>
std::vector<Field> Matrix<M, N, Field>::getRow(unsigned int index) const {
    return table[index];
}

template<unsigned int M, unsigned int N, typename Field>
std::vector<Field> Matrix<M, N, Field>::getColumn(unsigned int index) const {
    std::vector<Field> result(M);
    for (unsigned i = 0; i < M; ++i) {
        result[i] = table[i][index];
    }
    return result;
}

template<unsigned int M, unsigned int N, typename Field>
void Matrix<M, N, Field>::swapStrings(unsigned index1, unsigned index2) {
    std::swap(table[index1], table[index2]);
}

template<unsigned int M, unsigned int N, typename Field>
void Matrix<M, N, Field>::divideString(unsigned index, const Field divisor) {
    for (unsigned i = 0; i < N; ++i) {
        table[index][i] /= divisor;
    }
}

template<unsigned int M, unsigned int N, typename Field>
void Matrix<M, N, Field>::subtractStrings(unsigned indexMinuend, unsigned indexSubtrahend, const Field multiplier) {
    for (unsigned i = 0; i < N; ++i) {
        table[indexMinuend][i] -= table[indexSubtrahend][i] * multiplier;
    }
}

template<unsigned int M, unsigned int N, typename Field>
Field Matrix<M, N, Field>::makeSimple() {
    Field determinant(1);
    unsigned string = 0;
    for (unsigned i = 0; i < N; ++i) {
        unsigned found = M + 1;
        for (unsigned j = string; j < M; ++j) {
            if (table[j][i] != Field(0)) {
                found = j;
                break;
            }
        }
        if (found == M + 1){
            determinant = 0;
            continue;
        }
        swapStrings(string, found);
        determinant *= table[string][i];
        divideString(string, table[string][i]);
        for (unsigned j = 0; j < M; ++j) {
            if (j != string)
                subtractStrings(j, string, table[j][i]);
        }
        ++string;
    }
    return determinant;
}

template <unsigned N, typename Field = Rational>
using SquareMatrix = Matrix<N, N, Field>;