int64_t a, b, p;
void _main() {
cin >> a >> b >> p;
cout << (int64_t) ((__int128) a * b % p);
}
*使用 64 位浮点数
由取模的定义
ab \bmod p=ab-\lfloor \dfrac{ab}{p} \rfloor \times p
因为 p 是 64 位整数，所以结果可以对 2^{64} 取模，也就是 unsigned long long 自然溢出。现在 ab 和 \lfloor \dfrac{ab}{p} \rfloor \times p 都可以直接自然溢出解决。只要求出 \lfloor \dfrac{ab}{p} \rfloor 即可。使用 long double 算出 \dfrac{a}{p} 再乘上 b 即可。
int64_t a, b, p;
void _main() {
cin >> a >> b >> p;
uint64_t c = (uint64_t) a * b - (uint64_t) (1.0L * a / p * b + 0.5L) * p;
cout << (c < uint64_t(p) ? c : c + p);
}
所以把 p_i 减去 i 对 k 取模，然后用哈希表统计即可。注意删掉超出范围的点，以及负数取模。复杂度 O(n)。
const int N = 2e5 + 5;
long long n, k, a[N];
void _main() {
cin >> n >> k;
for (int i = 1; i <= n; i++) cin >> a[i], a[i] += a[i - 1];
long long cnt = 0;
for (int i = 0; i <= n; i++) a[i] = ((a[i] % k - i) % k + k) % k;
unordered_map<int, int> mp;
for (int i = 0; i <= n; i++) {
if (i >= k) mp[a[i - k]]--;
cnt += mp[a[i]], mp[a[i]]++;
} cout << cnt;
}
P1154 奶牛分厩
暴力的做法是枚举 k，然后 O(n) 验证，O(nk) 无法通过。
形式化一下题意，题意即为 \nexists a,b \in S, a \equiv b \pmod p。根据同余的差性质，等价为 \nexists i \in \mathbb{N}^+, pi \mid (a-b) 。所以考虑 O(n^2) 求出所有 s_i-s_j 并标记。
仍然暴力枚举 k，再枚举 k 的倍数 ik，只要判断 ik 是否标记即可。复杂度 O(n^2 + k \log k) 可以通过。
const int N = 5e3 + 5, M = 1e6 + 5;
int n, a[N];
bool tag[M << 1];
void _main() {
cin >> n;
for (int i = 1; i <= n; i++) cin >> a[i];
for (int i = 1; i <= n; i++) {
for (int j = i + 1; j <= n; j++) tag[abs(a[i] - a[j])] = true;
}
for (int i = 1; i < M; i++) {
if (tag[i]) continue;
bool flag = true;
for (int j = i; j < M; j += i) {
if (tag[j]) {flag = false; break;}
}
if (flag) return cout << i, void();
}
}
AT_abc357_d [ABC357D] 88888888
介绍三种做法。法一：设 d 为 n 的位数，则 v(n) 可以看作一个 10^d 进制数，根据等比数列求和：
using ull = unsigned long long;
ull n;
const ull mod = 998244353;
ull power(ull a, ull b) {
ull res = 1;
for (a %= mod; b; b >>= 1) {
if (b & 1) res = res * a % mod;
a = a * a % mod;
}
return res;
}
void _main() {
cin >> n;
ull b = 1;
while (b <= n) b *= 10;
cout << n % mod * (power(b, n) - 1) % mod * power(b - 1, mod - 2) % mod;
}
法二：如果 10^d-1 没有逆元，一个想法就是广义快速幂。定义 a \times a 表示将 a 和 a 首尾拼接所得的数，求出 a^n 即可。复杂度 O(\log^2 n)。
可以枚举 [1,\sqrt{n}] 中的整数，判定其是否为 p 的约数。因为 a \mid p 时 \dfrac{p}{a}\mid p，所以无需枚举到 n。复杂度为 O(\sqrt{n})。
这里给出一个 \dfrac{1}{6} 常数的实现，筛掉 2 和 3 的倍数：
inline bool isprime(int x) {
if (x <= 1) return false;
if (x == 2 || x == 3) return true;
if (x % 6 != 1 && x % 6 != 5) return false;
for (int i = 5; i * i <= x; i += 6) {
if (x % i == 0 || x % (i + 2) == 0) return false;
}
return true;
}
*2.2.2 Miller-Rabin 法
先引入二次探测定理，即：若 p 为奇质数，则 x^2 \equiv 1 \pmod p 的解为 x \equiv 1 或 x \equiv p-1。
接下来使用费马小定理：若 p 为质数且 \gcd(a,p)=1，则 a^{p-1} \equiv 1 \pmod p。这个结论的讲解可到文章下部。
设 p-1=u \times 2^t，随机一个 a 值并求得 v=a^u \bmod p，然后执行 t 次 v \gets v^2 \bmod p，检查是否满足 a^{p-1} \equiv 1 \pmod p。
在 long long 范围内，只需选取 a \in \{2,325,9375,28178,450775,9780504,1795265022\} 即可保证正确性。复杂度 O(\log n)。这里给出一个比较科技的写法。
inline long long power(long long a, long long b, long long p) {
long long res = 1; for (a %= p; b; b >>= 1) {
if (b & 1) res = (__int128) res * a % p;
a = (__int128) a * a % p;
} return res;
}
const long long BASE[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
inline bool miller_rabin(long long n) {
if (n < 2 || n % 6 % 4 != 1) return (n | 1) == 3;
long long s = __builtin_ctzll(n - 1), d = n >> s;
for (long long a : BASE) {
long long p = power(a, d, n), i = s;
while (p != 1 && p != n - 1 && a % n && i--) p = (__int128) p * p % n;
if (p != n - 1 && i != s) return false;
} return true;
}
2.3 威尔逊定理
若 p 为质数，则
(p-1)! \equiv -1 \pmod p
其逆定理也成立。即若 (p-1)! \equiv -1 \pmod p，则 p 为质数。
2.4 算术基本定理
任何一个合数 n 可以唯一分解成有限个质数的乘积。
存在性：用反证法，假设 n 是最小的不能被分解的合数，则存在 n=ab，若 a,b 都可分解，则 n 可以被分解；若 a,b 有不可分解的数，则 a,b 才是最小的数，矛盾。
int p[N], c[N];
int decompose(int n) {
int m = 0;
for (int i = 2; i * i <= n; i++) {
if (n % i) continue;
p[++m] = i, c[m] = 0;
while (n % i == 0) n /= i, c[m]++;
}
if (n > 1) p[++m] = n, c[m] = 1;
return m;
}
需要注意的是，如果知道 n 的值域，可以先用下面的质数筛打出一个质数表。根据质数分布规律，单次试除的复杂度降为 O(\dfrac{\sqrt{n}}{\log n})。
*2.4.3 除数函数
定义除数函数d(n) 为 n 的正约数数目。容易发现 d(n) \le 2\sqrt{n}。
设 n 可以质因数分解为 n=p_1^{c_1} p_2^{c_2} p_3^{c_3} \cdots p_m^{c_m}，根据推论 1 有
d(n)=\prod_{i=1}^{m} (c_i+1)
除数函数 d(n) 的级别远小于 O(\sqrt{n})。事实上，我们有
\sum_{n=1}^{x} \dfrac{d(n)}{n}=\dfrac{1}{2} \log^2 x + 2 \gamma \log x+C+\delta
证明参见这篇知乎讨论。其中 C \approx \dfrac{1}{2}，\delta 为无穷小。这是一个均值估计。
bitset<N> s, p;
void solve(int a, int b) { // 筛出[a, b)之间的质数，用 p[i - a] 判断i是否为质数
s.set(), p.set(), s[0] = s[1] = false;
if (a <= 1) p[1 - a] = false;
int z = (int) sqrt(b) + 1;
for (int i = 2; i < z; i++) {
for (int j = 2; i * j < z; j++) s[i * j] = false;
}
for (int i = 2; i * i <= b; i++) {
if (!s[i]) continue;
for (int j = max(i * i, (a + i - 1) / i * i); j < b; j += i) p[j - a] = false;
}
}
int n, m1, m2, a[N], tot;
int solve(int x) {
int res = 0;
for (int i = 1; i <= tot; i++) {
if (x % p[i]) return INT_MAX;
int cnt = 0;
while (x % p[i] == 0) x /= p[i], cnt++;
res = max(res, c[i] / cnt + (c[i] % cnt != 0));
}
return res;
}
void _main() {
cin >> n >> m1 >> m2;
for (int i = 1; i <= n; i++) cin >> a[i];
tot = decompose(m1);
for (int i = 1; i <= tot; i++) c[i] *= m2;
int res = INT_MAX;
for (int i = 1; i <= n; i++) res = min(res, solve(a[i]));
cout << (res == INT_MAX ? -1 : res);
}
P1865 A % B Problem
由 m \le 10^6 可以筛一下 m 以内的质数，询问用前缀和处理。远古代码太丑了不想放。
P7960 [NOIP2021] 报数
质数筛的思想延伸。把含有 7 的数先暴力找见，然后把它的倍数筛掉。
const int N = 1e7 + 1000; // 注意这里要开大点
int t, x, nxt[N];
bitset<N> dis;
inline bool check(int x) {
for (; x != 0; x /= 10) if (x % 10 == 7) return true;
return false;
}
void _main() {
int last = 0;
dis.reset();
for (int i = 1; i < N; i++) {
if (dis[i]) continue;
if (check(i)) {
for (int j = 1; j * i < N; j++) dis[j * i] = 1;
continue;
}
nxt[last] = i, last = i;
}
for (cin >> t; t--; ) {
cin >> x;
cout << (dis[x] ? -1 : nxt[x]) << '\n';
}
}
AT_abc172_d [ABC172D] Sum of Divisors
不会推公式，只会大力筛法。注意到在埃氏筛标记倍数的过程中可以顺便求出 d(n)，直接做就好了。
const int N = 1e7 + 5;
bitset<N> isprime;
int n, d[N];
void _main() {
cin >> n;
isprime.set(), isprime[1] = 0;
for (int i = 2; i <= n; i++) {
if (!isprime[i]) continue;
int mul = 2;
for (int j = i * 2; j <= n; j += i, mul++) isprime[j] = 0, d[j] += d[mul] + 1;
}
long long res = 0;
d[1] = -1;
for (int i = 1; i <= n; i++) res += 1LL * i * (d[i] + 2);
cout << res;
}
#define int long long
const int N = 1e7 + 5;
int l, r;
bitset<N> s, p;
void solve(int a, int b) {
s.set(), p.set();
int z = (int) sqrt(b) + 1;
for (int i = 2; i < z; i++) {
for (int j = 2; i * j < z; j++) s[i * j] = false;
}
for (int i = 2; i * i <= b; i++) {
if (!s[i]) continue;
for (int j = max(i * i, (a + i - 1) / i * i); j < b; j += i) p[j - a] = false;
}
for (int i = 2; i * i <= b; i++) {
if (!s[i]) continue;
int j = i * i;
while (j < a) j *= i;
for (; j < b; j *= i) p[j - a] = true;
}
}
void _main() {
cin >> l >> r;
if (l == r) return cout << 1, void();
solve(l, r + 1);
int cnt = 0;
for (int i = l; i <= r; i++) cnt += p[i - l];
if (!p[0]) cnt++;
cout << cnt;
}
[模拟赛] 舔狗的付出
给你一个十进制正整数 x，在 x 后面添加若干数字使它成为一个质数，可以不添加，求这个质数的最小值。
const int N = 6725; // d(10^12) <= 6720
long long n, x[N];
mint dp[N][N];
gp_hash_table<int, int> p;
void _main() {
x[0] = 0, p.clear();
cin >> n;
for (long long i = 1; i * i <= n; i++) {
if (n % i) continue;
x[++x[0]] = i;
if (i * i != n) x[++x[0]] = n / i;
}
sort(x + 1, x + x[0] + 1);
for (int i = 1; i <= x[0]; i++) p[x[i]] = i;
for (int i = 1; i <= x[0]; i++) {
dp[i][1] = (i == 1);
for (int j = 2; j <= x[0]; j++) {
dp[i][j] = dp[i][j - 1];
if (x[i] % x[j]) continue;
dp[i][j] += dp[p[x[i] / x[j]]][j - 1];
}
} cout << dp[x[0]][x[0]] - 1 << '\n';
}
CF1878F Vasilije Loves Number Theory
注意到 n=d(a)d(n)，故原问题等价于判断 d(n) \mid n。
用一个 std::map 维护当前分解质因数的结果 {p_1}^{c_1}{p_2}^{c_2}{p_3}^{c_3}\cdots。那么 d(n) 就是 \prod (c_i+1)。因为我们只需判断 n \bmod d(n)=0，所以算 n 的时候对 d(n) 取模即可。用一个快速幂来实现。
int n, q, opt, x;
long long mpow(long long a, long long b, long long p) {
long long res = 1; for (a %= p; b; a = a * a % p, b >>= 1) {
if (b & 1) res = res * a % p;
} return res;
}
void decompose(int x, map<int, int>& mp) {
for (int i = 2; i * i <= x; i++) {
while (x % i == 0) mp[i]++, x /= i;
}
if (x != 1) mp[x]++;
}
void _main() {
cin >> n >> q;
map<int, int> mp;
decompose(n, mp);
map<int, int> cur = mp;
while (q--) {
cin >> opt;
if (opt == 1) {
cin >> x;
decompose(x, cur);
long long a = 1, b = 1;
for (const auto& i : cur) b *= i.second + 1;
for (const auto& i : cur) a = a * mpow(i.first, i.second, b) % b;
cout << (a % b ? "NO\n" : "YES\n");
} else cur = mp;
} cout << '\n';
}
template <class T> constexpr T gcd(T a, T b) {
int az = ctz(a), bz = ctz(b), z = min(az, bz);
for (b >>= bz; a; ) {
a >>= az; int d = a - b;
az = ctz(a - b), b = min(a, b), a = abs(d);
} return b << z;
}
由上述结论得：\gcd(\dfrac{a_0}{a_1},\dfrac{x}{a_1})=1。又有 \operatorname{lcm}(x,b_0)=b_1，故 \gcd(\dfrac{b_1}{x},\dfrac{b_1}{b_0})=1。且 x\mid b_1，枚举 x 并判断即可，复杂度 O(\sqrt{b_1} \log V)。
int a0, a1, b0, b1;
void _main() {
cin >> a0 >> a1 >> b0 >> b1;
int cnt = 0;
for (int x = 1; x * x <= b1; x++) {
if (b1 % x) continue;
if (x % a1 == 0 && __gcd(x / a1, a0 / a1) == 1 && __gcd(b1 / x, b1 / b0) == 1) cnt++;
int y = b1 / x;
if (x == y) continue;
if (y % a1 == 0 && __gcd(y / a1, a0 / a1) == 1 && __gcd(b1 / y, b1 / b0) == 1) cnt++;
}
cout << cnt << '\n';
}
CF1499D The Number of Pairs
设 a=n\times \gcd(a,b),b=m \times \gcd(a,b)，大力推式子：
\begin{aligned}
c \times \operatorname{lcm(a,b)}-d\times \gcd(a,b)&=x\\
c \times \dfrac{ab}{\gcd(a,b)}-d\times \gcd(a,b)&=x\\
c nm \times \gcd(a,b)-d\times \gcd(a,b)&=x\\
\gcd(a,b)&=\dfrac{x}{cnm-d}
\end{aligned}
说明 (cnm-d) \mid x。O(\sqrt{x}) 地枚举 x 的因子 i。移项有 cnm=d+i，即 c \mid (d+i)。再将 \dfrac{d+i}{c} 分解质因数，对于每种质因子只有全给 n 和全给 m 两种选择，否则会使得 (cnm-d) \nmid x 。用线性筛先得到每个数的质因子个数，则贡献为 2^{cnt}。复杂度 O(n+T\sqrt{n})。
const int N = 2e7 + 5;
bitset<N> isprime;
int cnt[N], prime[N];
void init() {
isprime.set(), isprime[0] = isprime[1] = false;
for (int i = 2; i < N; i++) {
if (isprime[i]) prime[++prime[0]] = i, cnt[i] = 1;
for (int j = 1; j <= prime[0] && 1LL * i * prime[j] < N; j++) {
isprime[i * prime[j]] = false;
if (i % prime[j] == 0) {
cnt[i * prime[j]] = cnt[i];
break;
}
cnt[i * prime[j]] = cnt[i] + 1;
}
}
}
long long c, d, x;
long long f(int x) {
if ((x + d) % c) return 0;
return 1LL << cnt[(x + d) / c];
}
void _main() {
cin >> c >> d >> x;
long long res = 0;
for (long long i = 1; i * i <= x; i++) {
if (x % i) continue;
res += f(i);
if (i * i != x) res += f(x / i);
} cout << res << '\n';
}
以 m_1=0,m_2=1 为例，同加 b 再同除 2，化简得到 2 \gcd(a,b) \mid (x+b)。
剩下两种类似。四种可能或起来即可。
long long a, b, x, y;
void _main() {
cin >> a >> b >> x >> y;
long long g = __gcd(a, b);
if (x % g || y % g) return cout << "N\n", void();
g *= 2;
if (x % g == 0 && y % g == 0) return cout << "Y\n", void();
if ((x + b) % g == 0 && (y + a) % g == 0) return cout << "Y\n", void();
if ((x + a) % g == 0 && (y + b) % g == 0) return cout << "Y\n", void();
if ((x + a + b) % g == 0 && (y + a + b) % g == 0) return cout << "Y\n", void();
cout << "N\n";
}
*P3518 [POI 2011] SEJ-Strongbox
设密码集合为 S，则 \forall i,j \in S, (i+j) \bmod n \in S。考虑对于 i \in S，必然有 \forall k \in \mathbb{N}^+, ki \bmod n \in S，所以首先考虑 ki \bmod n 能取到哪些数。
首先讨论 \gcd(i,n)=1，此时 ki \bmod n 取遍 0,1,2,3,\cdots,n-1。因为同余方程 ki \equiv x \pmod n 必然有解。一般地，猜想 i \in S 使得 ki=x \pmod n 有解，即 \gcd(k,n)=1，反证法得到 ki \bmod n 取到了所有 \gcd(i,n) 的倍数。
对于 i \notin S，i 的所有因子也不是密码。根据这个判定方法，枚举 d \mid \gcd(m_k,n)，若合法，则对于 d \in S 的密码数量就是 \dfrac{n}{d}。复杂度 O(k \sqrt {\gcd(m_k,n)})，可以获得 76pts。
#define int long long
const int N = 2.5e5 + 5;
int n, k, a[N];
bool check(int x) {
for (int i = 1; i < k; i++) {
if (a[i] % x == 0) return false;
} return true;
}
void _main() {
cin >> n >> k;
for (int i = 1; i <= k; i++) cin >> a[i];
int x = __gcd(n, a[k]);
for (int i = 1; i * i <= x; i++) {
if (x % i) continue;
if (check(i)) return cout << n / i, void();
}
for (int i = sqrt(x) + 1; i >= 1; i--) {
if (x % i) continue;
if (check(x / i)) return cout << n / (x / i), void();
}
}
#define int long long
const int N = 2.5e5 + 5;
int n, k, a[N];
bool check(int x) {
for (int i = 1; i < k; i++) {
if (a[i] % x == 0) return false;
} return true;
}
vector<int> p;
void decompose(int x) {
for (int i = 2; i * i <= x; i++) {
if (x % i) continue;
p.emplace_back(i);
while (x % i == 0) x /= i;
}
if (x != 1) p.emplace_back(x);
}
unordered_set<int> f;
void dfs(int x) {
if (f.count(x)) return;
f.emplace(x);
for (int i : p) {
if (x % i == 0) dfs(x / i);
}
}
void _main() {
cin >> n >> k;
for (int i = 1; i <= k; i++) cin >> a[i];
int x = __gcd(n, a[k]);
decompose(x);
for (int i = 1; i < k; i++) dfs(__gcd(a[i], x));
for (int i = 1; i * i <= x; i++) {
if (x % i) continue;
if (!f.count(i)) return cout << n / i, void();
}
for (int i = sqrt(x) + 1; i >= 1; i--) {
if (x % i) continue;
if (!f.count(x / i)) return cout << n / (x / i), void();
}
}
P1082 [NOIP 2012 提高组] 同余方程
这东西我们在之前的乘法逆元讲过方法了，就是用 exGCD 解方程 ax-by=1。输入保证有解，意味着 \gcd(a,b)=1，所以这就是 exGCD 方程的标准形式。细节上，求出 x 以后要处理负数问题。
template <class T> void exgcd(T a, T b, T& x, T& y) {
if (b == 0) return x = 1, y = 0, void();
exgcd(b, a % b, y, x), y -= a / b * x;
}
long long a, b;
void _main() {
cin >> a >> b;
long long x, y; exgcd(a, b, x, y);
cout << (x % b + b) % b;
}
template <class T> void exgcd(T a, T b, T& x, T& y) {
if (b == 0) return x = 1, y = 0, void();
exgcd(b, a % b, y, x), y -= a / b * x;
}
long long x, y, m, n, L;
void _main() {
cin >> x >> y >> m >> n >> L;
long long a = n - m, b = L, c = x - y;
if (a < 0) a = -a, c = -c;
long long d = __gcd(a, b);
if (c % d) return cout << "Impossible", void();
a /= d, b /= d, c /= d;
long long x0, y0; exgcd(a, b, x0, y0);
long long x1 = c * x0;
cout << ((x1 > 0 && x1 % b != 0) ? x1 % b : x1 % b + b);
}
若 p \ge 0，则 a-b=p，即 a+b=2b+p，当 b=0 且 p 为最小非负整数解时，a+b 取得最小值。
若 p<0，则 b-a=-p，即 a+b=2a-p，当 a=0 且 p 为最大非负整数解的相反数时，a-b 取得最小值。
至此，套用模板题的方法，本题解决。
const int N = 2005;
#define int long long
void exgcd(int a, int b, int& x, int &y) {
if (b == 0) return x = 1, y = 0, void();
exgcd(b, a % b, y, x), y -= a / b * x;
}
int w, n, d, x[N];
void _main(int kase) {
cout << "Case #" << kase << ": ";
cin >> w >> n >> d;
for (int i = 1; i <= w; i++) cin >> x[i];
int g = __gcd(n, d), n0 = n / g, res = 0;
for (int i = 1; i <= w / 2; i++) {
if ((x[i] - x[w - i + 1]) % g) return cout << "IMPOSSIBLE\n", void();
int a, b; exgcd(d, n, a, b);
int p = a * ((x[w - i + 1] - x[i]) / g % n0) % n0;
if (p >= 0) p %= n0, res += min(p, n0 - p);
else {
int q = p + n0 * ceil(-1.0 * p / n0);
res += min(q, n0 - q);
}
} cout << res << '\n';
}
P4777 【模板】扩展中国剩余定理（EXCRT）
板子题。给个代码：
const int N = 1e5 + 5;
template <class T> void exgcd(T a, T b, T& x, T& y) {
if (!b) return x = 1, y = 0, void();
exgcd(b, a % b, y, x), y -= a / b * x;
}
template <class T> T crt(int n, const T* a, const T* m) {
T x = 0, y = 0, p = m[1], res = a[1];
for (int i = 2; i <= n; i++) {
T a0 = p, b0 = m[i], c = (a[i] - res % b0 + b0) % b0, g = __gcd(a0, b0);
if (c % g) return -1;
exgcd(a0, b0, x, y);
x = (__int128) x * c / g % b0, res += x * p, p = p / __gcd(p, b0) * b0, res = (res % p + p) % p;
}
return res;
}
int n;
long long a[N], m[N];
void _main() {
cin >> n;
for (int i = 1; i <= n; i++) cin >> m[i] >> a[i];
cout << crt(n, a, m);
}
细节上注意 a_i 可能为负，需要对 b_i 取模处理，另外这题会爆 long long，不过上面板子里已经开了 __int128 了。
const int N = 15;
int n;
long long a[N], b[N];
void _main() {
cin >> n;
for (int i = 1; i <= n; i++) cin >> a[i];
for (int i = 1; i <= n; i++) cin >> b[i];
for (int i = 1; i <= n; i++) a[i] = (a[i] % b[i] + b[i]) % b[i];
cout << crt(n, a, b);
}
const int N = 1e5 + 5;
int n, m;
long long a[N], b[N], p[N], c[N], x, x0[N];
multiset<long long> st;
void _main() {
st.clear();
cin >> n >> m;
for (int i = 1; i <= n; i++) cin >> a[i];
bool f1 = false, f2 = true;
for (int i = 1; i <= n; i++) {
cin >> p[i];
if (a[i] > p[i]) f1 = true;
if (a[i] != p[i]) f2 = false;
}
for (int i = 1; i <= n; i++) cin >> c[i];
for (int i = 1; i <= m; i++) cin >> x, st.emplace(x);
for (int i = 1; i <= n; i++) {
auto it = st.upper_bound(a[i]);
if (it != st.begin()) it--;
b[i] = *it, st.erase(it), st.emplace(c[i]);
debug(b[i]);
}
if (f1) {
long long res = 0;
for (int i = 1; i <= n; i++) res = max(res, (a[i] + b[i] - 1) / b[i]);
return cout << res << '\n', void();
}
if (f2) {
long long res = 1;
for (int i = 1; i <= n; i++) {
long long v = p[i] / __gcd(p[i], b[i]);
res = res / __gcd(res, v) * v;
} return cout << res << '\n', void();
}
for (int i = 1; i <= n; i++) {
if (b[i] % p[i] == 0) {
if (p[i] == a[i]) x0[i] = 0, p[i] = 1;
else return cout << -1 << '\n', void();
continue;
}
long long g = __gcd(b[i], p[i]);
if (a[i] % g) return cout << -1 << '\n', void();
long long y0;
exgcd(b[i], p[i], x0[i], y0);
p[i] /= g, x0[i] = (x0[i] % p[i] + p[i]) % p[i];
x0[i] = (int128) x0[i] * (a[i] / g) % p[i];
}
cout << crt(n, x0, p) << '\n';
}
4. 欧拉函数
4.1 定义
欧拉函数\varphi(x) 指 [1,n] 中与 n 互质的数的个数。
若由算术基本定理可将 n 分解为 p_1^{c_1} p_2^{c_2} p_3^{c_3} \cdots p_m^{c_m}，则
inline int phi(int n) {
if (n == 1) return 1;
int res = n;
for (int i = 2; i * i <= n; i++) {
if (n % i == 0) {
res = res / i * (i - 1);
while (n % i == 0) n /= i;
}
}
if (n > 1) res = res / n * (n - 1);
return res;
}
int n;
void _main() {
cin >> n;
if (n == 1) return cout << 0, 0;
int res = 1;
eular(40000);
for (int i = 1; i < n; i++) res += 2 * phi[i];
cout << res;
}
P2398 GCD SUM
这题做法很多，读者可以翻到下文 8.3 容斥例题学习容斥做法，这里写的是欧拉函数做法。
可以发现 \gcd(i,j) 的值只有 n 种，不妨考虑 gcd(i,j)=k 的数目。对于一对互质的 i,j，有 \gcd(ik,jk)=k，所以可以得到结果为 k 的数目为
const int N = 1e5 + 5;
int n;
long long pre[N];
void _main() {
cin >> n;
eular(n);
for (int i = 1; i <= n; i++) pre[i] = pre[i - 1] + phi[i];
long long res = 0;
for (int i = 1; i <= n; i++) res += 1LL * i * (2 * pre[n / i] - 1);
cout << res;
}
#define int long long
int n;
inline int phi(int n) {
if (n == 1) return 1;
int res = n;
for (int i = 2; i * i <= n; i++) {
if (n % i == 0) {
res = res / i * (i - 1);
while (n % i == 0) n /= i;
}
}
if (n > 1) res = res / n * (n - 1);
return res;
}
inline void _main() {
cin >> n;
long long res = 0;
for (int i = 1; i * i <= n; i++) {
if (n % i) continue;
res += i * phi(n / i);
if (i * i != n) res += n / i * phi(n / (n / i));
}
cout << res;
}
long long p;
long long power(long long a, long long b, long long p) {
long long res = 1; for (a %= p; b; a = a * a % p, b >>= 1) {
if (b & 1) res = res * a % p;
} return res;
}
long long f(long long p) {return p == 1 ? 0 : power(2, f(phi[p]) + phi[p], p);}
void _main() { // 这里已经筛过phi[n]了
cin >> p; cout << f(p) << '\n';
}
long long power(long long a, long long b, long long p) {
long long res = 1; for (a %= p; b; a = (__int128) a * a % p, b >>= 1) {
if (b & 1) res = (__int128) res * a % p;
} return res;
}
vector<long long> factors(long long n) {
vector<long long> res;
for (long long i = 1; i * i <= n; i++) {
if (n % i) continue;
res.emplace_back(i);
if (n / i != i) res.emplace_back(n / i);
}
sort(res.begin(), res.end());
return res;
}
long long n;
void _main() {
for (int kase = 1; ; kase++) {
cin >> n;
if (n == 0) break;
cout << "Case " << kase << ": ";
n = 9 * n / __gcd(n, 8LL);
if (__gcd(n, 10LL) != 1) {cout << "0\n"; continue;}
vector<long long> f = factors(phi(n));
for (long long x : f) {
if (power(10, x, n) == 1) {cout << x << '\n'; break;}
}
}
}
#define int long long
int n, g;
int power(int a, int b, int p) {
int res = 1; for (a %= p; b; b >>= 1) {
if (b & 1) res = res * a % p;
a = a * a % p;
} return res;
}
namespace Lucas {
int fac[N], ifac[N];
inline void init(int p) {
fac[0] = fac[1] = ifac[0] = ifac[1] = 1;
for (int i = 2; i <= p; i++) fac[i] = fac[i - 1] * i % p, ifac[i] = power(fac[i], p - 2, p);
} inline int C(int n, int m, int p) {
if (n < m) return 0;
return fac[n] * ifac[m] % p * ifac[n - m] % p;
} int lucas(int n, int m, int p) {
if (m == 0) return 1;
if (n < p && m < p) return C(n, m, p);
return lucas(n / p, m / p, p) * C(n % p, m % p, p) % p;
}
}
namespace CRT {
template <class T> void exgcd(T a, T b, T& x, T& y) {
if (!b) return x = 1, y = 0, void();
exgcd(b, a % b, y, x), y -= a / b * x;
}
template <class T> T crt(int n, const T* a, const T* m) {
T x = 0, y = 0, p = m[1], res = a[1];
for (int i = 2; i <= n; i++) {
T a0 = p, b0 = m[i], c = (a[i] - res % b0 + b0) % b0, g = __gcd(a0, b0);
if (c % g) return -1;
exgcd(a0, b0, x, y);
x = (__int128) x * c / g % b0, res += x * p, p = p / __gcd(p, b0) * b0, res = (res % p + p) % p;
}
return res;
}
}
const int P[] = {0, 2, 3, 4679, 35617};
int a[5];
void _main() {
cin >> n >> g;
if (__gcd(g, 999911659LL) != 1) return cout << 0, void(); // 注意判无解
for (int i = 1; i <= 4; i++) {
Lucas::init(P[i]);
for (int d = 1; d * d <= n; d++) {
if (n % d) continue;
a[i] = (a[i] + Lucas::lucas(n, d, P[i])) % P[i];
if (n / d != d) a[i] = (a[i] + Lucas::lucas(n, n / d, P[i])) % P[i];
}
}
int val = CRT::crt(4, a, P);
if (val == -1) return cout << 0, void(); // 注意判无解
cout << power(g, val, 999911659);
}
template <class F_t, class G_t>
long long sqrt_decomposition(long long n, long long k, F_t f, G_t g) {
long long res = 0;
for (long long l = 1, r = 0; l <= n; l = r + 1) {
r = (k / l) ? min(n, (k / (k / l))) : n;
res += f(r, l - 1) * g(k / l);
}
return res;
}
int n;
void _main() {
cin >> n;
cout << sqrt_decomposition(n, n,
[](long long x, long long y) {return x - y;}, // f(x) = 1
[](long long x) {return x;} // g(x) = x
) << '\n';
}
P2424 约数和
首先考虑求 f(x) 的前缀和 g(x)。有一个结论是 i 的约数在 n 以内有 \lfloor \dfrac{n}{i} \rfloor 个，所以
g(x) = \sum_{i=1}^{n} i \times \lfloor \dfrac{n}{i} \rfloor
然后套板子去求即可。
long long g(long long n) {
if (n <= 1) return n;
return sqrt_decomposition(n, n,
[](long long x, long long y) {return (x - y) * (x + y + 1) / 2;},
[](long long x) {return x;}
);
}
long long l, r;
void _main() {
cin >> l >> r;
cout << g(r) - g(l - 1);
}
P2261 [CQOI2007] 余数求和
根据取模的定义，有 a \bmod b=a-b \times \lfloor \dfrac{a}{b} \rfloor，然后推一波式子：
\begin{aligned}
G(n, k) &= \sum_{i = 1}^n k \bmod i \\
&= \sum_{i = 1}^n k -i \times \lfloor \dfrac{k}{b} \rfloor \\
&= nk-\sum_{i = 1}^n i \times \lfloor \dfrac{k}{b} \rfloor
\end{aligned}
后面这个东西可以数论分块来做。
long long n, k;
void _main() {
cin >> n >> k;
cout << n * k - sqrt_decomposition(n, k,
[](long long x, long long y) {return (x - y) * (x + y + 1) / 2;}, // f(x) = x
[](long long x) {return x;} // g(x) = x
);
}
long long n;
void _main() {
cin >> n;
long long res = 0;
for (int c = 1; c <= n; c++) {
long long cnt = 0;
if (c % 2 == 0 || c % 5 == 0) continue;
for (int i = 1; i <= n; i <<= 1) {
for (int j = 1; j <= n; j *= 5) {
if (1LL * c * i * j > n) break;
cnt++;
}
} res += cnt * (n / c);
} cout << res;
}
#define int long long
int n;
int cnt(int n) {
int res = 0;
for (int i = 1; i <= n; i <<= 1) {
for (int j = 1; j <= n; j *= 5) {
if (1LL * i * j > n) break;
res++;
}
} return res;
}
int g(int l, int r, int d) {return r / d - (l - 1) / d;}
void _main() {
cin >> n;
int res = 0;
for (int l = 1, r = 0; l <= n; l = r + 1) {
r = min(n, n / (n / l));
res += cnt(n / l) * (n / l) * (r - l + 1 - g(l, r, 2) - g(l, r, 5) + g(l, r, 10));
} cout << res;
}
const int N = 2e5 + 5;
int n, a[N], b[N], p[N], q[N];
long long calc(int len) {return 1LL * len * (len - 1) / 2;}
void _main() {
cin >> n;
for (int i = 1; i <= n; i++) cin >> a[i], p[a[i]] = i;
for (int i = 1; i <= n; i++) cin >> b[i], q[b[i]] = i;
int s = p[1], t = q[1];
if (s > t) swap(s, t);
long long res = 0;
if (1 <= s - 1) res += calc(s);
if (t + 1 <= n) res += calc(n - t + 1);
if (s < t) res += calc(t - s);
int l = s, r = t;
for (int m = 2; m <= n; m++) {
s = p[m], t = q[m];
if (s > t) swap(s, t);
if ((l <= s && s <= r) || (l <= t && t <= r)) {
l = min(l, s), r = max(r, t);
continue;
}
if (s < l && t < l) res += 1LL * (n - r + 1) * (l - t);
if (s > r && t > r) res += 1LL * l * (s - r);
if (s < l && r < t) res += 1LL * (l - s) * (t - r);
l = min(l, s), r = max(r, t);
} cout << res;
}
OpenJudge 9285 盒子与小球之三
推不出来公式，考虑计数 DP。设 dp_{i,j} 表示前 i 个盒子放 j 个球的方案数，由加法原理
const int N = 5005;
int n, m, k;
mint dp[N][N];
void _main() {
cin >> n >> m >> k;
for (int i = 0; i <= m; i++) dp[i][0] = 1;
for (int i = 1; i <= m; i++) {
mint sum = i;
for (int j = 1; j <= n; j++) {
dp[i][j] = sum, sum += dp[i - 1][j + 1];
if (j >= k) sum -= dp[i - 1][j - k];
}
} cout << dp[m][n];
}
P6146 [USACO20FEB] Help Yourself G
将所有线段按左端点排序。设 dp_i 表示前 i 条线段的子集复杂度之和。
考虑转移，计数 DP 非常常见的套路是从插入角度考虑第 i 条的贡献。显然不选的贡献为 dp_{i-1}，只要考虑加入第 i 条的贡献。
由于我们按左端点排好了序，原复杂度单调不降。加入这条线段，会使得某些子集中所有线段与之不交，从而复杂度增加 1。设 [1,i) 中有 x 条线段不与第 i 条相交，选择这 x 条的一个子集会使得复杂度加一。根据加法原理：
const int N = 1e5 + 5;
long long f[N];
int n, m;
void _main() {
cin >> n >> m;
for (int k = n; k >= 1; k--) {
f[k] = 1LL * (n / k) * (m / k);
for (int i = 2; i * k <= min(n, m); i++) f[k] -= f[i * k];
}
long long res = 0;
for (int i = 1; i <= n; i++) res += f[i] * i;
cout << res * 2 - 1LL * n * m;
}
AT_abc162_e [ABC162E] Sum of gcd of Tuples (Hard)
笔者数论训练赛の T5。和上面题的套路一样，我们发现，\gcd a_i=k 时当且仅当 a_1,a_2, \cdots,a_n 均为 k 的倍数，且这个倍数是互质的。考虑容斥 + DP，令 dp_x 表示 [1,x] 内选出 n 个互质数的方法数目。考虑用减法原理，总数目减去不互质的方案，枚举公约数 i，则根据约数性质可得
int n, k;
mint dp[N];
mint solve(int x) {
if (dp[x] != 0) return dp[x];
if (x == 1) return dp[x] = 1;
mint res = mint(x).pow(n);
for (int i = 2; i <= x; i++) {
res -= solve(x / i);
} return dp[x] = res;
}
void _main() {
cin >> n >> k;
mint res = 0;
for (int i = 1; i <= k; i++) res += solve(k / i) * i;
cout << res;
}
|A\cup B \cup C \cup D|=
|A|+|B|+|C|+|D|
-|A\cap B|-|A\cap C|-|A\cap D|-|B\cap C|-|B \cap D|-|C \cap D|
+|A\cap B\cap C|+|A\cap B\cap D|+|A\cap C\cap D|+|B\cap C \cap D|
-|A\cap B\cap C \cap D|
最后再用一个总的减法原理即可。
#define int long long
const int N = 1e5 + 5;
int a, b, c, d, s, c1, c2, c3, c4, q, dp[N];
void work(int x) {
for (int i = x; i < N; i++) dp[i] += dp[i - x];
}
int f(int x) {return x < 0 ? 0 : dp[x];}
void _main() {
cin >> c1 >> c2 >> c3 >> c4;
dp[0] = 1;
work(c1), work(c2), work(c3), work(c4);
for (cin >> q; q--; ) {
cin >> a >> b >> c >> d >> s;
a = (a + 1) * c1, b = (b + 1) * c2, c = (c + 1) * c3, d = (d + 1) * c4;
cout << f(s) - f(s - a) - f(s - b) - f(s - c) - f(s - d)
+ f(s - a - b) + f(s - a - c) + f(s - a - d) + f(s - b - c) + f(s - b - d) + f(s - c - d)
- f(s - a - b - c) - f(s - a - b - d) - f(s - a - c - d) - f(s - b - c - d)
+ f(s - a - b - c - d) << '\n';
}
}
*P6651 「SWTR-5」Chain
一步一步考虑，先想无修改怎么做，我们在拓扑序上做 DP，由加法原理
f_{i} \gets f_{i}+f_j
进一步地，由于 n\le 2\times 10^3，我们可以利用拓扑排序，在 O(n^2) 时间内预处理从 u 到 v 的链数，则
dp_{k,v} \gets dp_{k,v}+dp_{k,u}
类似 Floyd 那样，对于有向边 (u,v)，枚举中继点 k \in [1,n]，由加法原理合并。记入度为 ideg_i，出度为 odeg_i，则总链数
tot=\sum_{ideg_i=0} \sum_{odeg_j=0} dp_{i,j}
因为一条链的起终点一定满足上述条件。接着我们考虑 k=1 的做法，套路地，记 f_i 为链起点到点 i 的方案数，然后建反图，记 g_i 是反图上的 f_i，发现 g_i 等价于原图中点 i 到链终点的方案数。由 dp_{i,j} 的定义易得转移：
若 u,v 可以在同一条链上，那么 u 到终端的路径包含了 v。考虑容斥把重复的加回来，这里的重复是 f_ug_v dp_{u,v}，答案为 tot-f_ug_u-f_vg_v+f_ug_v dp_{u,v}=tot-f_ug_u-g_v(f_v-f_udp_{u,v})。
接着我们扩展到 k \le 15。我们先按拓扑序排序，对于每个点处理要容斥掉的数 h_i，则
h_i=f_i-\sum h_j dp_{j,i}
其中 j 的拓扑序在 i 之前。于是答案为
tot-\sum h_ig_i
至此，我们在预处理 O(n^2+m)，单次查询 O(k^2) 的复杂度下解决了此问题。
const int N = 2e3 + 5, M = 2e4 + 5;
int tot = 0, head[N];
struct Edge {
int next, to;
} edge[M];
inline void add_edge(int u, int v) {
edge[++tot].next = head[u], edge[tot].to = v, head[u] = tot;
}
int n, m, u, v, q, len, ideg[N], odeg[N];
int cnt, bfn[N];
mint dp[N][N], f[N], g[N], h[20];
void topo() {
queue<int> q;
for (int i = 1; i <= n; i++) {
if (ideg[i] == 0) q.emplace(i);
dp[i][i] = 1;
}
while (!q.empty()) {
int u = q.front(); q.pop();
bfn[u] = ++cnt;
for (int j = head[u]; j != 0; j = edge[j].next) {
int v = edge[j].to;
for (int k = 1; k <= n; k++) dp[k][v] += dp[k][u];
if (--ideg[v] == 0) q.emplace(v);
}
}
}
vector<int> st, ed;
void _main() {
cin >> n >> m;
for (int i = 1; i <= m; i++) {
cin >> u >> v;
add_edge(u, v);
odeg[u]++, ideg[v]++;
}
for (int i = 1; i <= n; i++) {
if (ideg[i] == 0) st.emplace_back(i);
if (odeg[i] == 0) ed.emplace_back(i);
} topo();
mint tot = 0;
for (int i : st) {
for (int j : ed) tot += dp[i][j];
}
for (int i = 1; i <= n; i++) {
for (int j : st) f[i] += dp[j][i];
for (int j : ed) g[i] += dp[i][j];
}
for (cin >> q; q--; ) {
fill(h, h + len, 0);
cin >> len;
vector<int> c;
for (int i = 0; i < len; i++) cin >> u, c.emplace_back(u);
sort(c.begin(), c.end(), [](int x, int y) -> bool {
return bfn[x] < bfn[y];
});
for (int i = 0; i < len; i++) {
h[i] = f[c[i]];
for (int j = 0; j < i; j++) h[i] -= dp[c[j]][c[i]] * h[j];
}
mint res = 0;
for (int i = 0; i < len; i++) res += h[i] * g[c[i]];
cout << tot - res << '\n';
}
}
int n, m, cur[N];
void dfs(int x) {
if (x > m) return ..., void(); // 这里cur就是一个组合，进行处理
for (int i = cur[x - 1] + 1; i <= n; i++) cur[x] = i, dfs(x + 1);
}
dfs(1);
它用于枚举 1 到 n 所有自然数选 m 个的组合。如果套上 next_permutation，还能实现排列的枚举。
9.2.2 二项式定理
(a+b)^n=\sum_{i=0}^{n} C_{n}^{i} a^{n-i} b^{i}
这个定理表明，二项式 (a+b)^n 展开项的系数与组合数有直接关系。
我们知道杨辉三角：
其每一个位置的数字通过左上 + 右上确定，每一行所对应的就是二项式展开的系数。
9.2.3 组合数的性质
对称性：
C_{n}^{m} = C_{n}^{n-m}
代数推导易证，组合意义就是把选出的集合取补集。
递推式：
C_{n}^{m}=C_{n-1}^{m}+C_{n-1}^{m-1}
代数推导略去，组合意义类似 dp 的思想可以证明。这个式子实际上就是杨辉三角的递推式。
二项式定理的特殊情况 1：
2^n=\sum_{i=0}^{n} C_n^i
也就是杨辉三角每一行的和。
斐波那契数列：
fib_{n+1}=\sum_{i=0}^{n} C_{n-i}^i
把杨辉三角每条斜 30 \degree 角的线取出来相加可以发现。
范德蒙恒等式：
C_{m+n}^{k}=\sum_{i=0}^{k} C_m^i C_n^{k-i}
假设有两堆物品，每堆分别有 m,n 个物品，总共取 k 个，则方案数可分解为：从第一堆取 i 个物品，第二堆取 k-i 个物品，且两种选择独立，故由乘法得到贡献。最后方案即为求和。
二项式定理的特殊情况 2：
\sum_{i=0}^{n} (-1)^i C_n^i = 0
特殊情况是 n=0 时上式的值为 1。
9.2.4 计算组合数
定义法
优点：写起来简单，不用预处理。
缺点：查询复杂度 O(m)。
根据定义直接计算，注意先乘后除。
inline long long C(int n, int m) {
if (m > n) return 0;
long long res = 1;
for (int i = 1; i <= m; i++) {
res = res * (n - i + 1) % mod * power(i, mod - 2) % mod;
}
return res;
}
递推法
优点：任意模数，写起来简单，查询 O(1)。
缺点：预处理时空复杂度均为 O(n^2)。
即利用组合数性质 2 预处理杨辉三角。
for (int i = 0; i < N; i++) c[i][0] = 1, c[i][i] = 1;
for (int i = 0; i < N; i++) {
for (int j = 1; j < i; j++) c[i][j] = c[i - 1][j] + c[i - 1][j - 1];
}
逆元法
优点：预处理可做到 O(n)，查询 O(1)。
缺点：只能在模数为质数时使用。
通过预处理阶乘和阶乘的逆元，直接用定义式计算组合数。
inline long long power(long long a, long long b) {
long long res = 1;
for (a %= mod; b; b >>= 1) {
if (b & 1) res = res * a % mod;
a = a * a % mod;
}
return res;
}
long long fac[N], ifac[N];
inline void pre() {
fac[0] = fac[1] = ifac[0] = ifac[1] = 1;
for (int i = 2; i < N; i++) fac[i] = fac[i - 1] * i % mod, ifac[i] = power(fac[i], mod - 2);
}
inline long long C(int n, int m) {
if (n < m) return 0;
return fac[n] * ifac[m] % mod * ifac[n - m] % mod;
}
inline long long power(long long a, long long b) {
long long res = 1;
for (a %= mod; b; b >>= 1) {
if (b & 1) res = res * a % mod;
a = a * a % mod;
}
return res;
}
long long fac[N], ifac[N];
inline void pre() {
fac[0] = fac[1] = ifac[0] = ifac[1] = 1;
for (int i = 2; i < N; i++) fac[i] = fac[i - 1] * i % mod, ifac[i] = power(fac[i], mod - 2);
}
inline long long C(int n, int m) {
if (n < m) return 0;
return fac[n] * ifac[m] % mod * ifac[n - m] % mod;
}
long long lucas(long long n, long long m) {
if (m == 0) return 1;
if (n < mod && m < mod) return C(n, m);
return lucas(n / mod, m / mod) * C(n % mod, m % mod) % mod;
}
P3807 【模板】卢卡斯定理/Lucas 定理。
Lucas 定理可以实现快速求组合数前缀和，可以看进阶部分例题中的 P4345 [SHOI2015] 超能粒子炮·改。
exLucas 法
在 10.2 介绍。
9.3 经典模型
这里介绍排列组合题里常用的转化思想，每种模型给出一些例题。
9.3.1 捆绑法
Q: 有 n+m 个不同元素要进行排列，其中 m 个元素必须连续，求方案数。
A: 将 m 个元素捆绑在一起，内部排列方案为 m! 种，这 m 个元素的整体视为一个元素，则外部排列方案为 (n+1)! 种，由乘法原理得方案数为 (n+1)!m!。
9.3.2 插板法
Q1: 现有 n 个相同元素，分为 k 组，每组至少有一个元素，求方案数。
（可以抽象为：求方程 x_1+x_2+\cdots+x_k=n 的正整数解数目）
int n, m;
BigInteger A(int m, int n) {
BigInteger res = 1;
for (int i = n - m + 1; i <= n; i++) res *= i;
return res;
}
void _main() {
cin >> n >> m;
cout << A(n + 2, n + 2) * A(m, n + 3) - A(2, 2) * A(n + 1, n + 1) * A(m, n + 2);
}
const int N = 2e5 + 5;
int n, m;
mint fac[N], ifac[N];
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];}
void _main() {
cin >> n >> m;
fac[0] = ifac[0] = 1;
for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
mint res = 1;
for (int i = 2; i * i <= m; i++) {
if (m % i) continue;
int c = 0;
while (m % i == 0) m /= i, c++;
res *= C(n + c - 1, n - 1);
}
if (m != 1) res *= n;
cout << res;
}
P11250 [GESP202409 八级] 手套配对
首先先从 n 对手套中拿出 k 对，方案数 C_n^k。又因为要恰好拿走 k 对手套，说明剩下的 m-2k 只手套中不能同时选走一对，从其中取出 m 对左或右手套，由乘法原理得答案为 2^{m-2k} \times C_{n-k}^{m-2k}，最后再用乘法原理合并答案。预处理组合数用杨辉三角即可。
const int N = 2005;
int n, m, k;
modint c[N][N];
void _main() {
cin >> n >> m >> k;
if (m < 2 * k) return cout << 0 << '\n', void();
cout << c[n][k] * c[n - k][m - 2 * k] * modint(2).pow(m - 2 * k) << '\n';
} signed main() {
for (int i = 0; i < N; i++) c[i][0] = 1, c[i][i] = 1;
for (int i = 0; i < N; i++) {
for (int j = 1; j < i; j++) c[i][j] = c[i - 1][j] + c[i - 1][j - 1];
}
int t = 1; for (cin >> t; t--; ) _main();
}
P6870 [COCI 2019/2020 #5] Zapina
用 dp_{i,j} 表示前 i 个人分配 j 道题的合法方案数。分类讨论：
让第 i 个人满意，剩下 j-i 道题随便分，由乘法原理得答案为 (i-1)^{j-i} C_i^j。
让第 i 个人不满意，枚举给他分几道题加起来即可。
由加法原理合并答案，复杂度 O(n^3)，注意一下边界。
const int N = 355;
int n;
mint dp[N][N], fac[N], ifac[N];
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];}
void _main() {
cin >> n;
fac[0] = ifac[0] = 1, dp[1][1] = 1;
for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
for (int i = 2; i <= n; i++) {
for (int j = 1; j <= n; j++) {
if (j >= i) dp[i][j] = mint(i - 1).pow(j - i) * C(j, i);
for (int k = 0; k <= j; k++) {
if (i != k) dp[i][j] += dp[i - 1][j - k] * C(j, k);
}
}
} cout << dp[n][n];
}
P1350 车的放置
线性做法。设 f(n,m,k) 表示在 n \times m 的网格中放置 k 个车的方案数。枚举在上方棋盘放置多少个车 i，根据加法原理答案为
\sum_{i=0}^k f(a,b,i) \times f(a+c-i,d,k-i)
考虑算出 f(n,m,k)。考虑 f(n,n,n) 就是将 n 行全排列，方案数为 n!。从 n\times m 的网格中选出一个 k \times k 的子图即可。故
f(n,m,k)=C_n^k C_m^k k!
问题解决。
const int N = 2e3 + 5;
int a, b, c, d, k;
mint fac[N], ifac[N];
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];}
mint f(int n, int m, int k) {return C(n, k) * C(m, k) * fac[k];}
void _main() {
fac[0] = ifac[0] = 1;
for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
cin >> a >> b >> c >> d >> k;
mint res = 0;
for (int i = 0; i <= k; i++) {
res += f(a, b, i) * f(a + c - i, d, k - i);
} cout << res;
}
CF1436C Binary Search
考虑二分查找的过程，由于 n 是固定的，则二分查找的路径是唯一的，因此我们可以求出一定小于等于 x 的数的个数，一定大于 x 的数的个数，分别记作 a,b。那么在小于等于 x 的 x-1 个数中，就要找 a-1 个数出来排列（因为还有一个数等于 x） ，而大于 x 的 n-x 个数中找 b 个来排列，剩下的任意排列，是全排列。
由乘法原理可得答案
A_{x-1}^{a-1} A_{n-x}^{b} A_{n-a-b}^{n-a-b}
计算排列数也可以用逆元法。
const int N = 1005;
int n, x, pos;
mint fac[N], ifac[N];
mint A(int m, int n) {
if (m > n) return 0;
return fac[n] * ifac[n - m];
}
void _main() {
cin >> n >> x >> pos;
int l = 0, r = n, a = 0, b = 0;
while (l < r) {
int mid = (l + r) >> 1;
if (mid <= pos) l = mid + 1, a++;
else r = mid, b++;
} fac[0] = 1, ifac[0] = 1;
for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
cout << A(a - 1, x - 1) * A(b, n - x) * A(n - a - b, n - a - b);
}
const int N = 5e5 + 5;
int n, p[N], q[N];
mint factorial(int x) {
mint res = 1;
for (int i = 2; i <= x; i++) res *= i;
return res;
}
void _main() {
cin >> n;
for (int i = 1; i <= n; i++) cin >> p[i];
for (int i = 1; i <= n; i++) cin >> q[i];
int i = 1, j = 1;
for (int x = 1; x <= n; x++) {
if (p[x] == n) i = x;
if (q[x] == n) j = x;
}
if (i == j) cout << "2 " << factorial(n - 1);
else cout << "3 " << factorial(n - 1) / (n - q[i]) + factorial(n - 1) / (n - p[j]);
}
再考虑不合法的方案。不妨设 c 为最长边，枚举 i \in [0,l] 分配给 c，则 l-i 个数分配给 a,b，且要求 a+b \le c，所以不合法分配为 u=\min(c+i-a-b, l-i)，方案数为 C^{2}_{u+2}=\dfrac{(u+1)(u+2)}{2}，加法原理合并方案。然后对于 a,b 为最长边同理。
#define int long long
int a, b, c, l;
int solve(int a, int b, int c) {
int res = 0;
for (int i = 0; i <= l; i++) {
int u = min(c + i - a - b, l - i);
if (u >= 0) res += (u + 1) * (u + 2) / 2;
} return res;
}
void _main() {
cin >> a >> b >> c >> l;
int res = 1;
for (int i = 1; i <= l; i++) res += (i + 1) * (i + 2) / 2;
res -= solve(a, b, c), res -= solve(a, c, b), res -= solve(b, c, a);
cout << res;
}
CF1929F Sasha and the Wedding Binary Search Tree
BST 题经典套路中序遍历，问题转化为给 -1 赋值使序列单调不降的方案数。进一步地，对于两个已经确定的点 i,j，则对 (i,j) 中的数产生限制 num \in [a_i,a_j]，于是我们只需解决这个问题：
给出一个长为 n 的序列，值域为 [l,r]，求让其单调不降的方案数。
很遗憾的是笔者赛时没有做出来。考虑先把第 i 个数加上 i，则单调不降转为严格递增。而递增序列直接选出 n 个不同的数即可，值域为 [l+1,r+n]，故答案就是 C_{r-l+n}^{n}。或者考虑插板法，n 块板子插入 r-l+n 个空也能得出来。
注意 c 很大，直接逆元法会 RE，但 \sum n \le 10^6，所以分子可以暴力，仍然处理一下逆元。
const int N = 2e6 + 5;
mint fac[N], ifac[N];
mint C(int n, int m) {
if (n < m) return 0;
//return fac[n] * ifac[m] * ifac[n - m];
mint res = 1;
for (int i = n - m + 1; i <= n; i++) res *= i;
return res * ifac[m];
} mint solve(int n, int l, int r) {
if (n <= 0 || l > r) return 1;
return C(r - l + n, n);
}
int n, c, a[N], ls[N], rs[N], v[N];
void dfs(int rt) {
if (rt == -1) return;
dfs(ls[rt]), a[++a[0]] = v[rt], dfs(rs[rt]);
}
void _main() {
cin >> n >> c;
for (int i = 1; i <= n; i++) cin >> ls[i] >> rs[i] >> v[i];
a[0] = 0, dfs(1);
mint res = 1;
int l = 0;
for (int i = 1; i <= n; i++) {
if (a[i] == -1) continue;
res *= solve(i - l - 1, l ? a[l] : 1, a[i]);
l = i;
}
res *= solve(n - l, l ? a[l] : 1, c);
cout << res << '\n';
}
*P5689 [CSP-S2019 江西] 多叉堆
看到题目中的操作形式，想到带权并查集。我们设当前操作是把 u 接到 v 上。那么 v 的位置只能填 0。
维护 a_x 表示答案，s_x 表示树的大小。填入 u 子树中的数字就有 C_{s_v-1}^{s_u} 种选择方案，剩下的数字填入 v 子树中。
const int N = 3e5 + 5;
mint fac[N], ifac[N];
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];}
int n, q, opt, x, y, last;
int fa[N], sz[N];
mint w[N];
int find(int x) {return x == fa[x] ? x : fa[x] = find(fa[x]);}
void unite(int x, int y) {
x = find(x), y = find(y);
if (x == y) return;
fa[x] = y, sz[y] += sz[x], w[y] *= w[x] * C(sz[y] - 1, sz[x]);
}
void _main() {
cin >> n >> q;
fac[0] = ifac[0] = 1;
for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
iota(fa + 1, fa + n + 1, 1), fill(sz + 1, sz + n + 1, 1), fill(w + 1, w + n + 1, 1);
while (q--) {
cin >> opt >> x;
if (opt == 1) {
cin >> y;
x = (x + last) % n + 1, y = (y + last) % n + 1;
unite(x, y);
} else {
x = (x + last) % n + 1;
cout << (last = w[find(x)].val) << '\n';
}
}
}
const int N = 4e5 + 5;
int n, a[N];
mint fac[N], ifac[N];
mint C(int n, int m) {
if (n < m) return 0;
return fac[n] * ifac[m] * ifac[n - m];
}
void _main() {
cin >> n;
for (int i = 0; i <= n; i++) cin >> a[i];
fac[0] = ifac[0] = 1;
for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
mint res = 0;
for (int i = 0; i <= n; i++) res += C(i + a[i], i + 1);
cout << res;
}
*P5505 [JSOI2011] 分特产
正难则反，记 f(x) 为至少 x 人没有分到的方案数，根据容斥原理答案为
\sum_{i=0}^{n-1} (-1)^i f(i)
需要求出 f(x)。考虑插板法，第 k 个特产没有分到的方案是 C_{a_k+n-x-1}^{n-x-1}，钦定 x 人没有分到的方案数为 C_n^x，根据乘法原理合并：
f(x)=C_n^x \prod_{k=1}^m C_{a_k+n-x-1}^{n-x-1}
逆元法预处理组合数即可，复杂度 O(nm)。
const int N = 2005;
mint fac[N], ifac[N];
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m]; }
int n, m, a[N];
mint f(int x) {
mint res = 1;
for (int k = 1; k <= m; k++) res *= C(a[k] + n - x - 1, n - x - 1);
return res * C(n, x);
}
void _main() {
fac[0] = ifac[0] = 1;
for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
cin >> n >> m;
for (int i = 1; i <= m; i++) cin >> a[i];
mint res = 0;
for (int i = 0; i < n; i++) {
if (i & 1) res -= f(i);
else res += f(i);
} cout << res;
}
10. 排列组合进阶
由于把这坨东西放到上面会导致目录阅读体验较差，所以把这些较难的内容单独开出来了。
10.1 排列进阶
*10.1.1 康托展开
康托展开解决的是求排列字典序的问题。先给出公式：
1+\sum_{i=1}^{n} rk_i \times (n-i)!
其中 rk_i 表示 [i,n] 中有多少个数小于 a_i，即 a_i 的后缀排名。理解一下这个公式，字典序就是比当前排列小的排列数目，贪心地枚举 i，使得第 i 位小于 a_i，而后置位与字典序无关。那么我们就需要把 i 后面小于 a_i 的数换到前面，且根据全排列的原理有 (n-i)! 种方案。
long long power(long long a, long long b, long long p) {
long long res = 1; for (a %= p; b; b >>= 1) {
if (b & 1) res = res * a % p;
a = a * a % p;
} return res;
}
void exgcd(long long a, long long b, long long& x, long long& y) {
b == 0 ? (x = 1, y = 0) : (exgcd(b, a % b, y, x), y -= a / b * x);
}
long long inverse(long long x, long long p) {
long long a, b;
return exgcd(x, p, a, b), (a % p + p) % p;
}
long long calc(long long n, long long x, long long p) {
if (n == 0) return 1;
long long s = 1;
for (long long i = 1; i <= p; i++) {
if (i % x) s = s * i % p;
} s = power(s, n / p, p);
for (long long i = n / p * p + 1; i <= n; i++) {
if (i % x) s = i % p * s % p;
} return s * calc(n / x, x, p) % p;
}
long long multilucas(long long m, long long n, long long x, long long p) {
int cnt = 0;
for (long long i = m; i; i /= x) cnt += i / x;
for (long long i = n; i; i /= x) cnt -= i / x;
for (long long i = m - n; i; i /= x) cnt -= i / x;
return power(x, cnt, p) % p * calc(m, x, p) % p * inverse(calc(n, x, p), p) % p
* inverse(calc(m - n, x, p), p) % p;
}
long long x[20];
int crt(int n, long long* a, long long* m) {
long long mod = 1, res = 0;
for (int i = 1; i <= n; i++) mod *= m[i];
for (int i = 1; i <= n; i++) {
x[i] = mod / m[i];
long long x0 = -1, y0 = -1;
exgcd(x[i], m[i], x0, y0);
res += a[i] * x[i] * (x0 >= 0 ? x0 : x0 + m[i]);
} return res % mod;
}
long long exlucas(long long m, long long n, long long p) {
int len = 0;
long long p0[20], a0[20];
for (long long i = 2; i * i <= p; i++) {
if (p % i) continue;
p0[++len] = 1;
while (p % i == 0) p0[len] *= i, p /= i;
a0[len] = multilucas(m, n, i, p0[len]);
}
if (p > 1) p0[++len] = p, a0[len] = multilucas(m, n, p, p);
return crt(len, a0, p0);
}
P4720 【模板】扩展卢卡斯定理/exLucas。
*10.3 二项式反演
如果某些计数问题的限制是“某些物品恰好若干个”，就可以考虑使用二项式反演。
10.3.1 原理
设 f_n 表示恰好使用 n 个不同元素形成特定结构的方案数，g_n 表示从这 n 个不同元素中选出若干个元素形成特定结构的方案数。
const int N = 4e5 + 5;
long long n, m, k;
mint fac[N], ifac[N];
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];}
mint f(int x, int y) {
if (x < 0 || y < 0) return 0;
return C(x + y + 1, x) * (y + 1) - C(x + y + 1, x + 1) * 2;
}
void _main() {
cin >> n >> m >> k;
if (k < n + m - 2) cout << 0 << '\n';
else if (k == n + m - 2) cout << C(n + m - 2, m - 1) << '\n';
else if (k == n + m - 1) {
mint x = n * (m - 1) + m * (n - 1) - (n + m - 2);
cout << x * C(n + m - 2, m - 1) << '\n';
} else {
mint y = 2 * n * m - n - m - k + 1;
cout << C(n + m - 2, m - 1) * y * (y + 1) / 2
- C(n + m - 4, n - 2) * (n + m - 3)
+ f(n, m - 3) + f(m, n - 3) << '\n';
}
}
其中 p 表示分数不变的方案数。现在我们只要算出 p 即可。设有 m 个位置使得 h_i 与 h_{i\bmod n+1} 不同，分数不变就是说移动时对了 i 个，错了 i 个，则 2i \le m，即 i \in [0,\lfloor \dfrac{m}{2} \rfloor]。枚举 i，考察 i 产生的贡献，有三部分：
对了 i 个，从 m 个位置选出 i 个，方案数 C^i_m；
错了 i 个，从剩下 m-i 个位置再选 i 个，然后考虑不同选项，一共是错 m-2i 个，且去掉正确的选项与自身，总方案是 C^{i}_{m-i} (k-2)^{m-2i}。
const int N = 2e5 + 5;
int n, k, a[N];
mint fac[N], ifac[N];
mint C(int n, int m) {
if (n < m) return 0;
return fac[n] * ifac[m] * ifac[n - m];
}
void _main() {
cin >> n >> k;
if (k == 1) return cout << 0, void();
fac[0] = ifac[0] = 1;
for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
int m = 0;
for (int i = 1; i <= n; i++) cin >> a[i];
for (int i = 1; i <= n; i++) {
if (a[i] != a[i % n + 1]) m++;
} mint res = 0;
for (int i = 0; i <= m / 2; i++) res += C(m, i) * C(m - i, i) * mint(k - 2).pow(m - 2 * i) * mint(k).pow(n - m);
cout << (mint(k).pow(n) - res) / 2;
}
P5367 【模板】康托展开
康托展开板子题，但你要是直接按公式求是 O(n^2) 的。显然枚举 i 不能删，瓶颈在计算 rk_i，使用权值树状数组维护即可。
const int N = 1e6 + 5;
int n, a[N];
int tr[N];
void add(int x, int c) {for (; x <= n; x += (x & -x)) tr[x] += c;}
int ask(int x) {
int res = 0;
for (; x != 0; x -= (x & -x)) res += tr[x];
return res;
}
mint fac[N];
void _main() {
cin >> n;
for (int i = 1; i <= n; i++) cin >> a[i];
fac[0] = 1;
for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i, add(i, 1);
mint res = 1; // 初值为1
for (int i = 1; i <= n; i++) res += fac[n - i] * (ask(a[i]) - 1), add(a[i], -1);
cout << res;
}
int n, q, a[N], ans[N];
long long x, fac[N], w[N];
char opt;
long long tr[N];
void add(int x, long long c) {for (; x <= n; x += (x & -x)) tr[x] += c;}
long long ask(int x) {
long long res = 0;
for (; x != 0; x -= (x & -x)) res += tr[x];
return res;
}
int kth(int k) {
long long pos = 0, x = 0;
for (int i = __lg(n); i >= 0; i--) {
x += (1 << i);
if (x >= n || pos + tr[x] >= k) x -= (1 << i);
else pos += tr[x];
} return x + 1;
}
void _main() {
cin >> n >> q;
fac[0] = 1;
for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i;
while (q--) {
memset(tr, 0, sizeof(tr));
cin >> opt;
if (opt == 'Q') {
for (int i = 1; i <= n; i++) cin >> a[i], add(i, 1);
long long res = 1;
for (int i = 1; i <= n; i++) res += fac[n - i] * (ask(a[i]) - 1), add(a[i], -1);
cout << res << '\n';
} else if (opt == 'P') {
cin >> x; x--;
for (int i = 1; i <= n; i++) w[n - i + 1] = x % i, x /= i, add(i, 1);
for (int i = 1; i <= n; i++) ans[i] = kth(w[i] + 1), add(ans[i], -1);
for (int i = 1; i <= n; i++) cout << ans[i] << ' ';
cout << '\n';
}
}
}
在这个问题中，对于位置 i 有且仅有一个 k 使得 p_i 右移 k 位以后满足 p_i=i。我们预处理出每个 k 能使得多少个 p_i 归位，根据抽屉原理，合法的 k 的数目不超过
\lfloor \dfrac{n}{n-2m} \rfloor
个。在本题中，这个数值不超过 3。加上这个剪枝即可通过。
const int N = 3e5 + 5;
int n, m, a[N], pos[N], cnt[N], fa[N];
int find(int x) {return x == fa[x] ? x : fa[x] = find(fa[x]);}
void _main() {
memset(cnt, 0, sizeof(int) * (n + 1));
cin >> n >> m;
for (int i = 1; i <= n; i++) cin >> a[i], pos[a[i]] = i, cnt[(i - a[i] + n) % n]++;
vector<int> res;
for (int k = 0; k < n; k++) {
if (cnt[k] < n - 2 * m) continue;
iota(fa + 1, fa + n + 1, 1);
int c = 0;
for (int i = 1; i <= n; i++) {
int x = pos[i], y = (i + k - 1) % n + 1;
x = find(x), y = find(y);
if (x != y) fa[x] = y, c++;
}
if (c <= m) res.emplace_back(k);
}
cout << res.size() << ' ';
for (int i : res) cout << i << ' ';
cout << '\n';
}
*AT_arc124_d [ARC124D] Yet Another Sorting Problem
看到排列对换问题，要想到图论建模，也就是把当前位置和目标位置连边。对于本题，我们把 i 和 p_i 之间连一条边。
设前 n 个点为红点，后 m 个点为白点。对 (x,y) 进行一次对换就是交换 x,y 的出点。以第二个样例为例：
如果交换 (2, 9)，那么新图为
如果交换 (5,9)，则新图为
我们可以发现两条结论：
如果交换两个不在同一连通块内的点，则等价于在一个环上插入一条链。
如果交换两个在同一连通块内的点，则原环在两个位置断开分裂为两个独立环。
最终的目标是做出 n 个自环。我们可以发现：对于一个大小为 n 的异色环，因为每次操作需要保留一对异色点，总次数为 n-1。
对于一个大小为 n 的同色环，需要借助外部力量完成交换，最小次数是 n+1。有一种做法是对于两个异色的同色环一起操作，先将两个环合并为一个大的异色环再消除，总次数为二者大小之和。于是我们先做完异色环，然后按大小排序贪心即可。复杂度 O(n \log n)。
const int N = 2e5 + 5;
int n, m, p[N], fa[N];
int find(int x) {return x == fa[x] ? x : fa[x] = find(fa[x]);}
vector<int> r[N], a, b;
void _main() {
cin >> n >> m;
iota(fa + 1, fa + n + m + 1, 1);
for (int i = 1; i <= n + m; i++) {
cin >> p[i];
fa[find(i)] = find(p[i]);
}
for (int i = 1; i <= n + m; i++) r[find(i)].emplace_back(i);
int res = 0;
for (int i = 1; i <= n + m; i++) {
if (find(i) != i) continue;
if (r[i].size() == 1) continue; // 自环
if (r[i].back() <= n) a.emplace_back(r[i].size()); // 红色环
else if (r[i][0] > n) b.emplace_back(r[i].size()); // 白色环
else res += r[i].size() - 1;
}
sort(a.begin(), a.end(), greater<int>()), sort(b.begin(), b.end(), greater<int>());
int la = a.size(), lb = b.size();
for (int i = 0; i < min(la, lb); i++) res += a[i] + b[i];
for (int i = min(la, lb); i < max(la, lb); i++) res += (i < la ? a[i] : b[i]) + 1;
cout << res;
}
*P4778 Counting swaps
排列对换问题，先把 i \to a_i 连边，最后的目标状态是 n 个自环。
容易证明将一个大小为 n 的环变成 n 个自环至少需要 n- 1 次对换。对于一个长度为 n 的环，达成最少次数的方法是将它拆成两个大小为 x,y 的环，满足 x+y=n。令 T(x,y) 表示分裂环的方案数，不难得到当且仅当 x=y 时 T(x,y)=x，否则 T(x,y)=x+y.
考虑一个 DP，令 f_n 表示大小为 n 的环以最优步数变成自环的方案数，根据多重集的排列数得
f_n=\sum_{x+y=n, x \le y} T(x,y) f_x f_y \dfrac{(n-2)!}{(x-1)!(y-1)!}
const int N = 1e5 + 5;
int n, a[N], fa[N], cnt[N];
mint f[N], fac[N];
int find(int x) {return x == fa[x] ? x : fa[x] = find(fa[x]);}
mint T(int x, int y) {return x == y ? x : x + y;}
void init() {
fac[0] = 1;
for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i;
f[1] = 1;
for (int i = 2; i < 1000; i++) {
for (int j = 1; j <= i / 2; j++) {
f[i] += T(j, i - j) * f[j] * f[i - j] * fac[i - 2] / fac[j - 1] / fac[i - j - 1];
}
}
}
void _main() {
cin >> n;
iota(fa + 1, fa + n + 1, 1), fill(cnt + 1, cnt + n + 1, 0);
for (int i = 1; i <= n; i++) cin >> a[i], fa[find(i)] = find(a[i]);
for (int i = 1; i <= n; i++) cnt[find(i)]++;
mint res = 1;
int k = 0;
for (int i = 1; i <= n; i++) {
if (find(i) != i) continue;
res *= f[cnt[i]] / fac[cnt[i] - 1], k++;
}
cout << res * fac[n - k] << '\n';
}
const int N = 2e5 + 5;
int n, a[N], v[N];
mint f[N], fac[N], ifac[N];
mint C(int n, int m) {
if (n < m) return 0;
return fac[n] * ifac[m] * ifac[n - m];
}
void _main() {
cin >> n;
for (int i = 1; i <= n; i++) cin >> a[i];
for (int i = 1; i <= n; i++) cin >> v[i];
fac[0] = ifac[0] = 1;
for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
mint sum = 1, res = 0;
for (int i = 1; i <= n; i++) f[i] = mint(2).pow(a[i]), sum *= f[i];
for (int i = 1; i <= n; i++) {
mint cur = 0;
res += (mint(v[i] + 1).pow(a[i]) - 1) * (sum / f[i]);
} cout << res;
}
long long n, m, l, r;
void _main() {
cin >> n >> m >> l >> r;
if (n * m % 2) return cout << mint(r - l + 1).pow(n * m), void();
long long a = r / 2 - (l - 1) / 2, b = r - l + 1 - a;
cout << (mint(a + b).pow(n * m) + mint(a - b).pow(n * m)) / 2;
}
#define popcount __builtin_popcount
const int N = 25;
int n, m;
long long y, a[N], g[N], c[N][N];
void _main() {
cin >> n >> m >> y;
for (int i = 1; i <= n; i++) cin >> a[i];
for (int i = 0; i <= n; i++) c[i][0] = 1, c[i][i] = 1;
for (int i = 0; i <= n; i++) {
for (int j = 1; j < i; j++) c[i][j] = c[i - 1][j] + c[i - 1][j - 1];
}
for (int s = 0; s < (1 << n); s++) {
__int128 x = 1;
for (int i = 1; i <= n; i++) {
if (s >> (i - 1) & 1) x = x / __gcd<__int128>(x, a[i]) * a[i];
if (x > y) break; // 注意这句
} g[popcount(s)] += y / x;
}
long long res = 0;
for (int i = m; i <= n; i++) {
if ((m - i) & 1) res -= c[i][m] * g[i];
else res += c[i][m] * g[i];
} cout << res;
}
P10596 BZOJ2839 集合计数
看到不好求的“恰好”，考虑二项式反演。设 f_i 表示交集大小恰好为 i 的方案数，g_i 表示钦定 i 个元素属于交集的方案数。
那么对于 g_i，先从 n 个元素中选出 i 个，方案数 C_n^i。剩下 n-i 个元素可选可不选，就是求大小为 2^{n-i} 的集合的非空子集数，答案为 2^{2^{n-i}}-1，乘法原理：
g_i=C_n^i (2^{2^{n-i}}-1)
上二项式反演：
f_k=\sum_{i=k}^{n}(-1)^{i-k} C_i^k g_i
复杂度 O(n \log n)。注意求 g_i 要用欧拉定理。
const int N = 1e6 + 5;
mint fac[N], ifac[N], g[N];
mint C(int n, int m) {
return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];
}
int n, k;
void _main() {
cin >> n >> k;
fac[0] = ifac[0] = 1;
for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
for (int i = 0; i <= n; i++) g[i] = C(n, i) * (mint(2).pow(mod32<1000000006>(2).pow(n - i).val) - 1);
mint res = 0;
for (int i = k; i <= n; i++) {
if ((i - k) & 1) res -= C(i, k) * g[i];
else res += C(i, k) * g[i];
} cout << res;
}
*P4859 已经没有什么好害怕的了
形式化题意：给出长为 n 的序列 a,b，将其两两配对使得 a_i>b_i 的数目减去 a_i<b_i 的数目恰为 k，求方案数。
const int N = 2005;
int n, m, k, a[N], b[N], last[N];
mint fac[N], ifac[N], dp[N][N], g[N];
mint C(int n, int m) {
if (n < m) return 0;
return fac[n] * ifac[m] * ifac[n - m];
}
void _main() {
cin >> n >> k;
if ((n + k) & 1) return cout << 0, void();
fac[0] = ifac[0] = 1;
for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
m = (n + k) / 2;
for (int i = 1; i <= n; i++) cin >> a[i];
for (int i = 1; i <= n; i++) cin >> b[i];
sort(a + 1, a + n + 1), sort(b + 1, b + n + 1);
for (int i = 1; i <= n; i++) last[i] = lower_bound(b + 1, b + n + 1, a[i]) - b - 1;
dp[0][0] = 1;
for (int i = 1; i <= n; i++) {
dp[i][0] = dp[i - 1][0];
for (int j = 1; j <= i; j++) {
dp[i][j] = dp[i - 1][j] + dp[i - 1][j - 1] * max(0, last[i] - j + 1);
}
}
for (int i = 0; i <= n; i++) g[i] = dp[n][i] * fac[n - i];
mint res = 0;
for (int i = m; i <= n; i++) {
if ((i - m) & 1) res -= C(i, m) * g[i];
else res += C(i, m) * g[i];
} cout << res;
}
*P6478 [NOI Online #2 提高组] 游戏
“恰好 k 次非平局”很难处理，不妨转为“钦定 k 次非平局”。设恰好 k 次非平局的方案数为 g(k)，钦定 k 次非平局的方案数为 f(k)，显然有 f(n)=\sum_{i=n}^m C_i^n g(i)，根据二项式反演：
g(n)=\sum_{i=n}^m (-1)^{i-n} C_i^n f(i)
考虑求 f(i) 即可。考虑树形 DP，设 dp_{u,i} 表示在以 u 为根的子树中钦定 i 个点且必须有胜负的方案数。
const int N = 2400, p = 2333;
long long q, n, k;
mint c[N][N], pre[N][N];
mint lucas(long long n, long long m) {
if (n < m) return 0;
if (n == m) return 1;
return c[n % p][m % p] * lucas(n / p, m / p);
}
mint f(long long n, long long k) {
if (k < 0) return 0;
if (n == 0 || k == 0) return 1;
if (n < p && k < p) return pre[n][k];
return f(n / p, k / p - 1) * pre[n % p][p - 1] + pre[n % p][k % p] * lucas(n / p, k / p);
}
void _main() {
for (int i = 0; i < N; i++) c[i][0] = 1, c[i][i] = 1;
for (int i = 0; i < N; i++) {
for (int j = 1; j < i; j++) c[i][j] = c[i - 1][j] + c[i - 1][j - 1];
}
for (int i = 0; i < N; i++) {
pre[i][0] = 1;
for (int j = 1; j < N; j++) pre[i][j] = pre[i][j - 1] + c[i][j];
}
for (cin >> q; q--; ) cin >> n >> k, cout << f(n, k) << '\n';
}
*CF785D Anton and School - 2
考虑枚举子序列的最后一个左括号，可以动态统计其左侧及自身的左括号数目 a，右括号数目 b，然后枚举再选 i 个括号，则其贡献为
const int N = 2e5 + 5;
int n, a[N], b[N];
char s[N];
mint fac[N], ifac[N];
mint C(int n, int m) {
if (n < m) return 0;
return fac[n] * ifac[m] * ifac[n - m];
}
void _main() {
cin >> (s + 1); n = strlen(s + 1);
fac[0] = ifac[0] = 1;
for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
for (int i = 1; i <= n; i++) a[i] = a[i - 1] + (s[i] == '(');
for (int i = n; i >= 1; i--) b[i] = b[i + 1] + (s[i] == ')');
mint res = 0;
for (int i = 1; i <= n; i++) {
if (s[i] == '(') res += C(a[i] + b[i] - 1, a[i]);
} cout << res;
}
CF451E Devu and Flowers
【模板】多重集的组合数。将一个盒子视作一组重复元素，根据公式计算即可。
实现上，通过二进制枚举集合 s，设 s 在二进制下第 i_1,i_2,\cdots,i_k 位为 1，则贡献为
(-1)^x C_{n+m-\sum _{j=1}^k a_{i_j}-(k+1)}^{n-1}
还有一个问题是 n \le 20 但是 m \le 10^{14}。需要把逆元法和定义法结合起来算组合数。复杂度 O(n2^n)。
#define int long long
const int N = 22;
int n, m, a[N];
mint inv[N];
mint C(int n, int m) {
if (m < 0 || n < 0 || n < m) return 0;
if (m == 0 || mint(n) == 0) return 1;
mint res = 1;
for (int i = 0; i < m; i++) res *= n - i;
for (int i = 1; i <= m; i++) res *= inv[i];
return res;
}
void _main() {
cin >> n >> m;
for (int i = 1; i <= n; i++) inv[i] = ~mint(i);
for (int i = 1; i <= n; i++) cin >> a[i];
mint res = C(n + m - 1, n - 1);
for (int s = 1; s < (1 << n); s++) {
int sum = n + m, k = 0;
for (int i = 1; i <= n; i++) {
if (s >> (i - 1) & 1) sum -= a[i], k++;
}
sum -= k + 1;
if (k & 1) res -= C(sum, n - 1);
else res += C(sum, n - 1);
} cout << res;
}
11. 错位排列
错排数D_i 是指将 1 到 n 的自然数排列重新排列为 P_i 后，满足 \forall i \in [1,n], P_i \ne i 的方案数。
const int N = 1e6 + 5;
const long long mod = 1e9 + 7;
long long d[N], fac[N], ifac[N];
long long power(long long a, long long b) {
long long res = 1;
for (a %= mod; b; b >>= 1) {
if (b & 1) res = res * a % mod;
a = a * a % mod;
}
return res;
}
long long t, n, m;
void _main() {
d[2] = 1;
for (int i = 3; i < N; i++) d[i] = 1LL * (i - 1) * ((d[i - 1] + d[i - 2]) % mod) % mod;
fac[0] = 1;
for (int i = 1; i < N; i++) {
fac[i] = fac[i - 1] * i % mod;
ifac[i] = power(fac[i], mod - 2);
}
for (cin >> t; t--; ) {
cin >> n >> m;
if (m == 0) {cout << d[n] << '\n'; continue;}
if (n == m) {cout << 1 << '\n'; continue;}
if (n == m + 1) {cout << 0 << '\n'; continue;}
cout << (fac[n] * ifac[m] % mod * ifac[n - m] % mod) * (n >= m ? d[n - m] : 1) % mod << '\n';
}
}
同样的思想可以解决 P8788 的问题 A。
*P4921 [MtOI2018] 情侣？给我烧了！
定义一个“广义错排数” f(x) 表示 x 对情侣都错开的方案数。方法和上面那题一样，区别在于情侣位置有序，且一对情侣的相对位置有两种选择，由乘法原理得
const int N = 1005;
int q, n;
mint fac[N], ifac[N], f[N];
mint A(int n, int m) {
if (n < m) return 0;
return fac[n] * ifac[n - m];
}
mint C(int n, int m) {
if (n < m) return 0;
return fac[n] * ifac[m] * ifac[n - m];
}
void _main() {
fac[0] = ifac[0] = 1, f[0] = 1, f[1] = 0;
for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
for (int i = 2; i < N; i++) f[i] = (f[i - 1] + f[i - 2] * 2 * (i - 1)) * 2 * i * (2 * i - 2);
for (cin >> q; q--; ) {
cin >> n;
for (int k = 0; k <= n; k++) cout << C(n, k) * A(n, k) * mint(2).pow(k) * f[n - k] << '\n';
}
}
const int N = 2e6 + 5;
int n, ls[N], rs[N], sz[N];
mint h[N];
void dfs1(int u) {
if (!u) return;
dfs1(ls[u]), dfs1(rs[u]), sz[u] = sz[ls[u]] + sz[rs[u]] + 1;
} mint dfs2(int u, mint x) {
if (!u) return 0;
mint res = 0;
for (int i = 0; i < sz[ls[u]]; i++) res += x * h[i] * h[sz[u] - i - 1];
res += dfs2(ls[u], x * h[sz[rs[u]]]);
res += dfs2(rs[u], x);
return res;
}
void _main() {
cin >> n;
for (int i = 1; i <= n; i++) cin >> ls[i] >> rs[i];
h[0] = 1;
for (int i = 1; i <= (n << 1); i++) h[i] = h[i - 1] * (4 * i - 2) / (i + 1);
mint res = 0;
for (int i = 1; i < n; i++) res += h[i];
dfs1(1), res += dfs2(1, 1);
cout << res;
}
mint dfs2(int u, mint x) {
if (!u) return 0;
mint res = 0;
if (sz[ls[u]] <= sz[rs[u]]) {
for (int i = 0; i < sz[ls[u]]; i++) res += x * h[i] * h[sz[u] - i - 1];
} else {
res += x * h[sz[u]];
for (int i = 0; i <= sz[rs[u]]; i++) res -= x * h[i] * h[sz[u] - i - 1];
}
res += dfs2(ls[u], x * h[sz[rs[u]]]);
res += dfs2(rs[u], x);
return res;
}
13. Stirling 数
13.1 第二类 Stirling 数
第二类 Stirling 数S(n,k) 表示将 n 个不同元素分为 k 个互不区分的非空子集的方案数。
因为 \operatorname{lcm}_{i=1}^{k}=\dfrac{\prod_{i=1}^{k} a_i}{\gcd_{i=1}^{k} a_i}，所以 \gcd_{i=1}^{k} a_i=1，而且又说了 a 单调不降，就是说 a 序列两两互质。对 n 分解质因数得 n=p_1^{a_1} p_2^{a_2} \cdots，发现若 a_x \ne 0，则这些 p_x 只能被分配到同一个 a_i 中。记有 c 个不同质因子，然后你发现这是把 c 个元素划成 k 个子集的方案数，即第二类 Stirling 数。
实现细节上，注意特判掉 c<k 的情况。
int n, k;
int decompose(int n) {
int m = 0;
for (int i = 2; i * i <= n; i++) {
if (n % i) continue;
m++;
while (n % i == 0) n /= i;
}
if (n > 1) m++;
return m;
}
long long s[15][15];
void _main() {
cin >> n >> k;
int m = decompose(n);
cout << (m < k ? 0 : s[m][k]) << '\n';
} signed main() {
s[0][0] = 1;
for (int i = 1; i < 15; i++) {
for (int j = 1; j <= i; j++) s[i][j] = s[i - 1][j - 1] + s[i - 1][j] * j;
int t = 1; for (cin >> t; t--; ) _main();
}
*P4609 [FJOI2016] 建筑师
单调栈经典模型。 首先高度为 n 的建筑肯定不会被挡，用它将建筑划分为两段，左边的看不到右边，右边的看不到左边。
推广一下，把 n 的建筑分成 a+b-1 个部分，把一个可以看到的和被它的分到一组去，一共是 a+b-2 组。我们发现，每一组除了最高的建筑都可以任意排列，而且这还是一个环形排列，所以方案数是第一类 Stirling 数 s(n-1,a+b-2)。
然后再思考每一组的放置方法，这是一个组合数 C_{a+b-2}^{a-1}。由乘法原理合并答案。
int q, n, a, b;
const int N = 205;
mint c[N][N], s[50005][N];
void _main() {
s[0][0] = 1;
for (int i = 1; i < 50005; i++) {
for (int j = 1; j < N; j++) s[i][j] = s[i - 1][j - 1] + s[i - 1][j] * (i - 1);
}
for (int i = 0; i < N; i++) c[i][0] = 1, c[i][i] = 1;
for (int i = 0; i < N; i++) {
for (int j = 1; j < i; j++) c[i][j] = c[i - 1][j] + c[i - 1][j - 1];
}
for (cin >> q; q--; ) {
cin >> n >> a >> b;
cout << s[n - 1][a + b - 2] * c[a + b - 2][a - 1] << '\n';
}
}
const int N = 1e6 + 5;
int n, k;
mint fac[N], ifac[N];
mint S(int n, int k) {
mint res = 0;
for (int i = 0; i <= k; i++) {
if ((k - i) & 1) res -= mint(i).pow(n) * ifac[i] * ifac[k - i];
else res += mint(i).pow(n) * ifac[i] * ifac[k - i];
} return res;
}
void _main() {
cin >> n >> k;
fac[0] = ifac[0] = 1;
for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
cout << S(n, n - k);
}
*P10591 BZOJ4671 异或图
好题。前置知识：14 章 Bell 数、21.3 高斯消元。建议先做完 P10499 开关问题。
看到 n \le 10 一眼状压，然后就做不下去了。设 f(x) 表示钦定 x 个连通块之间两两不连通的方案数，g(x) 表示恰好有 x 个连通块两两不连通的方案数。根据第二类 Stirling 数定义
f(x)=\sum_{i=x}^n S(i,x) \times g(i)
直接上 Stirling 反演：
g(x)=\sum_{i=x}^{n} (-1)^{i-x} s(i,x) \times f(i)
所求为
g(1)=\sum_{i=1}^n (-1)^{i-1} (i-1)! f(i)
只要求出 f(x)。因为 n \le 10 可以考虑直接爆搜。注意到 B_{10}=115975，考虑搜索哪些点被分到同一个集合。对于每个集合，对于端点所属集合不同的边必须选偶数次。
到这里考虑列方程。设 x_i 表示第 i 个图是否属于子集，那么可以列出若干形如 \bigoplus_{k \in S} x_k=0 的异或方程组。高斯消元以后可以确定一些 x 一定属于 / 不属于该子集，剩下的自由元任选，方案数 2^c，加法原理合并即可。总复杂度 O(n^4B_n)。
实现上，注意输入格式转换成图。同时 n \le 10，不用 bitset 优化，可以直接状压存下。
#define int long long
const int N = 105;
int n, m, a[N][20][20], x[N], id[N], f[N], fac[N];
char s[N * N];
int guass(int tot) {
int cnt = m;
for (int i = 1, cur = 1; i <= tot && cur <= m; cur++) {
for (int j = i; j <= tot; j++) {
if (x[j] >> cur & 1) {swap(x[i], x[j]); break;}
}
if (!(x[i] >> cur & 1)) continue;
for (int j = i + 1; j <= tot; j++) {
if (x[j] >> cur & 1) x[j] ^= x[i];
} cnt--, i++;
} return cnt;
}
void dfs(int step) {
if (step > n) {
int tot = 0;
for (int i = 1; i <= n; i++) {
for (int j = i + 1; j <= n; j++) {
if (id[i] == id[j]) continue;
x[++tot] = 0;
for (int k = 1; k <= m; k++) {
if (a[k][i][j]) x[tot] |= 1LL << k;
}
}
}
return f[id[0]] += 1LL << guass(tot), void();
}
id[step] = ++id[0], dfs(step + 1), id[0]--;
for (int i = 1; i <= id[0]; i++) id[step] = i, dfs(step + 1);
}
void _main() {
cin >> m;
for (int i = 1; i <= m; i++) {
cin >> (s + 1);
int len = strlen(s + 1), tot = 0;
for (int j = 1; ; j++) {
if (j * (j - 1) / 2 == len) {n = j; break;}
}
for (int j = 1; j <= n; j++) {
for (int k = j + 1; k <= n; k++) a[i][j][k] = s[++tot] ^ 48;
}
}
dfs(1);
fac[0] = 1;
for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i;
int res = 0;
for (int i = 1; i <= n; i++) {
if (i & 1) res += fac[i - 1] * f[i];
else res -= fac[i - 1] * f[i];
} cout << res;
}
const int N = 5005;
int n, k;
mint s[N][N], fac[N];
void _main() {
cin >> n >> k;
if (k == 0) return cout << 1, void();
s[0][0] = 1;
for (int i = 1; i <= k; i++) {
for (int j = 1; j <= i; j++) s[i][j] = s[i - 1][j - 1] + s[i - 1][j] * j;
}
mint res = 0, p = 1;
for (int i = 0; i <= min(n, k); i++) {
res += s[k][i] * p * mint(2).pow(n - i);
p *= n - i;
} cout << res;
}
逆用二项式定理，最终得到一个 O(m) 的式子。我们可以 O(m^2) 递推出第二类 Stirling 数和 n 的下降幂，此题解决。
const int N = 1005;
int n, x, p, m, a[N], b[N], s[N][N], low[N];
int madd(int x, int y) {return x += y, x >= p ? x -= p : x;}
int mpow(int a, int b) {
int res = 1; for (a %= p; b; a = 1LL * a * a % p, b >>= 1) {
if (b & 1) res = 1LL * res * a % p;
} return res;
}
void _main() {
cin >> n >> x >> p >> m;
for (int i = 0; i <= m; i++) cin >> a[i];
s[0][0] = 1;
for (int i = 1; i <= m; i++) {
for (int j = 1; j <= i; j++) s[i][j] = madd(s[i - 1][j - 1], 1LL * j * s[i - 1][j] % p);
} // 第二类 Stirling 数
for (int i = 0; i <= m; i++) {
for (int j = i; j <= m; j++) b[i] = madd(b[i], 1LL * a[j] * s[j][i] % p);
} // 下降幂系数
for (int i = 0; i <= m; i++) {
low[i] = 1;
for (int j = 0; j < i; j++) low[i] = 1LL * low[i] * (n - j) % p;
} // n 的下降幂
int res = 0;
for (int i = 0; i <= m; i++) res = madd(res, 1LL * b[i] * mpow(x, i) % p * low[i] % p * mpow(x + 1, n - i) % p);
cout << res;
}
类似套路的题：P6667，但需要用到笔者不会的多项式科技。
14. Bell 数
Bell 数B_n 表示大小为 n 的集合的划分方法数，参见 A000110。例如，对于集合 \{1,2,3\}，有 5 种划分方法：
设集合 A(x)=\{y \mid d(f(x),y)=1\}，则对于 A(x) 中的所有元素 y 有 A(x)=A(y)。把所有相等的 A(x) 视作一个集合，则一个神圣的 S 会导致这 n 位被划分为若干非空子集。
这样我们就得到了，一个合法的 S 与将 m 位划分为若干子集的方案一一对应。因此，S 的数目为 Bell 数 B_m。
回头看 T \subseteq S 这条限制，就是告诉你哪些 f(x) 不属于一个集合。将所有连通块状压预处理出来，由乘法原理，答案即 \prod B_{sz}。时间复杂度 O(m(m+n))。
const int N = 1005;
long long n, m, a[N];
char s[N];
mint b[N][N];
void _main() {
cin >> m >> n;
for (int i = 1; i <= n; i++) {
cin >> (s + 1);
for (int j = 1; j <= m; j++) {
if (s[j] == '1') a[j] |= 1LL << (i - 1);
}
}
map<long long, int> sz;
for (int i = 1; i <= m; i++) sz[a[i]]++;
b[0][0] = 1;
for (int i = 1; i <= m; i++) {
b[i][0] = b[i - 1][i - 1];
for (int j = 1; j <= i; j++) b[i][j] = b[i - 1][j - 1] + b[i][j - 1];
}
mint res = 1;
for (const auto& k : sz) res *= b[k.second][0];
cout << res;
}