This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub tko919/library
#include "Math/bigint.hpp"
#pragma once #include "Convolution/arbitrary.hpp" template <int D> struct bigint { using u128 = __uint128_t; static const int B = pow(10, D); int sign = 0; vector<int> v; static int get_D() { return D; } static int get_B() { return B; } bigint() {} bigint(const vector<int> &_v, bool _s = false) : sign(_s), v(_v) {} bigint(ll x) { if (x < 0) x *= -1, sign = 1; while (x) { v.push_back(x % B); x /= B; } } bigint(string s) { if (s[0] == '-') s.erase(s.begin()), sign = 1; int add = 0, cnt = 0, base = 1; while (s.size()) { if (cnt == D) { v.push_back(add); cnt = 0; add = 0; base = 1; } add = (s.back() - '0') * base + add; cnt++; base *= 10; s.pop_back(); } if (add) v.push_back(add); } bigint operator-() const { bigint res = *this; res.sign ^= 1; return res; } bigint abs() const { bigint res = *this; res.sign = 0; return res; } int &operator[](const int i) { return v[i]; } int size() const { return v.size(); } void norm() { rep(i, 0, v.size() - 1) { if (v[i] >= 0) { v[i + 1] += v[i] / B; v[i] %= B; } else { int c = (-v[i] + B - 1) / B; v[i] += c * B; v[i + 1] -= c; } } while (!v.empty() and v.back() >= B) { int c = v.back() / B; v.back() %= B; v.push_back(c); } while (!v.empty() and v.back() == 0) v.pop_back(); } string to_str() const { string res; if (v.empty()) return "0"; if (sign) res += '-'; res += to_string(v.back()); for (int i = v.size() - 2; i >= 0; i--) { string add; int w = v[i]; rep(_, 0, D) { add += ('0' + (w % 10)); w /= 10; } reverse(ALL(add)); res += add; } return res; } friend istream &operator>>(istream &is, bigint<D> &x) { string tmp; is >> tmp; x = bigint(tmp); return is; } friend ostream &operator<<(ostream &os, bigint<D> x) { os << x.to_str(); return os; } bigint &operator<<=(const int &x) { if (!v.empty()) { vector<int> add(x, 0); v.insert(v.begin(), ALL(add)); } return *this; } bigint &operator>>=(const int &x) { v = vector<int>(v.begin() + min(x, (int)v.size()), v.end()); return *this; } bigint &operator+=(const bigint &x) { if (sign != x.sign) { *this -= (-x); return *this; } if ((int)v.size() < (int)x.size()) v.resize(x.size(), 0); rep(i, 0, x.size()) { v[i] += x.v[i]; } norm(); return *this; } bigint &operator-=(const bigint &x) { if (sign != x.sign) { *this += (-x); return *this; } if (abs() < x.abs()) { *this = x - (*this); sign ^= 1; return *this; } rep(i, 0, x.size()) { v[i] -= x.v[i]; } norm(); return *this; } bigint &operator*=(const bigint &x) { sign ^= x.sign; auto v1 = ArbitraryMult<u128>(v, x.v); u128 add = 0; v.clear(); v.reserve(v1.size() + 3); for (int i = 0;; i++) { if (i >= int(v1.size()) and add == 0) break; if (i < int(v1.size())) add += v1[i]; v.push_back(add % B); add /= B; } norm(); return *this; } bigint div_naive(const bigint &a, const bigint &b) { if (SZ(b) == 1) return a.div(b.v.back()); if (a < b) return bigint(); int norm = B / (b.v.back() + 1); bigint x = a.mul(norm), y = b.mul(norm); int yb = y.v.back(); bigint quo, rem; quo.v.resize(x.size() - y.size() + 1); rem.v = {x.v.end() - y.size(), x.v.end()}; for (int i = x.size() - y.size(); i >= 0; i--) { if (rem.size() == y.size()) { if (rem >= y) { quo[i] = 1; rem -= y; } } else if (rem.size() > y.size()) { ll rb = ll(rem.v.back()) * B + rem[rem.size() - 2]; int q = rb / yb; bigint yq = y.mul(q); while (rem < yq) { q--; yq -= y; } rem -= yq; while (rem >= y) { q++; rem -= y; } quo[i] = q; } if (i) rem.v.insert(rem.v.begin(), x[i - 1]); } return quo; } bigint &operator/=(const bigint &x) { bigint a = abs(), b = x.abs(); sign ^= x.sign; if (a < b) return *this = bigint(); if (b.size() == 1) return *this = a.div(b.v.back()); int deg = a.size() - b.size() + 2, k = deg; while (k > 64) k = (k + 1) >> 1; bigint base(1); base <<= (b.size() + k); bigint inv(div_naive(base, b)); while (k < deg) { bigint y = inv.square(); y.v.insert(y.v.begin(), 0); int l = min(SZ(b), k * 2 + 1); bigint pref; pref.v = {b.v.end() - l, b.v.end()}; y *= pref; y >>= l; inv = inv + inv; inv <<= k + 1; inv -= y; inv.v.erase(inv.v.begin()); k <<= 1; } inv >>= (k - deg); (*this) = a * inv; (*this) >>= int(a.size() + 2); bigint mul = (*this) * b; while (mul + b <= a) { (*this) += bigint(1); mul += b; } while (mul > a) { (*this) -= bigint(1); mul -= b; } return *this; } bigint &operator%=(const bigint &x) { bigint div = (*this) / x; (*this) -= div * x; return *this; } bigint square() const { bigint ret; auto v1 = ArbitraryMult<u128>(v, v); u128 add = 0; ret.v.reserve(v1.size() + 3); for (int i = 0;; i++) { if (i >= int(v1.size()) and add == 0) break; if (i < int(v1.size())) add += v1[i]; ret.v.push_back(add % B); add /= B; } return ret; } bigint mul(ll x) const { bigint res; if (x < 0) res.sign ^= 1, x *= -1; u128 add = 0; res.v.reserve(v.size() + 3); for (int i = 0;; i++) { if (i >= int(v.size()) and add == 0) break; if (i < int(v.size())) add += ll(v[i]) * x; res.v.push_back(add % B); add /= B; } return res; } bigint div(ll x) const { bigint res = *this; if (x < 0) res.sign ^= 1, x *= -1; ll add = 0; for (int i = res.v.size() - 1; i >= 0; i--) { add = add * B + res.v[i]; int q = add / x, r = add % x; res.v[i] = q, add = r; } res.norm(); return res; } bigint operator<<(const int &x) const { return bigint(*this) <<= x; } bigint operator>>(const int &x) const { return bigint(*this) >>= x; } bigint operator+(const bigint &x) const { return bigint(*this) += x; } bigint operator-(const bigint &x) const { return bigint(*this) -= x; } bigint operator*(const bigint &x) const { return bigint(*this) *= x; } bigint operator/(const bigint &x) const { return bigint(*this) /= x; } bigint operator%(const bigint &x) const { return bigint(*this) %= x; } bool operator<(const bigint &x) const { if (sign != x.sign) return sign > x.sign; if ((int)v.size() != (int)x.size()) { if (sign) return (int)x.size() < (int)v.size(); else return (int)v.size() < (int)x.size(); } for (int i = v.size() - 1; i >= 0; i--) if (v[i] != x.v[i]) { if (sign) return x.v[i] < v[i]; else return v[i] < x.v[i]; } return false; } bool operator>(const bigint &x) const { return x < *this; } bool operator<=(const bigint &x) const { return !(*this > x); } bool operator>=(const bigint &x) const { return !(*this < x); } bool operator==(const bigint &x) const { return !(*this < x) and !(*this > x); } bool operator!=(const bigint &x) const { return !(*this == x); } }; typedef bigint<9> Bigint; struct Bigfloat { Bigint v; int p = 0; Bigfloat() {} Bigfloat(const ll &_v) { v = Bigint(_v); } Bigfloat(const Bigint &_v, int _p = 0) : v(_v), p(_p) {} void set(int _p) { if (p < _p) { v >>= (_p - p); } else { v <<= (p - _p); } p = _p; } Bigint to_int() const { if (p < 0) return v >> (-p); else return v << p; } Bigfloat &operator+=(const Bigfloat &x) { if (p > x.p) set(x.p), v += x.v; else v += x.v << (x.p - p); return *this; } Bigfloat &operator-=(const Bigfloat &x) { if (p > x.p) set(x.p), v -= x.v; else v -= x.v << (x.p - p); return *this; } Bigfloat square() { Bigfloat res = *this; res.p *= 2; res.v = res.v.square(); return res; } Bigfloat mul(ll x) { Bigfloat res = *this; res.v = v.mul(x); return res; } Bigfloat div(ll x) { Bigfloat res = *this; res.v = v.div(x); return res; } Bigfloat &operator*=(const Bigfloat &x) { p += x.p; v *= x.v; return *this; } Bigfloat &operator/=(const Bigfloat &x) { p -= x.p; v /= x.v; return *this; } Bigfloat operator+(const Bigfloat &x) const { return Bigfloat(*this) += x; } Bigfloat operator-(const Bigfloat &x) const { return Bigfloat(*this) -= x; } Bigfloat operator*(const Bigfloat &x) const { return Bigfloat(*this) *= x; } Bigfloat operator/(const Bigfloat &x) const { return Bigfloat(*this) /= x; } string to_str() { string res = v.abs().to_str(); int d = Bigint::get_D(); if (p * d > 0) res += string(p * d, '0'); else if (-p * d >= 1 and -p * d < (int)res.size()) { res = res.substr(0, (int)res.size() + p * d) + '.' + res.substr((int)res.size() + p * d); } else if (-p * d >= (int)res.size()) res = "0." + string(-p * d - (int)res.size(), '0') + res; if (v.sign) { res.insert(res.begin(), '-'); } return res; } friend ostream &operator<<(ostream &os, Bigfloat x) { os << x.to_str(); return os; } }; Bigfloat sqrt(ll n, int d) { Bigfloat res(Bigint((ll)sqrt(1LL * Bigint::get_B() * Bigint::get_B() / n)), -1), pre; int cur = 1; while (pre.v != res.v) { cur = min(cur << 1, d); pre = res; Bigfloat add = Bigfloat(1) - res.square().mul(n); add.set(-cur); res += (res * add).div(2); res.set(-cur); } return res.mul(n); } /** * @brief Big Integer(Float) */
#line 2 "Convolution/ntt.hpp" template <typename T> struct NTT { static constexpr int rank2 = __builtin_ctzll(T::get_mod() - 1); std::array<T, rank2 + 1> root; // root[i]^(2^i) == 1 std::array<T, rank2 + 1> iroot; // root[i] * iroot[i] == 1 std::array<T, std::max(0, rank2 - 2 + 1)> rate2; std::array<T, std::max(0, rank2 - 2 + 1)> irate2; std::array<T, std::max(0, rank2 - 3 + 1)> rate3; std::array<T, std::max(0, rank2 - 3 + 1)> irate3; NTT() { T g = 2; while (g.pow((T::get_mod() - 1) >> 1) == 1) { g += 1; } root[rank2] = g.pow((T::get_mod() - 1) >> rank2); iroot[rank2] = root[rank2].inv(); for (int i = rank2 - 1; i >= 0; i--) { root[i] = root[i + 1] * root[i + 1]; iroot[i] = iroot[i + 1] * iroot[i + 1]; } { T prod = 1, iprod = 1; for (int i = 0; i <= rank2 - 2; i++) { rate2[i] = root[i + 2] * prod; irate2[i] = iroot[i + 2] * iprod; prod *= iroot[i + 2]; iprod *= root[i + 2]; } } { T prod = 1, iprod = 1; for (int i = 0; i <= rank2 - 3; i++) { rate3[i] = root[i + 3] * prod; irate3[i] = iroot[i + 3] * iprod; prod *= iroot[i + 3]; iprod *= root[i + 3]; } } } void ntt(std::vector<T> &a, bool type = 0) { int n = int(a.size()); int h = __builtin_ctzll((unsigned int)n); a.resize(1 << h); if (type) { int len = h; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed while (len) { if (len == 1) { int p = 1 << (h - len); T irot = 1; for (int s = 0; s < (1 << (len - 1)); s++) { int offset = s << (h - len + 1); for (int i = 0; i < p; i++) { auto l = a[i + offset]; auto r = a[i + offset + p]; a[i + offset] = l + r; a[i + offset + p] = (unsigned long long)(T::get_mod() + l.v - r.v) * irot.v; ; } if (s + 1 != (1 << (len - 1))) irot *= irate2[__builtin_ctzll(~(unsigned int)(s))]; } len--; } else { // 4-base int p = 1 << (h - len); T irot = 1, iimag = iroot[2]; for (int s = 0; s < (1 << (len - 2)); s++) { T irot2 = irot * irot; T irot3 = irot2 * irot; int offset = s << (h - len + 2); for (int i = 0; i < p; i++) { auto a0 = 1ULL * a[i + offset + 0 * p].v; auto a1 = 1ULL * a[i + offset + 1 * p].v; auto a2 = 1ULL * a[i + offset + 2 * p].v; auto a3 = 1ULL * a[i + offset + 3 * p].v; auto a2na3iimag = 1ULL * T((T::get_mod() + a2 - a3) * iimag.v).v; a[i + offset] = a0 + a1 + a2 + a3; a[i + offset + 1 * p] = (a0 + (T::get_mod() - a1) + a2na3iimag) * irot.v; a[i + offset + 2 * p] = (a0 + a1 + (T::get_mod() - a2) + (T::get_mod() - a3)) * irot2.v; a[i + offset + 3 * p] = (a0 + (T::get_mod() - a1) + (T::get_mod() - a2na3iimag)) * irot3.v; } if (s + 1 != (1 << (len - 2))) irot *= irate3[__builtin_ctzll(~(unsigned int)(s))]; } len -= 2; } } T e = T(n).inv(); for (auto &x : a) x *= e; } else { int len = 0; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed while (len < h) { if (h - len == 1) { int p = 1 << (h - len - 1); T rot = 1; for (int s = 0; s < (1 << len); s++) { int offset = s << (h - len); for (int i = 0; i < p; i++) { auto l = a[i + offset]; auto r = a[i + offset + p] * rot; a[i + offset] = l + r; a[i + offset + p] = l - r; } if (s + 1 != (1 << len)) rot *= rate2[__builtin_ctzll(~(unsigned int)(s))]; } len++; } else { // 4-base int p = 1 << (h - len - 2); T rot = 1, imag = root[2]; for (int s = 0; s < (1 << len); s++) { T rot2 = rot * rot; T rot3 = rot2 * rot; int offset = s << (h - len); for (int i = 0; i < p; i++) { auto mod2 = 1ULL * T::get_mod() * T::get_mod(); auto a0 = 1ULL * a[i + offset].v; auto a1 = 1ULL * a[i + offset + p].v * rot.v; auto a2 = 1ULL * a[i + offset + 2 * p].v * rot2.v; auto a3 = 1ULL * a[i + offset + 3 * p].v * rot3.v; auto a1na3imag = 1ULL * T(a1 + mod2 - a3).v * imag.v; auto na2 = mod2 - a2; a[i + offset] = a0 + a2 + a1 + a3; a[i + offset + 1 * p] = a0 + a2 + (2 * mod2 - (a1 + a3)); a[i + offset + 2 * p] = a0 + na2 + a1na3imag; a[i + offset + 3 * p] = a0 + na2 + (mod2 - a1na3imag); } if (s + 1 != (1 << len)) rot *= rate3[__builtin_ctzll(~(unsigned int)(s))]; } len += 2; } } } } vector<T> mult(const vector<T> &a, const vector<T> &b) { if (a.empty() or b.empty()) return vector<T>(); int as = a.size(), bs = b.size(); int n = as + bs - 1; if (as <= 30 or bs <= 30) { if (as > 30) return mult(b, a); vector<T> res(n); rep(i, 0, as) rep(j, 0, bs) res[i + j] += a[i] * b[j]; return res; } int m = 1; while (m < n) m <<= 1; vector<T> res(m); rep(i, 0, as) res[i] = a[i]; ntt(res); if (a == b) rep(i, 0, m) res[i] *= res[i]; else { vector<T> c(m); rep(i, 0, bs) c[i] = b[i]; ntt(c); rep(i, 0, m) res[i] *= c[i]; } ntt(res, 1); res.resize(n); return res; } }; /** * @brief Number Theoretic Transform */ #line 2 "Math/modint.hpp" template <unsigned mod = 1000000007> struct fp { unsigned v; static constexpr int get_mod() { return mod; } constexpr unsigned inv() const { assert(v != 0); int x = v, y = mod, p = 1, q = 0, t = 0, tmp = 0; while (y > 0) { t = x / y; x -= t * y, p -= t * q; tmp = x, x = y, y = tmp; tmp = p, p = q, q = tmp; } if (p < 0) p += mod; return p; } constexpr fp(ll x = 0) : v(x >= 0 ? x % mod : (mod - (-x) % mod) % mod) {} fp operator-() const { return fp() - *this; } fp pow(ull t) { fp res = 1, b = *this; while (t) { if (t & 1) res *= b; b *= b; t >>= 1; } return res; } fp &operator+=(const fp &x) { if ((v += x.v) >= mod) v -= mod; return *this; } fp &operator-=(const fp &x) { if ((v += mod - x.v) >= mod) v -= mod; return *this; } fp &operator*=(const fp &x) { v = ull(v) * x.v % mod; return *this; } fp &operator/=(const fp &x) { v = ull(v) * x.inv() % mod; return *this; } fp operator+(const fp &x) const { return fp(*this) += x; } fp operator-(const fp &x) const { return fp(*this) -= x; } fp operator*(const fp &x) const { return fp(*this) *= x; } fp operator/(const fp &x) const { return fp(*this) /= x; } bool operator==(const fp &x) const { return v == x.v; } bool operator!=(const fp &x) const { return v != x.v; } friend istream &operator>>(istream &is, fp &x) { return is >> x.v; } friend ostream &operator<<(ostream &os, const fp &x) { return os << x.v; } }; template <unsigned mod> void rd(fp<mod> &x) { fastio::rd(x.v); } template <unsigned mod> void wt(fp<mod> x) { fastio::wt(x.v); } template <typename T> T Inv(ll n) { static const int md = T::get_mod(); static vector<T> buf({0, 1}); assert(n > 0); n %= md; while (SZ(buf) <= n) { int k = SZ(buf), q = (md + k - 1) / k; buf.push_back(buf[k * q - md] * q); } return buf[n]; } template <typename T> T Fact(ll n, bool inv = 0) { static const int md = T::get_mod(); static vector<T> buf({1, 1}), ibuf({1, 1}); assert(n >= 0 and n < md); while (SZ(buf) <= n) { buf.push_back(buf.back() * SZ(buf)); ibuf.push_back(ibuf.back() * Inv<T>(SZ(ibuf))); } return inv ? ibuf[n] : buf[n]; } template <typename T> T nPr(int n, int r, bool inv = 0) { if (n < 0 || n < r || r < 0) return 0; return Fact<T>(n, inv) * Fact<T>(n - r, inv ^ 1); } template <typename T> T nCr(int n, int r, bool inv = 0) { if (n < 0 || n < r || r < 0) return 0; return Fact<T>(n, inv) * Fact<T>(r, inv ^ 1) * Fact<T>(n - r, inv ^ 1); } template <typename T> T nHr(int n, int r, bool inv = 0) { return nCr<T>(n + r - 1, r, inv); } /** * @brief Modint */ #line 4 "Convolution/arbitrary.hpp" using M1 = fp<1045430273>; using M2 = fp<1051721729>; using M3 = fp<1053818881>; NTT<M1> N1; NTT<M2> N2; NTT<M3> N3; constexpr int r_12 = M2(M1::get_mod()).inv(); constexpr int r_13 = M3(M1::get_mod()).inv(); constexpr int r_23 = M3(M2::get_mod()).inv(); constexpr int r_1323 = M3(ll(r_13) * r_23).v; constexpr ll w1 = M1::get_mod(); constexpr ll w2 = ll(w1) * M2::get_mod(); template <typename T> vector<T> ArbitraryMult(const vector<int> &a, const vector<int> &b) { if (a.empty() or b.empty()) return vector<T>(); int n = a.size() + b.size() - 1; vector<T> res(n); if (min(a.size(), b.size()) <= 60) { rep(i, 0, a.size()) rep(j, 0, b.size()) res[i + j] += T(a[i]) * b[j]; return res; } vector<int> vals[3]; vector<M1> a1(ALL(a)), b1(ALL(b)), c1 = N1.mult(a1, b1); vector<M2> a2(ALL(a)), b2(ALL(b)), c2 = N2.mult(a2, b2); vector<M3> a3(ALL(a)), b3(ALL(b)), c3 = N3.mult(a3, b3); for (M1 x : c1) vals[0].push_back(x.v); for (M2 x : c2) vals[1].push_back(x.v); for (M3 x : c3) vals[2].push_back(x.v); rep(i, 0, n) { ll p = vals[0][i]; ll q = (vals[1][i] + M2::get_mod() - p) * r_12 % M2::get_mod(); ll r = ((vals[2][i] + M3::get_mod() - p) * r_1323 + (M3::get_mod() - q) * r_23) % M3::get_mod(); res[i] = (T(r) * w2 + q * w1 + p); } return res; } template <typename T> vector<T> ArbitraryMult(const vector<T> &a, const vector<T> &b) { vector<int> A, B; for (auto &x : a) A.push_back(x.v); for (auto &x : b) B.push_back(x.v); return ArbitraryMult<T>(A, B); } /** * @brief Arbitrary Mod Convolution */ #line 3 "Math/bigint.hpp" template <int D> struct bigint { using u128 = __uint128_t; static const int B = pow(10, D); int sign = 0; vector<int> v; static int get_D() { return D; } static int get_B() { return B; } bigint() {} bigint(const vector<int> &_v, bool _s = false) : sign(_s), v(_v) {} bigint(ll x) { if (x < 0) x *= -1, sign = 1; while (x) { v.push_back(x % B); x /= B; } } bigint(string s) { if (s[0] == '-') s.erase(s.begin()), sign = 1; int add = 0, cnt = 0, base = 1; while (s.size()) { if (cnt == D) { v.push_back(add); cnt = 0; add = 0; base = 1; } add = (s.back() - '0') * base + add; cnt++; base *= 10; s.pop_back(); } if (add) v.push_back(add); } bigint operator-() const { bigint res = *this; res.sign ^= 1; return res; } bigint abs() const { bigint res = *this; res.sign = 0; return res; } int &operator[](const int i) { return v[i]; } int size() const { return v.size(); } void norm() { rep(i, 0, v.size() - 1) { if (v[i] >= 0) { v[i + 1] += v[i] / B; v[i] %= B; } else { int c = (-v[i] + B - 1) / B; v[i] += c * B; v[i + 1] -= c; } } while (!v.empty() and v.back() >= B) { int c = v.back() / B; v.back() %= B; v.push_back(c); } while (!v.empty() and v.back() == 0) v.pop_back(); } string to_str() const { string res; if (v.empty()) return "0"; if (sign) res += '-'; res += to_string(v.back()); for (int i = v.size() - 2; i >= 0; i--) { string add; int w = v[i]; rep(_, 0, D) { add += ('0' + (w % 10)); w /= 10; } reverse(ALL(add)); res += add; } return res; } friend istream &operator>>(istream &is, bigint<D> &x) { string tmp; is >> tmp; x = bigint(tmp); return is; } friend ostream &operator<<(ostream &os, bigint<D> x) { os << x.to_str(); return os; } bigint &operator<<=(const int &x) { if (!v.empty()) { vector<int> add(x, 0); v.insert(v.begin(), ALL(add)); } return *this; } bigint &operator>>=(const int &x) { v = vector<int>(v.begin() + min(x, (int)v.size()), v.end()); return *this; } bigint &operator+=(const bigint &x) { if (sign != x.sign) { *this -= (-x); return *this; } if ((int)v.size() < (int)x.size()) v.resize(x.size(), 0); rep(i, 0, x.size()) { v[i] += x.v[i]; } norm(); return *this; } bigint &operator-=(const bigint &x) { if (sign != x.sign) { *this += (-x); return *this; } if (abs() < x.abs()) { *this = x - (*this); sign ^= 1; return *this; } rep(i, 0, x.size()) { v[i] -= x.v[i]; } norm(); return *this; } bigint &operator*=(const bigint &x) { sign ^= x.sign; auto v1 = ArbitraryMult<u128>(v, x.v); u128 add = 0; v.clear(); v.reserve(v1.size() + 3); for (int i = 0;; i++) { if (i >= int(v1.size()) and add == 0) break; if (i < int(v1.size())) add += v1[i]; v.push_back(add % B); add /= B; } norm(); return *this; } bigint div_naive(const bigint &a, const bigint &b) { if (SZ(b) == 1) return a.div(b.v.back()); if (a < b) return bigint(); int norm = B / (b.v.back() + 1); bigint x = a.mul(norm), y = b.mul(norm); int yb = y.v.back(); bigint quo, rem; quo.v.resize(x.size() - y.size() + 1); rem.v = {x.v.end() - y.size(), x.v.end()}; for (int i = x.size() - y.size(); i >= 0; i--) { if (rem.size() == y.size()) { if (rem >= y) { quo[i] = 1; rem -= y; } } else if (rem.size() > y.size()) { ll rb = ll(rem.v.back()) * B + rem[rem.size() - 2]; int q = rb / yb; bigint yq = y.mul(q); while (rem < yq) { q--; yq -= y; } rem -= yq; while (rem >= y) { q++; rem -= y; } quo[i] = q; } if (i) rem.v.insert(rem.v.begin(), x[i - 1]); } return quo; } bigint &operator/=(const bigint &x) { bigint a = abs(), b = x.abs(); sign ^= x.sign; if (a < b) return *this = bigint(); if (b.size() == 1) return *this = a.div(b.v.back()); int deg = a.size() - b.size() + 2, k = deg; while (k > 64) k = (k + 1) >> 1; bigint base(1); base <<= (b.size() + k); bigint inv(div_naive(base, b)); while (k < deg) { bigint y = inv.square(); y.v.insert(y.v.begin(), 0); int l = min(SZ(b), k * 2 + 1); bigint pref; pref.v = {b.v.end() - l, b.v.end()}; y *= pref; y >>= l; inv = inv + inv; inv <<= k + 1; inv -= y; inv.v.erase(inv.v.begin()); k <<= 1; } inv >>= (k - deg); (*this) = a * inv; (*this) >>= int(a.size() + 2); bigint mul = (*this) * b; while (mul + b <= a) { (*this) += bigint(1); mul += b; } while (mul > a) { (*this) -= bigint(1); mul -= b; } return *this; } bigint &operator%=(const bigint &x) { bigint div = (*this) / x; (*this) -= div * x; return *this; } bigint square() const { bigint ret; auto v1 = ArbitraryMult<u128>(v, v); u128 add = 0; ret.v.reserve(v1.size() + 3); for (int i = 0;; i++) { if (i >= int(v1.size()) and add == 0) break; if (i < int(v1.size())) add += v1[i]; ret.v.push_back(add % B); add /= B; } return ret; } bigint mul(ll x) const { bigint res; if (x < 0) res.sign ^= 1, x *= -1; u128 add = 0; res.v.reserve(v.size() + 3); for (int i = 0;; i++) { if (i >= int(v.size()) and add == 0) break; if (i < int(v.size())) add += ll(v[i]) * x; res.v.push_back(add % B); add /= B; } return res; } bigint div(ll x) const { bigint res = *this; if (x < 0) res.sign ^= 1, x *= -1; ll add = 0; for (int i = res.v.size() - 1; i >= 0; i--) { add = add * B + res.v[i]; int q = add / x, r = add % x; res.v[i] = q, add = r; } res.norm(); return res; } bigint operator<<(const int &x) const { return bigint(*this) <<= x; } bigint operator>>(const int &x) const { return bigint(*this) >>= x; } bigint operator+(const bigint &x) const { return bigint(*this) += x; } bigint operator-(const bigint &x) const { return bigint(*this) -= x; } bigint operator*(const bigint &x) const { return bigint(*this) *= x; } bigint operator/(const bigint &x) const { return bigint(*this) /= x; } bigint operator%(const bigint &x) const { return bigint(*this) %= x; } bool operator<(const bigint &x) const { if (sign != x.sign) return sign > x.sign; if ((int)v.size() != (int)x.size()) { if (sign) return (int)x.size() < (int)v.size(); else return (int)v.size() < (int)x.size(); } for (int i = v.size() - 1; i >= 0; i--) if (v[i] != x.v[i]) { if (sign) return x.v[i] < v[i]; else return v[i] < x.v[i]; } return false; } bool operator>(const bigint &x) const { return x < *this; } bool operator<=(const bigint &x) const { return !(*this > x); } bool operator>=(const bigint &x) const { return !(*this < x); } bool operator==(const bigint &x) const { return !(*this < x) and !(*this > x); } bool operator!=(const bigint &x) const { return !(*this == x); } }; typedef bigint<9> Bigint; struct Bigfloat { Bigint v; int p = 0; Bigfloat() {} Bigfloat(const ll &_v) { v = Bigint(_v); } Bigfloat(const Bigint &_v, int _p = 0) : v(_v), p(_p) {} void set(int _p) { if (p < _p) { v >>= (_p - p); } else { v <<= (p - _p); } p = _p; } Bigint to_int() const { if (p < 0) return v >> (-p); else return v << p; } Bigfloat &operator+=(const Bigfloat &x) { if (p > x.p) set(x.p), v += x.v; else v += x.v << (x.p - p); return *this; } Bigfloat &operator-=(const Bigfloat &x) { if (p > x.p) set(x.p), v -= x.v; else v -= x.v << (x.p - p); return *this; } Bigfloat square() { Bigfloat res = *this; res.p *= 2; res.v = res.v.square(); return res; } Bigfloat mul(ll x) { Bigfloat res = *this; res.v = v.mul(x); return res; } Bigfloat div(ll x) { Bigfloat res = *this; res.v = v.div(x); return res; } Bigfloat &operator*=(const Bigfloat &x) { p += x.p; v *= x.v; return *this; } Bigfloat &operator/=(const Bigfloat &x) { p -= x.p; v /= x.v; return *this; } Bigfloat operator+(const Bigfloat &x) const { return Bigfloat(*this) += x; } Bigfloat operator-(const Bigfloat &x) const { return Bigfloat(*this) -= x; } Bigfloat operator*(const Bigfloat &x) const { return Bigfloat(*this) *= x; } Bigfloat operator/(const Bigfloat &x) const { return Bigfloat(*this) /= x; } string to_str() { string res = v.abs().to_str(); int d = Bigint::get_D(); if (p * d > 0) res += string(p * d, '0'); else if (-p * d >= 1 and -p * d < (int)res.size()) { res = res.substr(0, (int)res.size() + p * d) + '.' + res.substr((int)res.size() + p * d); } else if (-p * d >= (int)res.size()) res = "0." + string(-p * d - (int)res.size(), '0') + res; if (v.sign) { res.insert(res.begin(), '-'); } return res; } friend ostream &operator<<(ostream &os, Bigfloat x) { os << x.to_str(); return os; } }; Bigfloat sqrt(ll n, int d) { Bigfloat res(Bigint((ll)sqrt(1LL * Bigint::get_B() * Bigint::get_B() / n)), -1), pre; int cur = 1; while (pre.v != res.v) { cur = min(cur << 1, d); pre = res; Bigfloat add = Bigfloat(1) - res.square().mul(n); add.set(-cur); res += (res * add).div(2); res.set(-cur); } return res.mul(n); } /** * @brief Big Integer(Float) */