「学习笔记」组合计数与中国剩余定理
点击查看目录
Warning:
本文内含有大量 LATEX 公式，可能会引起不适与页面卡顿.
#define 谔 二
知识点#
排列#
n 个元素里选出 m 个元素排成一列的方案数.
计算公式：
Amn=n(n−1)(n−2)···(n−m+1)=n!(n−m)!
还有一种特殊情况：全排列
Ann=n!
错排列#
求对于一个排列 1∼n，满足任意 i 都不在第 i 位上的排列有多少个.
递推公式为：
Dn=(n−1)(Dn−1+Dn−2)
证明点这里
我们先把 n 放在第 n 位，然后对于任意一个有 n−1 个数的排列，我们分情况讨论：
- 这个排列满足要求：随便找一个数和 n 换，方案数为 (n−1)Dn−1.
- 这个排列有且只有一位 k(1≤k≤n−1) 不满足要求：把 k 和 n 交换过来，方案数为 (n−1)Dn−2.
- 这个排列有 k(2≤k≤n−1) 位不满足要求：不可能一次换完，应该在计算 Dn−1 时就已经换完了.
合并一下，总方案数为：
Dn=(n−1)(Dn−1+Dn−2)
组合数#
式子#
n 个元素里选出一个 m 个元素的集合的方案数.
组合与排列的区别：组合没有顺序
计算公式：
(nm)=AmnAmm=n!m!(n−m)!
可以理解为排列数去掉顺序.
一些性质#
(nm)=(nn−m)=(n−1m)+(n−1m−1)
卢卡斯定理#
用于求较大的且模数 p∈P 的组合数.
公式：
(nm)modp=(⌊n/p⌋⌊m/p⌋)⋅(nmodpmmodp)modp
谔项式定理#
公式：
(a+b)n=n∑i=0(ni)aibn−i
点击查看证明
数学归纳法.
(a+b)1=a0b1+a1b0=a+b(a+b)n+1=(a+b)n∑i=0(ni)aibn−i=n∑i=0(ni)ai+1bn−i+n∑i=0(ni)aibn+1−i=n+1∑i=1(ni−1)aibn+1−i+n∑i=0(ni)aibn+1−i=an+1b0+a0bn+1+n∑i=1((ni−1)+(ni))aibn+1−i=(n+1n+1)an+1b0+(n+10)a0bn+1+n∑i=1(n+1i)aibn+1−i=n+1∑i=0(n+1i)aibn+1−i
组合意义
显然展开后每一项都是 n 次的.
那么考虑构成 a 为 i 次的项的方案数.
显然为从 n 个 a 里选出 i 个 a 的方案数，即 (ni).
那么这一项为 (ni)aibn−i.
全部展开后就是 ∑ni=0(ni)aibn−i.
谔项式反演#
相当有趣但相当长.
这里挂几个式子，以后有机会单独整理一下学习笔记并挂在这里.
形式零#
f(n)=n∑i=0(−1)i(ni)g(i)⟺g(n)=n∑i=0(−1)i(ni)f(i)
形式一#
f(n)=n∑i=0(ni)g(i)⟺g(n)=n∑i=0(−1)n−i(ni)f(i)
形式谔#
f(n)=m∑i=n(in)g(i)⟺g(n)=m∑i=n(−1)i−n(in)f(i)
小技巧：线性推阶乘逆元#
好像是 SoyTony 教的我，%%%
求完阶乘后求出最后一个的逆元，然后用一点阶乘的逆元的性质，反推回去即可.
在组合题中及其常用，可以线性预处理后直接 Θ(1) 求组合数.
点击查看代码
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;
}
中国剩余定理（CRT）#
解决如下形式的方程组：
{x≡a1(modm1)x≡a2(modm2)…x≡an(modmn)
其中，x,ai,mi 均为正整数，保证 mi 互质.
做法#
令 M=∏ni=1mi，ni=Mmi，n−1i 表示 ni 在膜 mi 意义下的逆元.
答案是：
n∑i=1ainin−1i(modM)
注意：mi 互质但 mi 不一定是质数，求逆元不能用费马小定理，只能用扩欧！
证明#
首先解这样一组方程：
{xi≡0(modm1)…xi≡ai(modmi)…xi≡0(modmn)
显然，xi=ainin−1i 是一组合法解：
-
对于第 i 组方程，由于在膜 mi 意义下 nin−1i=1，所以 ainin−1i≡ai(modmi).
-
对于第 j(j≠i) 组方程，由于 ni 为 mj 倍数，所以 ajnjn−1j≡0(modmj).
那么如何合并？可以发现 xi+xj 仍是原方程的解，因此答案就是 ∑ni=1xi(modM)=∑ni=1ainin−1i(modM).
EXCRT#
解决问题没有变，但 mi 不互质.
该算法主要运用合并的思想.
首先解这个方程：
{x≡a1(modm1)x≡a2(modm2)
其中 m1,m2 不互质.
转化为不定方程：
x=a1+p⋅m1=a2+q⋅m2
进行一个移项：
p⋅m1−q⋅m2=a2−a1
当 a2−a1 不能被 gcd(m1,m2) 整除时，原方程无解，否则一定可以解出来一组 p,q.
那么最后可以合并得到一个方程：
x≡p⋅m1+a1(modlcm(m1,m2))
按这种方法两两合并，最后得到答案.
ExLucas#
因为前置比较多，所以放到最后写.
问题#
求：
(nm)modP(P∈N∗)
拆为 CRT#
P 不一定是质数，那么我们考虑将它拆成质数解决.
设 P=pk11⋅pk22⋅pk33⋯pktt(pi∈P).
那么可以列出方程：
{x≡(nm)(modpk11)x≡(nm)(modpk22)…x≡(nm)(modpktt)
可以发现 x 就是最终结果，这里可以使用 CRT 解决.
那么现在把焦点放在这个方程上：
x≡(nm)(modpk)
构造余数#
我们现在需要解决的是计算 (nm)modpk，即 n!m!(n−m)!modpk.
然而 n!,m! 和 (n−m)! 可能与 pk 不互质，不能直接求逆元，那么我们就把它们的因数中的 p 提前去掉.
我们设 g(n) 表示 n 的因数里有多少个 p， f(n) 表示 n!pg(n).
那么原式就是：
f(n)f(m)f(n−m)×Pg(n)−g(m)−g(n−m)modpk
构造函数#
现在我们只需要解决函数 f(n) 和 g(n) 即可，那么如何计算？
首先我们拆分一下 n!:
n!=1⋅2⋅3⋯n(modpk)=(p⋅2p⋅3p⋯⌊np⌋p)(1⋅2⋅3⋯n)(modpk)=p⌊np⌋(⌊np⌋)!(n∏i=1,imodp≠0i)(modpk)=p⌊np⌋(⌊np⌋)!(pk∏i=1,imodp≠0i)(2pk∏i=pk+1,imodp≠0i)⋯(pk⌊npk⌋∏i=pk(⌊npk⌋−1)+1,imodp≠0i)(n∏i=pk⌊npk⌋+1,imodp≠0i)(modpk)=p⌊np⌋(⌊np⌋)!(pk∏i=1,imodp≠0i)(pk∏i=1,imodp≠0i)⋯(pk∏i=1,imodp≠0i)(n∏i=pk⌊npk⌋+1,imodp≠0i)(modpk)=p⌊np⌋(⌊np⌋)!(pk∏i=1,imodp≠0i)⌊npk⌋(nmodpk∏i=1,imodp≠0i)(modpk)
我们的目的是除去 p，因此 p⌊np⌋ 应去除.(⌊np⌋)! 里还可能会有 p，所以最后式子为：
f(n)=f(⌊np⌋)(pk∏i=1,imodp≠0i)⌊npk⌋(nmodpk∏i=1,imodp≠0i)(modpk)
边界为 f(0)=1.
从刚才的式子也可以看出来，每次递推会诞生 ⌊np⌋ 个 p，由于它还在往下递推，所以还会产生 g(⌊np⌋) 个 p.
那么递推式为：
g(n)=⌊np⌋+g(⌊np⌋)
代码#
点击查看代码
const ll N = 1e5 + 10, INF = 1ll << 40;
namespace MathBasic {
inline void GetFactor (ll x, std::vector <ll>& f1, std::vector <ll>& f2) {
_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 - 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 FDP (ll x, ll P, ll pk) { //FacDivP
if (x == 0) 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) % pk;
return FDP (x / P, P, pk) * ans % pk;
}
ll Index (ll x, ll P) {
if (x < P) return 0;
return (x / P) + Index (x / P, P);
}
ll a[N], md[N];
inline ll ExLucas (ll n, ll m, ll P) {
std::vector <ll> p, k;
p.push_back (0), k.push_back (0);
GetFactor (P, p, k);
ll len = p.size () - 1, ans = 0;
_for (i, 1, len) {
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 (n - m, p[i]) - Index (m, p[i]), md[i]) % md[i];
}
_for (i, 1, len) {
ll q = P / md[i], x, y;
ExGcd (q, md[i], x, y);
ans = (ans + a[i] * q % P * ((x % P + P) % P) % P) % P;
}
return ans;
}
}
例题#
排列组合#
排队#
题意
n 名男同学，m 名女同学和两名老师要排成一条直线，并且任意两名女同学不能相邻，两名老师也不能相邻，求一共有多少种排法？
n,m≤2000
思路
要分类讨论.
如果两名老师被男生隔开，则方案数为 男生的排列乘老师的排列乘女生的排列
即：
n!×A2n+1×Amn+3
如果两名老师被女生隔开，则 用来隔开老师的女生 乘上 老师的排列数 再乘上 这两名老师与这名女生可插的空 的方案数为 2m(n+1)，再乘上男生的排列数和剩余女生的排列数得到：
2m(n+1)×n!×Am−1n+2
合并一下可得：
n!(A2n+1×Amn+3+2m(n+1)×Am−1n+2)
但是由于 n,m≤2000，要手写高精才能过……
Code
点击查看代码
const ll N = 10000, k = 1000000000;
ll n, m;
class BigNum {
public:
ll len = 0, a[N] = {0};
BigNum () { memset(a, 0, sizeof(a));}
inline void In (ll num){
while (num) {
a[++ len] = num % k;
num /= k;
}
return;
}
inline void Out () {
for_(i, len, 1)
printf((i == len ? "%lld" : "%.9lld"), a[i]);
if (!len) puts("0");
puts("");
return;
}
void operator = (BigNum another) {
len = another.len;
_for(i, 1, len) a[i] = another.a[i];
return;
}
BigNum operator + (BigNum another) {
BigNum answer;
answer.len = max(len, another.len);
_for(i, 1, answer.len) {
answer.a[i] += a[i] + another.a[i];
answer.a[i + 1] += answer.a[i] / k;
answer.a[i] %= k;
}
while (answer.a[answer.len + 1]) ++ answer.len;
return answer;
}
BigNum operator * (BigNum another) {
BigNum answer;
answer.len = len + another.len - 1;
_for(i, 1, answer.len) {
_for(j, 1, another.len) {
answer.a[i + j - 1] += a[i] * another.a[j];
answer.a[i + j] += answer.a[i + j - 1] / k;
answer.a[i + j - 1] %= k;
}
}
while (answer.a[answer.len + 1]) ++ answer.len;
return answer;
}
} mm, nn, ans;
namespace SOLVE {
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();
nn.In(n), mm.In(m);
return;
}
inline void Solve () {
BigNum a, b, one;
one.In(1), a.In(2), b.In(1);
a = a * (nn + one) * mm;
for_(i, n + 2, n - m + 4) {
BigNum ii;
ii.In(i);
a = a * ii;
}
b = nn * (nn + one);
for_(i, n + 3, n - m + 4) {
BigNum ii;
ii.In(i);
b = b * ii;
}
ans = a + b;
_for(i, 2, n) {
BigNum ii;
ii.In(i);
ans = ans * ii;
}
return;
}
inline void Out () {
ans.Out();
return;
}
}
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 个数为稳定的，方案数为 (nm).
在剩下的 n−m 个位置里要保证每个数不稳定.
欸那不就是错排列 (3) 吗？
那么方案数就是 Dn−m.
总方案数就是 (nm)Dn−m，Θ(n) 预处理一下错排列，阶乘与逆元可用 Θ(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]排列计数#
思路
观察一下可以发现满足性质的序列是一个小根堆.
那么设 si 表示以 i 为根的堆的大小，fi 表示以 i 为根的堆的可行方案数（此时该子堆里的序号不是最终序号，而是在子堆内大小的排名，因为归并到父堆时要算分配给子堆不同序号的方案数）.
那么转移方程就是：
si=si∗2+si∗2+1+1fi=(si−1si∗2)fi∗2fi∗2+1
（自己必须是最小的所以只能从 si−1 个序号选 si∗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) 表示交集数量 ≥i 的方案数，g(i) 表示交集个数恰好为 i 个的方案数，那么答案为 g(k).
那么：
f(i)=(ni)(22n−i−1)
即先确定 i 个必选，包含这 i 个的集合数为 2n−k 个，每个集合都可以选或不选但不能一个不选，即 22n−i−1.
同时：
f(k)=n∑i=k(ik)g(i)
等一下这式子是不是在哪里见过?
这不是 (10) 吗？！
那么愉快的套一个谔项式反演：
g(k)=n∑i=k(−1)i−k(ik)f(i)=n∑i=k(−1)i−k(ik)(ni)(22n−i−1)
再加上一点预处理，就可以解决了.
代码
点击查看代码
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×(i−1) 只牡牛被固定住，此时剩下 w(i)=(n−i−k×(i−1))×[i>0]+n×[i=0] 只牡牛可以随便选位置.
观察一下，看上去是只有 k+1 个地方可以插空，然而两只牝牛之间可以放多只牡牛，如何解决这个问题？
既然可以重复放，那我们就把重复放的位置 new 出来！
即把空的个数改为 k+1+(w(i)−1)=k+i.
这样会不会导致选的全都是 new 出来的呢？不会，因为我们只 new 出来了 w(i)−1 个空，剩下的一只牛必然会被放在原有的位置.
那么答案就是：
n∑i=0[w(i)≥0](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 的序列的方案数为 (m+i−1i).
然后推式子：
n∑i=1(m+i−1i)=n∑i=1(m+i−1m−1)+(mm)−1=n∑i=2(m+i−1m−1)+(mm−1)+(mm)−1=n∑i=2(m+i−1m−1)+(m+1m)−1=n∑i=3(m+i−1m−1)+(m+2m)−1=n∑i=4(m+i−1m−1)+(m+3m)−1=⋯=n∑i=n(m+i−1m−1)+(m+n−1m)−1=(m+nm)−1
但 n,m 过大，需要用到 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]地精部落#
思路
fi,0 表示长度为 i 且第一段山为山谷的序列数量.
fi,1 表示长度为 i 且第一段山为山峰的序列数量.
fi,k=i∑j=1[jmod2=k](i−1j−1)fj−1,kfi−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]看电影#
思路
这个式子还是蛮有意思的，但要用高精.
首先答案 =合法情况数总情况数，总情况数显然是 kn，难点在于如何算出合法情况数.
首先我们在 k 后面新增一个座位 k+1 然后拉链为环，让没有座位的人从头开始往后坐，这样一定所有人都会有座位，那么这样的环一共会有 (k+1)n−1 个（即圆排列）.注意：这里暂时不考虑标号.
但是我们怎么判断做法是否合法呢？如果 k+1 这个位置最后有人，那么在没有环与新座位时就一定会有人站在这里，否则没人站着（即合法情况），也就是说 我们在环中找到空的位置当成 k+1 即可，空的位置一共有 k+1−n 个，所以最后答案为：
(k+1)n−1(k+1−n)kn
这就是拉链为环前莫名其妙地新增一个座位 k+1 的原因.
然后这个题非常恶心，要用高精，化简分数时要用高精除，这里考虑一种简单的方法：
显然 k+1 与 k 互质，只能化简 k+1−nkn.
这里 k+1−n 为低精，我们提前做一次 gcd 把 kn 膜 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;
}
}
礼物#
思路
式子显然是：
m∏i=1(n−sumi−1wi)modP
直接扩卢即可.
代码
点击查看代码
const ll N = 1e5 + 10, INF = 1ll << 40;
namespace MathBasic {
inline void GetFactor (ll x, std::vector <ll>& f1, std::vector <ll>& 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 <ll> 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∑n|k(nk)mod999911659
用一个费马小定理：
g∑n|k(nk)mod999911658mod999911659
那么现在的问题就是：如何求出 ∑n|k(nk)mod999911658？
仔细看两眼发现就是个扩卢.
然后就出来了.
（其实 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;
}
}
[POI2008]PER-Permutation#
Reference#
- 排列组合
——OI-Wiki- 学习笔记-组合数学——组合数基础，组合数取模，CRT
——Rolling_star- 浅谈Lucas定理应用及组合数建模
——BerryKanry- 谔项式反演
——GXZlegend- 中国剩余定理(CRT)的证明
——Pycr- 中国剩余定理
——OI-Wiki- 扩展Lucas定理
——HorizonWind
TheEnd