「学习笔记」组合计数与中国剩余定理
「学习笔记」组合计数与中国剩余定理
点击查看目录
目录
- Combination
思路
Lucas 定理 \((6)\) 板子题.
Code
点击查看代码
namespace SOLVE { const ll P = 1e4 + 7, N = 1e4 + 10; ll T, x, y, fac[N], inv[N]; inline ll rnt () { ll x = 0, w = 1; char c = getchar(); while (!isdigit(c)) { if (c == '-') w = -1; c = getchar();} while (isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar(); return x * w; } inline ll FastPow (ll a, ll b) { ll ans = 1; while (b) { if (b & 1) ans = ans * a % P; a = a * a % P, b >>= 1; } return ans; } inline void Pre () { fac[0] = 1; _for (i, 1, P) fac[i] = fac[i - 1] * i % P; inv[P - 1] = FastPow(fac[P - 1], P - 2); for_ (i, P - 2, 0) inv[i] = inv[i + 1] * (i + 1) % P; return ; } inline ll C (ll n, ll m) { if (m > n) return 0; return fac[n] * inv[n - m] % P * inv[m] % P; } inline ll Lucas (ll n, ll m) { if (m == 0) return 1; return C(n % P, m % P) * Lucas(n / P, m / P) % P; } inline void In () { x = rnt(), y = rnt(); return ; } inline void Out () { printf("%lld\n", Lucas(x, y)); return ; } }[SDOI2016]排列计数
思路
我们钦定 \(m\) 个数为稳定的,方案数为 \(\dbinom{n}{m}\).
在剩下的 \(n-m\) 个位置里要保证每个数不稳定.
欸那不就是错排列 \((3)\) 吗?
那么方案数就是 \(D_{n-m}\).
总方案数就是 \(\dbinom{n}{m}D_{n-m}\),\(\Theta(n)\) 预处理一下错排列,阶乘与逆元可用 \(\Theta(1)\) 求出单次询问.
代码
点击查看代码
namespace SOLVE { const ll P = 1e9 + 7, N = 1e6 + 10, M = 1e6; ll T, n, m, d[N], fac[N], inv[N]; inline ll rnt () { ll x = 0, w = 1; char c = getchar(); while (!isdigit(c)) { if (c == '-') w = -1; c = getchar();} while (isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar(); return x * w; } inline ll FastPow(ll a, ll b) { ll ans = 1; while (b) { if(b & 1) ans = ans * a % P; a = a * a % P, b >>= 1; } return ans; } inline void Pre () { d[0] = 1, d[1] = 0, fac[0] = 1; _for (i, 2, M) d[i] = (i - 1) * ((d[i - 1] + d[i - 2]) % P) % P; _for (i, 1, M) fac[i] = fac[i - 1] * i % P; inv[M] = FastPow(fac[M], P - 2); for_ (i, M - 1, 0) inv[i] = inv[i + 1] * (i + 1) % P; return; } inline ll C(ll n, ll m) { return fac[n] * inv[n - m] % P * inv[m] % P; } inline void In () { n = rnt(), m = rnt(); return ; } inline void Out () { printf("%lld\n", C(n, m) * d[n-m] % P); return ; } }[ZJOI2010]排列计数
思路
观察一下可以发现满足性质的序列是一个小根堆.
那么设 \(s_i\) 表示以 \(i\) 为根的堆的大小,\(f_i\) 表示以 \(i\) 为根的堆的可行方案数(此时该子堆里的序号不是最终序号,而是在子堆内大小的排名,因为归并到父堆时要算分配给子堆不同序号的方案数).
那么转移方程就是:
\[s_{i}=s_{i*2}+s_{i*2+1}+1\\ f_{i}=\dbinom{s_{i}-1}{s_{i*2}}f_{i*2}f_{i*2+1} \](自己必须是最小的所以只能从 \(s_{i}-1\) 个序号选 \(s_{i*2}\) 分配给左儿子,剩下的全给右儿子)
代码
点击查看代码
namespace SOLVE { const ll N = 4e6 + 10; ll T, n, P, sz[N], f[N], fac[N]; inline ll rnt () { ll x = 0, w = 1; char c = getchar(); while (!isdigit(c)) { if (c == '-') w = -1; c = getchar();} while (isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar(); return x * w; } inline ll FastPow (ll a, ll b) { ll ans = 1; while (b) { if(b & 1) ans = ans * a % P; a = a * a % P, b >>= 1; } return ans; } inline void Pre () { fac[0] = 1; _for (i, 1, std::min(P, n)) fac[i] = fac[i - 1] * i % P; _for (i, 1, n * 2 + 1) f[i] = 1; return; } inline ll Inv (ll n) { return FastPow(fac[n], P - 2); } inline ll C (ll n, ll m) { if(!n || !m) return 1; return fac[n] * Inv(n - m) % P * Inv(m) % P; } inline ll Lucas (ll n, ll m) { if(!n || !m) return 1; return C(n % P, m % P) * Lucas(n / P, m / P) % P; } inline void In () { n = rnt(), P = rnt(); return ; } inline void Solve () { for_ (i, n, 1) { sz[i] = sz[i << 1] + sz[(i << 1) + 1] + 1; f[i] = f[i << 1] * f[(i << 1) + 1] % P * Lucas(sz[i] - 1, sz[i << 1]) % P; } return ; } inline void Out () { printf("%lld\n", f[1]); return ; } }BZOJ2839 集合计数
思路
谔项式反演.
设 \(f(i)\) 表示交集数量 \(\ge i\) 的方案数,\(g(i)\) 表示交集个数恰好为 \(i\) 个的方案数,那么答案为 \(g(k)\).
那么:
\[f(i)=\dbinom{n}{i}(2^{2^{n-i}}-1) \]即先确定 \(i\) 个必选,包含这 \(i\) 个的集合数为 \(2^{n-k}\) 个,每个集合都可以选或不选但不能一个不选,即 \(2^{2^{n-i}}-1\).
同时:
\[f(k)=\sum_{i=k}^{n}\dbinom{i}{k}g(i) \]等一下这式子是不是在哪里见过?
这不是 \((10)\) 吗?!
那么愉快的套一个谔项式反演:
\[\begin{aligned} g(k) &=\sum_{i=k}^{n}(-1)^{i-k}\dbinom{i}{k}f(i)\\ &=\sum_{i=k}^{n}(-1)^{i-k}\dbinom{i}{k}\dbinom{n}{i}(2^{2^{n-i}}-1)\\ \end{aligned} \]再加上一点预处理,就可以解决了.
代码
点击查看代码
namespace SOLVE { typedef long double ldb; typedef long long ll; typedef double db; const ll N = 1e6 + 10, P = 1e9 + 7; ll T, n, k, er[N], fac[N], inv[N], ans; inline ll rnt () { ll x = 0, w = 1; char c = getchar (); while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); } while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar (); return x * w; } inline ll FastPow (ll a, ll b) { ll ans = 1; while (b) { if (b & 1) ans = ans * a % P; a = a * a % P, b >>= 1; } return ans; } inline void Pre () { fac[0] = 1, er[0] = 2; _for (i, 1, n) { fac[i] = fac[i - 1] * i % P; er[i] = (er[i - 1] * er[i - 1]) % P; } inv[n] = FastPow (fac[n], P - 2); for_ (i, n - 1, 0) inv[i] = inv[i + 1] * (i + 1) % P; return; } inline ll C (ll n, ll m) { if (!m) return 1; return fac[n] * inv[n - m] % P * inv[m] % P; } inline void In () { n = rnt (), k = rnt (); return; } inline void Solve () { _for (i, k, n) { ll w = ((i - k) & 1) ? -1 : 1; ans = (ans + (er[n - i] - 1 + P) % P * C (n, i) % P * C (i, k) % P * w + P) % P; } return; } inline void Out () { printf ("%lld\n", ans); return; } }牡牛和牝牛
思路
我们枚举牝牛的数量 \(i\),那么一定会有 \(k\times(i-1)\) 只牡牛被固定住,此时剩下 \(w(i)=(n-i-k\times(i-1))\times[i>0]+n\times[i=0]\) 只牡牛可以随便选位置.
观察一下,看上去是只有 \(k+1\) 个地方可以插空,然而两只牝牛之间可以放多只牡牛,如何解决这个问题?
既然可以重复放,那我们就把重复放的位置 \(\text{new}\) 出来!
即把空的个数改为 \(k+1+(w(i)-1)=k+i\).
这样会不会导致选的全都是 \(\text{new}\) 出来的呢?不会,因为我们只 \(\text{new}\) 出来了 \(w(i)-1\) 个空,剩下的一只牛必然会被放在原有的位置.
那么答案就是:
\[\sum_{i=0}^{n}[w(i)\ge0]\dbinom{i+w(i)}{w(i)} \]代码
点击查看代码
namespace SOLVE { typedef long double ldb; typedef long long ll; typedef double db; const ll N = 1e5 + 10, P = 5e6 + 11; ll n, k, fac[N], inv[N], ans; inline ll rnt () { ll x = 0, w = 1; char c = getchar(); while (!isdigit(c)) { if (c == '-') w = -1; c = getchar();} while (isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar(); return x * w; } inline ll FastPow (ll a, ll b) { ll ans = 1; while (b) { if (b & 1) ans = ans * a % P; a = a * a % P, b >>= 1; } return ans; } inline void Pre () { fac[0] = 1; _for (i, 1, n) fac[i] = fac[i - 1] * i % P; inv[n] = FastPow (fac[n], P - 2); for_ (i, n - 1, 0) inv[i] = inv[i + 1] * (i + 1) % P; return; } inline ll C (ll n, ll m) { return fac[n] * inv[n - m] % P * inv[m] % P; } inline void In () { n = rnt (), k = rnt (); return; } inline void Solve () { _for (i, 0, n) { ll w = i ? (n - i - k * (i - 1)) : n; if (w < 0) break; ans = (ans + C (i + w, w)) % P; } return ; } inline void Out () { printf ("%lld\n", ans); return ; } }序列统计
思路
本题和上一题有些类似,每个数也是可以重复选的.
那么设 \(m=r-l+1\),长度为 \(i\) 的序列的方案数为 \(\dbinom{m+i-1}{i}\).
然后推式子:
\[\begin{aligned} \sum_{i=1}^{n}\dbinom{m+i-1}{i} &=\sum_{i=1}^{n}\dbinom{m+i-1}{m-1}+\dbinom{m}{m}-1\\ &=\sum_{i=2}^{n}\dbinom{m+i-1}{m-1}+\dbinom{m}{m-1}+\dbinom{m}{m}-1\\ &=\sum_{i=2}^{n}\dbinom{m+i-1}{m-1}+\dbinom{m+1}{m}-1\\ &=\sum_{i=3}^{n}\dbinom{m+i-1}{m-1}+\dbinom{m+2}{m}-1\\ &=\sum_{i=4}^{n}\dbinom{m+i-1}{m-1}+\dbinom{m+3}{m}-1\\ &=\cdots\\ &=\sum_{i=n}^{n}\dbinom{m+i-1}{m-1}+\dbinom{m+n-1}{m}-1\\ &=\dbinom{m+n}{m}-1\\ \end{aligned} \]但 \(n,m\) 过大,需要用到 \(\text{Lucas}\) 定理 \((6)\).
代码
点击查看代码
namespace SOLVE { typedef long double ldb; typedef long long ll; typedef double db; const ll N = 1e6 + 10, P = 1e6 + 3; ll T, n, m, l, r, fac[N], inv[N]; inline ll rnt () { ll x = 0, w = 1; char c = getchar (); while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); } while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar (); return x * w; } inline ll FastPow (ll a, ll b) { ll ans = 1; while (b) { if (b & 1) ans = ans * a % P; a = a * a % P, b >>= 1; } return ans; } inline void Pre () { fac[0] = 1; _for (i, 1, P - 1) fac[i] = fac[i - 1] * i % P; inv[P - 1] = FastPow (fac[P - 1], P - 2); for_ (i, P - 2, 0) inv[i] = inv[i + 1] * (i + 1) % P; return; } inline ll C (ll n, ll m) { if (n < m) return 0; if (!n || !m) return 1; return fac[n] * inv[n - m] % P * inv[m] % P; } inline ll Lucas (ll n, ll m) { if (n < m) return 0; if (!n || !m) return 1; return C (n % P, m % P) * Lucas (n / P, m / P) % P; } inline void In () { n = rnt (), l = rnt (), r = rnt (); m = r - l + 1; return; } inline void Out () { printf ("%lld\n", (Lucas (m + n, m) + P - 1) % P); return; } }[SDOI2009] 虔诚的墓主人
思路
代码
感觉以前写的代码太丑了.
于是又写了一份.
点击查看代码
namespace SOLVE { typedef long double ldb; typedef long long ll; typedef double db; const ll N = 1e5 + 10, P = 2147483648; ll n, m, w, k, C[N][20], ans; ll cx[N], cy[N], nx[N], ny; class TREE { public: ll x, y; inline bool operator < (TREE another) { return (y == another.y) ? (x < another.x) : (y < another.y); } } tr[N]; class TreeArray { public: ll b[N]; inline ll lowbit (ll x) { return x & -x; } inline void Update (ll x, ll y) { while (x <= n) { b[x] = (b[x] + y) % P; x += lowbit (x); } return; } inline ll Query (ll x) { ll ans = 0; while (x) { ans = (ans + b[x]) % P; x -= lowbit (x); } return ans; } } ta; namespace LISAN { ll ls1[N], ls2[N]; inline void lisan () { _for (i, 1, w) ls1[i] = tr[i].x; _for (i, 1, w) ls2[i] = tr[i].y; std::sort (ls1 + 1, ls1 + w + 1); std::sort (ls2 + 1, ls2 + w + 1); n = std::unique (ls1 + 1, ls1 + w + 1) - ls1; m = std::unique (ls2 + 1, ls2 + w + 1) - ls2; _for (i, 1, w) { tr[i].x = std::lower_bound (ls1 + 1, ls1 + n + 1, tr[i].x) - ls1; tr[i].y = std::lower_bound (ls2 + 1, ls2 + m + 1, tr[i].y) - ls2; } return; } } inline void Pre () { C[0][0] = 1; _for (i, 1, w) { C[i][0] = 1; _for (j, 1, k) C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % P; } return; } inline ll rnt () { ll x = 0, w = 1; char c = getchar (); while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); } while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar (); return x * w; } inline void In () { n = rnt (), m = rnt (), w = rnt (); _for (i, 1, w) tr[i].x = rnt (), tr[i].y = rnt (); k = rnt (); return; } inline void Solve () { LISAN::lisan (); std::sort (tr + 1, tr + w + 1); _for (i, 1, w) ++cx[tr[i].x], ++cy[tr[i].y]; _for (i, 1, w - 1) { ++ny, ++nx[tr[i].x]; ll last = (ta.Query (tr[i].x) - ta.Query (tr[i].x - 1) + P) % P; if (tr[i].y == tr[i + 1].y) { ll up_down = C[ny][k] * C[cy[tr[i].y] - ny][k] % P; ll left_right = (ta.Query (tr[i + 1].x - 1) - ta.Query (tr[i].x) + P) % P; ans = (ans + up_down * left_right % P) % P; } else ny = 0; ta.Update (tr[i].x, (C [nx[tr[i].x]][k] * C [cx[tr[i].x] - nx[tr[i].x]][k] % P - last + P) % P); } return; } inline void Out () { printf ("%lld\n", ans); return; } }[SDOI2010]地精部落
思路
\(f_{i,0}\) 表示长度为 \(i\) 且第一段山为山谷的序列数量.
\(f_{i,1}\) 表示长度为 \(i\) 且第一段山为山峰的序列数量.
\[f_{i,k}= \sum_{j=1}^{i}[j\bmod{2}=k]\dbinom{i-1}{j-1}f_{j-1,k}f_{i-j,0} \]代码
点击查看代码
namespace SOLVE { typedef long double ldb; typedef long long ll; typedef double db; const ll N = 4200 + 10; int n, P, f[N][2], C[N][N], ans; inline ll rnt () { ll x = 0, w = 1; char c = getchar (); while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); } while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar (); return x * w; } inline void In () { n = rnt (), P = rnt (); return; } inline void Solve () { f[0][0] = f[0][1] = f[1][0] = 1, C[0][0] = 1; _for (i, 1, n) { C[i][0] = 1; _for (j, 1, i) { f[i][j & 1] = ((ll)(f[i][j & 1]) + (ll)(f[j - 1][j & 1]) * (ll)(f[i - j][0]) % P * C[i - 1][j - 1] % P) % P; C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % P; } } return; } inline void Out () { printf ("%lld\n", (f[n][0] + f[n][1]) % P); return; } }[ZJOI2011]看电影
思路
这个式子还是蛮有意思的,但要用高精.
首先答案 \(=\frac{合法情况数}{总情况数}\),总情况数显然是 \(k^n\),难点在于如何算出合法情况数.
首先我们在 \(k\) 后面新增一个座位 \(k+1\) 然后拉链为环,让没有座位的人从头开始往后坐,这样一定所有人都会有座位,那么这样的环一共会有 \((k+1)^{n-1}\) 个(即圆排列).注意:这里暂时不考虑标号.
但是我们怎么判断做法是否合法呢?如果 \(k+1\) 这个位置最后有人,那么在没有环与新座位时就一定会有人站在这里,否则没人站着(即合法情况),也就是说 我们在环中找到空的位置当成 \(k+1\) 即可,空的位置一共有 \(k+1-n\) 个,所以最后答案为:
\[\frac{(k+1)^{n-1}(k+1-n)}{k^n} \]这就是拉链为环前莫名其妙地新增一个座位 \(k+1\) 的原因.
然后这个题非常恶心,要用高精,化简分数时要用高精除,这里考虑一种简单的方法:
显然 \(k+1\) 与 \(k\) 互质,只能化简 \(\dfrac{k+1-n}{k^n}\).
这里 \(k+1-n\) 为低精,我们提前做一次 \(\gcd\) 把 \(k^n\) 膜 \(k+1-n\) 转换为低精,就可以低精求 \(\gcd\) 了.
也就是说,最后我们只需要高精乘,高精除低精和高精膜低精即可.
代码
点击查看代码
namespace SOLVE { typedef long double ldb; typedef long long ll; typedef double db; const ll N = 110, B = 10000000; // Base ll T, n, k; class BigNum { public: ll num[N]; inline void Print () { printf ("%lld", num[num[0]]); for_ (i, num[0] - 1, 1) printf ("%07lld", num[i]); return; } inline void Clear () { memset (num, 0, sizeof (num)); num[0] = 1; return; } inline void In (ll number) { num[1] = number; } BigNum operator * (ll ano) { BigNum ans; ans.Clear (); ans.num[0] = num[0]; _for (i, 1, num[0]) { ans.num[i] += num[i] * ano; ans.num[i + 1] += ans.num[i] / B; ans.num[i] %= B; } while (ans.num[ans.num[0] + 1]) ++ans.num[0]; return ans; } BigNum operator * (BigNum ano) { BigNum ans;ans.Clear (); ans.num[0] = num[0] + ano.num[0] - 1; _for (i, 1, num[0]) { _for (j, 1, ano.num[0]) { ans.num[i + j - 1] += num[i] * ano.num[j]; ans.num[i + j] += ans.num[i + j - 1] / B; ans.num[i + j - 1] %= B; } } while (ans.num[ans.num[0] + 1]) ++ans.num[0]; return ans; } inline BigNum operator / (ll ano) { BigNum ans (*this); for_ (i, num[0], 1) { if (i > 1) ans.num[i - 1] += (ans.num[i] % ano) * B; ans.num[i] /= ano; } while (!ans.num[ans.num[0]] && ans.num[0] > 1) --ans.num[0]; return ans; } inline ll operator % (ll ano) { ll ans = 0; for_ (i, num[0], 1) { ans = ans * B % ano; ans += num[i] % ano; } return ans; } } a, b; inline ll Gcd (ll a, ll b) { if (!b) return a; return Gcd (b, a % b); } inline BigNum FastPow (BigNum a, ll b) { BigNum ans;ans.Clear (); ans.num[0] = ans.num[1] = 1; while (b) { if (b & 1) ans = ans * a; a = a * a, b >>= 1; } return ans; } inline ll rnt () { ll x = 0, w = 1; char c = getchar (); while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); } while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar (); return x * w; } inline void In () { n = rnt (), k = rnt (); return; } inline void Solve () { a.Clear (), b.Clear (); a.num[1] = k + 1, b.num[1] = k; a = FastPow (a, n - 1) * (k + 1 - n); b = FastPow (b, n); ll c = b % (k + 1 - n); ll g = Gcd (k + 1 - n, c); a = a / g, b = b / g; return; } inline void Out () { a.Print (), putchar (' '); b.Print (), puts (""); return; } }中国剩余定理
【模板】中国剩余定理(CRT)/ 曹冲养猪
思路
模板题.
代码
点击查看代码
namespace SOLVE { typedef long double ldb; typedef long long ll; typedef double db; const ll N = 20; ll n, a[N], m[N], M = 1, ans; inline ll rnt () { ll x = 0, w = 1; char c = getchar (); while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); } while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar (); return x * w; } inline void exgcd (ll a, ll b, ll& x, ll& y) { if (!b) { x = 1, y = 0; return; } exgcd (b, a % b, x, y); ll _x = x; x = y, y = _x - (a / b) * y; return; } inline void In () { n = rnt (); _for (i, 1, n) { m[i] = rnt (), a[i] = rnt (); M *= m[i]; } return; } inline void Solve () { _for (i, 1, n) { ll Mi = M / m[i], inv, y; exgcd (Mi, m[i], inv, y); ans = (ans + a[i] * Mi % M * (inv + m[i]) % M) % M; } return; } inline void Out () { printf ("%lld\n", ans); return; } }Strange Way to Express Integers
思路
EXCRT 模板题.
代码
点击查看代码
namespace SOLVE { typedef long double ldb; typedef long long ll; typedef double db; const ll N = 1e5 + 10; ll n, a[N], m[N], b, M; inline ll rnt () { ll x = 0, w = 1; char c = getchar (); while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); } while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar (); return x * w; } ll Exgcd (ll a, ll b, ll& x, ll& y) { if (!b) { x = 1, y = 0; return a; } ll g = Exgcd (b, a % b, x, y), _x = x; x = y, y = _x - (a / b) * y; return g; } inline ll Lcm (ll a, ll b) { return a * b / std::__gcd (a, b); } inline ll FastMul (ll a, ll b, ll MOD) { ll ans = 0; while (b) { if (b & 1) ans = (ans + a) % MOD; a = (a + a) % MOD, b >>= 1; } return (ans + MOD) % MOD; } inline void In () { n = rnt (); _for (i, 1, n) m[i] = rnt (), a[i] = rnt (); return; } inline ll EXCRT () { b = a[1], M = m[1]; _for (i, 2, n) { ll x, y, num = (a[i] - b % m[i] + m[i]) % m[i]; ll g = Exgcd (M, m[i], x, y); if (num % g) return -1; b += M * FastMul(x, num / g, m[i] / g); M *= m[i] / g, b = (b % M + M) % M; } return (b % M + M) % M; } inline void Out () { printf ("%lld\n", EXCRT()); return; } }礼物
思路
式子显然是:
\[\prod_{i=1}^{m}\dbinom{n-sum_{i-1}}{w_i}\bmod{P} \]直接扩卢即可.
代码
点击查看代码
const ll N = 1e5 + 10, INF = 1ll << 40; namespace MathBasic { inline void GetFactor (ll x, std::vector& f1, std::vector & f2) { f1.push_back (0), f2.push_back (0); _for (i, 2, x) { if (!(x % i)) { f1.push_back (i), f2.push_back (0); while (!(x % i)) ++f2[f2.size () - 1], x /= i; } } return; } inline ll FastPow (ll a, ll b, ll Mod = INF) { ll ans = 1; while (b) { if (b & 1) ans = ans * a % Mod; a = a * a % Mod, b >>= 1; } return ans; } ll ExGcd (ll a, ll b, ll& x, ll& y) { if (!b) { x = 1, y = 0; return a; } ll g = ExGcd (b, a % b, x, y), _x = x; x = y, y = _x - (a / b) * y; return g; } inline ll Inv (ll a, ll P) { ll x, y; ExGcd (a, P, x, y); return (x % P + P) % P; } } namespace EXLUCAS { using namespace MathBasic; inline ll FDP (ll x, ll P, ll pk) { if (!x) return 1; ll ans = 1; _for (i, 1, pk) if (i % P) ans = ans * i % pk; ans = FastPow (ans, x / pk, pk); _for (i, 1, x % pk) if (i % P) ans = ans * i % pk; return ans * FDP (x / P, P, pk) % pk; } inline ll Index (ll x, ll P) { if (x < P) return 0; return (x / P) + Index (x / P, P); } ll a[N], md[N], P; std::vector p, k; inline void Pre (ll _P) { GetFactor (_P, p, k); P = _P; return; } inline ll ExLucas (ll n, ll m) { ll ans = 0, sz = p.size () - 1; _for (i, 1, sz) { md[i] = FastPow (p[i], k[i]); a[i] = FDP (n, p[i], md[i]) * Inv (FDP (m, p[i], md[i]), md[i]) % md[i] * Inv (FDP (n - m, p[i], md[i]), md[i]) % md[i]; a[i] = a[i] * FastPow (p[i], Index (n, p[i]) - Index (m, p[i]) - Index (n - m, p[i]), md[i]) % md[i]; ans = (ans + a[i] * (P / md[i]) % P * Inv (P / md[i], md[i]) % P) % P; } return ans; } } namespace SOLVE { ll P, n, m, w[10], sum[10], ans = 1; inline ll rnt () { ll x = 0, w = 1; char c = getchar (); while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); } while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar (); return x * w; } inline void In () { P = rnt (), n = rnt (), m = rnt (); _for (i, 1, m) { w[i] = rnt (); sum[i] = sum[i - 1] + w[i]; } return; } inline void Solve () { EXLUCAS::Pre (P); _for (i, 1, m) { if (n - sum[i - 1] < w[i]) { ans = -1; return; } ans = ans * EXLUCAS::ExLucas (n - sum[i - 1], w[i]) % P; } return; } inline void Out () { printf ((ans == -1) ? "Impossible\n" : "%lld\n", ans); return; } } [SDOI2010]古代猪文
思路
(为啥 SDOI2010 经典题这么多啊,猪国杀和 \(k\) 短路模板也是这里出的)
数论全家桶.
显然式子为:
\[g^{\sum_{n|k}\binom{n}{k}}\bmod{999911659} \]用一个费马小定理:
\[g^{\sum_{n|k}\binom{n}{k}\bmod{999911658}}\bmod{999911659} \]那么现在的问题就是:如何求出 \(\sum_{n|k}\binom{n}{k}\bmod{999911658}\)?
仔细看两眼发现就是个扩卢.
然后就出来了.
(其实 \(999911658\) 就是 \(2, 3, 4679, 35617\) 这几个质数之积,可以简化一点写.)
代码
点击查看代码
const ll N = 50000, INF = 1ll << 40; namespace MathBasic { inline ll FastPow (ll a, ll b, ll MOD = INF) { ll ans = 1; while (b) { if (b & 1) ans = ans * a % MOD; a = a * a % MOD, b >>= 1; } return ans; } inline ll ExGcd (ll a, ll b, ll& x, ll& y) { if (!b) { x = 1, y = 0; return a; } ll g = ExGcd (b, a % b, x, y), _x = x; x = y, y = _x - y * (a / b); return g; } inline ll Inv (ll a, ll P) { ll x, y; ExGcd (a, P, x, y); return (x % P + P) % P; } } namespace EXLUCAS { using namespace MathBasic; ll a[5], p[5] = { 0, 2, 3, 4679, 35617 }, fac[5][N], q[5]; ll FDP (ll x, ll P, ll qwq) { if (!x) return 1; ll ans = FastPow (fac[qwq][P - 1], x / P, P); ans = ans * fac[qwq][x % P] % P; return ans * FDP (x / P, P, qwq) % P; } ll Index (ll x, ll P) { if (x < P) return 0; return (x / P) + Index (x / P, P); } inline void Pre () { fac[1][0] = fac[2][0] = fac[3][0] = fac[4][0] = 1; _for (k, 1, 4) { _for (i, 1, 36000) fac[k][i] = fac[k][i - 1] * i % p[k]; q[k] = (999911658 / p[k]) * Inv (999911658 / p[k], p[k]); } return; } inline ll ExLucas (ll n, ll m, ll P) { ll ans = 0; _for (i, 1, 4) { a[i] = FDP (n, p[i], i) * Inv (FDP (m, p[i], i), p[i]) % P * Inv (FDP (n - m, p[i], i), p[i]) % p[i]; a[i] = a[i] * FastPow (p[i], Index (n, p[i]) - Index (m, p[i]) - Index (n - m, p[i])) % p[i]; ans = (ans + a[i] * q[i] % P) % P; } return ans; } } namespace SOLVE { ll n, g, idx, P = 999911659; inline ll rnt () { ll x = 0, w = 1; char c = getchar (); while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); } while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar (); return x * w; } inline void In () { n = rnt (), g = rnt (); EXLUCAS::Pre (); return; } inline void Solve () { for (ll i = 1; i * i <= n; ++i) { if (n % i) continue; idx = (idx + EXLUCAS::ExLucas (n, i, P - 1)) % (P - 1); if (i * i != n) idx = (idx + EXLUCAS::ExLucas (n, n / i, P - 1)) % (P - 1); } return; } inline void Out () { printf ("%lld\n", g == P ? 0 : MathBasic::FastPow (g, idx, P)); return; } }Reference
- 排列组合
——OI-Wiki
——Rolling_star- 浅谈Lucas定理应用及组合数建模
——BerryKanry
——GXZlegend
——Pycr- 中国剩余定理
——OI-Wiki- 扩展Lucas定理
——HorizonWind