1.Lucas定理
首先给出式子:(C_n^m\%p = C_{lfloorfrac{n}{p}
floor}^{lfloorfrac{m}{p}
floor} * 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_{lfloorfrac{n}{p}
floor}^{lfloorfrac{m}{p}
floor} * 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,其阶乘内部有(lfloorfrac{n}{p}
floor)个p的倍数,把他们全部提取出来,结果就是(p^{lfloorfrac{n}{p}
floor} * lfloorfrac{n}{p}
floor!),其中(lfloorfrac{n}{p}
floor!)就可以递归计算。
对于不是p的倍数的,每个(p^k)成一个循环节,在每个循环节里面是直接把乘积算出来,最后再套上(lfloorfrac{n}{p^k}
floor)的指数。最后会剩余几项,那些直接暴力乘起来就行。
举个例子。(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^{lfloorfrac{n}{p}
floor}),我们不在递归的时候计算,而是全部提出来,先上下消去,最后再做乘方。这个看代码实现就好。
这里有一道板子题国家集训队 礼物,答案显然是(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('
')
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
"),exit(0);
rep(i,1,m) ans *= exlucas(n,w[i]),ans %= p,n -= w[i];
printf("%lld
",ans);
return 0;
}