「學習筆記」組合計數與中國剩余定理
點擊查看目錄
目錄
- 「學習筆記」組合計數與中國剩余定理
- 知識點
- 排列
- 錯排列
- 組合數
- 式子
- 一些性質
- 盧卡斯定理
- 諤項式定理
- 諤項式反演
- 形式零
- 形式一
- 形式諤
- 小技巧：線性推階乘逆元
- 中國剩余定理（CRT）
- 做法
- 證明
- EXCRT
- ExLucas
- 問題
- 拆為 CRT
- 構造余數
- 建構式
- 代碼
- 例題
- 排列組合
- 排隊
- 題意
- 思路
- Code
- Combination
- 思路
- Code
- [SDOI2016]排列計數
- 思路
- 代碼
- [ZJOI2010]排列計數
- 思路
- 代碼
- BZOJ2839 集合計數
- 思路
- 代碼
- 牡牛和牝牛
- 思路
- 代碼
- 序列統計
- 思路
- 代碼
- [SDOI2009] 虔誠的墓主人
- 思路
- 代碼
- [SDOI2010]地精部落
- 思路
- 代碼
- [ZJOI2011]看電影
- 思路
- 代碼
- 中國剩余定理
- 【模板】中國剩余定理（CRT）/ 曹沖養豬
- 思路
- 代碼
- Strange Way to Express Integers
- 思路
- 代碼
- 禮物
- 思路
- 代碼
- [SDOI2010]古代豬文
- 思路
- 代碼
- Reference
Warning:
本文內含有大量 $\LaTeX$ 公式，可能會引起不適與頁面卡頓.
知識點
排列
\(n\) 個元素里選出 \(m\) 個元素排成一列的方案數.
計算公式：
\[A_{n}^{m}=n(n-1)(n-2)···(n-m+1)=\frac{n!}{(n-m)!}\tag{1} \]還有一種特殊情況：全排列
\[A_{n}^{n}=n!\tag{2} \]錯排列
求對于一個排列 \(1\sim n\)，滿足任意 \(i\) 都不在第 \(i\) 位上的排列有多少個.
遞推公式為：
\[D_{n}=(n-1)(D_{n-1}+D_{n-2})\tag{3} \]證明點這里
我們先把 \(n\) 放在第 \(n\) 位，然后對于任意一個有 \(n-1\) 個數的排列，我們分情況討論：
- 這個排列滿足要求：隨便找一個數和 \(n\) 換，方案數為 \((n-1)D_{n-1}\).
- 這個排列有且只有一位 \(k(1\le k\le n-1)\) 不滿足要求：把 \(k\) 和 \(n\) 交換過來，方案數為 \((n-1)D_{n-2}\).
- 這個排列有 \(k(2\le k\le n-1)\) 位不滿足要求：不可能一次換完，應該在計算 \(D_{n-1}\) 時就已經換完了.
合并一下，總方案數為：
\[D_{n}=(n-1)(D_{n-1}+D_{n-2}) \]組合數
式子
\(n\) 個元素里選出一個 \(m\) 個元素的集合的方案數.
組合與排列的區別：組合沒有順序
計算公式：
\[\dbinom{n}{m}=\frac{A^{m}_{n}}{A^{m}_{m}}=\frac{n!}{m!(n-m)!}\tag{4} \]可以理解為排列數去掉順序.
一些性質
\[\dbinom{n}{m}=\dbinom{n}{n-m}=\dbinom{n-1}{m}+\dbinom{n-1}{m-1}\tag{5} \]盧卡斯定理
用于求較大的且模數 \(p\in\mathbb{P}\) 的組合數.
公式：
\[\dbinom{n}{m}\bmod p=\dbinom{\lfloor n/p\rfloor}{\lfloor m/p\rfloor}\cdot\dbinom{n\bmod p}{m\bmod p}\bmod p\tag{6} \]諤項式定理
公式：
\[(a+b)^n=\sum_{i=0}^{n}\dbinom{n}{i}a^{i}b^{n-i}\tag{7} \]點擊查看證明
數學歸納法.
\[\begin{aligned} (a+b)^1&=a^0b^1+a^1b^0=a+b\\ (a+b)^{n+1} &=(a+b)\sum_{i=0}^{n}\dbinom{n}{i}a^{i}b^{n-i}\\ &=\sum_{i=0}^{n}\dbinom{n}{i}a^{i+1}b^{n-i}+\sum_{i=0}^{n}\dbinom{n}{i}a^{i}b^{n+1-i}\\ &=\sum_{i=1}^{n+1}\dbinom{n}{i-1}a^{i}b^{n+1-i}+\sum_{i=0}^{n}\dbinom{n}{i}a^{i}b^{n+1-i}\\ &=a^{n+1}b^{0}+a^{0}b^{n+1}+\sum_{i=1}^{n}\left(\dbinom{n}{i-1}+\dbinom{n}{i}\right)a^{i}b^{n+1-i}\\ &=\dbinom{n+1}{n+1}a^{n+1}b^{0}+\dbinom{n+1}{0}a^{0}b^{n+1}+\sum_{i=1}^{n}\dbinom{n+1}{i}a^{i}b^{n+1-i}\\ &=\sum_{i=0}^{n+1}\dbinom{n+1}{i}a^{i}b^{n+1-i}\\ \end{aligned} \]組合意義
顯然展開后每一項都是 \(n\) 次的.
那么考慮構成 \(a\) 為 \(i\) 次的項的方案數.
顯然為從 \(n\) 個 \(a\) 里選出 \(i\) 個 \(a\) 的方案數，即 \(\dbinom{n}{i}\).
那么這一項為 \(\dbinom{n}{i}a^{i}b^{n-i}\).
全部展開后就是 \(\sum_{i=0}^{n}\dbinom{n}{i}a^{i}b^{n-i}\).
諤項式反演
相當有趣但相當長.
這里掛幾個式子，以后有機會單獨整理一下學習筆記并掛在這里.
Good Reference
形式零
\[f(n)=\sum_{i=0}^{n}(-1)^i\dbinom{n}{i}g(i)\Longleftrightarrow g(n)=\sum_{i=0}^{n}(-1)^i\dbinom{n}{i}f(i)\tag{8} \]形式一
\[f(n)=\sum_{i=0}^{n}\dbinom{n}{i}g(i)\Longleftrightarrow g(n)=\sum_{i=0}^{n}(-1)^{n-i}\dbinom{n}{i}f(i)\tag{9} \]形式諤
\[f(n)=\sum_{i=n}^{m}\dbinom{i}{n}g(i)\Longleftrightarrow g(n)=\sum_{i=n}^{m}(-1)^{i-n}\dbinom{i}{n}f(i)\tag{10} \]小技巧：線性推階乘逆元
好像是 SoyTony 教的我，%%%
求完階乘后求出最后一個的逆元，然后用一點階乘的逆元的性質，反推回去即可.
在組合題中及其常用，可以線性預處理后直接 \(\Theta(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）
解決如下形式的方程組：
\[\begin{cases} x \equiv a_1 \pmod{m_1}\\ x \equiv a_2 \pmod{m_2}\\ \dots\\ x \equiv a_n \pmod{m_n}\\ \end{cases} \]其中，\(x,a_i,m_i\) 均為正整數，保證 \(m_i\) 互質.
做法
令 \(M=\prod_{i=1}^{n}m_i\)，\(n_i=\frac{M}{m_i}\)，\(n_i^{-1}\) 表示 \(n_i\) 在膜 \(m_i\) 意義下的逆元.
答案是：
\[\sum_{i=1}^{n}a_i n_i n_i^{-1}\pmod{M} \]注意：\(m_i\) 互質但 \(m_i\) 不一定是質數，求逆元不能用費馬小定理，只能用擴歐！
證明
首先解這樣一組方程：
\[\begin{cases} x_i \equiv 0 \pmod{m_1}\\ \dots\\ x_i \equiv a_i \pmod{m_i}\\ \dots\\ x_i \equiv 0 \pmod{m_n}\\ \end{cases} \]顯然，\(x_i=a_i n_i n_i^{-1}\) 是一組合法解：
-
對于第 \(i\) 組方程，由于在膜 \(m_i\) 意義下 \(n_i n_i^{-1}=1\)，所以 \(a_i n_i n_i^{-1} \equiv a_i \pmod{m_i}\).
-
對于第 \(j(j\neq i)\) 組方程，由于 \(n_i\) 為 \(m_j\) 倍數，所以 \(a_j n_j n_j^{-1} \equiv 0 \pmod{m_j}\).
那么如何合并？可以發現 \(x_i+x_j\) 仍是原方程的解，因此答案就是 \(\sum_{i=1}^{n}x_i\pmod{M}=\sum_{i=1}^{n}a_i n_i n_i^{-1}\pmod{M}\).
EXCRT
解決問題沒有變，但 \(m_i\) 不互質.
該演算法主要運用合并的思想.
首先解這個方程：
\[\begin{cases} x \equiv a_1 \pmod{m_1}\\ x \equiv a_2 \pmod{m_2}\\ \end{cases} \]其中 \(m_1,m_2\) 不互質.
轉化為不定方程：
\[\begin{aligned} x &= a_1 + p \cdot m_1\\ &= a_2 + q \cdot m_2\\ \end{aligned} \]進行一個移項：
\[p \cdot m_1 - q \cdot m_2 = a_2 - a_1 \]當 \(a_2-a_1\) 不能被 \(\gcd(m_1,m_2)\) 整除時，原方程無解，否則一定可以解出來一組 \(p,q\).
那么最后可以合并得到一個方程：
\[x \equiv p \cdot m_1 + a_1 \pmod{\operatorname{lcm}(m_1,m_2)} \]按這種方法兩兩合并，最后得到答案.
ExLucas
因為前置比較多，所以放到最后寫.
問題
求：
\[\dbinom{n}{m} \bmod{P} (P\in \mathbb{N^*}) \]拆為 CRT
\(P\) 不一定是質數，那么我們考慮將它拆成質數解決.
設 \(P=p_1^{k_1}\cdot p_2^{k_2}\cdot p_3^{k_3}\cdots p_t^{k_t}(p_i\in \mathbb{P})\).
那么可以列出方程：
\[\begin{cases} x \equiv \dbinom{n}{m} \pmod{p_1^{k_1}}\\ x \equiv \dbinom{n}{m} \pmod{p_2^{k_2}}\\ \dots\\ x \equiv \dbinom{n}{m} \pmod{p_t^{k_t}}\\ \end{cases} \]可以發現 \(x\) 就是最終結果，這里可以使用 CRT 解決.
那么現在把焦點放在這個方程上：
\[x \equiv \dbinom{n}{m} \pmod{p^k}\\ \]構造余數
我們現在需要解決的是計算 \(\dbinom{n}{m} \bmod{p^k}\)，即 \(\dfrac{n!}{m!(n-m)!} \bmod{p^k}\).
然而 \(n!,m!\) 和 \((n-m)!\) 可能與 \(p^k\) 不互質，不能直接求逆元，那么我們就把它們的因數中的 \(p\) 提前去掉.
我們設 \(g(n)\) 表示 \(n\) 的因數里有多少個 \(p\)， \(f(n)\) 表示 \(\dfrac{n!}{p^{g(n)}}\).
那么原式就是：
\[\dfrac{f(n)}{f(m)f(n-m)} \times P^{g(n)-g(m)-g(n-m)} \bmod{p^k} \]建構式
現在我們只需要解決函式 \(f(n)\) 和 \(g(n)\) 即可，那么如何計算？
首先我們拆分一下 \(n!\):
\[\begin{aligned} n! &=1\cdot2\cdot3\cdots n&\pmod{p^k}\\ &=(p\cdot2p\cdot3p\cdots\left\lfloor\frac{n}{p}\right\rfloor p)(1\cdot2\cdot3\cdots n)&\pmod{p^k}\\ &=p^{\left\lfloor\frac{n}{p}\right\rfloor}(\left\lfloor\frac{n}{p}\right\rfloor)!(\prod_{i=1,i\bmod p\neq0}^{n}i)&\pmod{p^k}\\ &=p^{\left\lfloor\frac{n}{p}\right\rfloor}(\left\lfloor\frac{n}{p}\right\rfloor)!(\prod_{i=1,i\bmod p\neq0}^{p^k}i)(\prod_{i=p^k+1,i\bmod p\neq0}^{2p^k}i)\cdots(\prod_{i=p^k(\left\lfloor\frac{n}{p^k}\right\rfloor-1)+1,i\bmod{p}\neq0}^{p^k\left\lfloor\frac{n}{p^k}\right\rfloor}i)(\prod_{i=p^k\left\lfloor\frac{n}{p^k}\right\rfloor+1,i\bmod{p}\neq0}^{n}i)&\pmod{p^k}\\ &=p^{\left\lfloor\frac{n}{p}\right\rfloor}(\left\lfloor\frac{n}{p}\right\rfloor)!(\prod_{i=1,i\bmod p\neq0}^{p^k}i)(\prod_{i=1,i\bmod{p}\neq0}^{p^k}i)\cdots(\prod_{i=1,i\bmod{p}\neq0}^{p^k}i)(\prod_{i=p^k\left\lfloor\frac{n}{p^k}\right\rfloor+1,i\bmod{p}\neq0}^{n}i)&\pmod{p^k}\\ &=p^{\left\lfloor\frac{n}{p}\right\rfloor}(\left\lfloor\frac{n}{p}\right\rfloor)!(\prod_{i=1,i\bmod{p}\neq0}^{p^k}i)^{\left\lfloor\frac{n}{p^k}\right\rfloor}(\prod_{i=1,i\bmod{p}\neq0}^{n\bmod{p^k}}i)&\pmod{p^k} \end{aligned} \]我們的目的是除去 \(p\)，因此 \(p^{\left\lfloor\frac{n}{p}\right\rfloor}\) 應去除.\((\left\lfloor\frac{n}{p}\right\rfloor)!\) 里還可能會有 \(p\)，所以最后式子為：
\[f(n)=f(\left\lfloor\frac{n}{p}\right\rfloor)(\prod_{i=1,i\bmod{p}\neq0}^{p^k}i)^{\left\lfloor\frac{n}{p^k}\right\rfloor}(\prod_{i=1,i\bmod{p}\neq0}^{n\bmod{p^k}}i)\pmod{p^k} \]邊界為 \(f(0)=1\).
從剛才的式子也可以看出來，每次遞推會誕生 \(\left\lfloor\frac{n}{p}\right\rfloor\) 個 \(p\)，由于它還在往下遞推，所以還會產生 \(g(\left\lfloor\frac{n}{p}\right\rfloor)\) 個 \(p\).
那么遞推式為：
\[g(n)=\left\lfloor\frac{n}{p}\right\rfloor+g(\left\lfloor\frac{n}{p}\right\rfloor) \]代碼
點擊查看代碼
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\le 2000\)
思路
要分類討論.
如果兩名老師被男生隔開，則方案數為 男生的排列乘老師的排列乘女生的排列
即：
\[n!\times A_{n+1}^{2}\times A_{n+3}^{m} \]如果兩名老師被女生隔開，則 用來隔開老師的女生 乘上 老師的排列數 再乘上 這兩名老師與這名女生可插的空 的方案數為 \(2m(n+1)\)，再乘上男生的排列數和剩余女生的排列數得到：
\[2m(n+1)\times n!\times A_{n+2}^{m-1} \]合并一下可得：
\[n!(A_{n+1}^{2}\times A_{n+3}^{m}+2m(n+1)\times A_{n+2}^{m-1}) \]但是由于 \(n,m\le 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\) 個數為穩定的，方案數為 \(\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 <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^{\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- 學習筆記-組合數學——組合數基礎，組合數取模，CRT
——Rolling_star- 淺談Lucas定理應用及組合數建模
——BerryKanry- 諤項式反演
——GXZlegend- 中國剩余定理(CRT)的證明
——Pycr- 中國剩余定理
——OI-Wiki- 擴展Lucas定理
——HorizonWind
本文來自博客園，作者：Keven-He，轉載請注明原文鏈接：https://www.cnblogs.com/Keven-He/p/CombinationAndCRT.html
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/527850.html
標籤：其他
上一篇：【小記】二八十十六，進制團團轉
下一篇：演算法題--用兩個堆疊實作佇列