Lucas定理和扩展Lucas定理
1.Lucas定理
首先给出式子：\(C_n^m\%p = C_{\lfloor\frac{n}{p}\rfloor}^{\lfloor\frac{m}{p}\rfloor} * C_{n\%p}^{m\%p}\% p\),其中p为质数。
这里给出证明……证明是我在luogu上看到的lance1ot大佬的证明，个人认为是写的很好的，在此还要做一下补充。
首先，对于质数p，可以保证\(C_p^i(1 <= i <= p-1) \equiv 0(mod\ p)\)，这个比较显然，因为组合数一定是整数，而质数p的因子只有自己和1，也就是说并没有某个数能整除它，所以答案必然是p的倍数。
根据二项式定理，\((1+x)^p = C_p^0x^0 + C_p^1x^1 + C_p^2x^2+…C_p^px^p\),结合上面的结论可以知道\((1+x)^p \equiv 1 + x^p (mod\ p)\)
之后我们要证明lucas定理 ，我们假设\(a = kp + j,b = sp + g\),那我们只要证明\(C_a^b = C_k^s * C_j^g\%p\)即可。
继续从二项式定理入手。\((1+x)^a \equiv (1+x)^{kp+j} \equiv (1+x)^{kp} * (1+x)^j \equiv (1+x^p)^k * (1+x)^j(mod\ p)\)
这时候我们观察二项式展开后第b+1项，就是\(C_a^bx^b\)，这项显然也是可以由\((1+x^p)^k\)中的一项和 \((1+x)^j\)中的一项相乘得到的。而且是唯一的两项，就是\(C_k^sx^{sp}\)和\(C_j^gx^g\) ,因为\(b = sp+g\),那么对于\((1+x^p)^k\)的展开式，显然x的指数是\((s+1)p\)或者更高次的指数比b大，而\((s-1)p\)或者更低次的，后面的指数会不够用，也无法匹配成b。所以\(C_a^bx^b \equiv C_k^sx^{sp}*C_j^gx^g(mod\ p)\),那么\(C_a^b \equiv C_k^sx^{sp}*C_j^gx^g(mod\ p)\),即\(C_n^m\%p = C_{\lfloor\frac{n}{p}\rfloor}^{\lfloor\frac{m}{p}\rfloor} * C_{n\%p}^{m\%p}\% p\)。
一开始的前提条件要求了p必须是质数。否则的话\(C_p^i(1 <= i <= p-1) \equiv 0(mod\ p)\)无法保证，因为p很有可能被自己的因子筛掉了。比如\(C_6^4\%6\)，答案就是3而不是0.
2.扩展Lucas定理。
当我们遇到p不是质数的时候怎么办呢……
如果p能分解成几个质数的乘积，而且这些质数的指数都是1的话，可以直接套用lucas然后用CRT合并。比如SDOI2010古代猪文
不过如果它可以分解成质数的k次幂的乘积，这样就不行了。
于是我们有了扩展Lucas。
首先我们把p唯一分解，假设我们现在用\(p_i^k\)做模数，把所有的计算出来以后还是可以用CRT合并的。
因为\(C_n^m = \frac{n!}{m!(n-m)!}\),所以问题变成了如何快速在模意义下算出阶乘。
大致的方法就是，首先我们先提取出所有p的倍数，对于n，其阶乘内部有\(\lfloor\frac{n}{p}\rfloor\)个p的倍数，把他们全部提取出来，结果就是\(p^{\lfloor\frac{n}{p}\rfloor} * \lfloor\frac{n}{p}\rfloor!\)，其中\(\lfloor\frac{n}{p}\rfloor!\)就可以递归计算。
对于不是p的倍数的，每个\(p^k\)成一个循环节，在每个循环节里面是直接把乘积算出来，最后再套上\(\lfloor\frac{n}{p^k}\rfloor\)的指数。最后会剩余几项，那些直接暴力乘起来就行。
举个例子。(luogu上的例子)
假设\(n = 19,p = 3,k = 2\),首先我们先把p的倍数提取出来，就变成\(1 * 2 * 4 * 5 * 7 * 8 * 10 * 11 * 13 * 14 * 16 * 17 * 19 * 3^6 * 6!\)
之后可以看出\(1 \equiv 10(mod\ 9)\),\(2 \equiv 11 (mod\ 9)\),所以原式就变成\((1 * 2 * 4 * 5 * 7 * 8)^2 * 19 * 3^6 * 6!\)
最后那个单个的19就是不在循环节里面的，但是\(1 \equiv 19(mod\ 9)\)嘛，所以暴力计算就好了啦。
这样递归下去计算就可以。代码实现中略微有些不同，就是对于里面每一次计算次方的\(p^{\lfloor\frac{n}{p}\rfloor}\),我们不在递归的时候计算，而是全部提出来，先上下消去，最后再做乘方。这个看代码实现就好。
这里有一道板子题国家集训队 礼物,答案显然是\(C_n^{w[1]}*C_{n-w[1]}^{w[2]}…\),直接套扩展lucas即可。
#include<bits/stdc++.h>
#define rep(i,a,n) for(ll i = a;i <= n;i++)
#define per(i,n,a) for(ll i = n;i >= a;i--)
#define enter putchar('\n')
using namespace std;
typedef long long ll;
const int M = 1000005;
const int INF = 1000000009;
const double eps = 1e-8;
ll read()
{
ll ans = 0,op = 1;char ch = getchar();
while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
return ans * op;
}
ll n,m,p,w[105],sum,ans = 1;
ll mul(ll a,ll b,ll t)
{
ll res = a * b - (ll)((long double)a / t * b + eps) * t;
return (res % t + t) % t;
}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
if(!b){x = 1,y = 0;return a;}
ll d = exgcd(b,a%b,y,x);
y -= a / b * x;
return d;
}
ll inv(ll a,ll b)
{
ll x,y;
exgcd(a,b,x,y);
return (x % b + b) % b;
}
ll CRT(ll b,ll t) {return mul(mul(b,inv(p/t,t),p),p/t,p);}
ll qpow(ll a,ll b,ll t)
{
ll p = 1;
while(b)
{
if(b & 1) p = mul(p,a,t);
a = mul(a,a,t),b >>= 1;
}
return p;
}
ll fac(ll n,ll pi,ll pk)
{
if(!n) return 1;
ll res = 1;
rep(i,2,pk) if(i % pi) res *= i,res %= pk;
res = qpow(res,n/pk,pk);
rep(i,2,n%pk) if(i % pi) res *= i,res %= pk;
return res * fac(n / pi,pi,pk) % pk;
}
ll C(ll n,ll m,ll pi,ll pk)
{
ll d = fac(n,pi,pk),d1 = fac(m,pi,pk),d2 = fac(n-m,pi,pk);
ll k = 0;
for(ll i = n;i;i /= pi) k += i / pi;
for(ll i = m;i;i /= pi) k -= i / pi;
for(ll i = n-m;i;i /= pi) k -= i / pi;
return mul(mul(d,inv(d1,pk),pk),mul(qpow(pi,k,pk),inv(d2,pk),pk),pk);
}
ll exlucas(ll n,ll m)
{
ll res = 0,tmp = p,pk;
ll lim = sqrt(p);
rep(i,2,lim) if(!(tmp % i))
{
pk = 1;
while(!(tmp%i)) pk *= i,tmp /= i;
res += CRT(C(n,m,i,pk),pk),res %= p;
}
if(tmp > 1) res += CRT(C(n,m,tmp,tmp),tmp),res %= p;
return res;
}
int main()
{
p = read();
n = read(),m = read();
rep(i,1,m) w[i] = read(),sum += w[i];
if(sum > n) printf("Impossible\n"),exit(0);
rep(i,1,m) ans *= exlucas(n,w[i]),ans %= p,n -= w[i];
printf("%lld\n",ans);
return 0;
}
