卢卡斯定理 学习笔记

卢卡斯定理

概述

卢卡斯定理是这样一个式子：\(C_n^m\equiv C_{n\bmod p}^{m\bmod p}C_{\lfloor \frac{n}{p}\rfloor}^{\lfloor \frac{m}{p}\rfloor}\pmod{p}\)，其中 \(p\) 是质数。

通过这个式子可以快速求组合数。对组合数使用卢卡斯定理后，\(C_{n\bmod p}^{m\bmod p}\) 小于 \(p\)，而对于 \(C_{\lfloor \frac{n}{p}\rfloor}^{\lfloor \frac{m}{p}\rfloor}\pmod{p}\) 递归求解。预处理 \(1\sim p\) 的阶乘与逆元可以 \(O(p)-O(\log_p m)\) 求组合数。

当 \(p\) 可能小于 \(m\) 时，处理阶乘更劣，并且逆元可能不存在，这时需要用卢卡斯定理求组合数。

证明

当 \(0<n<p\)，\(C_p^n=\frac{p!}{n!(p-n)!}\)。因为 \(p\) 为质数，有 \(p\nmid n!(p-n)!\)。所以 \(C_p^n\bmod p=0\)。故 \(C_p^n \bmod p=[n=0\vee n=p]\)。

\[\begin{aligned}(a+b)^p&\equiv\sum_{i=0}^{p}C_p^i a^{n-i}b^i\\&\equiv\sum_{i=0}^{p}[i=0\vee i=p] a^{n-i}b^i\\&\equiv a^p+b^p\pmod{p}\end{aligned} \]

则 \((1+x)^p\equiv 1+x^p\pmod{p}\)。

\[\begin{aligned}(1+x)^n&\equiv (1+x)^{p\lfloor\frac{n}{p}\rfloor}(1+x)^{n\bmod p}\\&\equiv(1+x^p)^{\lfloor\frac{n}{p}\rfloor}(1+x)^{n\bmod p}\pmod{p}\end{aligned} \]

\(C_n^m\bmod p\) 是 \((1+x)^n\) 的 \(x^m\) 次项系数。观察式子发现 \((1+x^p)^{\lfloor\frac{n}{p}\rfloor}\) 只贡献 \(p\) 的倍数，\((1+x)^{n\bmod p}\) 只贡献小于 \(p\) 的数。因此 \(x^m\) 项应该是左边的 \((x^p)^{\lfloor\frac{m}{p}\rfloor}\) 乘右边的 \(x^{m\bmod p}\)。系数就是 \(C_{n\bmod p}^{m\bmod p}C_{\lfloor \frac{n}{p}\rfloor}^{\lfloor \frac{m}{p}\rfloor}\bmod{p}\)，证毕。

例题

卢卡斯定理也可用于数学推导。

P3773

对于 \(C_n^m\bmod 2\)，不断使用卢卡斯定理，发现这就是一个二进制拆位的过程。这个式子就等于 \(n,m\) 的每一位求组合数的积。

\(C_n^m\bmod 2>0\) 相当于 \(n,m\) 的每一位求组合数都不为 \(0\)。\(C_1^0=1,C_1^1=1,C_0^0=1,C_0^1=0\)。因此 \(m\) 的一位为 \(1\) 时，\(n\) 的这一位也应为 \(1\)。\(m\) 的二进制表示是 \(n\) 的子集。

问题就变成求多少个不上升的子序列满足后一项是前一项的子集，可以 DP。设 \(f_i\) 表示以数字 \(i\) 结尾的合法子序列的个数。开一个桶存数字出现的位置，枚举子集转移。

#include<bits/stdc++.h>
using namespace std;
const int mod=1000000007;
int n,a[262144],pos[262144],f[262144],ans;
int main(){
  cin>>n;
  for(int i=1;i<=n;i++)cin>>a[i],pos[a[i]]=i;
  for(int i=n;i>=1;i--){
    f[i]=1;
    for(int j=a[i]&(a[i]-1);j;j=(j-1)&a[i])if(pos[j]>i)f[i]=(f[pos[j]]+f[i])%mod;
    ans=(ans+f[i])%mod;
  }
  cout<<ans<<'\n';
  return 0;
}

扩展卢卡斯定理

一种当模数不是质数时求组合数的方法。

