预处理结果
问题类型
n组查询：(1 leq n leq 10^4)
每次查询(C_a^b) : (1 leq b leq a leq 2000)
如果按照公式计算，最坏时间复杂度是O(n * b)，当n取(10^4)，a 和 b取 (10^5) 时，总运算次数为 (10^4 * 2 * 10^5 = 2 * 10^9)
算法原理及流程
根据公式 (C_a^b = c_{a - 1}^b + C_{a - 1}^{b - 1}) 进行递推
代码实现
#include <iostream>
using namespace std;
const int N = 2010, mod = 1e9 + 7;
int n;
int c[N][N];
void init()
{
for (int i = 0; i < N; ++ i)
for (int j = 0; j <= i; ++ j) // 注意这里要求j<=i，所以能够处理所有边界情况。从实际意义上来将，从i个物品中选，最多选择i个，同样能够表示j的范围的意义
if (!j) c[i][j] = 1; // c_i^0都等于0
else c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % mod;
}
int main()
{
init();
cin >> n;
while (n --)
{
int a, b;
cin >> a >> b;
cout << c[a][b] << endl;
}
return 0;
}
预处理阶乘和逆元
问题类型
n组查询：(1 leq n leq 10^4)
每次查询(C_a^b) : (1 leq b leq a leq 10^5)
之所以不能预处理结果，是因为a和b的范围过大，用数组下标已经存不下了
算法原理及流程
预处理阶乘和阶乘的逆元。由于根据公式计算涉及到除法，所以取模需要用到逆元
代码实现
解释一下代码中计算阶乘的逆元时的两种写法等价的原因：((k!)^{-1} = (k!)^{mod - 2} = ((k - 1)! * k)^{mod - 2} = ((k - 1)!)^{mod - 2} * (k)^{mod - 2} = ((k - 1)!)^{-1} * (k)^{-1})
#include <iostream>
using namespace std;
typedef long long LL;
const int N = 1e5 + 10, mod = 1e9 + 7;
int n;
LL fact[N], infact[N]; // fact[i]:i的阶层 infact[i]:i的阶乘的逆元
int qpow(int a, int b, int mod)
{
int res = 1;
while (b)
{
if (b & 1) res = (LL)res * a % mod;
a = (LL)a * a % mod;
b >>= 1;
}
return res;
}
void init()
{
fact[0] = infact[0] = 1;
for (int i = 1; i < N; ++ i)
{
fact[i] = (LL)fact[i - 1] * i % mod;
infact[i] = qpow(fact[i], mod - 2, mod); // 因为mod是一个质数，所以可以使用费马小定理计算逆元
// 或者 infact[i] = (LL)infact[i - 1] * qpow(i, mod - 2, mod); 原因见上
}
}
int main()
{
init();
cin >> n;
while (n --)
{
int a, b;
cin >> a >> b;
cout << ((LL)fact[a] * infact[b] % mod * infact[a - b] % mod) << endl;
// 两个大int可能会爆long long，所以需要先%一下
}
return 0;
}
卢卡斯优化
问题类型
n组查询：(1 leq n leq 20)
每次查询(C_a^b) : (1 leq b leq a leq 10^{18})
mod p: (1 leq p leq 10^5)
查询次数不多；该类型具有一个很明显的特征就是 (C_a^b)中的a和b都是大于模数p的，所以可以采用卢卡斯定理进行优化
算法原理及流程
卢卡斯(lucas)定理：(C_a^b equiv C_{a / p}^{b / p} * C_{a \% p}^{b \% p} pmod {p})
首先调用卢卡斯定理进行优化，然后直接使用公式进行计算，除法取模需要求逆元
代码实现
#include <iostream>
using namespace std;
typedef long long LL;
int n;
int qpow(int a, int b, int p)
{
int res = 1;
while (b)
{
if (b & 1) res = (LL)res * a % p;
a = (LL)a * a % p;
b >>= 1;
}
return res;
}
int c(int a, int b, int p) // 调用组合数时候就一定保证a和b是小于p的了，所以a和b定为int
{
if (b > a) return 0;
// 按照公式计算
int res = 1;
for (int i = a - b + 1, j = 1; j <= b; ++ i, ++ j)
{
res = (LL)res * i % p;
// 保证b是小于p的，因为p质数，所以1到b的所有数都一定与p互质，所以计算逆元可以采用费马小定理
res = (LL)res * qpow(j, p - 2, p) % p; // 理解原理需要用到“预处理中间值-代码实现”中提到的等式
}
return res;
}
int lucas(LL a, LL b, int p)
{
if (a < p && b < p) return c(a, b, p);
return (LL)c(a % p, b % p, p) * lucas(a / p, b / p, p) % p; // %p保证运算后的结果小于p，无法再使用lucas优化了，所以直接调用组合式函数即可，但/运算则无法保证
}
int main()
{
cin >> n;
while (n --)
{
int p;
LL a, b;
cin >> a >> b >> p;
cout << lucas(a, b, p) << endl;
}
return 0;
}
质因数分解
问题类型
计算结果不能取模，结果可能很大
对于(c_a^b)，我们只需筛选出(<=a)的所有质因子即可，因为根据下方介绍的计算公式，我们计算的是(1)到(n)中质因子(p)的个数，当(p>n)时，其个数一定为(0)，没有计算的必要了
算法原理及流程
原理
很明显的是，结果很大需要用高精度，如果单纯通过公式进行计算，(c_a^b = frac{(a - b + 1) * ... * a}{1 * ... * b})，需要用到高精度乘法和除法，整起来就有些复杂了，所以考虑将(C_a^b) 转化为 $ p_1^{k_1} * p_2^{k_2} * ... * p_n^{k_n} $，，从而将运算转变为单纯的乘法。
如果对公式 (c_a^b = frac{(a - b + 1) * ... * a}{1 * ... * b}) 中的每个数进行一次质因数分解，复杂度就是 (O(n sqrt[2]{n}))，估计也没什么大问题，但是代码写起来还是有点复杂的，所以这里利用一个很神奇的计算方法：
a! 中 质因子 p 的个数为 (lfloor frac{a}{p} floor + lfloor frac{a}{p^2} floor + ... + lfloor frac{a}{p^n} floor).
这里的(lfloor frac{a}{p} floor)具有一个特定的含义
1到n中p的倍数(或者说能被p整除)的数有 (lfloor frac{n}{p} floor) 个(该性质的详细描述请参照 倍数个数的计算)
根据(lfloor frac{n}{p}
floor)的含义，(lfloor frac{n}{p}
floor) 计算的为1到n中p的倍数个数； (lfloor frac{n}{p^2}
floor) 计算的为1到n中(p^2)的倍数个数。。。以此类推
将这些倍数个数相加，计算得出的就是1到n的所有数中质因子p所占据的个数。
不妨举个例子来理解这个方法，计算6！中质因子2的个数: (6! = 1 * 2 * 3 * 4 * 5 * 6), 其中2，4，6是2的倍数，4是(2^2)的倍数，而2，6对质因子2个数的贡献值为1，4的贡献值为2。4中有两个2，在(lfloor frac{a}{p}
floor)时考虑到的是它其中一个2，(lfloor frac{a}{p^2}
floor)考虑的则是另一个2。也就是说1到6的所有数，2的倍数有3个，(2^2)的倍数有1个，所以1到6中所有数中质因子2的个数为4个。
流程
(C_a^b = frac{a!}{b! * (a - b)!}),我们首先预处理出所有质因子，然后分别计算出所有质因子的个数，最后调用高精度乘法计算得出最后结果
代码实现
#include <iostream>
#include <vector>
using namespace std;
const int N = 5010;
int sum[N];
int primes[N], cnt;
bool st[N];
vector<int> mul(vector<int> &a, int b) // 高精度乘法
{
vector<int> c;
int t = 0;
for (int i = 0; i < a.size() || t; ++ i)
{
if (i < a.size()) t += a[i] * b;
c.push_back(t % 10);
t /= 10;
}
while (c.size() != 1 && c.back() == 0) c.pop_back();
return c;
}
void get_primes(int n) // 打素数表
{
st[0] = st[1] = true;
for (int i = 2; i <= n; ++ i)
{
if (!st[i]) primes[cnt ++] = i;
for (int j = 0; primes[j] <= n / i; ++ j)
{
st[i * primes[j]] = true;
if (i % primes[j] == 0) break;
}
}
}
int get(int a, int p) // 返回a!中质因子p的个数
{
int res = 0;
while (a)
{
res += a / p;
a /= p;
}
return res;
}
int main()
{
int a, b;
cin >> a >> b;
// 预处理质数表
get_primes(a);
// 计算公式中每个质数的个数
for (int i = 0; i < cnt; ++ i)
{
int p = primes[i];
sum[i] = get(a, p) - get(b, p) - get(a - b, p);
}
// 高精度乘法计算最终结果
vector<int> res;
res.push_back(1);
for (int i = 0; i < cnt; ++ i)
for (int j = 0; j < sum[i]; ++ j)
res = mul(res, primes[i]);
// 输出答案
for (int i = res.size() - 1; ~i; -- i) cout << res[i];
cout << endl;
return 0;
}
应用
计算卡特兰数(Catalan)
该类问题的难点在于如何将实际问题抽象为求解卡特兰数，而非代码实现
公式推导
下面以一个例子说明卡特兰数的公式是如何推导的
给定n个0和n个1，它们将按照某种顺序排成长度为2n的序列，求它们能排列成的所有序列中，能够满足任意前缀序列中0的个数都不少于1的个数的序列有多少个
问题分析
需要把问题变为一个几何问题，在一个网格图中，类似哈夫曼编码，我们规定 向右走为0，向上走为1
假设分别有6个0和1，任意一个序列如果映射到网格图中，就是一条从(0, 0)到(6, 6)的路径，0的数量代表着横坐标的值，1的数量代表纵坐标的值，所以满足要求的序列就是保证路径上任意一个点的横坐标>=纵坐标，也就是说从(0, 0)到(6, 6)的直线下方的路径才是合法路径，合法序列数量等价为合法路径数量，也就是我们需要求解的从(0, 0)到(6, 6)且不经过对角线的路径总数
求解
按照公式合法方案数 = 总方案数 - 非法方案数，总方案数 = (C_{12}^6)，问题的难点在于如何求解非法方案数。观看下图，我们首先画出一条非法的路径，交对角线于1号点，将非法路径中1号点之后的路径关于对角线做对称可以得到3号线到(5, 7)。很神奇的是不管1号点在哪里，3号线的终点都是(5, 7)(我没想明白为什么)。也就是说任意一条从(0, 0)到(6, 6)的非法路径都可以等价为从(0, 0)到(5, 7)的路径，所以非法路径数 = (C_{12}^5) 或 (C_{12}^7)。
将上述描述总结为公式就是：(C_{2n}^{n} - C_{2n}^{n - 1} = frac{C_{2n}^{n}}{n + 1}) 或 (C_{2n}^{n} - C_{2n}^{n + 1} = frac{C_{2n}^{n}}{n + 1})，但一般都记前面那个式子
代码实现
只需要按照公式计算即可，实际采用的是卢卡斯优化代码中组合数的求解方法
#include <iostream>
using namespace std;
typedef long long LL;
const int mod = 1e9 + 7;
int qpow(int a, int b, int mod)
{
int res = 1;
while (b)
{
if (b & 1) res = (LL)res * a % mod;
a = (LL)a * a % mod;
b >>= 1;
}
return res;
}
int main()
{
int n;
cin >> n;
int res = 1;
for (int i = n + 1, j = 1; j <= n; ++ i, ++ j)
{
res = (LL)res * i % mod;
res = (LL)res * qpow(j, mod - 2, mod) % mod;
}
res = (LL)res * qpow(n + 1, mod - 2, mod) % mod;
cout << res << endl;
return 0;
}