#include #include #include #include #include #include #include #include #include #include // M(n): O(n^1.47) #define HAVE_UINT128 using uint32 = uint32_t; using int64 = int64_t; using uint64 = uint64_t; #ifdef HAVE_UINT128 using uint128 = __uint128_t; using word_t = uint64; using dword_t = uint128; #else using word_t = uint32; using dword_t = uint64; #endif namespace bits { inline constexpr uint32 pop_count(uint32 n) { return __builtin_popcount(n); } inline constexpr uint64 pop_count(uint64 n) { return __builtin_popcountll(n); } inline constexpr uint32 ctz(uint32 n) { return __builtin_ctz(n); } inline constexpr uint64 ctz(uint64 n) { return __builtin_ctzll(n); } inline constexpr uint32 clz(uint32 n) { return __builtin_clz(n); } inline constexpr uint64 clz(uint64 n) { return __builtin_clzll(n); } template inline constexpr T ilog2(T n) { return n == 0 ? 0 : (sizeof(n) * 8 - 1) - clz(n); } #ifdef HAVE_UINT128 inline uint64 floor_div_small(uint128 a, uint64 b) { uint64 q, r; __asm__ ( "divq\t%4" : "=a"(q), "=d"(r) : "0"(uint64(a)), "1"(uint64(a >> 64)), "rm"(b) ); return q; } #else inline uint32 floor_div_small(uint64 a, uint32 b) { uint32 q, r; __asm__ ( "divl\t%4" : "=a"(q), "=d"(r) : "0"(uint32(a)), "1"(uint32(a >> 32)), "rm"(b) ); return q; } #endif inline word_t mul_inv(word_t n) { word_t x = n; while (1) { word_t t = n * x; if (t == 1) { break; } x *= 2 - t; } return x; } } struct FastDiv { FastDiv() {} FastDiv(word_t n) { if (n == 1) { shamt = 0; magic_num = 1; } else { shamt = (8 * sizeof(dword_t)) - 1 - bits::clz(n - 1); magic_num = bits::floor_div_small((dword_t(1) << shamt) + n - 1, n); } mod = n; } friend word_t operator / (word_t n, FastDiv d) { return (dword_t(n) * d.magic_num) >> d.shamt; } friend word_t& operator /= (word_t& n, FastDiv d) { return n = n / d; } friend word_t operator % (word_t n, FastDiv d) { return n - n / d * d.mod; } friend word_t& operator %= (word_t& n, FastDiv d) { return n = n % d; } word_t magic_num; word_t shamt; word_t mod; }; class FastDiv21 { public: FastDiv21(word_t n) { shamt_ = bits::clz(n); d_ = n << shamt_; v_ = reciprocal(d_); }; FastDiv21() {} word_t divisor() const { return d_ >> shamt_; } // divl in x86. static void divmod( dword_t u, FastDiv21 fd, word_t& ret_q, word_t& ret_r) { u <<= fd.shamt_; word_t u_hi = u >> (8 * sizeof(word_t)); word_t u_lo = word_t(u); dword_t q = dword_t(u_hi) * fd.v_ + u; word_t q_hi = (q >> (8 * sizeof(word_t))) + 1; word_t r = u_lo - q_hi * fd.d_; if (r > word_t(q)) { q_hi -= 1; r += fd.d_; } if (r >= fd.d_) { q_hi += 1; r -= fd.d_; } ret_q = q_hi; ret_r = r >> fd.shamt_; } friend word_t operator / (dword_t n, FastDiv21 fd) { return FastDiv21::div(n, fd); } friend dword_t& operator /= (dword_t& n, FastDiv21 fd) { return n = n / fd; } friend word_t operator % (dword_t n, FastDiv21 fd) { return FastDiv21::mod(n, fd); } friend dword_t& operator %= (dword_t& n, FastDiv21 fd) { return n = n % fd; } private: static word_t div(dword_t u, FastDiv21 fd) { word_t q, dummy; divmod(u, fd, q, dummy); return q; } static word_t mod(dword_t u, FastDiv21 fd) { word_t dummy, r; divmod(u, fd, dummy, r); return r; } word_t reciprocal(word_t n) const { return bits::floor_div_small(~(dword_t(n) << (sizeof(word_t) * 8)), n); } word_t d_; word_t shamt_; word_t v_; }; class BigInt { public: enum { W_BITS = sizeof(word_t) * 8, #ifdef HAVE_UINT128 MUL_TOOM22_THRESHOLD = 24, // words MUL_TOOM33_THRESHOLD = 120, // >= 4 words SQ_TOOM22_THRESHOLD = 32, // words SQ_TOOM33_THRESHOLD = 160, // >= 4 words INT_THRESHOLD = 64, // words STR_THRESHOLD = 1 << 12, // bits DIVMOD_THRESHOLD = 1 << 13, // bits #else MUL_TOOM22_THRESHOLD = 16, // words MUL_TOOM33_THRESHOLD = 80, // >= 4 words SQ_TOOM22_THRESHOLD = 28, // words SQ_TOOM33_THRESHOLD = 140, // >= 4 words INT_THRESHOLD = 64, // words STR_THRESHOLD = 1 << 9, // bits DIVMOD_THRESHOLD = 1 << 12, // bits #endif }; static_assert(MUL_TOOM33_THRESHOLD >= 4, "mul_toom33"); static_assert(SQ_TOOM33_THRESHOLD >= 4, "sq_toom33"); public: struct BarrettReduction { ~BarrettReduction() { delete divisor; delete reciprocal; } BarrettReduction() : divisor(nullptr), reciprocal(nullptr), shamt(0) {} BarrettReduction(BigInt* d, BigInt* r, word_t s, word_t t) : divisor(d), reciprocal(r), shamt(s), two_e(t) {} static void divmod(const BigInt& n, const BarrettReduction& d, BigInt& q, BigInt& r) { if (d.divisor->is_one()) { q = n >> d.two_e; r = n.subinteger(0, d.two_e); return; } // not precomputed if (d.reciprocal == nullptr) { if (d.two_e == 0) { return BigInt::divmod(n, *d.divisor, q, r); } else { BigInt::divmod(n >> d.two_e, *d.divisor, q, r); r <<= d.two_e; r |= n.subinteger(0, d.two_e); return; } } // (base & 1) != 0 if (d.two_e == 0) { return _divmod(n, d, q, r); } word_t k = d.shamt + d.two_e; if (n.bit_length() < k) { q = BigInt(); r = BigInt(n); return; } // (base & 1) == 0 word_t nblock = 1 + (n.bit_length() - d.two_e) / d.shamt; word_t q_bit_len = n.bit_length() - k + 1; q.set_zero(); q.resize(bits_to_words(q_bit_len)); r = n.subinteger(d.two_e + (nblock - 1) * d.shamt, d.two_e + nblock * d.shamt); BigInt tmp_q, tmp_r; for (int i = nblock - 2; i >= 0; --i) { r <<= d.shamt; r |= n.subinteger(d.two_e + i * d.shamt, d.two_e + (i + 1) * d.shamt); _divmod(r, d, tmp_q, tmp_r); _lshifted_lor(q, tmp_q, i * d.shamt); r = tmp_r; } q.normalize(); r <<= d.two_e; r |= n.subinteger(0, d.two_e); } BigInt* divisor; BigInt* reciprocal; word_t shamt; word_t two_e; private: static void _divmod(const BigInt& n, const BarrettReduction& d, BigInt& q, BigInt& r) { q = ((n >> (d.shamt - 1)) * *d.reciprocal) >> (d.shamt + 1); r = n - q * *d.divisor; while (r >= *d.divisor) { r -= *d.divisor; q += 1; } } }; struct BaseData { BaseData(word_t base, bool lower) : base_s(base), lower(lower) { fd = FastDiv(base); word_t lim = (word_t(1) << (W_BITS - 1)) - 1; base_l = 1; e = 0; while (lim >= base_s) { base_l *= base_s; lim /= fd; e += 1; } fd21 = FastDiv21(base_l); } FastDiv fd; FastDiv21 fd21; word_t base_s; word_t base_l; word_t e; bool lower; private: BaseData() {} }; using container_t = std::vector; private: struct _bigint_t { _bigint_t(int s, const word_t* d, word_t sz) : sign(s), data(const_cast(d)), size(sz) {} _bigint_t(word_t* d) : sign(0), data(d), size(0) {} _bigint_t() {} int sign; word_t* data; word_t size; }; public: BigInt() : sign_(0) { n_ = container_t(1, 0); } BigInt(int n) : sign_(n < 0) { n_ = container_t(1, uint32(std::abs(n))); // fixed } BigInt(uint32 n) : sign_(0) { n_ = container_t(1, n); } #ifdef HAVE_UINT128 BigInt(int64 n) : sign_(n < 0) { n_ = container_t(1, std::abs(n)); } BigInt(uint64 n) : sign_(0) { n_ = container_t(1, n); } inline operator uint64() const { uint64 ret = (*this)[0]; if (is_neg()) { ret = -ret; } return ret; } #else BigInt(int64 n) : sign_(n < 0) { n_ = container_t(2); (*this)[0] = n; (*this)[1] = std::abs(int(n >> W_BITS)); } BigInt(uint64 n) : sign_(0) { n_ = container_t(2); (*this)[0] = n; (*this)[1] = n >> W_BITS; } inline operator uint64() const { uint64 ret = (size() >= 2 ? ((uint64((*this)[1]) << 32) | (*this)[0]) : (*this)[0]); if (is_neg()) { ret = -ret; } return ret; } #endif BigInt(const char* cstr, word_t base=10) { _int(cstr, std::strlen(cstr), base, *this); } BigInt(const std::string& str, word_t base=10) { _int(str.c_str(), str.size(), base, *this); } BigInt(word_t size, int n) : sign_(n < 0) { assert(size > 0); n_ = container_t(size, 0); n_[0] = n; } BigInt(const BigInt& n) { this->sign_ = n.sign_; this->n_ = n.n_; } // container inline word_t size() const { return n_.size(); } void push_back(word_t w) { n_.push_back(w); } void resize(word_t size) { if (size == 0) { set_zero(); } else { n_.resize(size, 0); } } inline word_t operator [] (word_t idx) const { return n_[idx]; } inline word_t& operator [] (word_t idx) { return n_[idx]; } inline word_t* data() { return n_.data(); } inline const word_t* data() const { return n_.data(); } // basic void change_sign() { if (!is_zero()) { sign_ ^= 1; } } void set_sign(int s) { if (!is_zero()) { sign_ = s; } } void clear_sign() { sign_ = 0; } int sign() const { return sign_; } bool is_neg() const { return sign_ == 1; } bool is_zero() const { return size() == 1 && (*this)[0] == 0; } bool is_one() const { // unsigned return size() == 1 && (*this)[0] == 1; } BigInt abs() const { BigInt ret = BigInt(*this); ret.clear_sign(); return ret; } void set_zero() { resize(1); clear_sign(); (*this)[0] = 0; } void normalize() { for (int i = size() - 1; i >= 0; --i) { if ((*this)[i]) { return resize(i + 1); } } return set_zero(); } word_t bit_length() const { if (is_zero()) { return 0; } word_t len = size(); return W_BITS * len - bits::clz((*this)[len - 1]); } static word_t bits_to_words(word_t bits) { return (bits + W_BITS - 1) / W_BITS; } word_t pop_count() const { word_t ret = 0; for (word_t i = 0; i < this->size(); ++i) { ret += bits::pop_count((*this)[i]); } return ret; } static BigInt mask(word_t bit_length) { if (bit_length == 0) { return BigInt(); } BigInt ret = BigInt(); word_t ret_size = 1 + (bit_length - 1) / W_BITS; ret.resize(ret_size); for (word_t i = 0; i < ret_size; ++i) { ret[i] = ~word_t(0); } ret[ret_size - 1] &= ~word_t(0) >> ((W_BITS - bit_length) & (W_BITS - 1)); return ret; } // string std::string str(word_t base=10, bool lower=true) const { if (is_zero()) { return "0"; } BaseData bdata = BaseData(base, lower); return _fast_str(bdata); } // relational bool operator == (const BigInt& rhs) const { if (this->is_neg() ^ rhs.is_neg()) { return false; } if (size() != rhs.size()) { return false; } for (word_t i = 0; i < rhs.size(); ++i) { if ((*this)[i] != rhs[i]) { return false; } } return true; } bool operator != (const BigInt& rhs) const { return !(*this == rhs); } bool operator < (const BigInt& rhs) const { if (this->is_neg()) { if (!rhs.is_neg()) { return true; } else { return _ucomp(rhs, *this, false); } } else { if (rhs.is_neg()) { return false; } else { return _ucomp(*this, rhs, false); } } } /* slow */ bool operator < (const int rhs) const { return *this < BigInt(rhs); } #ifndef HAVE_UINT128 /* tmp */ bool operator < (const uint64 rhs) const { return *this < BigInt(rhs); } #endif bool operator < (const word_t rhs) const { if (is_neg()) { return true; } return !(size() >= 2) && (*this)[0] < rhs; } bool operator <= (const BigInt& rhs) const { if (this->is_neg()) { if (!rhs.is_neg()) { return true; } else { return _ucomp(rhs, *this, true); } } else { if (rhs.is_neg()) { return false; } else { return _ucomp(*this, rhs, true); } } } bool operator <= (const word_t rhs) const { if (is_neg()) { return true; } return !(size() >= 2) && (*this)[0] <= rhs; } bool operator >= (const BigInt& rhs) const { return rhs <= *this; } bool operator >= (const word_t rhs) const { return !(*this < rhs); } bool operator > (const BigInt& rhs) const { return rhs < *this; } bool operator > (const word_t rhs) const { return !(*this <= rhs); } // unary operator bool operator ! () const { return is_zero(); } BigInt operator + () const { return BigInt(*this); } BigInt operator - () const { BigInt ret = BigInt(*this); ret.change_sign(); return ret; } // arith BigInt operator + (const BigInt& rhs) const { BigInt res; _add(*this, rhs, false, res); return res; } /* slow */ BigInt operator + (const int rhs) const { BigInt res; _add(*this, BigInt(rhs), false, res); return res; } BigInt& operator += (const BigInt& rhs) { _add(*this, rhs, false, *this); return *this; } BigInt operator - (const BigInt& rhs) const { BigInt res; _add(*this, rhs, true, res); return res; } /* slow */ BigInt operator - (const int rhs) const { BigInt res; _add(*this, BigInt(rhs), true, res); return res; } BigInt& operator -= (const BigInt& rhs) { _add(*this, rhs, true, *this); return *this; } BigInt operator * (const BigInt& rhs) const { BigInt ret; _mul(*this, rhs, ret); return ret; } BigInt operator * (const word_t rhs) const { BigInt ret; _mul1(*this, rhs, ret); return ret; } BigInt& operator *= (const BigInt& rhs) { return *this = *this * rhs; } BigInt& operator *= (const word_t rhs) { _mul1(*this, rhs, *this); return *this; } BigInt operator / (const BigInt& rhs) const { BigInt q, r; divmod(*this, rhs, q, r); return q; } BigInt operator / (const word_t rhs) const { BigInt ret; word_t dummy; divmod_n1(*this, rhs, ret, dummy); return ret; } BigInt operator % (const BigInt& rhs) const { BigInt q, r; divmod(*this, rhs, q, r); return r; } word_t operator % (const word_t rhs) const { BigInt dummy; word_t ret; divmod_n1(*this, rhs, dummy, ret); return ret; } static void divmod(const BigInt& a, const BigInt& b, BigInt& q, BigInt& r) { /* sign */ _fast_udivmod(a, b, q, r); } static void divmod_n1(const BigInt& n, const word_t d, BigInt& qs, word_t& r) { /* sign */ auto fd = FastDiv21(d); return divmod_n1(n, fd, qs, r); } static void divmod_n1(const BigInt& n, const FastDiv21 fd, BigInt& qs, word_t& r) { const word_t d = fd.divisor(); word_t size = n.size(); r = n[size - 1]; word_t q = 0; if (r >= d) { FastDiv21::divmod(r, fd, q, r); qs.resize(size); qs[size - 1] = q; } else { qs.resize(size - 1); } for (int i = size - 2; i >= 0; --i) { FastDiv21::divmod((dword_t(r) << W_BITS) | n[i], fd, q, r); qs[i] = q; } } // logical [unsigned] BigInt operator & (const BigInt& rhs) const { if (size() > rhs.size()) { return rhs & *this; } BigInt ret = BigInt(*this); _land(ret.data(), rhs.data(), ret.size()); ret.normalize(); return ret; } BigInt operator & (const int rhs) const { return BigInt((*this)[0] & rhs); } BigInt& operator &= (const BigInt& rhs) { if (size() > rhs.size()) { resize(rhs.size()); } _land(this->data(), rhs.data(), this->size()); this->normalize(); return *this; } BigInt operator ^ (const BigInt& rhs) const { if (size() < rhs.size()) { return rhs ^ *this; } BigInt ret = BigInt(*this); _lxor(ret.data(), rhs.data(), rhs.size()); ret.normalize(); return ret; } BigInt& operator ^= (const BigInt& rhs) { if (size() < rhs.size()) { resize(rhs.size()); } _lxor(this->data(), rhs.data(), this->size()); this->normalize(); return *this; } BigInt operator | (const BigInt& rhs) const { if (size() < rhs.size()) { return rhs | *this; } BigInt ret = BigInt(*this); _lor(ret.data(), rhs.data(), rhs.size()); return ret; } BigInt& operator |= (const BigInt& rhs) { if (size() < rhs.size()) { resize(rhs.size()); } _lor(this->data(), rhs.data(), rhs.size()); return *this; } // shift [unsigned] BigInt operator << (int shamt) const { if (shamt < 0) { return *this >> -shamt; } if (shamt == 0) { return BigInt(*this); } if (is_zero()) { return BigInt(); } BigInt ret = BigInt(); word_t bit_len = bit_length(); ret.resize(bits_to_words(bit_len + shamt)); _lshift(this->data(), this->size(), shamt, ret.data(), ret.size()); return ret; } BigInt operator << (word_t shamt) const { return *this << int(shamt); } BigInt& operator <<= (int shamt) { if (shamt < 0) { return *this >>= -shamt; } if (shamt == 0) { return *this; } if (is_zero()) { return *this; } word_t bit_len = bit_length(); word_t size = this->size(); this->resize(bits_to_words(bit_len + shamt)); _lshift(this->data(), size, shamt, this->data(), this->size()); return *this; } BigInt& operator <<= (word_t shamt) { return *this <<= int(shamt); } BigInt operator >> (int shamt) const { if (shamt < 0) { return *this << -shamt; } if (shamt == 0) { return BigInt(*this); } if (is_zero()) { return BigInt(); } word_t bit_len = bit_length(); if (bit_len <= word_t(shamt)) { return BigInt(); } BigInt ret = BigInt(); ret.resize(bits_to_words(bit_len - shamt)); _rshift(this->data(), this->size(), shamt, ret.data(), ret.size()); return ret; } BigInt operator >> (word_t shamt) const { return *this >> int(shamt); } BigInt& operator >>= (int shamt) { if (shamt < 0) { return *this <<= -shamt; } if (shamt == 0) { return *this; } if (is_zero()) { return *this; } word_t bit_len = bit_length(); if (bit_len <= word_t(shamt)) { this->set_zero(); return *this; } word_t ret_size = bits_to_words(bit_len - shamt); _rshift(this->data(), this->size(), shamt, this->data(), ret_size); this->resize(ret_size); return *this; } BigInt& operator >>= (word_t shamt) { return *this >>= int(shamt); } // subinteger ? strip ? BigInt subinteger(const word_t bit_beg, word_t bit_end) const { word_t bit_len = bit_length(); if (bit_end > bit_len) { bit_end = bit_len; } if (bit_beg >= bit_end) { return BigInt(0); } word_t bit_size = bit_end - bit_beg; word_t ret_size = bits_to_words(bit_size); word_t in_size = bits_to_words(bit_end); BigInt ret = BigInt(ret_size, 0); _rshift(this->data(), in_size, bit_beg, ret.data(), ret_size); word_t mask = ~word_t(0) >> ((W_BITS - bit_size) & (W_BITS - 1)); ret[ret_size - 1] &= mask; ret.normalize(); return ret; } // functions BigInt isqrt() const { if (is_neg()) { assert(0); } BigInt s, r; _isqrt(*this, s, r); return s; } BigInt sqrt() const { return isqrt(); } BigInt pow(word_t e) const { BigInt ret = BigInt(1); if (e == 0) { return ret; } word_t mask = word_t(1) << bits::ilog2(e); while (mask) { if (mask & e) { ret *= *this; } mask >>= 1; if (!mask) { break; } ret *= ret; } return ret; } static BigInt fib(word_t n) { BigInt a, b; _fib(n, a, b); return a; } static void fib(word_t n, BigInt& a, BigInt& b) { _fib(n, a, b); } // output friend std::ostream & operator << (std::ostream& os, const BigInt& n) { return os << n.str(); } private: static BigInt _int_small(word_t size, const word_t* smalls, const word_t base) { BigInt ret; for (int i = size - 1; i >= 0; --i) { word_t c = _umul1_add(ret.data(), ret.size(), base, ret.data(), smalls[i]); if (c) { ret.push_back(c); } } return ret; } static BigInt _int_rec(word_t size, word_t e, word_t base_l, const word_t* smalls, const std::vector& large_bases) { if (size <= INT_THRESHOLD) { return _int_small(size, smalls, base_l); } else { word_t d = word_t(1) << e; if (size > d) { return _int_rec(d, e - 1, base_l, smalls , large_bases) + _int_rec(size - d, e - 1, base_l, smalls + d, large_bases) * large_bases[e]; } else { return _int_rec(size, e - 1, base_l, smalls , large_bases); } } } static void _int(const char* cstr, word_t size, word_t base, BigInt& res) { static const word_t char_to_word[128] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 0, 0, 0, 0, 0, 0, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 0, 0, 0, 0, 0, 0, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 0, 0, 0, 0, 0 }; res.set_zero(); if (size == 0) { return; } bool neg = false; word_t pos = 0; if (cstr[pos] == '-') { neg = true; ++pos; } word_t lim = ~word_t(0); word_t base_l = 1; auto fd = FastDiv(base); word_t e = 0; while (lim >= base) { lim /= fd; e += 1; base_l *= base; } // cstr => word_size word_t size_s = (size - pos + e - 1) / e; auto smalls = std::vector(size_s, 0); for (word_t i = 0; i < size_s; ++i) { int end = size - i * e; int beg = std::max(int(pos), end - int(e)); word_t t = 0; for (int j = beg; j < end; ++j) { t = t * base + char_to_word[cstr[j] & 0x7f]; } smalls[i] = t; } if (size_s == 1) { res = BigInt(smalls[0]); return; } word_t large_e = bits::ilog2(size_s - 1); auto larges = std::vector(large_e + 1); BigInt large_base = BigInt(base_l); for (word_t i = 0; i < large_e; ++i) { larges[i] = large_base; large_base *= large_base; } larges[large_e] = large_base; res = _int_rec(size_s, large_e, base_l, smalls.data(), larges); if (neg) { res.change_sign(); } } std::string _ustr(const BaseData& base, word_t zpad=0) const { static const char* chars = base.lower ? "0123456789abcdefghijklmnopqrstuvwxyz" : "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; std::string ret; BigInt qs = BigInt(*this); qs.clear_sign(); auto fd = base.fd; word_t r = 0; while (qs >= base.base_l) { divmod_n1(qs, base.fd21, qs, r); for (word_t i = 0; i < base.e; ++i) { word_t s = r / fd; ret += chars[r - s * fd.mod]; r = s; } } r = qs[0]; while (r > 0) { word_t s = r / fd; ret += chars[r - s * fd.mod]; r = s; } if (zpad && ret.size() < zpad) { ret.append(zpad - ret.size(), '0'); } std::reverse(ret.begin(), ret.end()); return ret; } static void _fast_str_rec( const BigInt& n, word_t e, const BaseData& base, const std::vector& large_bases, std::string& ret, bool zpad=false) { if (e + 1 <= 10) { ret += n._ustr(base, zpad ? (1 << (e + 1)) : 0); } else { BigInt q, r; BarrettReduction::divmod(n, large_bases[e], q, r); if (zpad || !q.is_zero()) { _fast_str_rec(q, e - 1, base, large_bases, ret, zpad); _fast_str_rec(r, e - 1, base, large_bases, ret, true); } else { _fast_str_rec(r, e - 1, base, large_bases, ret, zpad); } } } std::string _fast_str(const BaseData& base) const { std::string ret; if (is_neg()) { ret += "-"; } word_t bit_len = bit_length(); if (bit_len <= STR_THRESHOLD) { ret += this->_ustr(base); return ret; } auto large_bases = std::vector(50); word_t base_two_e = bits::ctz(base.base_s); BigInt base_pow = BigInt(base.base_s >> base_two_e); word_t e = 0; while (1) { large_bases[e].divisor = new BigInt(base_pow); large_bases[e].two_e = base_two_e; if ((base_pow.bit_length() + base_two_e) * 2 - 1 > bit_len) { break; } e += 1; base_pow *= base_pow; base_two_e <<= 1; } BigInt two = BigInt(1); word_t prev_k = 0; for (int i = 0; i < int(e) - 2; ++i) { word_t k = large_bases[i].divisor->bit_length(); two <<= 2 * (k - prev_k); large_bases[i].reciprocal = new BigInt(two / *large_bases[i].divisor); large_bases[i].shamt = k; prev_k = k; } _fast_str_rec(*this, e, base, large_bases, ret); return ret; } // relational static bool _ucomp(const BigInt& a, const BigInt& b, bool case_eq) { if (a.size() < b.size()) { return true; } else if (a.size() > b.size()) { return false; } for (int i = a.size() - 1; i >= 0; --i) { if (a[i] != b[i]) { return a[i] < b[i]; } } return case_eq; } static bool _shifted_ucomp(const BigInt& a, word_t ofs, const BigInt& b, bool case_eq) { // assert(a >= b); for (int i = b.size() - 1; i >= 0; --i) { if (a[i + ofs] != b[i]) { return a[i + ofs] < b[i]; } } return case_eq; } // add sub static void _add(const BigInt& a, const BigInt& b, bool is_sub, BigInt& res) { word_t a_size = a.size(); word_t b_size = b.size(); res.resize(std::max(a_size, b_size) + 1); // a.data() (and b.data()) might be changed. _bigint_t a_ = _bigint_t(a.sign(), a.data(), a_size); _bigint_t b_ = _bigint_t(b.sign(), b.data(), b_size); _bigint_t r_ = _bigint_t(res.data()); _add_ar(a_, b_, is_sub, r_); res.set_sign(r_.sign); res.normalize(); } static void _add_ar(const _bigint_t& a, const _bigint_t& b, bool is_sub, _bigint_t& res) { // normalize word_t b_size = b.size; if (is_sub) { while (b_size > a.size && b.data[b_size - 1] == 0) { b_size -= 1; } } if (a.size < b_size) { _add_ar(b, a, is_sub, res); if (is_sub) { res.sign ^= 1; } } else { int res_sign = a.sign; res.size = a.size; if (a.sign ^ b.sign ^ is_sub) { res_sign ^= _abs_sub(a.data, a.size, b.data, b_size, res.data); } else { word_t c = _uadd_core(a.data, a.size, b.data, b_size, res.data); if (c) { res.data[res.size++] = c; } } res.sign = res_sign; } } static word_t _uadd_core(const word_t* a, word_t a_size, const word_t* b, word_t b_size, word_t* res) { // assert(a_size >= b_size); word_t i = 0; dword_t c = 0; for (; i < b_size; ++i) { c = dword_t(a[i]) + b[i] + word_t(c); res[i] = word_t(c); c >>= W_BITS; } while (c && i < a_size) { c += a[i]; res[i++] = word_t(c); c >>= W_BITS; } if (res != a) { for ( ; i < a_size; ++i) { res[i] = a[i]; } } return c; } static void _usub(BigInt& a, const BigInt& b) { // assert(a >= b); _usub_core(a.data(), a.size(), b.data(), b.size(), a.data()); a.normalize(); } static void _usub_unnorm(BigInt& a, const BigInt& b) { _usub_core(a.data(), a.size(), b.data(), b.size(), a.data()); } static void _shifted_usub_unnorm(BigInt& a, word_t ofs, const BigInt& b) { _usub_core(a.data() + ofs, a.size() - ofs, b.data(), b.size(), a.data() + ofs); } static void _usub_core(const word_t* a, word_t a_size, const word_t* b, word_t b_size, word_t* res) { word_t i = 0; dword_t c = 0; for (; i < b_size; ++i) { c = dword_t(a[i]) - b[i] - word_t(c); res[i] = word_t(c); c = ((c >> W_BITS) >> (W_BITS - 1)); } while (c && i < a_size) { c = dword_t(a[i]) - word_t(c); res[i++] = word_t(c); c = ((c >> W_BITS) >> (W_BITS - 1)); } if (res != a) { for (; i < a_size; ++i) { res[i] = a[i]; } } } static int _abs_sub(const word_t* a, word_t a_size, const word_t* b, word_t b_size, word_t* res) { // assert(a_size >= b_size); for (int i = a_size - 1; i >= int(b_size); --i) { if (a[i]) { _usub_core(a, i + 1, b, b_size, res); return 0; } res[i] = 0; } for (int i = b_size - 1; i >= 0; --i) { if (a[i] != b[i]) { if (a[i] < b[i]) { _usub_core(b, i + 1, a, i + 1, res); return 1; } else { _usub_core(a, i + 1, b, i + 1, res); return 0; } } res[i] = 0; } return 0; } // mul static void _mul1(const BigInt& a, const word_t b, BigInt& res) { int sign = a.sign(); _umul1(a, b, res); res.set_sign(sign); } static void _umul1(const BigInt& a, const word_t b, BigInt& res) { if (a.is_zero() || b == 0) { res = BigInt(); return; } if (b == 1) { res = BigInt(a); return; } res.resize(a.size()); word_t c = _umul1_add(a.data(), a.size(), b, res.data(), 0); if (c) { res.push_back(c); } } static word_t _umul1_add(const word_t* a, word_t a_size, const word_t b, word_t* res, const word_t gamma) { dword_t c = gamma; for (word_t i = 0; i < a_size; ++i) { c += dword_t(a[i]) * b; res[i] = c; c >>= W_BITS; } return c; } static void _mul(const BigInt& a, const BigInt& b, BigInt &res) { if (&a == &b) { _square(a, res); } else { int sign = a.sign() ^ b.sign(); _umul(a, b, res); res.set_sign(sign); } } static void _umul(const BigInt& a, const BigInt& b, BigInt &res) { word_t a_bit_len = a.bit_length(); word_t b_bit_len = b.bit_length(); if (a_bit_len < b_bit_len) { return _umul(b, a, res); } if (b.size() == 1) { return _umul1(a, b[0], res); } if (a.is_zero() || b.is_zero()) { res = BigInt(); return; } if (a.is_one()) { res = BigInt(b); return; } if (b.is_one()) { res = BigInt(a); return; } word_t a_size = a.size(); word_t b_size = b.size(); if (b.size() <= MUL_TOOM22_THRESHOLD) { res.resize(a_size + b_size); _umul_basecase(a.data(), a_size, b.data(), b_size, res.data()); } else { word_t block_size = _umul_toom33_calc_block_size(a_size, b_size); word_t nblock = (a_size + block_size - 1) / block_size; res.resize(block_size * (nblock + 1)); word_t* work = new word_t[(block_size * 16 / 3) + 3 * 128]; BigInt c = BigInt(b); if (block_size > b_size) { c.resize(block_size); } if (nblock > 1) { BigInt res_s = BigInt(); res_s.resize(2 * block_size); for (word_t i = 0; i < nblock - 1; ++i) { _umul_toom33(a.data() + i * block_size, c.data(), block_size, res_s.data(), work); (void) _uadd_core( res.data() + i * block_size, 2 * block_size, res_s.data(), 2 * block_size, res.data() + i * block_size); } // ... :( BigInt last_a = a.subinteger((nblock - 1) * block_size * W_BITS, a.bit_length()); last_a.resize(block_size); _umul_toom33(last_a.data(), c.data(), block_size, res_s.data(), work); (void) _uadd_core( res.data() + (nblock - 1) * block_size, res.size() - (nblock - 1) * block_size, res_s.data(), 2 * block_size, res.data() + (nblock - 1) * block_size); } else { _umul_toom33(a.data(), c.data(), block_size, res.data(), work); } delete [] work; } res.normalize(); } static word_t _umul_toom33_calc_block_size(word_t a_size, word_t b_size) { double ratio = double(a_size) / b_size; // >= 1.0 word_t l = std::floor(ratio); word_t h = std::ceil(ratio); word_t block_size_a = (a_size + l - 1) / l; double cost1 = l * _umul_toom33_cost(block_size_a); double cost2 = h * _umul_toom33_cost(b_size); if (cost1 < cost2) { return block_size_a; } else { return b_size; } } static double _umul_toom33_cost(word_t block_size) { return std::pow(double(block_size), 1.47); } static int _umul_toom33_d1_d3( const word_t* a0, const word_t* a1, const word_t* a2, word_t size_l, word_t size_h, word_t* d1, word_t* d3) { // d1 = a2 + a1 + a1 // d3 = a2 - a1 + a0 (sign) word_t d1_size = size_l; word_t c = _uadd_core(a0, size_l, a2, size_h, d1); if (c) { d1[d1_size++] = c; } word_t d3_size = d1_size; int sign = _abs_sub(d1, d1_size, a1, size_l, d3); c = _uadd_core(d1, d1_size, a1, size_l, d1); if (c) { d1[d1_size++] = c; } if (d3_size == size_l) { if (d1_size == size_l) { d1[d1_size++] = 0; } d3[d3_size++] = 0; } return sign; } static void _umul_toom33_d2( const word_t* a0, const word_t* a1, const word_t* a2, word_t size_l, word_t size_h, word_t* d2) { // d2 = 4 * a2 + 2 * a1 + a0 word_t i = 0; dword_t c = 0; for ( ; i < size_h; ++i) { c = (((dword_t(a2[i]) << 1) + a1[i]) << 1) + a0[i] + word_t(c); d2[i] = c; c >>= W_BITS; } for ( ; i < size_l; ++i) { c = (dword_t(a1[i]) << 1) + a0[i] + word_t(c); d2[i] = c; c >>= W_BITS; } d2[i++] = c; } static void _umul_toom33( const word_t* a, const word_t* b, word_t size, word_t* res, word_t*& work) { while (size > 1 && (a[size - 1] | b[size - 1]) == 0) { size -= 1; res[size * 2] = res[size * 2 + 1] = 0; } if (size < MUL_TOOM33_THRESHOLD) { return _umul_toom22(a, b, size, res, work); } const word_t size_l = (size + 2) / 3; const word_t size_h = size - size_l * 2; const word_t work_size = 8 * (size_l + 1); auto* a0 = a; auto* a1 = a0 + size_l; auto* a2 = a1 + size_l; auto* b0 = b; auto* b1 = b0 + size_l; auto* b2 = b1 + size_l; auto* ad1 = work; auto* bd1 = ad1 + (size_l + 1); auto* ad3 = bd1 + (size_l + 1); auto* bd3 = ad3 + (size_l + 1); auto* ad2 = ad1; auto* bd2 = bd1; auto* s0 = res; auto* s1 = work + 6 * (size_l + 1); auto* s3 = work + 4 * (size_l + 1); auto* s2 = work + 2 * (size_l + 1); auto* s4 = res + 4 * size_l; work += work_size; int sign_ad3 = _umul_toom33_d1_d3(a0, a1, a2, size_l, size_h, ad1, ad3); int sign_bd3 = _umul_toom33_d1_d3(b0, b1, b2, size_l, size_h, bd1, bd3); _umul_toom33(ad1, bd1, size_l + 1, s1, work); _umul_toom33(ad3, bd3, size_l + 1, s3, work); _umul_toom33_d2(a0, a1, a2, size_l, size_h, ad2); _umul_toom33_d2(b0, b1, b2, size_l, size_h, bd2); _umul_toom33(ad2, bd2, size_l + 1, s2, work); _umul_toom33(a0, b0, size_l, s0, work); _umul_toom33(a2, b2, size_h, s4, work); auto s0_ = _bigint_t(0, s0, 2 * size_l); auto s1_ = _bigint_t(0, s1, 2 * (size_l + 1)); auto s2_ = _bigint_t(0, s2, 2 * (size_l + 1)); auto s3_ = _bigint_t(sign_ad3 ^ sign_bd3, s3, 2 * (size_l + 1)); auto s4_ = _bigint_t(0, s4, 2 * size_h); auto r_ = _bigint_t(0, res + size_l, size * 2 - size_l); _add_ar(s2_, s3_, 1, s2_); _exact_div(s2, s2_.size, 3); _add_ar(s1_, s3_, 1, s3_); _rshift(s3, s3_.size, 1, s3, s3_.size); _add_ar(s1_, s0_, 1, s1_); _add_ar(s2_, s1_, 1, s2_); _rshift(s2, s2_.size, 1, s2, s2_.size); _add_ar(s1_, s3_, 1, s1_); _add_ar(s1_, s4_, 1, s1_); _add_ar(s2_, s4_, 1, s2_); _add_ar(s2_, s4_, 1, s2_); _add_ar(s3_, s2_, 1, s3_); std::fill(res + size_l * 2, res + size_l * 4, 0); s2_.size = std::min(s2_.size, r_.size - size_l * 2); _add_ar(r_, s3_, 0, r_); r_.data += size_l; r_.size -= size_l; _add_ar(r_, s1_, 0, r_); r_.data += size_l; r_.size -= size_l; _add_ar(r_, s2_, 0, r_); work -= work_size; } static void _umul_toom22( const word_t* a, const word_t* b, word_t size, word_t* res, word_t*& work) { if (size <= MUL_TOOM22_THRESHOLD) { return _umul_basecase(a, size, b, size, res); } const word_t size_l = (size + 1) >> 1; const word_t size_h = size - size_l; const word_t work_size = size_l * 2 + 1; word_t* d1 = res; word_t* d2 = res + size_l; word_t* c2 = work; work += work_size; int sign1 = _abs_sub(a, size_l, a + size_l, size_h, d1); int sign2 = _abs_sub(b, size_l, b + size_l, size_h, d2); _umul_toom22(d1, d2, size_l, c2, work); _umul_toom22(a, b, size_l, res, work); // c0 _umul_toom22(a + size_l, b + size_l, size_h, res + size_l * 2, work); // c1 word_t c2_size = size_l * 2; _bigint_t c2_ = _bigint_t(sign1 ^ sign2 ^ 1, c2, c2_size); _bigint_t c1_ = _bigint_t(0, res + size_l * 2, size_h * 2); _bigint_t c0_ = _bigint_t(0, res, size_l * 2); _bigint_t r_ = _bigint_t(0, res + size_l , size_l + size_h * 2); _add_ar(c2_, c0_, 0, c2_); _add_ar(c2_, c1_, 0, c2_); _add_ar(r_, c2_, 0, r_); work -= work_size; } static void _umul_basecase( const word_t* a, word_t a_size, const word_t* b, word_t b_size, word_t* res) { std::fill(res, res + a_size + b_size, 0); for (word_t i = 0; i < b_size; ++i) { if (b[i] == 0) { continue; } dword_t c = 0; word_t beta = b[i]; for (word_t j = 0; j < a_size; ++j) { c = dword_t(beta) * a[j] + res[i + j] + word_t(c); res[i + j] = c; c >>= W_BITS; } if (c) { res[i + a_size] = c; } } } // square static void _square(const BigInt& a, BigInt& res) { if (a.is_zero()) { res = BigInt(); return; } if (a.is_one()) { res = BigInt(1); return; } word_t a_size = a.size(); res.resize(a_size * 2); if (a_size <= SQ_TOOM22_THRESHOLD) { _square_basecase(a.data(), a_size, res.data()); } else if (a_size <= SQ_TOOM33_THRESHOLD) { word_t* work = new word_t[a_size * 2 + 3 * 32]; _square_toom22(a.data(), a_size, res.data(), work); delete [] work; } else { word_t* work = new word_t[a_size * 4 + 6 * 96]; _square_toom33(a.data(), a_size, res.data(), work); delete [] work; } res.normalize(); return; } static void _square_toom33( const word_t* a, word_t size, word_t* res, word_t*& work) { while (size > 1 && a[size - 1] == 0) { size -= 1; res[size * 2] = res[size * 2 + 1] = 0; } if (size <= SQ_TOOM33_THRESHOLD) { return _square_toom22(a, size, res, work); } const word_t size_l = (size + 2) / 3; const word_t size_h = size - size_l * 2; const word_t work_size = 6 * (size_l + 1); auto* a0 = a; auto* a1 = a0 + size_l; auto* a2 = a1 + size_l; auto* d1 = work; auto* d2 = d1 + size_l + 1; auto* s0 = res; auto* s1 = d1 + (size_l + 1) * 2; auto* s2 = s1 + (size_l + 1) * 2; auto* s3 = d1; auto* s4 = res + size_l * 4; work += work_size; // d1 = a2 + a1 + a0 // d2 = a2 - a1 + a0 (void) _umul_toom33_d1_d3(a0, a1, a2, size_l, size_h, d1, d2); std::copy(a2, a2 + size_h, res); std::fill(res + size_h, res + size_l, 0); auto* a2_ = res; _square_toom33(d2, size_l + 1, s2, work); // s2 = d2 ** 2 [w4, w6) _square_toom33(d1, size_l + 1, s1, work); // s1 = d1 ** 2 [w2, w4) _umul_toom33(a1, a2_, size_l, s3, work); // s3 = 2 * a1 * a2 [w0, w2) _lshift(s3, size_l * 2, 1, s3, size_l * 2 + 1); _square_toom33(a0, size_l, s0, work); // s0 = a0 ** 2 [r0, r2) _square_toom33(a2, size_h, s4, work); // s4 = a2 ** 2 [r4, r6) auto s0_ = _bigint_t(0, s0, 2 * size_l); auto s1_ = _bigint_t(0, s1, 2 * (size_l + 1)); auto s2_ = _bigint_t(0, s2, 2 * (size_l + 1)); auto s3_ = _bigint_t(0, s3, 2 * size_l + 1); auto s4_ = _bigint_t(0, s4, 2 * size_h); auto r_ = _bigint_t(0, res + size_l, size * 2 - size_l); _add_ar(s2_, s1_, 0, s2_); _rshift(s2, s2_.size, 1, s2, s2_.size); _add_ar(s1_, s2_, 1, s1_); _add_ar(s1_, s3_, 1, s1_); _add_ar(s2_, s4_, 1, s2_); _add_ar(s2_, s0_, 1, s2_); std::fill(res + size_l * 2, res + size_l * 4, 0); s3_.size = std::min(s3_.size, r_.size - size_l * 2); _add_ar(r_, s1_, 0, r_); r_.data += size_l; r_.size -= size_l; _add_ar(r_, s2_, 0, r_); r_.data += size_l; r_.size -= size_l; _add_ar(r_, s3_, 0, r_); work -= work_size; } static void _square_toom22( const word_t* a, word_t size, word_t* res, word_t*& work) { if (size <= SQ_TOOM22_THRESHOLD) { return _square_basecase(a, size, res); } const word_t size_l = (size + 1) >> 1; const word_t size_h = size - size_l; const word_t work_size = size_l * 2 + 1; word_t* d = res; word_t* c2 = work; work += work_size; (void) _abs_sub(a, size_l, a + size_l, size_h, d); _square_toom22(d, size_l, c2, work); _square_toom22(a, size_l, res, work); // c0 _square_toom22(a + size_l, size_h, res + size_l * 2, work); // c1 word_t c2_size = size_l * 2; _bigint_t c2_ = _bigint_t(1, c2, c2_size); _bigint_t c1_ = _bigint_t(0, res + size_l * 2, size_h * 2); _bigint_t c0_ = _bigint_t(0, res, size_l * 2); _bigint_t r_ = _bigint_t(0, res + size_l, size_l + size_h * 2); _add_ar(c2_, c0_, 0, c2_); _add_ar(c2_, c1_, 0, c2_); _add_ar(r_, c2_, 0, r_); work -= work_size; } static void _square_basecase( const word_t* a, word_t a_size, word_t* res) { std::fill(res, res + 2 * a_size, word_t(0)); for (word_t i = 0; i < a_size; ++i) { word_t alpha = a[i]; if (alpha == 0) { continue; } dword_t c = 0; for (word_t j = i + 1; j < a_size; ++j) { dword_t s = dword_t(alpha) * a[j] + res[i + j] + word_t(c); res[i + j] = s; c = s >> W_BITS; } if (c) { res[i + a_size] = c; } } word_t c = 0; for (word_t i = 0; i < a_size; ++i) { dword_t a2 = dword_t(a[i]) * a[i] + c; dword_t t = (dword_t(res[2 * i]) << 1) + word_t(a2); res[2 * i] = t; t = (dword_t(res[2 * i + 1]) << 1) + word_t(t >> W_BITS) + word_t(a2 >> W_BITS); res[2 * i + 1] = t; c = t >> W_BITS; } } // div-mod static void _udivmod(const BigInt& a, const BigInt& b, BigInt& qs, BigInt& rs) { // Modern computer arithmetic if (_ucomp(a, b, false)) { // fixed qs = BigInt(); rs = BigInt(a); return; } if (b.is_zero()) { assert(0); } if (b.is_one()) { qs = BigInt(a); rs = BigInt(); return; } word_t b_bit_len = b.bit_length(); word_t shamt = (W_BITS - b_bit_len) & (W_BITS - 1); BigInt numer = a << shamt; BigInt denom = b << shamt; word_t n_size = numer.size(); word_t d_size = denom.size(); word_t ofs = n_size - d_size; if (!_shifted_ucomp(numer, ofs, denom, false)) { qs.resize(ofs + 1); qs[ofs] = 1; _shifted_usub_unnorm(numer, ofs, denom); } else { qs.resize(ofs); } word_t d = denom[d_size - 1]; auto fd = FastDiv21(d); BigInt tmp = BigInt(d_size + 1, 0); for (int i = ofs - 1; i >= 0; --i) { word_t q; if (numer[i + d_size] >= d) { q = ~word_t(0); } else { q = ((dword_t(numer[i + d_size]) << W_BITS) | numer[i + d_size - 1]) / fd; } if (q > 0) { tmp[d_size] = _umul1_add(denom.data(), denom.size(), q, tmp.data(), 0); while (_shifted_ucomp(numer, i, tmp, false)) { _usub_unnorm(tmp, denom); q -= 1; } _shifted_usub_unnorm(numer, i, tmp); } qs[i] = q; numer.resize(numer.size() - 1); } numer.normalize(); rs = BigInt(numer); if (shamt) { rs >>= shamt; } } static void _fast_div32( const BigInt& a21, const BigInt& a0, const BigInt& b10, const BigInt& b1, const BigInt& b0, word_t n, BigInt& q, BigInt& r) { BigInt a2 = a21 >> n; if (_ucomp(a2, b1, false)) { // fixed _fast_div21(a21, b1, n, q, r); } else { q = mask(n); r = a21; r += b1; r -= b1 << n; } r <<= n; r |= a0; r -= b0 * q; // bottleneck while (r < 0) { q -= 1; r += b10; } } static void _fast_div21(const BigInt& a, const BigInt& b, word_t n, BigInt& q, BigInt& r) { if (_ucomp(a, b, false)) { q = BigInt(); r = BigInt(a); return; } if (n <= DIVMOD_THRESHOLD) { return _udivmod(a, b, q, r); } word_t ofs = n & 1; word_t nh = (n + ofs) >> 1; BigInt b10 = BigInt(b); BigInt b1 = b >> (nh - ofs); BigInt b0 = b.subinteger(0, nh - ofs); BigInt a32 = a >> n; BigInt a1 = a.subinteger(nh - ofs, n); BigInt a0 = a.subinteger(0, nh - ofs); if (ofs) { b0 <<= 1; a0 <<= 1; b10 <<= 1; } BigInt tmp_r; BigInt q0; _fast_div32(a32, a1, b10, b1, b0, nh, q, tmp_r); // fixed _fast_div32(tmp_r, a0, b10, b1, b0, nh, q0, r); q <<= nh; q |= q0; if (ofs) { r >>= 1; } } static void _fast_udivmod(const BigInt& a, const BigInt& b, BigInt& q, BigInt& r) { if (_ucomp(a, b, false)) { q = BigInt(); r = BigInt(a); return; } word_t m = a.bit_length(); word_t n = b.bit_length(); if (n <= DIVMOD_THRESHOLD) { return _udivmod(a, b, q, r); } word_t q_bit_len = m - n + 1; word_t nblock = 1 + m / n; q.set_zero(); // fixed q.resize(bits_to_words(q_bit_len)); BigInt tmp_q, tmp_r; r = a.subinteger((nblock - 1) * n, nblock * n); for (int i = nblock - 2; i >= 0; --i) { r <<= n; r |= a.subinteger(i * n, (i + 1) * n); _fast_div21(r, b, n, tmp_q, tmp_r); // fixed _lshifted_lor(q, tmp_q, i * n); r = tmp_r; } q.normalize(); } static void _exact_div(word_t* a, word_t size, word_t d) { assert(d & 1); word_t inv = bits::mul_inv(d); word_t c = 0; for (word_t i = 0; i < size; ++i) { a[i] = inv * (a[i] - c); c = (dword_t(a[i]) * d) >> W_BITS; } } // logical static void _land(word_t* a, const word_t* b, word_t b_size) { // assert(a_size == b_size); for (word_t i = 0; i < b_size; ++i) { a[i] &= b[i]; } } static void _lxor(word_t* a, const word_t* b, word_t b_size) { // assert(a_size >= b.size); for (word_t i = 0; i < b_size; ++i) { a[i] ^= b[i]; } } static void _lor(word_t* a, const word_t* b, word_t b_size) { // assert(a_size >= b_size); for (word_t i = 0; i < b_size; ++i) { a[i] |= b[i]; } } static void _lshifted_lor(BigInt& lhs, const BigInt& rhs, const word_t shamt) { word_t q = shamt / W_BITS; word_t r = shamt % W_BITS; BigInt shifted_rhs = rhs << r; _lor(lhs.data() + q, shifted_rhs.data(), shifted_rhs.size()); } static void _lshift(const word_t* n, const word_t size, word_t shamt, word_t* res, const word_t res_size) { word_t q = shamt / W_BITS; word_t r = shamt % W_BITS; if (r) { word_t c = n[size - 1]; if (res_size > size + q) { res[res_size - 1] = c >> (W_BITS - r); } c <<= r; for (int i = size - 1; i >= 1; --i) { word_t t = n[i - 1]; res[i + q] = c | (t >> (W_BITS - r)); c = t << r; } res[q] = c; } else { for (int i = size - 1; i >= 0; --i) { res[i + q] = n[i]; } } for (word_t i = 0; i < q; ++i) { res[i] = 0; } } static void _rshift(const word_t* n, const word_t size, word_t shamt, word_t* res, const word_t res_size) { word_t q = shamt / W_BITS; word_t r = shamt % W_BITS; if (r) { word_t c = n[q] >> r; for (word_t i = 0; i < res_size - 1; ++i) { word_t t = n[i + q + 1]; res[i] = c | (t << (W_BITS - r)); c = t >> r; } if (res_size + q < size) { c |= n[size - 1] << (W_BITS - r); } res[res_size - 1] = c; } else { for (word_t i = 0; i < size - q; ++i) { res[i] = n[i + q]; } } } // functions static void _fib(word_t n, BigInt& a, BigInt& b) { if (n <= 2) { a = BigInt((n + 1) >> 1); b = BigInt((n + 2) >> 1); return; } _fib((n >> 1) - 1, a, b); a *= a; b *= b; BigInt c = b + a; BigInt d = (b << 2) - a + ((n >> 1) & 1 ? -2 : 2); if (n & 1) { a = d; b = (d << 1) - c; } else { a = d - c; b = d; } } static void _isqrt(const BigInt& a, BigInt& s, BigInt& r) { if (a < (uint64(1) << 52)) { s = BigInt(word_t(std::sqrt(uint64(a)))); r = a - s * s; return; } word_t n = a.bit_length(); word_t ofs = (n - 1) & 2; if (ofs > 0) { n += 2; } word_t nq = (n + 1) >> 2; word_t nh = nq << 1; BigInt ah = a >> (nh - ofs); BigInt al = a.subinteger(0, nh - ofs); BigInt a1 = al >> (nq - ofs); BigInt a0 = al.subinteger(0, nq - ofs); if (ofs > 0) { a0 <<= ofs; } BigInt s1, r1; _isqrt(ah, s1, r1); BigInt q, u; divmod((r1 << nq) + a1, s1 << 1, q, u); s = (s1 << nq) + q; r = (u << nq) + a0 - q * q; if (r < 0) { r += (s << 1); r -= 1; s -= 1; } if (ofs > 0) { r >>= 2; if (s & 1) { s >>= 1; r += s; r += 1; } else { s >>= 1; } } } int sign_; container_t n_; }; int main() { using namespace std; uint32 n; while (~scanf("%u", &n)) { if (n == 2) { cout << "3\nINF" << endl; } else { BigInt ans; if (n & 1) { ans = BigInt::fib(n); } else { BigInt a, b; BigInt::fib(n / 2 - 1, a, b); ans = a * b; ans <<= 1; } cout << n << endl; cout << ans << endl; } } return 0; }