#line 2 ".lib/template.hpp" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define allof(obj) (obj).begin(), (obj).end() #define range(i, l, r) for(int i=l;i>1)|y_bit)) #define bit_kth(i, k) ((i >> k)&1) #define bit_highest(i) (i?63-__builtin_clzll(i):-1) #define bit_lowest(i) (i?__builtin_ctzll(i):-1) using ll = long long; using ld = long double; using ul = uint64_t; using pi = std::pair; using pl = std::pair; template using vec = std::vector; using namespace std; template std::ostream &operator<<(std::ostream &dest, const std::pair &p){ dest << p.first << ' ' << p.second; return dest; } template std::ostream &operator<<(std::ostream &dest, const std::vector> &v){ int sz = v.size(); if(sz==0) return dest; for(int i=0;i std::ostream &operator<<(std::ostream &dest, const std::vector &v){ int sz = v.size(); if(sz==0) return dest; for(int i=0;i std::ostream &operator<<(std::ostream &dest, const std::array &v){ if(sz==0) return dest; for(int i=0;i std::ostream &operator<<(std::ostream &dest, const std::set &v){ for(auto itr=v.begin();itr!=v.end();){ dest << *itr; itr++; if(itr!=v.end()) dest << ' '; } return dest; } template std::ostream &operator<<(std::ostream &dest, const std::map &v){ for(auto itr=v.begin();itr!=v.end();){ dest << '(' << itr->first << ", " << itr->second << ')'; itr++; if(itr!=v.end()) dest << '\n'; } return dest; } template vector make_vec(size_t sz, T val){return std::vector(sz, val);} template auto make_vec(size_t sz, Tail ...tail){ return std::vector(tail...))>(sz, make_vec(tail...)); } template vector read_vec(size_t sz){ std::vector v(sz); for(int i=0;i> v[i]; return v; } template auto read_vec(size_t sz, Tail ...tail){ auto v = std::vector(tail...))>(sz); for(int i=0;i(tail...); return v; } void io_init(){ std::cin.tie(nullptr); std::ios::sync_with_stdio(false); } #line 3 ".lib/math/integer.hpp" #include #line 8 ".lib/math/integer.hpp" #include #line 13 ".lib/math/integer.hpp" // a^k >= xとなる最小のa^k uint64_t ceil_pow(uint64_t x, uint64_t a){ assert(a > 1); if(x == 0) return 1; static uint64_t INF = std::numeric_limits::max(); if(a == 2){ uint64_t ret = uint64_t(1) << (63 - __builtin_clzll(x)); if(ret == x) return x; assert(ret <= (INF >> 1)); return ret << 1; } if(a > 10){ uint64_t ret = 1; while(true){ if(ret > x / a) break; ret *= a; } if(ret == x) return ret; assert(ret <= INF / a); return ret * a; } static std::vector> pow_table(11); if(pow_table[a].empty()){ uint64_t tmp = 1; pow_table[a].push_back(1); while(true){ if(tmp > INF / a) break; pow_table[a].push_back(tmp *= a); } } auto itr = std::lower_bound(pow_table[a].begin(), pow_table[a].end(), x); assert(itr != pow_table[a].end()); return *itr; } // 2^k >= xとなる最小の2^k uint64_t ceil_2_pow(uint64_t x){ static uint64_t INF = std::numeric_limits::max(); uint64_t ret = uint64_t(1) << (63 - __builtin_clzll(x)); if(ret == x) return x; assert(ret <= (INF >> 1)); return ret << 1; } // a^k <= xを満たす最大のa^k uint64_t floor_pow(uint64_t x, uint64_t a){ assert(x > 0 && a > 1); if(a == 2) return uint64_t(1) << (63 - __builtin_clzll(x)); if(a > 10){ uint64_t ret = 1; while(true){ if(ret > x / a) return ret; ret *= a; } } static std::vector> pow_table(11); static uint64_t INF = std::numeric_limits::max(); if(pow_table[a].empty()){ uint64_t tmp = 1; pow_table[a].push_back(1); while(true){ if(tmp > INF / a) break; pow_table[a].push_back(tmp *= a); } } return *--std::upper_bound(pow_table[a].begin(), pow_table[a].end(), x); } // 10 = 10 * 1 + 0 // 10 = 9 * 1 + 1 // 10 = 8 * 1 + 2 // 10 = 7 * 1 + 3 // 10 = 6 * 1 + 4 // 10 = 5 * 2 = 0 // 10 = 4 * 2 + 2 // 10 = 3 * 3 = 1 // 10 = 2 * 5 + 0 // 10 = 1 * 10 + 10 // 商としてありえる数は高々 2 * √x 通り // また, [dmin, dmax]が存在して, この区間の数で割った商は全て等しく, 区間に含まれない任意の数で割った商とは異なる // 可能な組{商, dmin, dmax}を列挙 // 商が[√x, x]の時はdmin = dmax std::vector> divrange(long long x){ if(x == 1) return {{1, 1, 1}}; std::vector> l{{x, 1, 1}}, r{{1, x, x}}; int ls = 0, rs = 0; long long i = 2; for(; i * i <= x; i++){ long long d = x / i; if(l[ls][0] != d) l.push_back({d, i, i}), ls++; else l[ls][1] = i; if(i * i == x) continue; if(r[rs][0] != i) r[rs][1] = d + 1, r.push_back({i, d, d}), rs++; else r[rs][2] = x; } if(l[ls][0] == r[rs][0]) l[ls][2] = r[rs][2], r.pop_back(); std::reverse(r.begin(), r.end()); r[0][1] = i; l.insert(l.end(), r.begin(), r.end()); return l; } // a ^ k <= xを満たす最大のa uint64_t kth_root_stable(uint64_t x, uint64_t k){ if(!x) return 0; if(k == 1 || x == 1) return x; if(k >= 64) return 1; uint64_t l = 1, r = x; const static uint64_t threshold = std::numeric_limits::max(); while(r - l > 1){ bool f = false; uint64_t mid = l + ((r - l) >> 1), z = 1; uint64_t lim = threshold / mid; for(int i = 0; i < k; i++){ if(z > lim){ f = true; break; } z *= mid; } if(f | (z > x)) r = mid; else l = mid; } return l; } // a^k <= x を満たす最大のa uint64_t kth_root_fast(uint64_t x, uint64_t k){ if(x <= 1) return (!k ? 1 : x); if(k <= 2) return (k <=1 ? !k ? 1 : x : sqrtl(x)); if(k >= 64||!(x >> k)) return 1; const static int sz[8] = {40, 31, 27, 24, 22, 21, 20, 19}; static std::vector> table; if(table.empty()){ table.resize(40); for(int i = 0; i < 40; i++){ for(uint64_t j = 0; j < 8 && sz[j] > i; j++){ table[i].push_back((!i ? 1 : table[i - 1][j]) *(j + 3)); } } } if(k >= 19){ int ans = 10; for(;ans > 2; ans--){ if(sz[ans - 3] x) continue; return ans; } return 2; } uint64_t r = (k != 3 ? pow(x, (1.0 + 1e-6) / k) : cbrt(x) + 1); const static uint64_t threshold = std::numeric_limits::max(); while(true){ uint64_t lim = threshold / r, z = 1; for(int i = 0; i < k; i++, z *= r) if(z > lim) goto upper; if(z > x) upper: r--; else break; } return r; } namespace prime_sieve{ std::vector primes, min_factor;// 素数, 各数を割り切る最小の素数 // O(MAX_N loglog MAX_N) void init(int MAX_N){ min_factor.resize(MAX_N + 1, -1); for(int i = 2; i <= MAX_N; i++){ if(min_factor[i] == -1){ primes.push_back(i); min_factor[i] = i; } for(int p : primes){ if((long long)p * i > MAX_N || p > min_factor[i]) break; min_factor[p * i] = p; } } } bool is_prime(int n){ assert(n < min_factor.size()); return n == min_factor[n]; } // {{素因数, 数}}, O(log n) std::vector> factorize(int n){ assert(n < min_factor.size()); std::vector> ret; while(n > 1){ int cnt = 0, f = min_factor[n]; while(n % f == 0){ n /= f; cnt++; } ret.push_back({f, cnt}); } return ret; } // 約数列挙, O(√n) std::vector divisor(int n){ auto p = factorize(n); std::vector> x; for(int i = 0; i < p.size(); i++){ x.push_back(std::vector{1}); for(int j = 0; j < p[i].second; j++) x[i].push_back(x[i][j] * p[i].first); } int l = 0, r = 1; std::vector ret{1}; for(int i = 0; i < x.size(); i++){ for(auto e : x[i]){ for(int j = l; j < r; j++){ ret.push_back(ret[j] * e); } } l = r; r = ret.size(); } return std::vector(ret.begin() + l, ret.end()); } }; std::vector v{1,2,4,6,12,24,36,48,60,120,180,240,360,720,840,1260,1680,2520,5040,7560,10080,15120,20160,25200,27720,45360,50400,55440,83160,110880,166320,221760,277200,332640,498960,554400,665280,720720,1081080,1441440,2162160,2882880,3603600,4324320,6486480,7207200,8648640,10810800,14414400,17297280,21621600,32432400,36756720,43243200,61261200,73513440,110270160,122522400,147026880,183783600,245044800,294053760,367567200,551350800,698377680,735134400,1102701600,1396755360,2095133040,2205403200,2327925600,2793510720,3491888400,4655851200,5587021440,6983776800,10475665200,13967553600,20951330400,27935107200,41902660800,48886437600,64250746560,73329656400,80313433200,97772875200,128501493120,146659312800,160626866400,240940299600,293318625600,321253732800,481880599200,642507465600,963761198400,1124388064800,1606268664000,1686582097200,1927522396800,2248776129600,3212537328000,3373164194400,4497552259200,6746328388800,8995104518400,9316358251200,13492656777600,18632716502400,26985313555200,27949074753600,32607253879200,46581791256000,48910880818800,55898149507200,65214507758400,93163582512000,97821761637600,130429015516800,195643523275200,260858031033600,288807105787200,391287046550400,577614211574400,782574093100800,866421317361600,1010824870255200,1444035528936000,1516237305382800,1732842634723200,2021649740510400,2888071057872000,3032474610765600,4043299481020800,6064949221531200,8086598962041600,10108248702552000,1212898443062400,18194847664593600,20216497405104000,24259796886124800,30324746107656000,36389695329187200,48519593772249600,60649492215312000,72779390658374400,74801040398884800,106858629141264000,112201560598327200,149602080797769600,224403121196654400,299204161595539200,374005201994424000,448806242393308800,673209363589963200,748010403988848000,897612484786617600}; std::vector ans{1,2,3,4,6,8,9,10,12,16,18,20,24,30,32,36,40,48,60,64,72,80,84,90,96,100,108,120,128,144,160,168,180,192,200,216,224,240,256,288,320,336,360,384,400,432,448,480,504,512,576,600,640,672,720,768,800,864,896,960,1008,1024,1152,1200,1280,1344,1440,1536,1600,1680,1728,1792,1920,2016,2048,2304,2400,2688,2880,3072,3360,3456,3584,3600,3840,4032,4096,4320,4608,4800,5040,5376,5760,6144,6720,6912,7168,7200,7680,8064,8192,8640,9216,10080,10368,10752,11520,12288,12960,13440,13824,14336,14400,15360,16128,16384,17280,18432,20160,20736,21504,23040,24576,25920,26880,27648,28672,28800,30720,32256,32768,34560,36864,40320,41472,43008,46080,48384,49152,51840,53760,55296,57600,61440,62208,64512,65536,69120,73728,80640,82944,86016,92160,96768,98304,103680 }; // 高度合成数(N以下の数の中で最も多い約数を持つ数) // {N以下の高度合成数, その約数} std::pair highly_composite_number(long long N){ assert(0 < N && N <= 1000000000000000000); int idx = upper_bound(v.begin(), v.end(), N) - v.begin() - 1; assert(idx != 0); return {v[idx], ans[idx]}; } long long llpow(long long a, long long b){ long long ret = 1, mul = a; while(b){ if(b & 1) ret *= mul; mul *= mul; b >>= 1; } return ret; } long long gcd(long long _a, long long _b) { uint64_t a = abs(_a), b = abs(_b); if(a == 0) return b; if(b == 0) return a; int shift = __builtin_ctzll(a | b); a >>= __builtin_ctzll(a); do{ b >>= __builtin_ctzll(b); if(a > b) std::swap(a, b); b -= a; }while(b); return (a << shift); } // 64bitに収まらない可能性あり long long lcm(long long a, long long b){ return __int128_t(a * b) / gcd(a, b); } std::tuple extgcd(long long a, long long b){ long long x, y; for(long long u = y = 1, v = x = 0; a;){ long long q = b / a; std::swap(x -= q * u, u); std::swap(y -= q * v, v); std::swap(b -= q * a, a); } return {x, y, b};//return {x, y, gcd(a, b)} s.t. ax + by = gcd(a, b) } // ak + b のl <= k < rにおける和 template T arithmetic_sum(T a, T b, T l, T r){ return a * (r * (r - 1) - l * (l - 1)) / 2 + b * (r - l); } // ! struct barrett_reduction{ unsigned int mod; unsigned long long m; barrett_reduction(unsigned int _mod) : mod(_mod){ m = ((__uint128_t)1 << 64) / mod; } unsigned int reduce(unsigned int a){ unsigned long long q = ((__uint128_t)a * m) >> 64; a -= q * mod; // 0 <= a < 2 * mod // return a; return a >= mod ? a - mod : a; } unsigned int mul(unsigned int a, unsigned int b){ return reduce((unsigned long long)a * b); } // {gcd(a, mod), x}, such that a * x ≡ gcd(a, mod) std::pair inv(unsigned int a){ if(a >= mod) a = reduce(a); if(a == 0) return {mod, 0}; unsigned int s = mod, t = a; long long m0 = 0, m1 = 1; while(t){ int u = s / t; s -= t * u; m0 -= m1 * u; std::swap(m0, m1); std::swap(s, t); } if(m0 < 0) m0 += mod / s; return {s, m0}; } }; // 64bit mod対応 struct montgomery_reduction_64bit{ private: // [0, 2 * MOD) inline uint64_t reduce_unsafe(__uint128_t x) const{ x = (x + ((uint64_t)x * (uint64_t)NEG_INV) * MOD) >> 64; return x; } void _set_mod(uint64_t mod){ assert(mod < (1ULL << 63)); MOD = mod; NEG_INV = 0; __uint128_t s = 1, t = 0; for(int i = 0; i < 64; i++){ if (~t & 1) { t += MOD; NEG_INV += s; } t >>= 1; s <<= 1; } R2 = ((__uint128_t)1 << 64) % MOD; R2 = R2 * R2 % MOD; } __uint128_t MOD, NEG_INV, R2; public: montgomery_reduction_64bit(){} montgomery_reduction_64bit(uint64_t mod){_set_mod(mod);} void set_mod(uint64_t mod){ _set_mod(mod); } uint64_t mod()const{ return MOD; } inline uint64_t generate(uint64_t x)const{ return reduce((__uint128_t)x * R2); } inline uint64_t reduce(__uint128_t x)const{ x = (x + ((uint64_t)x * (uint64_t)NEG_INV) * MOD) >> 64; return x < MOD ? x : x - MOD; } // [0, 2MOD) inline uint64_t mul(uint64_t mx, uint64_t my)const{ return reduce_unsafe((__uint128_t)mx * my); } inline uint64_t mul_safe(uint64_t mx, uint64_t my)const{ return reduce((__uint128_t)mx * my); } // [0, 2MOD) inline uint64_t add(uint64_t mx, uint64_t my)const{ return (mx >= MOD ? mx - MOD : mx) + (my >= MOD ? my - MOD : my); } inline uint64_t add_safe(uint64_t mx, uint64_t my)const{ uint64_t res = (mx >= MOD ? mx - MOD : mx) + (my >= MOD ? my - MOD : my); return res >= MOD ? res - MOD : res; } // [0, 2MOD) inline uint64_t sub(uint64_t mx, uint64_t my)const{ if(my >= MOD) my -= MOD; return mx >= my ? mx - my : mx + MOD - my; } inline uint64_t sub_safe(uint64_t mx, uint64_t my)const{ if(my >= MOD) my -= MOD; uint64_t res = mx >= my ? mx - my : mx + MOD - my; return res >= MOD ? res - MOD : res; } }; // maと戻り値はモンゴメリ表現 unsigned long long _mod_pow_mr(unsigned long long ma, unsigned long long b, const montgomery_reduction_64bit &mr){ unsigned long long res = mr.generate(1), mul = ma; while(b){ if(b & 1) res = mr.mul(res, mul); mul = mr.mul(mul, mul); b >>= 1; } if(res >= mr.mod()) return res - mr.mod(); return res; } unsigned long long mod_pow_mr(unsigned long long a, unsigned long long b, unsigned long long m){ montgomery_reduction_64bit mr(m); return mr.reduce(_mod_pow_mr(mr.generate(a), b, mr)); } static std::vector generate_random_a(int k){ static std::random_device seed_gen; static std::mt19937_64 engine(seed_gen()); std::vector res(k); for(int i = 0; i < k; i++) res[i] = engine(); return res; } bool _miller_rabin_mr(unsigned long long n, const montgomery_reduction_64bit &mr){ static constexpr int k = 10; static std::vector small_p{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47}; static std::vector A = generate_random_a(k); static std::vector B{2, 3, 5, 7};// https://t5k.org/prove/prove2_3.html if(n <= 1) return false; if(n <= 50){ for(int l = n < 20 ? 0 : 8, r = n < 20 ? 8 : 15; l < r; l++) if(small_p[l] == n) return true; return false; } if(!(n & 1)) return false; unsigned long long d = n - 1; unsigned long long one = mr.generate(1), mone = mr.generate(n - 1); d >>= __builtin_ctzll(d); for(unsigned long long a : (n < 118670087467 ? B : A)){ if(a % n == 0) continue; unsigned long long d2s = d; // d * 2^s, 0 <= s <= (n - 1)が2で割れる回数 unsigned long long y = _mod_pow_mr(mr.generate(a), d, mr); while(d2s != n - 1 && y != one && y != mone){ y = mr.mul_safe(y, y); d2s <<= 1; } if(y != mone && !(d2s & 1)) return false; } return true; } bool miller_rabin_mr(unsigned long long n){ montgomery_reduction_64bit mr(n); return _miller_rabin_mr(n, mr); } // https://en.wikipedia.org/wiki/Binary_GCD_algorithm unsigned long long binary_gcd(unsigned long long a, unsigned long long b){ if(!a || !b) return !a ? b : a; int shift = __builtin_ctzll(a | b); // [1] gcd(2a', 2b') = 2 * gcd(a', b') a >>= __builtin_ctzll(a); do{ // if b is odd // gcd(2a', b) = gcd(a', b), if a = 2a'(a is even) // gcd(a, b) = gcd(|a - b|, min(a, b)), if a is odd b >>= __builtin_ctzll(b); // make b odd if(a > b) std::swap(a, b); b -= a; }while(b); // gcd(a, 0) = a return a << shift; // [1] } unsigned long long generate_random_prime(unsigned long long min_n = 2, unsigned long long max_n = ~0){ static std::random_device seed_gen; static std::mt19937_64 engine(seed_gen()); unsigned long long len = max_n - min_n + 1; // about 40 tries per iteration. https://en.wikipedia.org/wiki/Prime_number_theorem // It brokes when π(max_n) == π(min_n - 1) while(true){ unsigned long long a = engine() % len + min_n; if(miller_rabin_mr(a)){ return a; } } } namespace rho_factorization{ unsigned long long rho(unsigned long long n){ static std::vector small_p{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47}; for(int sp : small_p) if(n % sp == 0) return sp; // n < 50 montgomery_reduction_64bit mr(n); if(_miller_rabin_mr(n, mr)) return n; auto try_factorize = [n, mr](unsigned long long c){ c = mr.generate(c); auto f = [mr, c](unsigned long long mx){ return mr.add(mr.mul(mx, mx), c); }; unsigned long long m = 1LL << ((64 - __builtin_clzll(n)) / 8); unsigned long long y = n, r = 1, q = 1; unsigned long long x, g, k, ys; do{ x = y; y = mr.generate(y); for(int i = 0; i < r; i++) y = f(y); y = mr.reduce(y); k = 0; while(k < r && g == 1){ q = mr.generate(q); y = mr.generate(y); ys = y; for(int i = 0; i < std::min(m, r - k); i++){ y = f(y); unsigned long long z = mr.reduce(y); q = mr.mul(q, mr.generate(x > z ? x - z : z - x)); } y = mr.reduce(y); g = binary_gcd(mr.reduce(q), n); k += m; } r <<= 1; }while(g == 1); if(g == n){ do{ ys = f(ys); unsigned long long z = mr.reduce(ys); g = binary_gcd(x > z ? x - z : z - x, n); }while(g == 1); } return g; // g == n when failure }; unsigned long long c = 1, res = n; do{ res = try_factorize(c); // c = generate_random_prime(2, n - 1); c = (c + 1) % n; }while(res == n); return res; } std::vector factorize(unsigned long long n){ if(n <= 1) return {}; unsigned long long x = rho(n); if(x == n) return {x}; auto l = factorize(x); auto r = factorize(n / x); l.insert(l.end(), r.begin(), r.end()); return l; } // 素因数の集合(重複なし, ソート済)を返す std::vector prime_factor(unsigned long long n){ auto p = factorize(n); sort(p.begin(), p.end()); p.erase(std::unique(p.begin(), p.end()), p.end()); return p; } // 10^18以下の高度合成数 897612484786617600の約数が103680個なので全列挙して良さそう std::vector divisor(unsigned long long n){ auto p = factorize(n); sort(p.begin(), p.end()); std::vector> x; for(int i = 0; i < p.size(); i++){ if(!i || p[i] != p[i - 1]) x.push_back({p[i], 1}); else x.back().second++; } int sz = 1; for(auto [p_cur, cnt] : x) sz *= cnt + 1; std::vector res(sz); res[0] = 1; int r_prev = 1, r = 1; for(auto [p_cur, cnt] : x){ unsigned long long ppow = 1; for(int c = 0; c < cnt; c++){ ppow *= p_cur; for(int i = 0; i < r_prev; i++){ res[r++] = res[i] * ppow; } } r_prev = r; } return res; } } //---tips--- //φ(n) := [1, n]の中でnと互いに素な正整数の数 //φ(n) = n * π{p:prime factor of n}(1 - 1/p) //aとnが互いに素な場合, a^φ(n) ≡ 1 (mod n) unsigned long long totient(unsigned long long n){ unsigned long long res = n; auto prims = rho_factorization::prime_factor(n); for(auto p : prims) res -= res / p; return res; } // n以下のtotient_sum std::vector totient_sum_table(unsigned long long n){ std::vector res(n + 1); std::iota(res.begin() + 1, res.end(), 0); res[1] = 1; for(int i = 2; i <= n; i++){ // prime if(res[i] == i - 1){ for(int j = 2 * i; j <= n; j += i){ res[j] -= res[j] / i; } } } for(int i = 2; i <= n; i++) res[i] += res[i - 1]; return res; } // [1] // totient_sum(n) = (i, j){1 <= i <= n, 1 <= j <= i}となるペアの数 - そのうちgcd(i, j) != 1であるものの数 // = n * (n + 1) / 2 - ∑{2 <= g <= n} gcd(i, j) == gのペア // = n * (n + 1) / 2 - ∑{2 <= g <= n} totient_sum(n / g) // この演算をメモ化再帰で行った時の n / g (切り捨て)の種類数について考える // [2] // floor(x / ab) = floor(floor(x / a) / b) // 証明: // x = qab + r (0 <= r < ab)とすると // 左辺 = floor(q + r / ab) = q // 右辺 = floor((qb + r') / b) (0 <= r' < b) = q // [3] // [2]よりnを正整数で0回以上除算した時に得られる商はO(√n)種類 __uint128_t totient_sum(unsigned long long n){ static std::vector low; if(low.empty()) low = totient_sum_table(std::min(n, (unsigned long long)4000000)); std::unordered_map mp; unsigned long long memo_max = 0; auto tsum = [&](auto &&tsum, unsigned long long m)->__uint128_t{ if(m < low.size()) return low[m]; if(m <= memo_max) return mp.at(m); __uint128_t res = (__uint128_t)m * (m + 1) / 2; auto d = divrange(m); std::reverse(d.begin(), d.end()); for(auto [q, a, b] : d){ if(q == m) continue; res -= (b - a + 1) * tsum(tsum, q); } mp.emplace(m, res); memo_max = m; return res; }; return tsum(tsum, n); } // p: 素数 unsigned long long is_primitive_root(unsigned long long p){ static std::random_device seed_gen; static std::mt19937_64 engine(seed_gen()); //assert(miller_rabin_mr(p)); auto primes = rho_factorization::prime_factor(p - 1); while(true){ bool f = true; unsigned long long a = engine() % (p - 1) + 1; for(unsigned long long pk : primes){ // a ^ (p - 1) / pk ≡ 1 (mod p) -> no if(mod_pow_mr(a, (p - 1) / pk, p) == 1){ f = false; break; } } if(f) return a; } } // ∑{i, 0 <= i < n} floor((ai + b) / m) // n, m, a, b <= 10^9 long long floor_sum(long long n, long long m, long long a, long long b){ long long ans = 0; if(a >= m){ ans += (n - 1) * n * (a / m) / 2; a %= m; } if(b >= m){ ans += n * (b / m); b %= m; } // 変形: // https://rsk0315.hatenablog.com/entry/2020/12/13/231307 // https://atcoder.jp/contests/practice2/editorial/579 long long y_max = a * n + b; if(y_max >= m) ans += floor_sum(y_max / m, a, m, y_max % m); return ans; } // min((ax + b) % mod) x in[l, r) , 0 <= l < r long long min_linear_mod(long long a, long long b, long long l, long long r, long long mod){ assert(l <= 0 && l < r); a %= mod; b = (b + a * l) % mod; // [0, r - l)に変換 r -= l; if(a == 0) return b; long long upper = floor_sum(r, mod, a, b); long long L = 0, R = b + 1; // 定数倍高速化 long long tmp_val = b, tmp_idx = 0; for(int i = 0; i < 10; i++){ long long sa = mod - tmp_val; tmp_idx += (sa + a - 1) / a; tmp_val = (tmp_idx * a + b) % mod; if(tmp_idx >= r) break; R = std::min(R, tmp_val + 1); } while(R - L > 1){ long long mid = (L + R) / 2; if(upper == floor_sum(r, mod, a, b - mid)) L = mid; else R = mid; } return L; } /// constexpr long long safe_mod(long long x, long long m) { x %= m; if (x < 0) x += m; return x; } // @param b `1 <= b` // @return pair(g, x) s.t. g = gcd(a, b), xa = g (mod b), 0 <= x < b/g constexpr std::pair inv_gcd(long long a, long long b) { a = safe_mod(a, b); if (a == 0) return {b, 0}; // Contracts: // [1] s - m0 * a = 0 (mod b) // [2] t - m1 * a = 0 (mod b) // [3] s * |m1| + t * |m0| <= b long long s = b, t = a; long long m0 = 0, m1 = 1; while (t) { long long u = s / t; s -= t * u; m0 -= m1 * u; // |m1 * u| <= |m1| * s <= b // [3]: // (s - t * u) * |m1| + t * |m0 - m1 * u| // <= s * |m1| - t * u * |m1| + t * (|m0| + |m1| * u) // = s * |m1| + t * |m0| <= b auto tmp = s; s = t; t = tmp; tmp = m0; m0 = m1; m1 = tmp; } // by [3]: |m0| <= b/g // by g != b: |m0| < b/g if (m0 < 0) m0 += b / s; return {s, m0}; } ///// // m_iで割るとr_i余る情報から(全ての条件を満たすlcm(m_i)以下の数)%modを返す // mが互いに素でないと壊れる // O(n^2) #line 754 ".lib/math/integer.hpp" long long garner_coprime(std::vector r, std::vector m, long long mod){ int n = r.size(); assert(n == m.size()); // ans = r_1 + x_1 * m_1 + x_2 * m_1 * m_2 ... // x_i = (r_i - (i - 1項までの和)) / (m_1 * m_2 ... m_i-1) mod m_i // M_i := (m_1 * m_2 ... m_i - 1) mod m_i std::vector accum(n + 1, 0), M(n + 1, 1); m.push_back(mod); for(int i = 0; i < n; i++){ auto [g, inv] = inv_gcd(M[i], m[i]); assert(g == 1); long long x = (((__int128_t)r[i] - accum[i]) * inv) % m[i]; if(x < 0) x += m[i]; for(int j = i + 1; j <= n; j++){ accum[j] = (accum[j] + (__int128_t)M[j] * x) % m[j]; M[j] = ((__int128_t)M[j] * m[i]) % m[j]; } } return accum[n]; } // m_iで割るとr_i余る情報から(全ての条件を満たすlcm(m_i)以下の数)%modを返す // 解が無い場合は-1 long long garner(std::vector &r, std::vector &m, long long mod){ int n = r.size(); assert(n == m.size()); // mを互いに素にする for(int i = 1; i < n; i++){ for(int j = 0; j < i; j++){ long long x = m[i], y = m[j], g = gcd(x, y); while(g != 1){ // g2 := gcd(g, x / g) != 1なら g2に含まれる素因数はy / gには含まれない. g3も同じ. // g2, g3はgの約数 long long g2 = gcd(g, x / g); if(g2 == 1){ x /= g; break; } long long g3 = gcd(g, y / g); if(g3 == 1){ y /= g; break; } x /= g3, y /= g2; g /= g2 * g3; } m[i] = x, r[i] %= x; m[j] = y, r[j] %= y; } } return garner_coprime(r, m, mod); } #line 3 "a.cpp" int main(){ io_init(); int n; std::cin >> n; vec r(n), m(n); range(i, 0, n) std::cin >> r[i] >> m[i]; std::cout << garner(r, m, 1000000007) << '\n'; }