首先将模数 \(p\) 分解。设 \(p=\prod p_i^{k_i}\)， 求单个 \(C_n^m\bmod p^k\)，用中国剩余定理合并答案。

\(C_n^m\bmod p^k=\frac{n!}{m!(n-m)!}\bmod p^k\)。此时分子的逆元不一定存在，需要提出 \(p\)，等于 \(\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}p^{x-y-z}\bmod p^k\)，其中 \(x,y,z\) 分别为 \(n!,m!,(n-m)!\) 的因数 \(p\) 个数。

那么求出单个 \(\frac{n!}{p^x}\bmod p^k\) 即可。

先想 \(n!\bmod p^k\)。把 \(n!\) 中 \(p\) 的倍数提出来，比如 \(n=22,p=3\)。

\[\begin{aligned}22!\bmod3^2&=(3\times6\times9\times12\times15\times18\times21)\times(1\times2\times4\times5\times7\times8\times10\times11\times13\times14\times16\times17\times19\times 20\times22)\bmod 3^2\\&=3^7\times7!\times(1\times2\times4\times5\times7\times8)\times(10\times11\times13\times14\times16\times17)\times(19\times20\times22)\end{aligned} \]

此时式子分为三个部分。第一部分是 \(p^{\lfloor\frac{n}{p}\rfloor}\)，第二部分是 \(\lfloor\frac{n}{p}\rfloor!\)。第三部分是一个循环的乘积，即 \(\prod_{i,p\nmid i}^{n}i\)。即：

\[n!\equiv p^{\left\lfloor\frac{n}{p}\right\rfloor}\cdot\left(\left\lfloor\frac{n}{p}\right\rfloor\right)!\cdot{\left(\prod_{i,p\nmid i}^{p^k}i\right)}^{\left\lfloor\frac{n}{p^k}\right\rfloor}\cdot\left(\prod_{i,p\nmid i}^{n\bmod p^k}i\right)\pmod{p^k} \]

第一部分直接快速幂，第二个部分也是阶乘，递归求解。第三个部分先暴力求单个循环节 \(\prod_{i,p\nmid i}^{p^k}i\)，剩下的循环节取模后相同，取 \(\lfloor\frac{n}{p^k}\rfloor\) 次方，暴力乘剩余部分。

对两边同除 \(p^x\)。因为 \(p^x\) 表示 \(n!\) 所有的因数 \(p\)，要去掉前两项中的因数 \(p\)。第一项直接没了，第二项变成 \(\frac{\left(\left\lfloor\frac{n}{p}\right\rfloor\right)!}{p^{x'}}\)。其中 \(p^{x'}\) 表示 \(\left(\left\lfloor\frac{n}{p}\right\rfloor\right)!\) 中的因数 \(p\)。这个式子的形式与 \(\frac{n!}{p^x}\) 一样，可以递归。

template<typename T>T f(T n,T q,T qk,T r=1){
  if(!n)return 1;
  for(T i=1;i<=qk;i++)if(i%q)r=r*(i%qk)%qk;
  r=qpow(r,n/qk,qk);
  for(T i=n/qk*qk+1;i<=n;i++)if(i%q)r=r*(i%qk)%qk;
  return r*f(n/q,q,qk)%qk;
}
template<typename T>T getb(T n,T m,T q,T qk,T cnt=0){
  for(T i=n;i;i/=q)cnt+=i/q;
  for(T i=m;i;i/=q)cnt-=i/q;
  for(T i=n-m;i;i/=q)cnt-=i/q;
  return qpow(q,cnt,qk)*f(n,q,qk)%qk*inv(f(m,q,qk),qk)%qk*inv(f(n-m,q,qk),qk)%qk;
}
template<typename T>T C(T n,T m,T p,T ans=0){
  vector<T>a,b,c;
  T temp=p;
  for(T i=2;i*i<=p;i++){
    if(p%i==0){
      T t=1;
      while(p%i==0)p/=i,t*=i;
      a.push_back(t),c.push_back(i);
    }
  }
  if(p>1)a.push_back(p),c.push_back(p);
  for(int i=0;i<a.size();i++)b.push_back(getb(n,m,c[i],a[i]));
  for(int i=0;i<a.size();i++)ans=(ans+b[i]*(temp/a[i])%temp*inv(temp/a[i],a[i])%temp)%temp;
  return ans;
}

[[数学]]