组合计数
Q4.2.2.2. 求组合数·土
考虑 \(\bmod p^k\) 的情况
递归计算阶乘,由于要计算逆元,需要除掉所有p的倍数
将阶乘中所有p的倍数拿出来, 递归计算...
点击查看代码
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <utility>
#include <queue>
using namespace std;
typedef long long LL;
LL qpow(LL b, LL p, LL mod) {
LL res = 1;
while(p) {
if(p & 1) res = (LL)res * b % mod;
b = (LL)b * b % mod, p >>= 1;
}
return res;
}
LL n, m, mod, prime, power;
LL fac(LL n) { // 求 1到n 的阶乘除去所有prime模power的结果
if(n == 0) return 1;
LL res = 1;
for(int i = 1; i <= power; i ++)
if(i % prime) (res *= i) %= power;
res = qpow(res, n / power, power);
for(int i = 1; i <= n % power; i ++)
if(i % prime) (res *= i) %= power;
return res * fac(n / prime) % power;
}
LL calc(LL Prime, LL Power) {
prime = Prime, power = Power;
LL k = 0, t, tmp = power - power / prime - 1; // k 表示 prime 的指数, tmp 为欧拉定理求逆元的指数
t = n; while(t) t /= prime, k += t;
t = m; while(t) t /= prime, k -= t;
t=n-m; while(t) t /= prime, k -= t;
return fac(n) * qpow(fac(m), tmp, power) % power * qpow(fac(n-m), tmp, power) % power * qpow(prime, k, power) % power;
}
LL exgcd(LL a, LL b, LL &x, LL &y) {
if(!b) return x = 1, y = 0, a;
LL d = exgcd(b, a % b, y, x);
return y -= a / b * x, d;
}
LL A[100], M[100];
LL a1, m1, a2, m2, d, p, q, x, c;
int tot;
bool solve() { // CRT
a1 = A[1], m1 = M[1];
for(int i = 2; i <= tot; i ++) {
a2 = A[i], m2 = M[i];
c = a1 - a2;
d = exgcd(m1, m2, p, q);
if(c % d) return true;
x = q * c / d;
m1 = m1 / d * m2;
a1 = (x % m1 * m2 % m1 + a2 + m1) % m1;
}
return false;
}
int main() {
int T;
scanf("%d", &T);
while(T --) {
tot = 0;
scanf("%lld%lld%lld", &n, &m, &mod);
for(int i = 2; (LL)i * i <= mod; i ++)
if(mod % i == 0) {
int pw = 1;
while(mod % i == 0) mod /= i, pw *= i;
A[++ tot] = calc(i, pw), M[tot] = pw;
}
if(mod > 1) A[++ tot] = calc(mod, mod), M[tot] = mod;
solve();
printf("%lld\n", a1 % m1);
}
return 0;
}
P2183 [国家集训队]礼物
计算组合数
\(ans=\frac{n!}{((n-\sum w)!w_1!w_2!\cdots w_m!}\mod p\)
求阶乘: exlucas的代码（类似上题）
点击查看代码
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <utility>
#include <queue>
using namespace std;
typedef long long LL;
LL qpow(LL b, LL p, LL mod) {
LL res = 1;
while(p) {
if(p & 1) res = (LL)res * b % mod;
b = (LL)b * b % mod, p >>= 1;
}
return res;
}
LL n, m, w[10], sum, mod, prime, power;
LL fac(LL n) { // 求 1到n 的阶乘除去所有prime模power的结果
if(n == 0) return 1;
LL res = 1;
for(LL i = 1; i <= power; i ++)
if(i % prime) (res *= i) %= power;
res = qpow(res, n / power, power);
for(LL i = 1; i <= n % power; i ++)
if(i % prime) (res *= i) %= power;
return res * fac(n / prime) % power;
}
LL calc(LL Prime, LL Power) {
prime = Prime, power = Power;
LL k = 0, t, tmp = power - power / prime - 1; // k 表示 prime 的指数, tmp 为欧拉定理求逆元的指数
t = n;
while(t) t /= prime, k += t;
for(LL i = 0; i <= m; i ++) {
t = w[i];
while(t) t /= prime, k -= t;
}
LL res = (__int128_t)fac(n) * qpow(prime, k, power) % power;
for(LL i = 0; i <= m; i ++) res = (__int128_t)res * qpow(fac(w[i]), tmp, power) % power;
return res;
}
LL exgcd(LL a, LL b, LL &x, LL &y) {
if(!b) return x = 1, y = 0, a;
LL d = exgcd(b, a % b, y, x);
return y -= a / b * x, d;
}
LL A[100], M[100];
LL a1, m1, a2, m2, d, p, q, x, c;
LL tot;
bool solve() { // CRT
a1 = A[1], m1 = M[1];
for(LL i = 2; i <= tot; i ++) {
a2 = A[i], m2 = M[i];
c = a1 - a2;
d = exgcd(m1, m2, p, q);
if(c % d) return true;
x = q * c / d;
m1 = m1 / d * m2;
a1 = (x % m1 * m2 % m1 + a2 + m1) % m1;
}
return false;
}
int main() {
scanf("%lld%lld%lld", &mod, &n, &m);
for(LL i = 1; i <= m; i ++) scanf("%lld", w + i), sum += w[i];
w[0] = n - sum;
if(sum > n) return puts("Impossible"), 0;
for(LL i = 2; (LL)i * i <= mod; i ++)
if(mod % i == 0) {
LL pw = 1;
while(mod % i == 0) mod /= i, pw *= i;
A[++ tot] = calc(i, pw), M[tot] = pw;
}
if(mod > 1) A[++ tot] = calc(mod, mod), M[tot] = mod;
solve();
printf("%lld\n", a1 % m1);
return 0;
}
P1350 车的放置
将棋盘分为上下两部分求解
答案为 \(\sum\limits_{i=0}^{k}C_b^iA_a^iC_d^{k-i}A_{a+c-i}^{k-i}\)
其中 i 枚举的是 上半部分的车的个数
点击查看代码
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;
typedef long long LL;
const int N = 2005, mod = 100003;
int a, b, c, d, k;
int fac[N], invfac[N];
int qpow(int b, int p, int mod) {
int res = 1;
while(p) {
if(p & 1) res = (LL)res * b % mod;
b = (LL)b * b % mod, p >>= 1;
}
return res;
}
int C(int a, int b) {
if(a < b) return 0;
return (LL)fac[a] * invfac[b] % mod * invfac[a - b] % mod;
}
int A(int a, int b) {
if(a < b) return 0;
return (LL)fac[a] * invfac[a - b] % mod;
}
int main() {
fac[0] = invfac[0] = 1;
for(int i = 1; i < N; i ++) {
fac[i] = (LL)fac[i - 1] * i % mod;
invfac[i] = qpow(fac[i], mod - 2, mod);
}
scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
int res = 0;
for(int i = 0; i <= k; i ++) {
res = (res + (LL)C(b, i) * A(a, i) % mod * C(d, k - i) % mod * A(a + c - i, k - i) % mod) % mod;
}
printf("%d\n", res);
return 0;
}
P2480 [SDOI2010] 古代猪文
给定 \(n,q\) 求 \(q^{\sum\limits_{d|n}C_n^d}\mod质数999911659\)
只需求 \(\sum\limits_{d|n}C_n^d \mod 999911658 = 2*3*4679*35617\)
使用Lucas求解, 再用CRT即可
点击查看代码
#include <stdio.h>
typedef long long LL; typedef __int128_t LLL;
const LL mod = 999911659, primes[] = {2, 3, 4679, 35617};
inline LL crt(LL a, LL b, LL c, LL d) // CRT
{ return (d * (LLL)877424796 + c * (LLL)289138806 + b * (LLL)333303886 + a * (LLL)499955829) % (mod - 1); }
LL n, q, p, val[4], n2;
LL qpow(LL b, LL k, LL m) {
LL res = 1;
while(k) {
if(k & 1) res = res * b % m;
b = b * b % m, k >>= 1;
}
return res;
}
LL calc(LL a, LL b) { // C_{a}^{b} 的朴素求法
if(a < b) return 0;
LL A = 1, B = 1;
for(int i = a, j = 1; j <= b; i --, j ++)
A = A * i % p, B = B * j % p;
return A * qpow(B, p - 2, p) % p;
}
LL Lucas(LL a, LL b) { // Lucas 求 C_{a}^{b}
if(a < p && b < p) return calc(a, b);
return Lucas(a / p, b / p) * calc(a % p, b % p) % p;
}
void modify(LL d) {
for(int i = 0; i < 4; i ++)
p = primes[i], (val[i] += Lucas(n2, d)) %= p;
}
int main() {
scanf("%lld%lld", &n, &q), n2 = n;
if(q % mod == 0) return puts("0"), 0; // 特判
for(int i = 1; (LL)i * i <= n; i ++)
if(n % i == 0) modify(i), i * i != n && (modify(n / i),1);
printf("%lld\n", qpow(q, crt(val[0], val[1], val[2], val[3]), mod));
return 0;
}
P4345 [SHOI2015]超能粒子炮·改
点击查看代码
#include <stdio.h>
#include <string.h>
typedef long long LL;
int C[2345][2345];
int S[2345][2345]; // S[i][j] 表示 C[i][0]+C[i][1]+C[i][2]+...+C[i][j] 的结果
int N[7], M[7], f[7][2]; // f : 数位DP, 0/1 表示这一位有没有上限
LL pow2333[7];
int main() {
for(int i = 0; i <= 2333; i ++) {
C[i][0] = S[i][0] = 1;
for(int j = 1; j <= i; j ++) C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % 2333;
for(int j = 1; j <= 2333; j ++) S[i][j] = (S[i][j - 1] + C[i][j]) % 2333;
}
for(int i = *pow2333 = 1; i < 6; i ++) pow2333[i] = pow2333[i - 1] * 2333;
int T; LL n, m;
scanf("%d", &T);
while(T --) {
scanf("%lld%lld", &n, &m);
for(int i = 5; i >= 0; i --) N[i] = n / pow2333[i] % 2333, M[i] = m / pow2333[i] % 2333;
memset(f, 0, sizeof(f)); // 数位DP, 0/1 表示这一位有没有上限
f[6][0] = 0, f[6][1] = 1;
for(int i = 5; i >= 0; i --) {
f[i][0] = (f[i+1][0] * S[N[i]][2332] % 2333 + (M[i] ? f[i+1][1] * S[N[i]][M[i]-1] % 2333 : 0)) % 2333;
f[i][1] = f[i+1][1] * C[N[i]][M[i]] % 2333;
}
printf("%d\n", (f[0][0] + f[0][1]) % 2333);
}
return 0;
}
P3166 [CQOI2014]数三角形
n 表示横向有 n 个点, m 表示纵向有 m 个点
答案为 \(c_{nm}^3\) 减去不合法即三点共线的点
不合法的点:
- 斜率为 \(0~~\) : \(nc_m^3\)
- 斜率为 \(\infin\) : \(mc_n^3\)
- 斜率 >0
先枚举左下角的点, 再枚举右上角的点, 确定中间的点
=> 先枚举左下角的点与右下角的点构成的直角三角形的长和宽, 设为 \(i(\in[1,n])\) 和 \(j(\in[1,m])\)
则这个直角三角形的斜边上有 \(\mathbb{gcd}(i,j)-1\) 个整点(不包含端点)
这个直角三角形再矩形中有 \((n-i)(m-j)\) 中可能 \(\Rightarrow\) 方案数为 \((\mathbb{gcd}(i,j)-1)(n-i)(m-j)\) - 斜率 <0 同 3
点击查看代码
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;
typedef long long LL;
LL gcd(LL x, LL y) {
while(x) {
y -= (y / x) * x;
x ^= y ^= x ^= y;
}
return y;
}
// c_n^3
LL C3(LL n) {
return ((n * (n - 1)) >> 1) * (n - 2) / 3;
}
int main() {
LL n, m;
scanf("%lld%lld", &n, &m);
n ++, m ++; // 原题为格子的个数
LL res = C3(n * m) - n * C3(m) - m * C3(n);
for(LL i = 1; i <= n; i ++)
for(LL j = 1; j <= m; j ++)
res -= 2 * (gcd(i, j) - 1) * (n - i) * (m - j);
printf("%lld\n", res);
return 0;
}
BZOJ 2982 Combination
此题使用Lucas定理求组合数
Lucas定理：设 \(p\) 为质数，则 \(C_a^b\equiv C_{a/p}^{b/p} C_{a\bmod p}^{b\bmod p}\pmod p\)
点击查看代码
#include <stdio.h>
typedef long long LL;
LL mod;
LL qpow(LL b, LL p) {
b %= mod;
LL res = 1;
while(p) {
if(p & 1) res = res * b % mod;
b = b * b % mod, p >>= 1;
}
return res;
}
LL C(LL a, LL b) {
LL A = 1, B = 1;
for(int i = a, j = 1; j <= b; i --, j ++)
A = A * i % mod, B = B * j % mod;
return A * qpow(B, mod - 2) % mod;
}
LL Lucas(LL a, LL b) {
if(a < mod || b < mod) return C(a, b);
return Lucas(a / mod, b / mod) * C(a % mod, b % mod) % mod;
}
int main() {
mod = 10007;
int T;
LL a, b;
scanf("%d", &T);
while(T --) {
scanf("%lld%lld", &a, &b);
printf("%lld\n", Lucas(a, b));
}
return 0;
}
BZOJ 3907 网格
使用排除法解决问题，用“折线法”：将第一个越过 \(y=x+1\) 的点及之前的按 \(y=x+1\) 翻转即可
答案为 \(C_{n+m}^{n}-C{n+m}^{n+1}\)
点击查看代码
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <vector>
using namespace std;
const int N = 10005;
int primes[N], cnt;
bool st[N];
void get_primes(int n) {
for(int i = 2; i <= n; i ++) {
if(!st[i]) primes[cnt ++] = i;
for(int j = 0; primes[j] * i <= n; j ++) {
st[i * primes[j]] = true;
if(i % primes[j] == 0) break;
}
}
}
vector<int> mul(const vector<int> &a, const vector<int> &b) {
int sa = a.size(), sb = b.size(), sc = sa + sb;
vector<int> c(sc + 1);
for(int i = 0; i < sa; i ++)
for(int j = 0; j < sb; j ++) c[i + j] += a[i] * b[j];
for(int i = 0; i < sc; i ++) c[i + 1] += c[i] / 10, c[i] %= 10;
while(c.size() > 1U && !c.back()) c.pop_back();
return c;
}
int n, m;
int get(int p, int t) {
int res = 0;
while(t) {
t /= p;
res += t;
}
return res;
}
vector<int> qpow(vector<int> b, int p) {
vector<int> res({1});
while(p) {
if(p & 1) res = mul(res, b);
b = mul(b, b), p >>= 1;
}
return res;
}
int main() {
scanf("%d%d", &n, &m);
get_primes(n + m);
vector<int> res({1});
for(int i = 0; i < cnt; i ++) {
int x = get(primes[i], n + m) - get(primes[i], n) - get(primes[i], m);
int t = n - m + 1;
while(t % primes[i] == 0) t /= primes[i], x ++;
t = n + 1;
while(t % primes[i] == 0) t /= primes[i], x --;
t = primes[i];
vector<int> tmp;
while(t) {
tmp.push_back(t % 10);
t /= 10;
}
tmp = qpow(tmp, x);
res = mul(res, tmp);
}
for(int i = res.size() - 1; i >= 0; i --) putchar(res[i] ^ 48);
return 0;
}
P3200 [HNOI2009] 有趣的数列
方法与上一题类似
卡特兰数 \(C_{2n}^{n}\over{n+1}\)
点击查看代码
#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;
typedef long long LL;
const int N = 2e6 + 5;
bool st[N];
vector<int> primes;
LL n, p;
LL qpow(LL b, LL k) {
LL res = 1;
while(k) {
if(k & 1) res = res * b % p;
b = b * b % p, k >>= 1;
}
return res;
}
void getprimes() {
int tag = 2 * n;
for(int i = 2; i <= tag; i ++) {
if(!st[i]) primes.push_back(i);
for(int j = 0; primes[j] <= tag / i; j ++) {
st[primes[j] * i] = true;
if(i % primes[j] == 0) break;
}
}
}
LL numberof(LL x, LL y) {
LL res = 0, t = x;
while(true) {
if(y < t) return res;
res += y / t, t *= x;
}
}
int main() {
cin >> n >> p;
getprimes();
LL res = 1;
for(int i = 0; i < primes.size(); i ++) {
LL pp = (LL)primes[i];
LL t = numberof(pp, 2 * n) - 2 * numberof(pp, n);
int __n = n + 1;
while(__n % pp == 0) t --, __n /= pp;
res = res * qpow(pp, t) % p;
}
cout << res << endl;
return 0;
}
P2532 [AHOI2012] 树屋阶梯
DP：枚举占了左下角的钢材的大小，剩下的就是子问题了
递推：（就是卡特兰数）\(f_0=1,f_n=f_0f_{n-1}+f_1f_{n-2}+\cdots+ f_{n-1}f_0\)
点击查看代码
f = [1]
n = int(input())
for i in range(1, n + 1):
now = 0
for j in range(0, i):
now += f[j] * f[i - 1 - j]
f.append(now)
print(f[n])
P1655 小朋友的球
第2类Stirling数：将 \(n\) 个不同的球放入 \(m\) 个相同的盒子中，无空盒的方案数
f[i][1] = 1, f[i][j] = f[i-1][j-1] + j * f[i-1][j]
容斥原理
枚举空盒子的个数，答案为 \(\sum\limits_{i=0}^{m}(-1)^iC_{m}^{i}(m-i)^n\)
点击查看代码
#include <stdio.h>
#include <vector>
using std::vector;
vector<int> f[105][105];
vector<int> operator + (const vector<int> &a, const vector<int> &b) {
size_t sa = a.size(), sb = b.size(), sc = (sa < sb ? sb : sa);
vector<int> c(sc + 1);
for(size_t i = 0; i < sc; i ++) {
if(i < sa) c[i] += a[i];
if(i < sb) c[i] += b[i];
if(c[i] >= 10) c[i] -= 10, c[i + 1] ++;
}
while(c.size() > 1U && !c.back()) c.pop_back();
return c;
}
vector<int> operator * (const vector<int> &a, int b) {
vector<int> c(a.size() + 4U);
for(size_t i = 0; i < a.size(); i ++) c[i] = a[i] * b;
for(size_t i = 0; i + 1 < c.size(); i ++) c[i + 1] += c[i] / 10, c[i] %= 10;
while(c.size() > 1U && !c.back()) c.pop_back();
return c;
}
int main() {
for(int i = 1; i <= 100; i ++) {
f[i][1] = {1};
if(i >= 2)
for(int j = 1; j <= i; j ++)
f[i][j] = f[i - 1][j - 1] + f[i - 1][j] * j;
}
int n, m;
while(scanf("%d%d", &n, &m) == 2) {
if(f[n][m].size()) for(int i = f[n][m].size() - 1; i >= 0; i --) putchar(f[n][m][i] ^ 48);
else putchar(48);
putchar(10);
}
return 0;
}
P4071 [SDOI2016] 排列计数
错排问题：递推公式 \(d_1=0,d_2=1,d_n=(n-1)(d_{n-1}+d_{n-2})\)
点击查看代码
#include <stdio.h>
typedef long long LL;
const int N = 1e6, mod = 1e9 + 7;
int d[N + 5], n = N, m, T;
int fac[N + 5], inv[N + 5];
int qpow(int b, int p) {
int res = 1;
while(p) {
if(p & 1) res = (LL)res * b % mod;
b = (LL)b * b % mod, p >>= 1;
}
return res;
}
int main() {
fac[0] = inv[0] = 1;
for(int i = 1; i <= n; i ++)
fac[i] = (LL)fac[i - 1] * i % mod, inv[i] = qpow(fac[i], mod - 2);
d[0] = 1, d[2] = 1;
for(int i = 3; i <= n; i ++) d[i] = (LL)(i - 1) * ((LL)d[i - 1] + d[i - 2]) % mod;
scanf("%d", &T);
while(T --) {
scanf("%d%d", &n, &m);
printf("%lld\n", (LL)fac[n] * inv[m] % mod * inv[n - m] % mod * d[n - m] % mod);
}
return 0;
}
P6475[NOI Online #2 入门组] 建设城市
分类讨论：地标是否在同一侧；枚举相同时的高度
使用阶乘预处理可以 O(1) 求出组合数
点击查看代码
#include <stdio.h>
const int N = 2e5 + 5, mod = 998244353;
typedef long long LL;
inline int qpow(int b, int p) {
int res = 1;
while(p) {
if(p & 1) res = (LL)res * b % mod;
b = (LL)b * b % mod, p >>= 1;
}
return res;
}
int fac[N], inv[N];
inline int C(int n, int m) { // C_{n}^{m}
return (LL)fac[n] * inv[m] % mod * inv[n - m] % mod;
}
int n, m, x, y, res;
int main() {
scanf("%d%d%d%d", &m, &n, &x, &y);
fac[0] = inv[0] = 1;
for(int i = 1; i <= n + m; i ++)
fac[i] = (LL)fac[i - 1] * i % mod,
inv[i] = qpow(fac[i], mod - 2);
if(x <= n && y > n) {
y = 2 * n + 1 - y;
for(int z = 1; z <= m; z ++)
(res += (LL)C(x - 2 + z, x - 1) * C(n - x + m - z, n - x) % mod * C(y - 2 + z, y - 1) % mod * C(n - y + m - z, n - y) % mod) %= mod;
} else {
res = (LL)C(n + m - 1, m - 1) * C(n + x - y + m - 1, m - 1) % mod;
}
printf("%d\n", res);
return 0;
}
