#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace std; auto getTime() { return (long long)std::chrono::steady_clock::now().time_since_epoch().count(); //return 0LL; } bool isDivided = false; namespace FastFourierTransform { typedef long double ld; typedef std::complex cld; typedef cld* ComplexPolynom; const ld PI = acosl(-1); std::vector Muliplication(const std::vector& a, const std::vector& b, std::vector& result); ComplexPolynom ToComplex(const std::vector& a, size_t m); void FastFourierTransform_(ComplexPolynom a, size_t n, bool invert); std::vector Muliplication(const std::vector& a, const std::vector& b, std::vector& result) { size_t n = std::max(a.size(), b.size()); size_t m = 1; while (m < n) m <<= 1; m <<= 1; ComplexPolynom ca = ToComplex(a, m); ComplexPolynom cb = ToComplex(b, m); FastFourierTransform_(ca, m, false); FastFourierTransform_(cb, m, false); for (size_t i = 0; i < m; ++i) { ca[i] *= cb[i]; } FastFourierTransform_(ca, m, true); std::vector v(m); for (size_t i = 0; i < m; ++i) { v[i] = round(ca[i].real()); } delete[] ca; delete[] cb; v.push_back(0); v.push_back(0); v.push_back(0); for (size_t i = 0; i < m; ++i) { v[i + 1] += v[i] / 1000; v[i] %= 1000; assert(v[i] >= 0); while (v[i] < 0) v[i] += 10000000, v[i + 1] -= 1; } result.resize(v.size()); std::copy(v.begin(), v.end(), result.begin()); return result; } ComplexPolynom ToComplex(const std::vector& a, size_t m) { ComplexPolynom res = new cld[m]; for (size_t i = 0; i < a.size(); ++i) res[i] = a[i]; for (size_t i = a.size(); i < m; ++i) res[i] = 0; return res; } size_t rev(size_t num, size_t lg_n) { int res = 0; for (size_t i = 0; i < lg_n; ++i) if (num & (1 << i)) res |= 1 << (lg_n - 1 - i); return res; } void FastFourierTransform_(ComplexPolynom a, size_t n, bool invert) { size_t lg_n = 0; while ((1u << lg_n) < n) ++lg_n; for (size_t i = 0; i < n; ++i) if (i < rev(i, lg_n)) swap(a[i], a[rev(i, lg_n)]); for (size_t len = 2; len <= n; len <<= 1) { double ang = 2 * PI / len * (invert ? -1 : 1); cld wlen(cos(ang), sin(ang)); for (size_t i = 0; i < n; i += len) { cld w(1); for (auto it1 = a + i, it2 = a + i + len / 2; it1 < a + i + len / 2; ++it1, ++it2) { cld u = *it1, v = *it2 * w; *it1 = u + v; *it2 = u - v; w *= wlen; } /* for (size_t j = 0; j < len / 2; ++j) { cld u = a[i + j], v = a[i + j + len / 2] * w; a[i + j] = u + v; a[i + j + len / 2] = u - v; w *= wlen; } */ } } if (invert) for (size_t i = 0; i < n; ++i) a[i] /= n; } }; // namespace FastFourierTransform class BigInteger; BigInteger operator-(const BigInteger& a, const BigInteger& b); BigInteger operator*(const BigInteger& a, const BigInteger& b); BigInteger operator/(const BigInteger& a, const BigInteger& b); BigInteger operator+(const BigInteger& a, const BigInteger& b); bool operator>(const BigInteger& a, const BigInteger& b); bool operator<=(const BigInteger& a, const BigInteger& b); bool operator>=(const BigInteger& a, const BigInteger& b); bool operator==(const BigInteger& a, const BigInteger& b); bool operator!=(const BigInteger& a, const BigInteger& b); BigInteger operator%(const BigInteger& a, const BigInteger& b); long long mul = 0; long long DIV = 0; long long PPLUS = 0; class BigInteger { public: BigInteger() {} BigInteger(int64_t x) { if (x < 0) { is_negative_ = true, x *= -1; } t_[0] = x % BASE; x /= BASE; while (x > 0) { t_.push_back(x % BASE); x /= BASE; } Normalize(); } BigInteger(const BigInteger& x) : t_(x.t_), is_negative_(x.is_negative_) {} void swap(BigInteger& x) { std::swap(t_, x.t_); std::swap(is_negative_, x.is_negative_); } BigInteger& operator=(const BigInteger& x) { BigInteger tmp = x; swap(tmp); return *this; } friend bool operator<(const BigInteger& a, const BigInteger& b) { if (a.is_negative_ ^ b.is_negative_) { return a.is_negative_; } if (a.t_.size() != b.t_.size()) { return (a.t_.size() < b.t_.size()) ^ a.is_negative_; } for (size_t i = a.t_.size(); i > 0; --i) { if (a.t_[i - 1] != b.t_[i - 1]) { return (a.t_[i - 1] < b.t_[i - 1]) ^ a.is_negative_; } } return a.is_negative_; } BigInteger& operator++() { *this += 1; Normalize(); return *this; } BigInteger operator++(int) { BigInteger tmp = *this; ++(*this); return tmp; } void addingDifferentSigns(const BigInteger& x) { int16_t carry[2] = {0, 0}; // Вычитаем, если *this < x, то последнее число будет отрицательным и // потом это пофиксим for (size_t i = 0; i < x.t_.size() || (carry[i & 1] && i < t_.size()); ++i) { if (i < x.t_.size()) { t_[i] -= x.t_[i]; } t_[i] += carry[i & 1]; if (t_[i] < 0 && i + 1 < t_.size()) { t_[i] += BASE; carry[(i & 1) ^ 1] = -1; } else { carry[(i & 1) ^ 1] = 0; } } // Убераем нули в начале и начинаем фиксить отрицательное число Normalize(); if (t_.back() < 0) { t_.pop_back(); for (size_t i = 0; i < t_.size(); ++i) { t_[i] = BASE - 1 - t_[i]; } ++t_[0]; while (t_.size() > 1 && t_.back() == 0) { t_.pop_back(); } is_negative_ ^= true; } // Когда мы пофиксили отрицательные числа, у нас в ячейках могли // образоваться числа >= BASE, фиксим это carry[0] = 0; carry[1] = 0; if (t_[0] >= BASE) { carry[1] = 1; t_[0] -= BASE; } for (size_t i = 1; i < t_.size() && carry[i & 1]; ++i) { t_[i] += carry[i & 1]; carry[i & 1] = 0; if (t_[i] >= BASE) { carry[(i & 1) ^ 1] = 1; t_[i] -= BASE; } } // У нас может произойти перенос в ячейку, которой не существует, добавим // ее) if (carry[t_.size() & 1]) { t_.push_back(1); } } BigInteger& operator+=(const BigInteger& x) { PPLUS -= getTime(); t_.resize(std::max(t_.size(), x.t_.size()) + 1); // Если числа одинаковых знаков - просто сложение if (!is_negative_ ^ x.is_negative_) { int16_t carry[2] = {0, 0}; for (size_t i = 0; i < x.t_.size() || carry[i & 1]; ++i) { if (i < x.t_.size()) { t_[i] += x.t_[i]; } t_[i] += carry[i & 1]; if (t_[i] >= BASE) { t_[i] -= BASE; carry[(i & 1) ^ 1] = 1; } else { carry[(i & 1) ^ 1] = 0; } } Normalize(); } else { addingDifferentSigns(x); } PPLUS += getTime(); return *this; } BigInteger& operator-=(const BigInteger& x) { // a - b = -(-a + b) is_negative_ ^= 1; *this += x; is_negative_ ^= 1; Normalize(); return *this; } BigInteger& operator*=(int64_t x) { // Это умножение на маленькое число, работает за O(n) и используется в // делении cout << "*= int" << endl; if (x < 0) { is_negative_ = !is_negative_; x *= -1; } if (x >= BASE) { cout << "!!!!!!!!!!!" << x << endl; return *this *= BigInteger(x); } t_.resize(t_.size() + 3); int64_t carry[2] = {0, 0}; for (size_t i = 0; i < t_.size(); ++i) { carry[!(i & 1)] = ((long long)(t_[i]) * x + carry[i & 1]) / BASE; t_[i] = ((long long)(t_[i]) * x + carry[i & 1]) % BASE; } while (t_.size() > 1 && t_.back() == 0) t_.pop_back(); Normalize(); return *this; } BigInteger& operator*=(const BigInteger& x) { // Если число на которое нам надо умножить небольшое, выгоднее умножать за // квадрат assert(!isDivided); cout << "*= Big" << endl; if (std::min(t_.size(), x.t_.size()) < 16 || true) { mul -= getTime(); SquareMultiplication(x); mul += getTime(); return *this; } assert(0); /* mul -= getTime(); FastFourierTransform::Muliplication(t_, x.t_, t_); is_negative_ ^= x.is_negative_; Normalize(); mul += getTime(); return *this; */ } // Деление работает за O(N^2), мы вычисляем все значащие биты в делении, // которых O(n), с помощью бинпоиска(константа, так как log 10000), и // проверкой в бинпоиске за O(n), потому что используется умножение за O(n), и // в итоге O(N^2) BigInteger& operator/=(const BigInteger& x) { cout << "/= " << endl; isDivided = true; DIV -= getTime(); bool is_negative = is_negative_ ^ x.is_negative_; is_negative_ = false; std::vector res(t_.size() - std::min(x.t_.size(), t_.size()) + 4); int64_t d = 1; BigInteger xx(x); if (xx < 0) { d = -1; xx *= -1; } BigInteger cur(0); for (int i = t_.size(); i > 0; --i) { int64_t l = 0, r = BASE - 1, m; cur.MuliplicationDegree10(1); cur += t_[i - 1]; if (cur >= xx) { if (cur.t_.size() == xx.t_.size()) { cout << "!!!" << cur.t_.size() << endl; l = cur.t_.back() / (xx.t_.back() + 1); r = cur.t_.back() / (xx.t_.back()) + 10; } else { cout << "!!!! " << cur.t_.size() << endl; l = ((int64_t)cur.t_.back() * BASE + cur.t_[cur.t_.size() - 2]) / (xx.t_.back() + 1); r = ((int64_t)cur.t_.back() * BASE + cur.t_[cur.t_.size() - 2]) / (xx.t_.back()) + 10; } l = max(l, (int64_t)0); assert(l < BASE); r = min(r, BASE - 1); while (l < r) { m = (l + r + 1) / 2; BigInteger tmp(xx); tmp *= d * m; if (tmp <= cur) { l = m; } else { r = m - 1; } } } else { continue; } res[i - 1] = l; BigInteger tmp(xx); tmp *= d * res[i - 1]; cur -= tmp; } t_ = res; is_negative_ = is_negative; Normalize(); DIV += getTime(); isDivided = false; return *this; } BigInteger& operator%=(const BigInteger& x) { *this -= (*this / x) * x; Normalize(); return *this; } std::string toString() const { std::ostringstream string; if (is_negative_) { string << "-"; } string << t_.back(); for (auto it = ++t_.rbegin(); it != t_.rend(); ++it) { string << std::setw(LOG) << std::setfill('0') << *it; } return string.str(); } explicit operator bool() const { return t_.size() > 1 || t_[0]; } friend std::istream& operator>>(std::istream& in, BigInteger& x) { assert(0); std::string s; in >> s; std::reverse(s.begin(), s.end()); if (s.back() == '-') { x.is_negative_ = true; s.pop_back(); } else { x.is_negative_ = false; } while (s.size() % 4) { s.push_back('0'); } x.t_.resize(s.size() / 1); for (size_t i = 0; i < s.size(); i += 1) { x.t_[i / 1] = s[i] - '0'; } x.Normalize(); return in; } friend std::ostream& operator<<(std::ostream& out, const BigInteger& x) { out << x.toString(); return out; } BigInteger operator-() { BigInteger tmp(*this); tmp.is_negative_ ^= true; tmp.Normalize(); return tmp; } private: void Normalize() { while (t_.size() > 1 && t_.back() == 0) { t_.pop_back(); } if (t_.size() == 1 && t_[0] == 0) { is_negative_ = false; } } // Это умножение на 10^degree BigInteger& MuliplicationDegree10(int degree) { size_t n = t_.size(); t_.resize(t_.size() + degree); std::rotate(t_.begin(), t_.begin() + n, t_.end()); Normalize(); return *this; } // Это используется, если число на которое мы умножаем не очень большое, // потому что FFT работает достаточно медленно BigInteger& SquareMultiplication(const BigInteger& x) { std::vector t(t_.size() + x.t_.size() + 1); for (size_t i = 0; i < t_.size(); ++i) { for (size_t j = 0; j < x.t_.size(); ++j) { t[i + j + 1] += (static_cast(t_[i]) * x.t_[j]) / BASE; t[i + j] += (static_cast(t_[i]) * x.t_[j]) % BASE; } } t_.resize(t_.size() + x.t_.size() + 1); for (size_t i = 0; i < t.size(); ++i) { t[i + 1] += t[i] / BASE; t_[i] = t[i] % BASE; } while (t_.size() > 1 && t_.back() == 0) t_.pop_back(); is_negative_ ^= x.is_negative_; Normalize(); return *this; } static constexpr int64_t LOG = 9; static constexpr int64_t BASE = 1000000000; std::vector t_ = {0}; bool is_negative_ = false; }; BigInteger operator-(const BigInteger& a, const BigInteger& b) { return BigInteger(a) -= b; } BigInteger operator*(const BigInteger& a, const BigInteger& b) { return BigInteger(a) *= b; } BigInteger operator/(const BigInteger& a, const BigInteger& b) { BigInteger tmp(a); tmp /= b; return tmp; } BigInteger operator%(const BigInteger& a, const BigInteger& b) { BigInteger tmp(a); tmp %= b; return tmp; } BigInteger operator+(const BigInteger& a, const BigInteger& b) { return BigInteger(a) += b; } bool operator>(const BigInteger& a, const BigInteger& b) { return b < a; } bool operator<=(const BigInteger& a, const BigInteger& b) { return !(a > b); } bool operator>=(const BigInteger& a, const BigInteger& b) { return b <= a; } bool operator==(const BigInteger& a, const BigInteger& b) { return !(a < b) && !(b < a); } bool operator!=(const BigInteger& a, const BigInteger& b) { return !(a == b); } class Rational; Rational operator+(const Rational& a, const Rational& b); Rational operator-(const Rational& a, const Rational& b); Rational operator*(const Rational& a, const Rational& b); Rational operator/(const Rational& a, const Rational& b); bool operator>(const Rational& a, const Rational& b); bool operator<=(const Rational& a, const Rational& b); bool operator>=(const Rational& a, const Rational& b); bool operator==(const Rational& a, const Rational& b); bool operator!=(const Rational& a, const Rational& b); class Rational { public: Rational() = default; Rational(const BigInteger& a) { numerator_ = a; } Rational(int a) { numerator_ = a; } Rational(const Rational& x) : numerator_(x.numerator_), denominator_(x.denominator_) {} void swap(Rational& x) { std::swap(numerator_, x.numerator_); std::swap(denominator_, x.denominator_); } Rational& operator=(const Rational& x) { Rational tmp(x); swap(tmp); return *this; } Rational& operator+=(const Rational& a) { if (this == &a) { return *this += Rational(a); } numerator_ *= a.denominator_; numerator_ += a.numerator_ * denominator_; denominator_ *= a.denominator_; Normalize(); return *this; } Rational& operator-=(const Rational& a) { if (this == &a) { return *this -= Rational(a); } numerator_ *= a.denominator_; numerator_ -= a.numerator_ * denominator_; denominator_ *= a.denominator_; Normalize(); return *this; } Rational& operator*=(const Rational& a) { numerator_ *= a.numerator_; denominator_ *= a.denominator_; Normalize(); return *this; } Rational& operator/=(const Rational& a) { if (this == &a) { return *this /= Rational(a); } numerator_ *= a.denominator_; denominator_ *= a.numerator_; Normalize(); return *this; } Rational operator-() const { Rational tmp(*this); tmp.numerator_ *= -1; return tmp; } bool operator<(const Rational& a) const { return numerator_ * a.denominator_ < denominator_ * a.numerator_; } std::string toString() { std::ostringstream str; str << numerator_; if (denominator_ != 1) { str << "/" << denominator_; } return str.str(); } std::string asDecimal(size_t precision) { BigInteger t = 1; for (size_t i = 0; i < precision; ++i) { t *= 10; } std::ostringstream str; BigInteger result = (numerator_ * t) / denominator_; if (result < 0) { str << "-", result *= -1; } str << result / t << "." << std::setw(precision) << std::setfill('0') << result % t; return str.str(); } explicit operator double() { return std::stod(asDecimal(20)); //while(true); return 1.0; } BigInteger gcd(const BigInteger& a, const BigInteger& b) { return b ? gcd(b, a % b) : a; } void Normalize(bool f = false) { f = true; if (f) { if (numerator_ == 0) { denominator_ = 1; return; } BigInteger d = gcd(numerator_, denominator_); numerator_ /= d; denominator_ /= d; if (denominator_ < 0) { numerator_ *= -1; denominator_ *= -1; } } } BigInteger numerator_ = 0; BigInteger denominator_ = 1; }; Rational operator+(const Rational& a, const Rational& b) { Rational tmp(a); tmp += b; return tmp; } Rational operator-(const Rational& a, const Rational& b) { Rational tmp(a); tmp -= b; return tmp; } Rational operator*(const Rational& a, const Rational& b) { Rational tmp(a); tmp *= b; return tmp; } Rational operator/(const Rational& a, const Rational& b) { Rational tmp(a); tmp /= b; return tmp; } bool operator>(const Rational& a, const Rational& b) { return b < a; } bool operator<=(const Rational& a, const Rational& b) { return a < b || a == b; } bool operator>=(const Rational& a, const Rational& b) { return a > b || a == b; } bool operator==(const Rational& a, const Rational& b) { return !(b < a) && !(a < b); } bool operator!=(const Rational& a, const Rational& b) { return !(a == b); } template struct is_prime_ { static const bool value = !(N % i == 0) && is_prime_::value; }; template struct is_prime_ { static const bool value = false; }; template struct is_prime_ { static const bool value = true; }; template struct is_prime { static const bool value = (is_prime_::value && (N != 1)) || (N == 2); }; template static const bool is_prime_v = is_prime::value; template struct is_power_of_prime_ { static const bool value = is_power_of_prime_= N ? 1 : N)), N % i == 0 ? i : i + 2, N % i == 0>::value; }; template struct is_power_of_prime_<1, i, find_DIVider> { static const bool value = true; }; template struct is_power_of_prime_<0, i, find_DIVider> { static const bool value = false; }; template struct is_power_of_prime { static const bool value = is_power_of_prime_::value; }; template struct has_primitive_root { static const bool value = is_power_of_prime<(N == 4 || N == 2 ? 1 : (N % 4 == 0 ? 0 : (N % 2 == 0 ? N / 2 : N)))>::value; }; template static const bool has_primitive_root_v = has_primitive_root::value; template void my_static_assert() { int* a = new int [T ? 1 : -1]; delete[] a; } unsigned phi(unsigned n) { unsigned result = n; for (int i = 2; i * i <= n; ++i) { if (n % i == 0) { while (n % i == 0) { n /= i; } result -= result / i; } } if (n > 1) result -= result / n; return result; } template class Residue { public: Residue() = default; Residue(const Residue& another) = default; explicit Residue(long long x) : x(x) { Normilize(); } explicit operator int() const { return x; } bool operator==(const Residue& another) const { return x == another.x; } bool operator!=(const Residue& another) { return !(*this == another); } Residue& operator+=(const Residue& another); Residue& operator-=(const Residue& another); Residue& operator*=(const Residue& another); Residue& operator/=(const Residue& another); Residue getInverse() const { my_static_assert>(); return Residue(pow(x, Modulus - 2)); } unsigned order() const { unsigned phi = ::phi(Modulus); unsigned copy_phi = phi; std::vector first = {1}, second = {phi}; for (unsigned i = 2; i * i <= phi; ++i) { if (copy_phi % i == 0) { first.push_back(i); second.push_back(phi / i); } } std::reverse(second.begin(), second.end()); first.insert(first.end(), second.begin(), second.end()); for (unsigned i: first) { if (pow(i) == Residue(1)) { return i; } } return 0; } static Residue getPrimitiveRoot() { my_static_assert>(); if (Modulus == 2) { return Residue(1); } else if (Modulus == 4) { return Residue(3); } unsigned phi = ::phi(Modulus); unsigned copy_phi = phi; std::vector deviders; for (unsigned i = 2; i * i <= phi; ++i) { if (copy_phi % i == 0) { deviders.push_back(i); deviders.push_back(phi / i); while (copy_phi % i == 0) copy_phi /= i; } } if (copy_phi > 1) deviders.push_back(copy_phi); for (unsigned res = 2; res < phi; ++res) { bool ok = pow(res, phi) == 1; for (unsigned i: deviders) { ok &= pow(res, i) != 1; } if (ok) { return Residue(res); } } return Residue(phi); } Residue pow(unsigned degree) const; private: static int pow(int x, unsigned degree); void Normilize() { x %= static_cast(Modulus); if (x < 0) x += Modulus; } long long x; }; template Residue& Residue::operator+=(const Residue& another) { x += another.x; Normilize(); return *this; } template Residue& Residue::operator-=(const Residue& another) { x -= another.x; Normilize(); return *this; } template Residue& Residue::operator*=(const Residue& another) { x *= another.x; Normilize(); return *this; } template Residue& Residue::operator/=(const Residue& another) { *this *= another.getInverse(); return *this; } template Residue operator+(const Residue& a, const Residue& b) { Residue tmp(a); tmp += b; return tmp; } template Residue operator-(const Residue& a, const Residue& b) { Residue tmp(a); tmp -= b; return tmp; } template Residue operator*(const Residue& a, const Residue& b) { Residue tmp(a); tmp *= b; return tmp; } template Residue operator/(const Residue& a, const Residue& b) { Residue tmp(a); tmp /= b; return tmp; } template Residue Residue::pow(unsigned degree) const { return Residue(pow(x, degree)); } template int Residue::pow(int x, unsigned degree) { if (degree == 0) { return 1; } else { long long t = pow(x, degree / 2); return (((t * t) % static_cast(Modulus)) * (degree % 2 == 0 ? 1 : x)) % static_cast(Modulus); } } std::istream& operator>>(std::istream& in, Rational& x) { std::string s; in >> s; x = 0; bool p = false; BigInteger a(0), b(1); for (size_t i = 0; i < s.size(); ++i) { if (s[i] == '-') { b *= -1; } else if (s[i] == '.') { p = true; } else { a = 10 * a + s[i] - '0'; if (p) { b *= 10; } } } x = a; x /= b; return in; } std::ostream& operator<<(std::ostream& out, Rational& x) { assert(0); while (true); out << x.toString(); return out; } template class Matrix { template friend class Matrix; public: Matrix() : t(N, vector(M)) {} template Matrix(std::initializer_list> a) { for (auto i: a) { t.push_back({}); for (auto j: i) { t.back().push_back(Field(j)); } } } static Matrix identityMatrix() { Matrix result; for (size_t i = 0; i < N; ++i) { result.t[i][i] = Field(1); } return result; } bool operator==(const Matrix& a) const { bool ok = true; for (size_t i = 0; i < N; ++i) { for (size_t j = 0; j < M; ++j) { ok &= t[i][j] == a.t[i][j]; } } return ok; } Matrix& operator+=(const Matrix& a) { for (size_t i = 0; i < N; ++i) { for (size_t j = 0; j < M; ++j) { t[i][j] += a.t[i][j]; } } return *this; } Matrix& operator-=(const Matrix& a) { for (size_t i = 0; i < N; ++i) { for (size_t j = 0; j < M; ++j) { t[i][j] -= a.t[i][j]; } } return *this; } Matrix& operator*=(const Matrix& a) { static_assert(N == M); *this = (*this * a); return *this; } bool operator!=(const Matrix& a) const { return !(*this == a); } vector& operator[](size_t i) { return t[i]; } const vector& operator[](size_t i) const { return t[i]; } Matrix& operator+=(Matrix& another) { for (size_t i = 0; i < N; ++i) { for (size_t j = 0; j < M; ++j) { t[i][j] += another[i][j]; } } } Matrix& operator-=(Matrix& another) { for (size_t i = 0; i < N; ++i) { for (size_t j = 0; j < M; ++j) { t[i][j] -= another[i][j]; } } } Matrix& operator*=(Field x) { for (size_t i = 0; i < N; ++i) { for (size_t j = 0; j < M; ++j) { t[i][j] *= x; } } return *this; } Matrix transposed() const { Matrix result; for (size_t i = 0; i < N; ++i) { for (size_t j = 0; j < M; ++j) { result[j][i] = t[i][j]; } } return result; } bool isZeroRow(size_t i) const { bool res = true; for (size_t j = 0; j < M; ++j) { res &= t[i][j] == Field(0); } return res; } void multiplyRow(size_t i, Field x) { for (size_t j = 0; j < M; ++j) { t[i][j] *= x; } } void subtractRows(size_t i, size_t j) { for (size_t k = 0; k < M; ++k) { t[i][k] -= t[j][k]; } } Field simplifiedView() { Field result = Field(1); for (size_t i = 0; i < std::min(N, M); ++i) { size_t found = M; for (size_t j = i; j < N && found == M; ++j) { if (t[i][j] != Field(0)) found = j; } if (found == M) { result = Field(0); continue; } swap(t[i], t[found]); result *= t[i][i]; multiplyRow(i, Field(1) / t[i][i]); for (size_t j = 0; j < N; ++j) { if (j != i) { if (t[j][i] != Field(0)) { vector tt = t[i]; multiplyRow(i, t[j][i]); subtractRows(j, i); swap(t[i], tt); } } } cout << clock() / CLOCKS_PER_SEC << "Div: " << endl << DIV << endl << " Mul: " << endl << mul << endl << " Plus: " << endl << PPLUS << endl; for (int ii = 0; ii < getCountRow(); ++ii) { for (int jj = 0; jj < getCountColumn(); ++jj) { cout << t[ii][jj].toString().size() << " "; } cout << endl; } cout << endl; cout << result.toString() << endl; } cerr << clock() / CLOCKS_PER_SEC << "Div: " << endl << DIV << endl << " Mul: " << endl << mul << endl << " Plus: " << endl << PPLUS << endl; return result; } Field det() const { static_assert(N == M); Matrix a(*this); return a.simplifiedView(); } template Matrix glueTogether(const Matrix& another) const { Matrix result; for (size_t i = 0; i < N; ++i) { result.t[i] = t[i]; result.t[i].insert(result.t[i].end(), another.t[i].begin(), another.t[i].end()); } return result; } template std::pair, Matrix> split() const { std::pair, Matrix> result; for (size_t i = 0; i < N; ++i) { result.first.t[i] = vector(t[i].begin(), t[i].begin() + K); result.second.t[i] = vector(t[i].begin() + K, t[i].end()); } return result; } Matrix inverted() const { // bool x = !(std::is_same::value && N == 20); // assert(x); static_assert(N == M); Matrix a = this->glueTogether(Matrix::identityMatrix()); a.simplifiedView(); return a.template split().second; } void invert() { *this = this->inverted(); } size_t getCountRow() const { return N; } size_t getCountColumn() const { return M; } size_t rank() const { Matrix a(*this); a.simplifiedView(); size_t res = 0; for (size_t i = 0; i < N; ++i) { res += !a.isZeroRow(i); } return res; } Field trace() const { static_assert(N == M); Field res = static_cast(0); for (size_t i = 0; i < N; ++i) res += t[i][i]; return res; } vector getRow(unsigned i) const { return t[i]; } vector getColumn(unsigned i) const { vector res(N); for (size_t j = 0; j < N; ++j) res[j] = t[j][i]; return res; } private: vector> t; }; template Matrix operator*(const Matrix& a, const Matrix& b) { bool ppp = is_same::value; if (ppp) { return Matrix::identityMatrix(); } // assert(!(N == 20 && ppp)); Matrix result; for (size_t i = 0; i < N; ++i) { for (size_t j = 0; j < K; ++j) { for (size_t k = 0; k < M; ++k) { result[i][j] += a[i][k] * b[k][j]; } } } return result; } template Matrix operator+(const Matrix& a, const Matrix& b) { Matrix tmp(a); tmp += b; return tmp; } template Matrix operator-(const Matrix& a, const Matrix& b) { Matrix tmp(a); tmp -= b; return tmp; } template Matrix operator*(const Matrix& a, Field x) { Matrix tmp(a); tmp *= x; return tmp; } template Matrix operator*(Field x, const Matrix& a) { Matrix tmp(a); tmp *= x; return tmp; } template using SquareMatrix = Matrix;