一个有趣的问题
题目大意
给定$n, m, p$，求$C_{n}^{m}$除以$p$后的余数。
Subtask#1 $0leqslant mleqslant n leqslant 2 imes 10^{3}$
直接杨辉恒等式$C_{n}^{m} = C_{n - 1}^{m - 1} + C_{n - 1}^{m}$递推。
时间复杂度$O(n^{2})$。
Subtask#2 $0leqslant mleqslant n leqslant 2 imes 10^{6}，p = 10^{9} + 7$
用式子$C_{n}^{m} = frac{n!}{m!(n - m)!}$来计算。
直接算阶乘，用快速幂或扩欧求逆元即可。
时间复杂度$O(n + log n)$
Subtask#3 $0leqslant mleqslant n leqslant 10^{9},pleqslant 10^{6}$且$p$为质数。
好吧，现在引入正题。
定理1（Lucas定理） 当$p$是一个质数时，$n = k_{1}p + r_{1}, m = k_{2}p + r_{2} left ( 0leqslant r_{1},r_{2} < p ight )$，则有$C_{n}^{m} equiv C_{k_{1}}^{k_{2}}C_{r_{1}}^{r_{2}} pmod{p}$
考虑怎么使用它。
对于每次操作的$r_{1}$和$r_{2}$都是小于$p$的，这个可以直接用Subtask 2的方法来做。
然后剩下的$C_{k_{1}}^{k_{2}}$就递归处理。
这很简单，就不给代码了。
现在来考虑Lucas定理的证明。
首先你需要一个二项式定理。
定理2（二项式定理） 当$n$是一个非负整数时，则有，$left (a + b ight ) ^{n} = sum_{i = 0} ^{n}C_{n}^{i}a^{i}b^{n - i}$
证明 考虑使用数学归纳法来证明。
当$n = 1$的时候，显然成立。
假设当$n = k - 1$的时候成立，考虑当$n = k$的时候
$left (a + b ight )^{k} = left ( a + b ight )^{k - 1}left ( a + b ight )\=left ( sum_{i = 0}^{k - 1}C_{k-1}^{i}a^{i}b^{k - i - 1} ight )(a + b)\=sum_{i = 1}^{k}C_{k - 1}^{i - 1}a^{i}b^{k - i} + sum_{i = 0}^{k - 1}C_{k - 1}^{i}a^{i}b^{k - i}\=C_{k - 1}^{k - 1}a^{n} + sum_{i = 1}^{k - 1}left ( C_{k - 1}^{i - 1} + C_{k - 1}^{i} ight )a^{i}b^{k - i} + C_{k - 1}^{0}b^{n}\=sum_{i = 0}^{k} C_{k}^{i} a^{i}b^{k - i}$
因此，定理得证。
定理3 若$p$为质数，若有整数$k$满足$1leqslant k < p$，则有$p mid C_{p}^{k}$
根据式子$C_{p}^{k} = frac{p!}{k!left( p - k ight)!}$和素数的定义易证。
推论4 若$p$为质数，则$left ( 1 + x ight )^{p}equiv 1 + x^{p} pmod{p}$。
根据定理3和二项式定理易证。
现在来考虑Lucas定理的证明。
Lucas定理的证明
$left ( 1 + x ight )^{n}\equiv left ( 1 + x ight )^{k_{1}p}left( 1 + x ight )^{r_{1}} \equiv left ( 1 + x^{p} ight )^{k_{1}}left( 1 + x ight )^{r_{1}}\equiv left ( sum_{i = 0}^{k_{1}}C_{k_{1}}^{i}x^{pi} ight )left ( sum_{i = 0}^{r_{1}}C_{r_{1}}^{i}x^{i} ight )\equiv sum_{i = 0}^{k_{1}}sum_{j = 0}^{r_{1}}C_{k_{1}}^{i}C_{r_{1}}^{j}x^{pi + j} pmod{p}$
然后考虑$x^{m}$的系数，同余式左边显然是$C_{n}^{m}$，右边是$C_{k_{1}}^{k_{2}}C_{r_{1}}^{r_{2}}$，因此定理得证。
扩展Lucas算法
继续回到上面的问题。
Subtask#4 $0 leqslant m leqslant n leqslant 10^{7},pleqslant 10^{9}$
继续考虑组合数的计算公式$C_{n}^{m} = frac{n!}{m!(n - m)!}$，考虑将$n!, m!, (n - m)!$质因数分解，这样直接就可以指数相加减，最后乘起来就好了。
考虑怎么数$n!$中质因子$p$的指数。
我们知道$n! = 1 imes 2 imes 3 imes cdots imes n$，其中是$p$的倍数的是$p, 2p, 3p, cdots$，考虑将这些$p$除掉，然后它会变成$1, 2, 3, cdots, left lfloor frac{n}{p} ight floor$。这是一个规模更小的子问题，递归处理即可。时间复杂度$O(log_{p}n)$。
可以用线性筛先预处理$10^{7}$以内的素数，然后考虑$n!$的时候对小于$n$的每个素数去跑一次这个算法，因为$n$以内约有$frac{n}{log n}$个素数，所以时间复杂度为$O(n)$。
Subtask#5 $0 leqslant m leqslant n leqslant 10^{9},pleqslant 10^{9}$，若将$p$质因数分解得到：$p = p_{1}^{a_{1}}cdots p_{k}^{a_{k}}$,，则对于任意$i = 1, 2, 3,cdots, k$都有$p_{i}^{a_{i}}leqslant 10^{5}$
数据范围怎么越来越不友好了？话说后面那一大句又是什么奇怪的性质？
继续考虑直接暴力，上面通过质因数分解的方法已经不可行。
以上NB的方法是阶乘加逆元对吧？因为模数不一定是质数，所以和$n$之类的不一定互质，也就是说逆元可能不存在。
设答案为$x$，那么可以将原问题拆成$k$个子问题：
$left{egin{matrix}xequiv C_{n}^{m} pmod{p_{1}^{a_{1}}} \ xequiv C_{n}^{m} pmod{p_{2}^{a_{2}}}\ vdots \ xequiv C_{n}^{m} pmod{p_{k}^{a_{k}}}end{matrix} ight.$
这有什么卵用？
假如能够把$p_{i}$抽出去单独计算，剩下就可以"快速阶乘"，因为此时和模数互质，所以可以直接用逆元在模意义下做除法。最后用中国剩余定理合并。
把$p_{i}$抽出去可以用上面的方法进行计算。
现在考虑如何快速计算模$p_{i}^{a_{i}}$下，被抽取所有有关$p_{i}$的因子的$n!$。
我们直接暴力枚举$1$ ~ $p_{i}^{a_{i}}$中间的数，然后暴力将不是$p_{i}$的倍数的数乘起来，对于中间在模意义下是一样的，快速幂乘一乘，最后末尾的暴力做一下。
然后把是$p_{i}$的倍数的除掉一个$p_{i}$，这样问题的规模会变小，递归处理即可
举个例子有助于说明：
假设$p_{i} = 3, a_{i} = 2, n = 19$，那么要计算的是：$1 imes 2 imes 3 imes 4 imes 5 imes 6 imes 7 imes 8 imes 9 imes 10 imes 11 imes 12 imes 13 imes 14 imes 15 imes 16 imes 17 imes 18 imes 19$
显然它同余于$left (1 imes 2 imes 4 imes 5 imes 7 imes 8 ight ) imes left (1 imes 2 imes 4 imes 5 imes 7 imes 8 ight ) imes 1 imes 3 imes left ( 1 imes 2 imes 3 imes 4 imes 5 imes 6 ight )$
前面的用快速幂，末尾余下的暴力算，后面的递归处理。
时间复杂度O（能过）
没错，这就是扩展Lucas，我很好奇它和Lucas有什么关系。这明明就是依赖特殊性质的暴力qwq。
最后附上代码。
Code
1 /** 2 * bzoj 3 * Problem#2142 4 * Accepted 5 * Time: 372ms 6 * Memory: 1292k 7 */ 8 #include <bits/stdc++.h> 9 #ifndef WIN32 10 #define Auto "%lld" 11 #else 12 #define Auto "%I64d" 13 #endif 14 using namespace std; 15 typedef bool boolean; 16 17 #define ll long long 18 19 int qpow(int a, int pos, int m) { 20 int pa = a, rt = 1; 21 for ( ; pos; pos >>= 1, pa = pa * 1ll * pa % m) 22 if (pos & 1) 23 rt = rt * 1ll * pa % m; 24 return rt; 25 } 26 27 void exgcd(ll a, ll b, ll& d, ll& x, ll& y) { 28 if (!b) 29 d = a, x = 1, y = 0; 30 else { 31 exgcd(b, a % b, d, y, x); 32 y -= (a / b) * x; 33 } 34 } 35 36 ll inv(ll a, ll n) { 37 ll d, x, y; 38 exgcd(a, n, d, x, y); 39 return (x < 0) ? (x + n) : (x); 40 } 41 42 ll p; 43 int n, m, k; 44 int top = 0; 45 int w[8]; 46 int fac[65]; 47 int val[65]; 48 49 inline void init() { 50 scanf(Auto"%d%d", &p, &n, &m); 51 ll ws = 0; 52 for (int i = 1; i <= m; i++) { 53 scanf("%d", w + i); 54 ws += w[i]; 55 } 56 if (ws > n) { 57 puts("Impossible"); 58 exit(0); 59 } 60 k = ws; 61 } 62 63 void getFactors(ll x) { 64 for (int i = 2; x != 1; i++) { 65 if (!(x % i)) { 66 fac[++top] = i, val[top] = 1; 67 while (!(x % i)) val[top] *= i, x /= i; 68 } 69 } 70 } 71 72 int getPos(int n, int p) { 73 int rt = 0; 74 while (n) rt += (n /= p); 75 return rt; 76 } 77 78 int calc(int n, int p, int pr) { 79 if (!n) return 1; 80 int rt = 1; 81 for (int i = 1; i < pr; i++) 82 if (i % p) 83 rt = rt * 1ll * i % pr; 84 rt = qpow(rt, n / pr, pr); 85 for (int i = 1; i <= n % pr; i++) 86 if (i % p) 87 rt = rt * 1ll * i % pr; 88 return rt * 1ll * calc(n / p, p, pr) % pr; 89 } 90 91 ll exLucas(int n, int m, ll ap) { 92 ll rt = 0, C; 93 for (int i = 1, t1, t2, t3, r1, r2, r3, p, pr; i <= top; i++) { 94 p = fac[i], pr = val[i]; 95 r1 = getPos(n, p), r2 = getPos(m, p), r3 = getPos(n - m, p); 96 t1 = calc(n, p, pr), t2 = calc(m, p, pr), t3 = calc(n - m, p, pr); 97 C = ((t1 * inv(t2, pr) % pr) * inv(t3, pr) % pr * 1ll * qpow(p, r1 - r2 - r3, pr)) % pr; 98 rt = (rt + (((ap / pr) * inv(ap / pr, pr) % ap) * C % ap)) % ap; 99 } 100 return rt; 101 } 102 103 inline void solve() { 104 getFactors(p); 105 ll C = exLucas(n, k, p); 106 for (int i = 1; i <= m; i++) { 107 C = C * exLucas(k, w[i], p) % p; 108 k -= w[i]; 109 } 110 printf(Auto, C); 111 } 112 113 int main() { 114 init(); 115 solve(); 116 return 0; 117 }