This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub tko919/library
#include "Math/binomquery.hpp"
#pragma once #include "Math/fastdiv.hpp" #include "Math/primitive.hpp" struct BinomialQuery{ struct X{ // for m=p^q int p,q,m,delta; vector<int> fact,ifac; FastDiv im,ip; X(){} X(int _p,int _q):p(_p),q(_q),ip(p){ m=1; while(_q--)m*=p; im=FastDiv(m); fact.resize(m); ifac.resize(m); fact[0]=fact[1]=ifac[0]=ifac[1]=1; rep(i,2,m){ if(i%ip==0)fact[i]=fact[i-1]; else fact[i]=(ll(fact[i-1])*i)%im; } ifac[m-1]=minv(fact[m-1],m); for(int i=m-2;i>1;i--){ if(i%ip==p-1)ifac[i]=ifac[i+1]; else ifac[i]=(ll(ifac[i+1])*(i+1))%im; } delta=(p==2 and q>=3?0:1); } int nCr(ll n,ll r){ if(n<0 or r<0 or n<r)return 0; ll s=n-r,ret=1; int e0=0,em=0; for(int i=1;n;i++){ ret=(ret*fact[n%m])%im; ret=(ret*ifac[r%m])%im; ret=(ret*ifac[s%m])%im; n=n/ip,r=r/ip,s=s/ip; int add=n-r-s; e0+=add; if(e0>=q)return 0; if(i>=q)em^=add; } if(delta and em&1)ret=m-ret; ret=(ret*mpow(p,e0,m))%im; return ret; } }; int mod; vector<ll> ms; vector<X> Xs; BinomialQuery(int M):mod(M){ for(int p=2;p*p<=M;p++)if(M%p==0){ int e=0,base=1; while(M%p==0)M/=p,e++,base*=p; ms.push_back(base); Xs.push_back(X(p,e)); } if(M!=1){ ms.push_back(M); Xs.push_back(X(M,1)); } } int nCr(ll n,ll r){ if(mod==1)return 0; vector<ll> ret; for(auto& buf:Xs)ret.push_back(buf.nCr(n,r)); return crt(ret,ms).first; } }; /** * @brief Binomial Coefficient for query */
#line 2 "Math/fastdiv.hpp" struct FastDiv{ using u64=uint64_t; using u128=__uint128_t; constexpr FastDiv():m(),s(),x(){} constexpr FastDiv(int _m) :m(_m),s(__lg(m-1)),x(((u128(1)<<(s+64))+m-1)/m){} constexpr int get(){return m;} constexpr friend u64 operator/(u64 n,const FastDiv& d){ return (u128(n)*d.x>>d.s)>>64; } constexpr friend int operator%(u64 n,const FastDiv& d){ return n-n/d*d.m; } constexpr pair<u64,int> divmod(u64 n)const{ u64 q=n/(*this); return {q,n-q*m}; } int m,s; u64 x; }; /** * @brief Fast Division */ #line 2 "Math/miller.hpp" struct m64 { using i64 = int64_t; using u64 = uint64_t; using u128 = __uint128_t; static u64 mod; static u64 r; static u64 n2; static u64 get_r() { u64 ret = mod; rep(_,0,5) ret *= 2 - mod * ret; return ret; } static void set_mod(u64 m) { assert(m < (1LL << 62)); assert((m & 1) == 1); mod = m; n2 = -u128(m) % m; r = get_r(); assert(r * mod == 1); } static u64 get_mod() { return mod; } u64 a; m64() : a(0) {} m64(const int64_t &b) : a(reduce((u128(b) + mod) * n2)){}; static u64 reduce(const u128 &b) { return (b + u128(u64(b) * u64(-r)) * mod) >> 64; } u64 get() const { u64 ret = reduce(a); return ret >= mod ? ret - mod : ret; } m64 &operator*=(const m64 &b) { a = reduce(u128(a) * b.a); return *this; } m64 operator*(const m64 &b) const { return m64(*this) *= b; } bool operator==(const m64 &b) const { return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a); } bool operator!=(const m64 &b) const { return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a); } m64 pow(u128 n) const { m64 ret(1), mul(*this); while (n > 0) { if (n & 1) ret *= mul; mul *= mul; n >>= 1; } return ret; } }; typename m64::u64 m64::mod, m64::r, m64::n2; bool Miller(ll n){ if(n<2 or (n&1)==0)return (n==2); m64::set_mod(n); ll d=n-1; while((d&1)==0)d>>=1; vector<ll> seeds; if(n<(1<<30))seeds={2, 7, 61}; else seeds={2, 325, 9375, 28178, 450775, 9780504}; for(auto& x:seeds){ if(n<=x)break; ll t=d; m64 y=m64(x).pow(t); while(t!=n-1 and y!=1 and y!=n-1){ y*=y; t<<=1; } if(y!=n-1 and (t&1)==0)return 0; } return 1; } /** * @brief Miller-Rabin */ #line 2 "Utility/random.hpp" namespace Random { mt19937_64 randgen(chrono::steady_clock::now().time_since_epoch().count()); using u64 = unsigned long long; u64 get() { return randgen(); } template <typename T> T get(T L) { // [0,L] return get() % (L + 1); } template <typename T> T get(T L, T R) { // [L,R] return get(R - L) + L; } double uniform() { return double(get(1000000000)) / 1000000000; } string str(int n) { string ret; rep(i, 0, n) ret += get('a', 'z'); return ret; } template <typename Iter> void shuffle(Iter first, Iter last) { if (first == last) return; int len = 1; for (auto it = first + 1; it != last; it++) { len++; int j = get(0, len - 1); if (j != len - 1) iter_swap(it, first + j); } } template <typename T> vector<T> select(int n, T L, T R) { // [L,R] if (n * 2 >= R - L + 1) { vector<T> ret(R - L + 1); iota(ALL(ret), L); shuffle(ALL(ret)); ret.resize(n); return ret; } else { unordered_set<T> used; vector<T> ret; while (SZ(used) < n) { T x = get(L, R); if (!used.count(x)) { used.insert(x); ret.push_back(x); } } return ret; } } void relabel(int n, vector<pair<int, int>> &es) { shuffle(ALL(es)); vector<int> ord(n); iota(ALL(ord), 0); shuffle(ALL(ord)); for (auto &[u, v] : es) u = ord[u], v = ord[v]; } template <bool directed, bool simple> vector<pair<int, int>> genGraph(int n) { vector<pair<int, int>> cand, es; rep(u, 0, n) rep(v, 0, n) { if (simple and u == v) continue; if (!directed and u > v) continue; cand.push_back({u, v}); } int m = get(SZ(cand)); vector<int> ord; if (simple) ord = select(m, 0, SZ(cand) - 1); else { rep(_, 0, m) ord.push_back(get(SZ(cand) - 1)); } for (auto &i : ord) es.push_back(cand[i]); relabel(n, es); return es; } vector<pair<int, int>> genTree(int n) { vector<pair<int, int>> es; rep(i, 1, n) es.push_back({get(i - 1), i}); relabel(n, es); return es; } }; // namespace Random /** * @brief Random */ #line 4 "Math/pollard.hpp" vector<ll> Pollard(ll n) { if (n <= 1) return {}; if (Miller(n)) return {n}; if ((n & 1) == 0) { vector<ll> v = Pollard(n >> 1); v.push_back(2); return v; } for (ll x = 2, y = 2, d;;) { ll c = Random::get(2LL, n - 1); do { x = (__int128_t(x) * x + c) % n; y = (__int128_t(y) * y + c) % n; y = (__int128_t(y) * y + c) % n; d = __gcd(x - y + n, n); } while (d == 1); if (d < n) { vector<ll> lb = Pollard(d), rb = Pollard(n / d); lb.insert(lb.end(), ALL(rb)); return lb; } } } /** * @brief Pollard-Rho */ #line 4 "Math/primitive.hpp" ll mpow(ll a, ll t, ll m) { ll res = 1; FastDiv im(m); while (t) { if (t & 1) res = __int128_t(res) * a % im; a = __int128_t(a) * a % im; t >>= 1; } return res; } ll minv(ll a, ll m) { ll b = m, u = 1, v = 0; while (b) { ll t = a / b; a -= t * b; swap(a, b); u -= t * v; swap(u, v); } u = (u % m + m) % m; return u; } ll getPrimitiveRoot(ll p) { vector<ll> ps = Pollard(p - 1); sort(ALL(ps)); rep(x, 1, inf) { for (auto &q : ps) { if (mpow(x, (p - 1) / q, p) == 1) goto fail; } return x; fail:; } assert(0); } ll extgcd(ll a, ll b, ll &p, ll &q) { if (b == 0) { p = 1; q = 0; return a; } ll d = extgcd(b, a % b, q, p); q -= a / b * p; return d; } pair<ll, ll> crt(const vector<ll> &vs, const vector<ll> &ms) { ll V = vs[0], M = ms[0]; rep(i, 1, vs.size()) { ll p, q, v = vs[i], m = ms[i]; if (M < m) swap(M, m), swap(V, v); ll d = extgcd(M, m, p, q); if ((v - V) % d != 0) return {0, -1}; ll md = m / d, tmp = (v - V) / d % md * p % md; V += M * tmp; M *= md; } V = (V % M + M) % M; return {V, M}; } ll ModLog(ll a, ll b, ll p) { ll g = 1; for (ll t = p; t; t >>= 1) g = g * a % p; g = __gcd(g, p); ll t = 1, c = 0; for (; t % g; c++) { if (t == b) return c; t = t * a % p; } if (b % g) return -1; t /= g, b /= g; ll n = p / g, h = 0, gs = 1; for (; h * h < n; h++) gs = gs * a % n; unordered_map<ll, ll> bs; for (ll s = 0, e = b; s < h; bs[e] = ++s) e = e * a % n; for (ll s = 0, e = t; s < n;) { e = e * gs % n, s += h; if (bs.count(e)) { return c + s - bs[e]; } } return -1; } ll mod_root(ll k, ll a, ll m) { if (a == 0) return k ? 0 : -1; if (m == 2) return a & 1; k %= m - 1; ll g = gcd(k, m - 1); if (mpow(a, (m - 1) / g, m) != 1) return -1; a = mpow(a, minv(k / g, (m - 1) / g), m); FastDiv im(m); auto _subroot = [&](ll p, int e, ll a) -> ll { // x^(p^e)==a(mod m) ll q = m - 1; int s = 0; while (q % p == 0) { q /= p; s++; } int d = s - e; ll pe = mpow(p, e, m), res = mpow(a, ((pe - 1) * minv(q, pe) % pe * q + 1) / pe, m), c = 1; while (mpow(c, (m - 1) / p, m) == 1) c++; c = mpow(c, q, m); map<ll, ll> mp; ll v = 1, block = sqrt(d * p) + 1, bs = mpow(c, mpow(p, s - 1, m - 1) * block % (m - 1), m); rep(i, 0, block + 1) mp[v] = i, v = v * bs % im; ll gs = minv(mpow(c, mpow(p, s - 1, m - 1), m), m); rep(i, 0, d) { ll err = a * minv(mpow(res, pe, m), m) % im; ll pos = mpow(err, mpow(p, d - 1 - i, m - 1), m); rep(j, 0, block + 1) { if (mp.count(pos)) { res = res * mpow(c, (block * mp[pos] + j) * mpow(p, i, m - 1) % (m - 1), m) % im; break; } pos = pos * gs % im; } } return res; }; for (ll d = 2; d * d <= g; d++) if (g % d == 0) { int sz = 0; while (g % d == 0) { g /= d; sz++; } a = _subroot(d, sz, a); } if (g > 1) a = _subroot(g, 1, a); return a; } ull floor_root(ull a, ull k) { if (a <= 1 or k == 1) return a; if (k >= 64) return 1; if (k == 2) return sqrtl(a); constexpr ull LIM = -1; if (a == LIM) a--; auto mul = [&](ull &x, const ull &y) { if (x <= LIM / y) x *= y; else x = LIM; }; auto pw = [&](ull x, ull t) -> ull { ull y = 1; while (t) { if (t & 1) mul(y, x); mul(x, x); t >>= 1; } return y; }; ull ret = (k == 3 ? cbrt(a) - 1 : pow(a, nextafter(1 / double(k), 0))); while (pw(ret + 1, k) <= a) ret++; return ret; } /** * @brief Primitive Function */ #line 4 "Math/binomquery.hpp" struct BinomialQuery{ struct X{ // for m=p^q int p,q,m,delta; vector<int> fact,ifac; FastDiv im,ip; X(){} X(int _p,int _q):p(_p),q(_q),ip(p){ m=1; while(_q--)m*=p; im=FastDiv(m); fact.resize(m); ifac.resize(m); fact[0]=fact[1]=ifac[0]=ifac[1]=1; rep(i,2,m){ if(i%ip==0)fact[i]=fact[i-1]; else fact[i]=(ll(fact[i-1])*i)%im; } ifac[m-1]=minv(fact[m-1],m); for(int i=m-2;i>1;i--){ if(i%ip==p-1)ifac[i]=ifac[i+1]; else ifac[i]=(ll(ifac[i+1])*(i+1))%im; } delta=(p==2 and q>=3?0:1); } int nCr(ll n,ll r){ if(n<0 or r<0 or n<r)return 0; ll s=n-r,ret=1; int e0=0,em=0; for(int i=1;n;i++){ ret=(ret*fact[n%m])%im; ret=(ret*ifac[r%m])%im; ret=(ret*ifac[s%m])%im; n=n/ip,r=r/ip,s=s/ip; int add=n-r-s; e0+=add; if(e0>=q)return 0; if(i>=q)em^=add; } if(delta and em&1)ret=m-ret; ret=(ret*mpow(p,e0,m))%im; return ret; } }; int mod; vector<ll> ms; vector<X> Xs; BinomialQuery(int M):mod(M){ for(int p=2;p*p<=M;p++)if(M%p==0){ int e=0,base=1; while(M%p==0)M/=p,e++,base*=p; ms.push_back(base); Xs.push_back(X(p,e)); } if(M!=1){ ms.push_back(M); Xs.push_back(X(M,1)); } } int nCr(ll n,ll r){ if(mod==1)return 0; vector<ll> ret; for(auto& buf:Xs)ret.push_back(buf.nCr(n,r)); return crt(ret,ms).first; } }; /** * @brief Binomial Coefficient for query */