题解 P4735 【【模板】扩展卢卡斯】

· · 题解

扩展卢卡斯定理。

首先，得确定先会扩展欧几里得算法和扩展中国剩余定理。至于卢卡斯定理，那并不重要。

设p=\prod p_i^{k_i}，则:

假如你已经求出了\binom{n}{m}\bmod{p_i^{k_i}}，那么明显是可以利用CRT来合并答案的。

那么问题转换为如何求出\binom{n}{m}\bmod{p^k}(p是质数)。

可以知道:

\binom{n}{m}=\frac{n!}{m!(n-m)!}

那么问题即转化为求出几个阶乘和阶乘的逆元。

现在，问题归为如何快速求出阶乘。

为了便于统计出现了多少个p的次幂，先将阶乘中所有p的倍数提出来。可以简单算出，一共有\displaystyle\lfloor\frac{n}{p}\rfloor个。这中间每一项都除去p，可以得到\displaystyle\lfloor\frac{n}{p}\rfloor!。该部分可以选择递归求解。

那么接下来只剩下非p的倍数的几项了。通过简单观察可以知道，剩余几项具有循环节，循环节长度小于p^k。原因是剩余项关于p具有循环节，而a+p^k\equiv a\pmod{p^k}，所以可以一起计算。结果会剩下几项凑不齐一个循环节，但是这几项长度已经小于一个循环节了，可以选择暴力求解。

为了更好地理解上方几条，可以举个栗子:

当n=19,p=3,k=2时:

n!=1*2*3*\cdots*19 =(1*2*4*5*7*8*10*11*13*14*16*17*19)*3^{6}*6! \equiv(1*2*4*5*7*8)^2*19*3^6*6!

这样就可以在可承受的时间复杂度内求出阶乘。但是，为了达到除法的效果，我们需要考虑p的次幂一共出现了多少次。根据前面的计算，可以知道只除去一个p时，n!内包括了\displaystyle\lfloor\frac{n}{p}\rfloor个p!。剩下的数字如果要可能存在p的次幂也可以归为一个阶乘，即\displaystyle\lfloor\frac{n}{p}\rfloor!。设f(n)表示n!中有多少个p的因数，那么，我们就可以得到一个递推式:

f(n)=f(\displaystyle\lfloor\frac{n}{p}\rfloor)+\displaystyle\lfloor\frac{n}{p}\rfloor

边界为:

f(x)=0(x<p)

开始计算时间复杂度。

首先是中国剩余定理的复杂度，可以简单查到是O(plogP)的。(其中P是模数，p是最大质因子)。

然后是求阶乘复杂度。当求n!\bmod p^k时，时间复杂度为O(p^klogn(log_pn-k))(玄学)

最后是求逆元的复杂度。利用扩展欧几里得可以做到O(logp)(可忽略)。

因此，总复杂度为

O(\sum p^klogn(log_pn-k)+plogP)

这个复杂度和PlogP同级。

ext:

如果要同时计算多个组合数膜同一个合数时，可以选择对每个质因子 p ，对每个 n\le p^k 预处理出小于等于它且不整除 p 的数字之积。这样的话，求阶乘的复杂度会变成 O(\log_{p^k}n) ，总复杂度就可以降到

O(\sum p^k +T\log_{p^k} n)

代码(ext部分没有贴代码):

#include<bits/stdc++.h>

using namespace std;

typedef long long ll;

ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(!b){x=1;y=0;return a;}
    ll res=exgcd(b,a%b,x,y),t;
    t=x;x=y;y=t-a/b*y;
    return res;
}

ll p;

inline ll power(ll a,ll b,ll mod)
{
    ll sm;
    for(sm=1;b;b>>=1,a=a*a%mod)if(b&1)
        sm=sm*a%mod;
    return sm;
}

ll fac(ll n,ll pi,ll pk)
{
    if(!n)return 1;
    ll res=1;
    for(register ll i=2;i<=pk;++i)
        if(i%pi)(res*=i)%=pk;
    res=power(res,n/pk,pk);
    for(register ll i=2;i<=n%pk;++i)
        if(i%pi)(res*=i)%=pk;
    return res*fac(n/pi,pi,pk)%pk;
}

inline ll inv(ll n,ll mod)
{
    ll x,y;
    exgcd(n,mod,x,y);
    return (x+=mod)>mod?x-mod:x;
}

inline ll CRT(ll b,ll mod){return b*inv(p/mod,mod)%p*(p/mod)%p;}

const int MAXN=11;

static ll n,m;

static ll w[MAXN];

inline ll C(ll n,ll m,ll pi,ll pk)
{
    ll up=fac(n,pi,pk),d1=fac(m,pi,pk),d2=fac(n-m,pi,pk);
    ll k=0;
    for(register ll i=n;i;i/=pi)k+=i/pi;
    for(register ll i=m;i;i/=pi)k-=i/pi;
    for(register ll i=n-m;i;i/=pi)k-=i/pi;
    return up*inv(d1,pk)%pk*inv(d2,pk)%pk*power(pi,k,pk)%pk;
}

inline ll exlucus(ll n,ll m)
{
    ll res=0,tmp=p,pk;
    static int lim=sqrt(p)+5;
    for(register int i=2;i<=lim;++i)if(tmp%i==0)
    {
        pk=1;while(tmp%i==0)pk*=i,tmp/=i;
        (res+=CRT(C(n,m,i,pk),pk))%=p;
    }
    if(tmp>1)(res+=CRT(C(n,m,tmp,tmp),tmp))%=p;
    return res;
}

int main()
{
    scanf("%lld%lld%d",&n,&m,&p);
    printf("%d\n",exlucus(n,m));
    return 0;
}