卢卡斯定理

卢卡斯定理

[largeinom{n}{m}equivinom{n mod p}{mmod p}inom{frac{n}{p}}{frac{m}{p}} pmod{p} ]

其中 (p) 为质数。即将 (n,m) 用 (p) 进制来表示。预处理阶乘后复杂度为 (O(log_p n))。

扩展卢卡斯定理

用来解决 (p) 不是质数的情况。设 (p=prod p_i^{k_i})，得：

[largeegin{cases} xequiv inom{n}{m} &pmod{p_1^{k_1}} \ xequiv inom{n}{m} &pmod{p_2^{k_2}} \ &cdots \ xequiv inom{n}{m} &pmod{p_m^{k_m}} \ end{cases} ]

求值后用中国剩余定理合并即可。

因为 (p) 不一定为质数，所以不能直接求逆元，考虑进行转化：

[largeegin{aligned} &inom{n}{m} pmod{p^k}\ =& frac{n!}{m!(n-m)!} pmod{p^k}\ =& frac{frac{n!}{p^x}}{frac{m!}{p^y}frac{(n-m)!}{p^z}}p^{x-y-z} pmod{p^k}\ end{aligned} ]

其中 (x,y,z) 都是对应的数质因数分解后 (p) 的指数。

先考虑如何求 (p^{x-y-z})，设 (g(n)) 为 (n!) 质因数分解后 (p) 的指数，因为有：

[large n!=p^{leftlfloor frac{n}{p} ight floor}left(leftlfloor frac{n}{p} ight floor ight)!prod_{i=1,p otmid i}^n i ]

所以得：

[large g(n)=leftlfloor frac{n}{p} ight floor+gleft(leftlfloor frac{n}{p} ight floor ight) ]

再考虑如何求 (frac{n!}{p^x} mod p^k)，设 (f(n)=frac{n!}{p^x})，得：

[largeegin{aligned} &f(n) pmod{p^k}\ =&frac{p^{leftlfloor frac{n}{p} ight floor}}{p^x}left(leftlfloor frac{n}{p} ight floor ight)!prod_{i=1,p otmid i}^n i pmod{p^k}\ =&fleft(leftlfloor frac{n}{p} ight floor ight)left( prod_{i=1,p otmid i}^{p^k} i ight)^{leftlfloor frac{n}{p^k} ight floor} prod_{i=1,p otmid i}^{n mod p^k} i pmod{p^k}\ end{aligned} ]

边界为 (f(0)=1)。

预处理复杂度为 (O(sqrt p + sum p_i^{k_i}))，询问一次复杂度为 (O(log_p nlog n))。

ll f(ll n,int x)
{
    if(!n) return 1;
    ll k=pk[x];
    return f(n/pri[x],x)*qp(a[x][k],n/k,k)%k*a[x][n%k]%k;
}
ll C(ll n,ll m,int x)
{
    ll v=0,p=pri[x],k=pk[x];
    for(ll i=n;i;i/=p) v+=i/p;
    for(ll i=m;i;i/=p) v-=i/p;
    for(ll i=n-m;i;i/=p) v-=i/p;
    return f(n,x)*inv(f(m,x),k)%k*inv(f(n-m,x),k)%k*qp(p,v,k)%k;
}
ll crt(ll x,ll m,ll p)
{
    return x*(m/p)*inv(m/p,p);
}
ll exlucas(ll n,ll m)
{
    ll v=0;
    for(int i=1;i<=cnt;++i) v=(v+crt(C(n,m,i),p,pk[i]))%p;
    return v;
}
void init(int x)
{
    a[x][0]=1;
    for(int i=1;i<=pk[x];++i)
    {
        a[x][i]=a[x][i-1];
        if(i%pri[x]) a[x][i]=a[x][i]*i%pk[x];
    }
}
void pre(ll x)
{
    for(int i=2;i*i<=x;++i)
    {
        if(x%i) continue;
        ll v=1;
        while(x%i==0) x/=i,v*=i;
        pri[++cnt]=i,pk[cnt]=v;
    }
    if(x!=1) pri[++cnt]=x,pk[cnt]=x;
    for(int i=1;i<=cnt;++i) init(i);
}