#include #include #include #include #include #include #ifndef __GNUC__ #include #endif #if !defined(MY_LOCAL_RUN) && !defined(NDEBUG) #define NDEBUG #endif #include #include #include #include using namespace std; #if defined(_M_X64) || defined(__amd64__) #define MULTIPREC_64BIT #endif #if defined(MY_LOCAL_RUN) #define MP_TEST #endif namespace uint_util { #if defined(_MSC_VER) && !defined(__clang__) #pragma section(".mycode", execute) __declspec(allocate(".mycode")) const unsigned char udiv128Data[] = {0x48, 0x89, 0xD0, 0x48, 0x89, 0xCA, 0x49, 0xF7, 0xF0, 0xC3}; uint64_t(__fastcall *udiv128)(uint64_t numhi, uint64_t numlo, uint64_t den) = (uint64_t(__fastcall *)(uint64_t, uint64_t, uint64_t))(void*)(udiv128Data); #endif template struct Utils {}; #ifdef MULTIPREC_64BIT template<> struct Utils { #if defined(_MSC_VER) && !defined(__clang__) static void umul_full(uint64_t a, uint64_t b, uint64_t *hi, uint64_t *lo) { *lo = _umul128(a, b, hi); } static uint64_t udiv_full(uint64_t num_hi, uint64_t num_lo, uint64_t den) { return udiv128(num_hi, num_lo, den); } #else static void umul_full(uint64_t a, uint64_t b, uint64_t *hi, uint64_t *lo) { const auto c = (unsigned __int128)a * b; *lo = (uint64_t)c; *hi = (uint64_t)(c >> 64); } static uint64_t udiv_full(uint64_t num_hi, uint64_t num_lo, uint64_t den) { uint64_t quot; asm("div %3":"=a"(quot):"a"(num_lo),"d"(num_hi),"r"(den)); return quot; } #endif static uint64_t div_preinv_invert(uint64_t b) { assert(b >> 63 & 1); return udiv_full(~0ULL - b, ~0ULL, b); } static uint64_t addcarry(uint64_t a, uint64_t b, uint64_t *c) { uint64_t t = a + b; bool ct = t < a; uint64_t u = t + *c; bool cu = u < t; *c = ct + cu; return u; } static uint64_t subborrow(uint64_t a, uint64_t b, uint64_t *c) { uint64_t t = a - b; bool ct = t > a; uint64_t u = t - *c; bool cu = u > t; *c = ct + cu; return u; } static void div_preinv(uint64_t num_hi, uint64_t num_lo, uint64_t den, uint64_t inv, uint64_t *quot, uint64_t *rem) { uint64_t qh, ql; umul_full(num_hi, inv, &qh, &ql); ql += num_lo; qh += ql < num_lo; qh += num_hi + 1; uint64_t r = num_lo - qh * den; uint64_t mask = -(r > ql); qh += mask; r += mask & den; if(r >= den) { qh += 1; r -= den; } *quot = qh; *rem = r; } static int count_leading_zeros(uint64_t x) { #ifndef __GNUC__ unsigned long pos; _BitScanReverse64(&pos, x); return 63 - pos; #else return __builtin_clzll(x); #endif } static int count_trailing_zeros(uint64_t x) { #ifndef __GNUC__ unsigned long pos; _BitScanForward64(&pos, x); return pos; #else return __builtin_ctzll(x); #endif } }; #endif template<> struct Utils { }; //for tests template<> struct Utils { static void umul_full(uint8_t a, uint8_t b, uint8_t *hi, uint8_t *lo) { auto c = (unsigned)a * b; *hi = (uint8_t)(c >> 8), *lo = (uint8_t)c; } static uint8_t udiv_full(uint8_t num_hi, uint8_t num_lo, uint8_t den) { assert(((unsigned)num_hi << 8 | num_lo) / den < 0x100); return (uint8_t)(((unsigned)num_hi << 8 | num_lo) / den); } static uint8_t div_preinv_invert(uint8_t b) { assert(b >> 7 & 1); return udiv_full(0xff - b, 0xff, b); } static uint8_t addcarry(uint8_t a, uint8_t b, uint8_t *c) { unsigned t = (unsigned)a + b + *c; *c = (uint8_t)(t >> 8); return (uint8_t)t; } static uint8_t subborrow(uint8_t a, uint8_t b, uint8_t *c) { unsigned t = (unsigned)a - b - *c; *c = -(uint8_t)(t >> 8); return (uint8_t)t; } static void div_preinv(uint8_t num_hi, uint8_t num_lo, uint8_t den, uint8_t inv, uint8_t *quot, uint8_t *rem) { unsigned num = ((unsigned)num_hi << 8 | num_lo); assert(num / den < 0x100); *quot = (uint8_t)(num / den); *rem = (uint8_t)(num % den); } static int count_leading_zeros(uint8_t x) { assert(x != 0); int r = 1; while((unsigned)x >> r) ++ r; return 8 - r; } static int count_trailing_zeros(uint64_t x) { assert(x != 0); int r = 0; while(x >> r & 1) ++ r; return r; } }; } namespace multiprec { class BigInt { public: #ifdef MP_TEST typedef uint8_t W; #else typedef uint64_t W; #endif typedef uint_util::Utils Util; enum { W_BITS = sizeof(W) * 8 }; enum : W { W_MAX = W(-1) }; private: std::vector _words; bool _sign; public: BigInt(): _words(), _sign(false) {} BigInt(W x) { set(x); } BigInt(int x) { set(x); } void clear() { _sign = false; _words.clear(); } bool zero() const { return _words.empty(); } bool negative() const { return _sign; } private: W *data() { return _words.data(); } const W *data() const { return _words.data(); } int size() const { return static_cast(_words.size()); } void resize(int n) { _words.resize(n); } bool normalized() const { return _words.empty() ? !_sign : _words.back() != 0; } void normalize_size() { while(!_words.empty() && _words.back() == 0) _words.pop_back(); } public: void set(W w) { if(w == 0) { clear(); } else { _sign = false; _words.assign({w}); } } void set(int s) { if(s > 0) { _sign = false; _words.assign({static_cast(s)}); } else if(s < 0) { _sign = true; _words.assign({static_cast(-s)}); } else { clear(); } } void set(const BigInt &that) { *this = that; } template void set_words(It begin_it, It end_it) { _words.assign(begin_it, end_it); normalize_size(); _sign = false; } bool operator==(const BigInt &that) const { return compare_to(that) == 0; } bool operator!=(const BigInt &that) const { return compare_to(that) != 0; } bool operator< (const BigInt &that) const { return compare_to(that) < 0; } bool operator> (const BigInt &that) const { return compare_to(that) > 0; } bool operator<=(const BigInt &that) const { return compare_to(that) <= 0; } bool operator>=(const BigInt &that) const { return compare_to(that) >= 0; } int compare_to(const BigInt &that) const { bool asign = _sign, bsign = that._sign; if(asign == bsign) { if(!asign) return compare_unsigned(*this, that); else return compare_unsigned(that, *this); } else { if(!asign) return +1; else return -1; } } static int compare_unsigned(const BigInt &a, const BigInt &b) { assert(a.normalized() && b.normalized()); int an = a.size(), bn = b.size(); if(an > bn) { return +1; } else if(an < bn) { return -1; } else { return _compare(a.data(), b.data(), an); } } BigInt &operator+=(const BigInt &that) { add_signed(*this, *this, _sign, that, that._sign); return *this; } BigInt operator+(const BigInt &that) const { BigInt res; add_signed(res, *this, _sign, that, that._sign); return res; } BigInt &operator-=(const BigInt &that) { add_signed(*this, *this, _sign, that, !that._sign); return *this; } BigInt operator-(const BigInt &that) const { BigInt res; add_signed(res, *this, _sign, that, !that._sign); return res; } BigInt operator-() const { BigInt res = *this; if(!zero()) res._sign = !res._sign; return res; } private: static void add_signed(BigInt &res, const BigInt &a, bool a_sign, const BigInt &b, bool b_sign) { if(a_sign == b_sign) { add_unsigned(res, a, b); res._sign = a_sign && !res.zero(); } else if(b_sign) { res._sign = subtract_unsigned(res, a, b); } else { res._sign = subtract_unsigned(res, b, a); } assert(res.normalized()); } static void add_unsigned(BigInt &res, const BigInt &a, const BigInt &b) { int an = a.size(), bn = b.size(); if(an < bn) return add_unsigned(res, b, a); res.resize(an + 1); res._words[an] = _add(res.data(), a.data(), an, b.data(), bn); res.normalize_size(); } static bool subtract_unsigned(BigInt &res, const BigInt &a, const BigInt &b) { int an = a.size(), bn = b.size(); if(compare_unsigned(a, b) >= 0) { subtract_unsigned_guarded(res, a, b); return false; } else { subtract_unsigned_guarded(res, b, a); return true; } } //a >= b static void subtract_unsigned_guarded(BigInt &res, const BigInt &a, const BigInt &b) { int an = a.size(), bn = b.size(); assert(an >= bn); res.resize(an); W borrow = _subtract(res.data(), a.data(), an, b.data(), bn); assert(borrow == 0); res.normalize_size(); } public: BigInt &operator*=(const BigInt &that) { multiply(*this, *this, that); return *this; } BigInt operator*(const BigInt &that) const { BigInt r; multiply(r, *this, that); return r; } static void multiply(BigInt &res, const BigInt &a, const BigInt &b) { int an = a.size(), bn = b.size(); if(an < bn) return multiply(res, b, a); if(&res == &a || &res == &b) { BigInt ws; multiply(ws, a, b); res = ws; return; } if(bn == 0) { res.clear(); } else { res.resize(an + bn); _multiply(res.data(), a.data(), an, b.data(), bn); res.normalize_size(); res._sign = a._sign != b._sign; } } public: BigInt &operator/=(const BigInt &that) { BigInt dummy; divide_trunc(*this, dummy, *this, that); return *this; } BigInt operator/(const BigInt &that) const { BigInt res, dummy; divide_trunc(res, dummy, *this, that); return res; } BigInt &operator%=(const BigInt &that) { BigInt dummy; divide_trunc(dummy, *this, *this, that); return *this; } BigInt operator%(const BigInt &that) const { BigInt res, dummy; divide_trunc(dummy, res, *this, that); return res; } static void divide_trunc(BigInt ", BigInt &rem, const BigInt &a, const BigInt &b) { bool asign = a._sign, bsign = b._sign; divide_floor_unsigned(quot, rem, a, b); quot._sign = asign != bsign && !quot.zero(); rem._sign = asign && !rem.zero(); } private: static void divide_floor_unsigned(BigInt ", BigInt &rem, const BigInt &a, const BigInt &b) { if(" == &rem || " == &b) { BigInt ws; divide_floor_unsigned(ws, rem, a, b); quot = ws; return; } if(&rem == &b) { BigInt ws; divide_floor_unsigned(quot, ws, a, b); rem = ws; return; } assert(a.normalized() && b.normalized()); int an = a.size(), bn = b.size(); if(an < bn) { rem = a; quot = BigInt(); return; } if(bn == 0) { std::cerr << "divide by 0" << std::endl; abort(); } rem = a; quot.resize(an - bn + 1); _divide(quot.data(), rem.data(), an, b.data(), bn); quot.normalize_size(); rem.resize(bn); rem.normalize_size(); } private: struct BaseData { W base; W bigbase_shifted; W invbigbase; int shift; int bigbase_size; double log_base_ratio; std::vector powers; //bigbase_shifted^(2^k) }; public: std::string to_string() const { assert(normalized()); static BaseData baseinfo_10 = _get_basedata(10); int n = size(); if(n == 0) return "0"; std::string str(_sign + _estimate_size_in_base(data(), n, baseinfo_10), '?'); if(_sign) str[0] = '-'; unique_ptr ws(new W[n]); _copy(ws.get(), data(), n); const auto begin_it = str.begin() + (_sign ? 1 : 0); const auto end_it = _from_binary(begin_it, ws.get(), n, baseinfo_10); for(auto it = begin_it; it != end_it; ++ it) *it += '0'; str.resize(end_it - str.begin()); return str; } friend std::ostream &operator<<(std::ostream &o, const BigInt &a) { return o << a.to_string(); } private: static void _copy(W *res, const W *a, int n); static void _fill_zero(W *res, int n); static int _compare(const W *a, const W *b, int n); static W _add_1(W *a, int an, W b); static W _subtract_1(W *a, int an, W b); static W _add(W *res, const W *a, int an, const W *b, int bn); static W _subtract(W *res, const W *a, int an, const W *b, int bn); static W _bitshift_left(W *res, const W *a, int n, int shift); static W _bitshift_right(W *res, const W *a, int n, int shift); static W _multiply_1(W *res, const W *a, int n, W b); static W _multiply_add_1(W *a, const W *b, int n, W c); static W _multiply_subtract_1(W *a, const W *b, int n, W c); enum { #ifdef MP_TEST KARATSUBA_THRESHOLD = 3 #else KARATSUBA_THRESHOLD = 8 #endif }; static void _multiply(W *res, const W *a, int an, const W *b, int bn); static void _square(W *res, const W *a, int an); static void _multiply_basecase(W *res, const W *a, int an, const W *b, int bn); static void _multiply_karatsuba(W *res, const W *a, int an, const W *b, int bn, W *ws); static bool _compare_subtract(W *res, int resn, const W *a, int an, const W *b, int bn); static W _divide_1(W *res_q, const W *a, int n, W b); static W _divide_1_preinv_noshift(W *res_q, W a_n, const W *a, int n, W b, W invb); static W _divide_1_preinv_shift(W *res_q, const W *a, int n, W b, W invb, int shift); static void _divide(W *res_q, W *a, int an, const W *b, int bn); static void _divide_preinv(W *res_q, int qn, W *a, int an, const W *b, int bn, W invbhi); static void _divide_basecase_preinv(W *res_q, int qn, W *a, int an, const W *b, int bn, W invbhi); static void _divide_blockwise(W *res_q, int qn, W *a, int an, const W *b, int bn, W invbhi, W *ws); static void _divide_preinv_invert(W *res, int resn, const W *a, int an); static BaseData _get_basedata(W base); static void _precompute_base_powers(int n, BaseData &data); using CharT = char; using StrIt = std::string::iterator; static StrIt _from_binary(StrIt res, W *a, int n, BaseData &data); static StrIt _from_binary_basecase(StrIt res, W *a, int n, const BaseData &data); static StrIt _from_binary_rec(StrIt res, W *a, int n, const BaseData &data); static size_t _estimate_size_in_base(const W *a, int n, const BaseData &data); }; void BigInt::_copy(W *res, const W *a, int n) { for(int i = 0; i < n; ++ i) res[i] = a[i]; } void BigInt::_fill_zero(W *res, int n) { for(int i = 0; i < n; ++ i) res[i] = 0; } int BigInt::_compare(const W *a, const W *b, int n) { for(int i = n - 1; i >= 0; -- i) { if(a[i] != b[i]) return a[i] > b[i] ? 1 : -1; } return 0; } BigInt::W BigInt::_add_1(W *a, int an, W b) { assert(an > 0); W c = a[0] + b; a[0] = c; if(c >= b) return 0; for(int i = 1; i < an; ++ i) { if((++ a[i]) != 0) return 0; } return 1; } BigInt::W BigInt::_subtract_1(W *a, int an, W b) { assert(an > 0); W c = a[0] - b; a[0] = c; if(c < (W)-b) return 0; for(int i = 1; i < an; ++ i) { if((-- a[i]) != W_MAX) return 0; } return 1; } BigInt::W BigInt::_add(W *res, const W *a, int an, const W *b, int bn) { assert(an >= bn); W carry = 0; for(int i = 0; i < bn; ++ i) res[i] = Util::addcarry(a[i], b[i], &carry); if(res == a) { if(carry == 0 || an == bn) return carry; else return _add_1(res + bn, an - bn, 1); } else { for(int i = bn; i < an; ++ i) res[i] = Util::addcarry(a[i], 0, &carry); return carry; } } BigInt::W BigInt::_subtract(W *res, const W *a, int an, const W *b, int bn) { assert(an >= bn); W borrow = 0; for(int i = 0; i < bn; ++ i) res[i] = Util::subborrow(a[i], b[i], &borrow); if(res == a) { if(borrow == 0 || an == bn) return borrow; else return _subtract_1(res + bn, an - bn, 1); } else { for(int i = bn; i < an; ++ i) res[i] = Util::subborrow(a[i], 0, &borrow); } return borrow; } BigInt::W BigInt::_bitshift_left(W *res, const W *a, int n, int shift) { assert(n > 0 && 0 < shift && shift < W_BITS); W carry = 0; for(int i = 0; i < n; ++ i) { W a_i = a[i]; res[i] = a_i << shift | carry; carry = a_i >> (W_BITS - shift); } return carry; } BigInt::W BigInt::_bitshift_right(W *res, const W *a, int n, int shift) { assert(n > 0 && 0 < shift && shift < W_BITS); W carry = 0; for(int i = n - 1; i >= 0; -- i) { W a_i = a[i]; res[i] = a_i >> shift | carry; carry = a_i << (W_BITS - shift); } return carry; } BigInt::W BigInt::_multiply_1(W *res, const W *a, int an, W b) { W carry = 0; for(int i = 0; i < an; ++ i) { W lo, hi; Util::umul_full(a[i], b, &hi, &lo); W c = lo + carry; res[i] = c; carry = hi + (c < lo); } return carry; } BigInt::W BigInt::_multiply_add_1(W *a, const W *b, int n, W c) { W carry = 0; for(int i = 0; i < n; ++ i) { W lo, hi; Util::umul_full(b[i], c, &hi, &lo); a[i] = Util::addcarry(a[i], lo, &carry); carry += hi; } return carry; } BigInt::W BigInt::_multiply_subtract_1(W *a, const W *b, int n, W c) { W borrow = 0; for(int i = 0; i < n; ++ i) { W lo, hi; Util::umul_full(b[i], c, &hi, &lo); a[i] = Util::subborrow(a[i], lo, &borrow); borrow += hi; } return borrow; } void BigInt::_multiply(W *res, const W *a, int an, const W *b, int bn) { assert(res != a && res != b); assert(an >= bn && bn > 0); if(a == b) { assert(an == bn); _square(res, a, an); return; } if(bn == 1) { res[an] = _multiply_1(res, a, an, b[0]); } else if(bn < KARATSUBA_THRESHOLD) { _multiply_basecase(res, a, an, b, bn); } else { unique_ptr ws(new W[an * 5]); _multiply_karatsuba(res, a, an, b, bn, ws.get()); } } void BigInt::_square(W *res, const W *a, int an) { if(an == 1) { Util::umul_full(a[0], a[0], &res[1], &res[0]); } else if(an < KARATSUBA_THRESHOLD) { //todo _multiply_basecase(res, a, an, a, an); } else { //todo unique_ptr ws(new W[an * 5]); _multiply_karatsuba(res, a, an, a, an, ws.get()); } } void BigInt::_multiply_basecase(W *res, const W *a, int an, const W *b, int bn) { res[an] = _multiply_1(res, a, an, b[0]); for(int i = 1; i < bn; ++ i) res[i + an] = _multiply_add_1(res + i, a, an, b[i]); } void BigInt::_multiply_karatsuba(W *res, const W *a, int an, const W *b, int bn, W *ws) { if(bn < KARATSUBA_THRESHOLD || (an + 1) / 2 >= bn) { //todo: unbalanced case _multiply_basecase(res, a, an, b, bn); return; } int lo = (an + 1) / 2, ha = an - lo, hb = bn - lo; assert(an >= bn && 0 < hb); W *am1 = ws, *bm1 = ws + lo, *rm1 = bm1 + lo, *nws = rm1 + lo * 2; W *r0 = res, *rinf = res + lo * 2; bool rm1_sign = true; rm1_sign ^= _compare_subtract(am1, lo, a, lo, a + lo, ha); rm1_sign ^= _compare_subtract(bm1, lo, b, lo, b + lo, hb); _multiply_karatsuba(r0, a, lo, b, lo, nws); _multiply_karatsuba(rinf, a + lo, ha, b + lo, hb, nws); _multiply_karatsuba(rm1, am1, lo, bm1, lo, nws); W rm1_carry = 0; if(!rm1_sign) rm1_carry += _add(rm1, rm1, lo * 2, r0, lo * 2); else rm1_sign = _compare_subtract(rm1, lo * 2, r0, lo * 2, rm1, lo * 2); if(!rm1_sign) rm1_carry += _add(rm1, rm1, lo * 2, rinf, ha + hb); else rm1_sign = _compare_subtract(rm1, lo * 2, rinf, ha + hb, rm1, lo * 2); if(rm1_sign) { _subtract(res + lo, res + lo, an + bn - lo, rm1, lo * 2); } else { _add(res + lo, res + lo, an + bn - lo, rm1, lo * 2); if(rm1_carry > 0) _add_1(res + lo * 3, ha + hb - lo, rm1_carry); } } bool BigInt::_compare_subtract(W *res, int resn, const W *a, int an, const W *b, int bn) { int aan = an; while(aan > 0 && a[aan - 1] == 0) -- aan; int abn = bn; while(abn > 0 && b[abn - 1] == 0) -- abn; if(aan < abn || (aan == abn && _compare(a, b, aan) < 0)) { _subtract(res, b, abn, a, aan); _fill_zero(res + abn, resn - abn); return true; } else { _subtract(res, a, aan, b, abn); _fill_zero(res + aan, resn - aan); return false; } } BigInt::W BigInt::_divide_1(W *res_q, const W *a, int n, W b) { assert(n > 0); int shift = Util::count_leading_zeros(b); W invb = Util::div_preinv_invert(b << shift); W rem = _divide_1_preinv_shift(res_q, a, n, b << shift, invb, shift); return rem >> shift; } BigInt::W BigInt::_divide_1_preinv_noshift(W *res_q, W a_n, const W *a, int n, W b, W invb) { assert(a_n < b); W rem = a_n; for(int i = n - 1; i >= 0; -- i) Util::div_preinv(rem, a[i], b, invb, &res_q[i], &rem); return rem; } BigInt::W BigInt::_divide_1_preinv_shift(W *res_q, const W *a, int n, W b, W invb, int shift) { if(shift == 0) return _divide_1_preinv_noshift(res_q, 0, a, n, b, invb); assert(n > 0); W a_i = a[n - 1], rem = a_i >> (W_BITS - shift); for(int i = n - 1; i > 0; -- i) { W a_im1 = a[i - 1]; W a_shifted = a_i << shift | a_im1 >> (W_BITS - shift); Util::div_preinv(rem, a_shifted, b, invb, &res_q[i], &rem); a_i = a_im1; } Util::div_preinv(rem, a_i << shift, b, invb, &res_q[0], &rem); return rem; } void BigInt::_divide(W *res_q, W *a, int an, const W *b, int bn) { assert(res_q != a && res_q != b && a != b && an >= bn && bn > 0); if(bn == 1) { a[0] = _divide_1(res_q, a, an, b[0]); _fill_zero(a + 1, an - 1); } else { int qn = an - bn + 1; int shift = Util::count_leading_zeros(b[bn - 1]); if(shift == 0) { W invbhi = Util::div_preinv_invert(b[bn - 1] << shift); _divide_preinv(res_q, qn, a, an, b, bn, invbhi); } else { unique_ptr ws(new W[an + 1 + bn]); W *shifted_a = ws.get(), *shifted_b = shifted_a + (an + 1); _bitshift_left(shifted_b, b, bn, shift); shifted_a[an] = _bitshift_left(shifted_a, a, an, shift); W invbhi = Util::div_preinv_invert(shifted_b[bn - 1]); _divide_preinv(res_q, qn, shifted_a, an + 1, shifted_b, bn, invbhi); // assert(shifted_a[an] == 0); _bitshift_right(a, shifted_a, an, shift); } } } void BigInt::_divide_preinv(W *res_q, int qn, W *a, int an, const W *b, int bn, W invbhi) { unique_ptr ws(new W[bn * 5]); _divide_blockwise(res_q, qn, a, an, b, bn, invbhi, ws.get()); } void BigInt::_divide_basecase_preinv(W *res_q, int qn, W *a, int an, const W *b, int bn, W invbhi) { assert(res_q != a && b != res_q && b != a && an >= bn && bn > 0 && (b[bn - 1] >> (W_BITS - 1) & 1)); assert(qn <= an - bn + 1); W bhi = b[bn - 1]; for(int i = an - bn; i >= 0; -- i) { W numhi = i == an - bn ? 0 : a[i + bn], q; if(numhi >= bhi) { q = W_MAX; } else { W dummy; Util::div_preinv(numhi, a[i + bn - 1], bhi, invbhi, &q, &dummy); } //[a / b] <= q <= [a / b] + 2 if(q != 0) { W borrow = _multiply_subtract_1(a + i, b, bn, q); if(numhi < borrow) { numhi -= borrow; -- q; W carry = _add(a + i, a + i, bn, b, bn); if((numhi += carry) >= carry) { -- q; numhi += _add(a + i, a + i, bn, b, bn); } } else { numhi -= borrow; } } assert(numhi == 0); if(i < an - bn) a[i + bn] = 0; if(i < qn) res_q[i] = q; } } void BigInt::_divide_blockwise(W *res_q, int qn, W *a, int an, const W *b, int bn, W invbhi, W *ws) { assert(an >= bn && bn > 0); assert(b[bn - 1] >> (W_BITS - 1) & 1); assert(qn <= an - bn + 1); if(bn < KARATSUBA_THRESHOLD) { _divide_basecase_preinv(res_q, qn, a, an, b, bn, invbhi); return; } int block = (bn + 1) / 2; W *q = ws, *prod = q + (block + 1), *num = prod; W *nws = ws + block * 2 + 1 + bn; int initi = 0; while(initi + block < an - bn + 1) initi += block; for(int i = initi; i >= 0; i -= block) { int hn; if(i == initi) { hn = an - (i + bn); }else { hn = block; } _copy(num, a + (i + bn - block), block + hn); _divide_blockwise(q, hn + 1, num, block + hn, b + (bn - block), block, invbhi, nws); if(q[hn] > 0) { if(hn == block) { for(int j = 0; j < hn; ++ j) q[j] = W_MAX; q[hn] = 0; } else { ++ hn; } } if(hn > 0) { _multiply_karatsuba(prod, b, bn, q, hn, nws); int prodn = bn + hn; if(i + prodn > an) { while(prod[prodn - 1] > 0) { _subtract_1(q, hn, 1); _subtract(prod, prod, prodn, b, bn); } -- prodn; } if(_subtract(a + i, a + i, prodn, prod, prodn) > 0) { _subtract_1(q, hn, 1); if(_add(a + i, a + i, prodn, b, bn) == 0) { _subtract_1(q, hn, 1); _add(a + i, a + i, prodn, b, bn); } } } if(qn > i) _copy(res_q + i, q, std::min(block, qn - i)); } } //todo void BigInt::_divide_preinv_invert(W *res, int resn, const W *a, int an) { /* assert(an > 0 && (a[an - 1] >> (W_BITS - 1) & 1)); unique_ptr buf(new W[resn * 4]); W *tmp1 = buf.get(), *tmp2 = tmp1 + resn * 2; W inv_a0 = Util::div_preinv_invert(a[0]); int curn = 1; while(curn < resn) { //res / W^{(an - 1) + curn} int nextn = std::min(resn, curn * 2); _square(tmp1, res, curn); int htmp1 = curn * 2, ha = std::min(nextn, an); _multiply(tmp2, tmp1 + (curn * 2 - htmp1), htmp1, a + (an - ha), ha); int diff = nextn - curn; _bitshift_left(res + diff, res, curn, 1); _subtract(res, res, nextn, tmp2, nextn); } */ } BigInt::StrIt BigInt::_from_binary(StrIt res, W *a, int n, BaseData &data) { _precompute_base_powers(n, data); StrIt it = _from_binary_rec(res, a, n, data); std::reverse(res, it); return it; } BigInt::StrIt BigInt::_from_binary_rec(StrIt res, W *a, int n, const BaseData &data) { if(n == 0) return res; if(n < KARATSUBA_THRESHOLD) return _from_binary_basecase(res, a, n, data); int j = 0; while(j + 1 < (int)data.powers.size() && data.powers[j + 1].size() < n) ++ j; BigInt x, quot, rem; x.set_words(a, a + n); const BigInt &pow = data.powers[j]; size_t pow_size = (size_t)data.bigbase_size << j; BigInt::divide_floor_unsigned(quot, rem, x, pow); StrIt mid = res + pow_size; StrIt it = _from_binary_rec(res, rem.data(), rem.size(), data); for(; it != mid; ++ it) *it = static_cast(0); return _from_binary_rec(it, quot.data(), quot.size(), data); } BigInt::StrIt BigInt::_from_binary_basecase(StrIt res, W *a, int n, const BaseData &data) { const W base = data.base; const W bigbase_shifted = data.bigbase_shifted; const W invbigbase = data.invbigbase; const int shift = data.shift; const int bigbase_size = data.bigbase_size; assert(n > 0); auto p = res; int remn = n; while(remn > 1 || a[0] >= bigbase_shifted >> shift) { W rem, frac, dummy, digit; rem = _divide_1_preinv_shift(a, a, n, bigbase_shifted, invbigbase, shift); Util::div_preinv(rem, 0, bigbase_shifted, invbigbase, &frac, &dummy); ++ frac; for(int i = bigbase_size - 1; i >= 0; -- i) { Util::umul_full(frac, base, &digit, &frac); p[i] = static_cast(digit); } p += bigbase_size; remn -= a[remn - 1] == 0; } W frac, dummy, digit; Util::div_preinv(a[0] << shift, 0, bigbase_shifted, invbigbase, &frac, &dummy); ++ frac; int last_digits = bigbase_size; for(;; -- last_digits) { Util::umul_full(frac, base, &digit, &frac); if(digit != 0) break; } p[last_digits - 1] = static_cast(digit); for(int i = last_digits - 2; i >= 0; -- i) { Util::umul_full(frac, base, &digit, &frac); p[i] = static_cast(digit); } p += last_digits; return p; } BigInt::BaseData BigInt::_get_basedata(W base) { assert(base > 1); W bigbase_shifted = base; int bigbase_size = 1; const W limit = W_MAX / base; while(bigbase_shifted <= limit) bigbase_shifted *= base, ++ bigbase_size; int shift = Util::count_leading_zeros(bigbase_shifted); bigbase_shifted <<= shift; const W invbigbase = Util::div_preinv_invert(bigbase_shifted); BaseData res; res.base = base; res.bigbase_shifted = bigbase_shifted; res.invbigbase = invbigbase; res.shift = shift; res.bigbase_size = bigbase_size; if((base & (base - 1)) == 0) res.log_base_ratio = Util::count_trailing_zeros(base); else res.log_base_ratio = log(2.0) / std::log(double(base)); return res; } void BigInt::_precompute_base_powers(int n, BaseData &data) { std::vector &powers = data.powers; if(powers.empty()) powers.push_back(BigInt(data.bigbase_shifted >> data.shift)); while(powers.back().size() * 2 < n) { const auto &a = powers.back(); powers.push_back(a * a); } } size_t BigInt::_estimate_size_in_base(const W *a, int n, const BaseData &data) { if(n == 0) return 1; size_t bit_size = static_cast(n) * W_BITS - Util::count_leading_zeros(a[n - 1]); return static_cast(ceil(static_cast(bit_size) * data.log_base_ratio)); } } using multiprec::BigInt; pair fib_pair(int n) { if(n <= 2) return make_pair(BigInt((n + 1) / 2), BigInt((n + 2) / 2)); BigInt a, b; tie(a, b) = fib_pair(n / 2 - 1); a *= a; b *= b; BigInt c = b + a; BigInt d = b * 4 - a + (n >> 1 & 1 ? -2 : 2); if(n & 1) return make_pair(d, d * 2 - c); else return make_pair(d - c, d); } struct Xor128 { unsigned x, y, z, w; Xor128(): x(123456789), y(362436069), z(521288629), w(88675123) { } //[0, 2^32) unsigned operator()() { unsigned t = x ^ (x << 11); x = y; y = z; z = w; return w = w ^ (w >> 19) ^ (t ^ (t >> 8)); } //[0, 2^64) unsigned long long nextLL() { unsigned x = (*this)(); unsigned y = (*this)(); return (unsigned long long)x << 32 | y; } //[0, n) unsigned operator()(unsigned n) { unsigned mask = calculateMask(n - 1), x; do { x = (*this)() & mask; }while(x >= n); return x; } //[0, n) signed int operator()(signed int n) { return (*this)((unsigned int)n); } //[L, U] signed int operator()(signed int L, signed int U) { return L + (*this)(U - L + 1); } //[0, n) unsigned long long operator()(unsigned long long n) { unsigned long long mask = calculateMask(n - 1), x; do { x = (*this).nextLL() & mask; }while(x >= n); return x; } //[0, n) signed long long operator()(signed long long n) { return (*this)((unsigned long long)n); } //[L, U] signed long long operator()(signed long long L, signed long long U) { return L + (*this)(U - L + 1); } private: static unsigned calculateMask(unsigned v) { v |= v >> 1; v |= v >> 2; v |= v >> 4; v |= v >> 8; v |= v >> 16; return v; } static unsigned long long calculateMask(unsigned long long v) { v |= v >> 1; v |= v >> 2; v |= v >> 4; v |= v >> 8; v |= v >> 16; v |= v >> 32; return v; } }; BigInt genBigInt(int maxN, Xor128 &xor128, bool sign = false) { int n = xor128(1, maxN); vector v(n); int P = xor128(100), X = xor128(2); for(int i = 0; i < n; ++ i) { if(xor128(100) < P) v[i] = (BigInt::W)xor128(-X, X); else v[i] = (BigInt::W)xor128.nextLL(); } BigInt a; a.set_words(v.begin(), v.end()); if(sign && xor128(2) == 0) a = -a; return a; } void test_add(int N, Xor128 &xor128) { BigInt a, b; a = genBigInt(N, xor128, true); b = genBigInt(N, xor128, true); BigInt sum = a + b; int P = 127; if((((a % P) + (b % P)) % P + P) % P != (sum % P + P) % P) { cerr << a << " + " << b << " != " << sum << endl; } } void test_multiply(int N, Xor128 &xor128) { BigInt a, b; a = genBigInt(N, xor128); b = genBigInt(N, xor128); BigInt prod = a * b; int P = 127; if((a % P) * (b % P) % P != prod % P) { cerr << a << " * " << b << " != " << prod << endl; } } void test_divide(int N, Xor128 &xor128) { BigInt a, b; a = genBigInt(N, xor128); b = genBigInt(N, xor128); if(b.zero()) return; BigInt q, r; BigInt::divide_trunc(q, r, a, b); if(a - b * q != r || !(BigInt() <= r && r < b)) { cerr << a << " `quotRem` " << b << endl; cerr << q << ", " << r << endl; cerr << "aaa" << endl; } } void test_all() { #ifdef NDEBUG cerr << "warning: NDEBUG is enabled" << endl; #endif Xor128 xor128; const int TT = 1000000; for(int k = 3; k < 10; ++ k) { const int N = k == 0 ? 2 : k == 1 ? 4 : k == 2 ? 6 : 1 << k, T = TT / max(1, N >> 3); cerr << "test N = " << N << ", T = " << T << endl; for(int ii = 0; ii < T; ++ ii) { // test_add(N, xor128); test_multiply(N, xor128); test_divide(N, xor128); } } cerr << "end" << endl; } int main() { #ifdef MY_LOCAL_RUN test_all(); #endif int n; while(cin >> n) { if(n == 2) { puts("3\nINF"); } else { BigInt ans; if(n % 2 == 1) { ans = fib_pair(n).first; } else { BigInt a, b; tie(a, b) = fib_pair(n / 2 - 1); ans = a * b * 2; } printf("%d\n", n); std::string s = ans.to_string(); puts(s.c_str()); } } }