Lucas定理和扩展Lucas定理

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;
}

原文地址:https://www.cnblogs.com/captain1/p/10382238.html