扩展卢卡斯

-------------------------------------------------------------------------------------------

这是蒟蒻对扩展卢卡斯的一些见解
如有错误欢迎指出，不胜感激

普通卢卡斯                                                                                                                      

快速求出$ C(n,m)$ ($ mod$ p)

约束条件：p为质数

 

考虑扩展                                                                                                                               

预备知识：中国剩余定理(可以参考我的前一篇博客)

我们可以对模数p分解质因数，使得p成为$ {p1}^{k1}{p2}^{k2}... {pn}^{kn}$的形式

列出n个同余方程：

ans≡x1 ($ mod$ $ {p1}^{k1}$)

ans≡x2 ($ mod$ $ {p2}^{k2}$)

......

 ans≡xn ($ mod$ $ {pn}^{kn}$)

由于模数两两互质，因此只要得到每一组的解x1...xn, 就可以通过中国剩余定理合并得到最终的解

 

如何得到每一组的解                                                                                                  

考虑每一组所求的式子为：

$ C(n,m)$ $ mod$ $ {pi}^{ki}$

= $\frac{n!}{m!(n-m)!}$ $ mod$ $ {pi}^{ki}$

由于模数只有一个质因数pi，因此提出先阶乘中的所有pi因子

 

如何求n！中剔除掉pi因子之后数的乘积 $ mod$ $ {pi}^{ki}$                          

下面以$ 22!$ $ mod$ $ 3^2$为例

$ (22!)=(1*2*3*4*5*6*7*8*9)*(10*11*12*13*14*15*16*17*18)*(19*20*21*22)$

剔除3因子之后

$ (22!)=(1*2*4*5*7*8)*(10*11*13*14*16*17)*(19*20*22)*3^6*(1*2*3*4*5*6*7)$

发现按照如下分组之后前$ \lfloor \frac{n}{pi^{ki}} \rfloor$ 组在$ mod$ $ {pi}^{ki}$ 后结果相同，因此只需要做一组然后快速幂即可

对于$ (19*20*22)$这是冗余，暴力计算即可

对于最后的$ (1*2*3*4*5*6*7)$递归计算即可

3因子被剔除暂不考虑

 

代码：                                                                                              

int Mul(const int n,const int pi,const int pk)//求n! % pi^ki pk=pi^ki
{
    if(n<=1)return 1;ll ans=1;
    if(n>=pk)
    {
        for(rt i=2;i<pk;i++)if(i%pi)ans=ans*i%pk;//计算一组
        ans=ksm(ans,n/pk,pk);//快速幂
    }
    for(rt i=2;i<=n%pk;i++)if(i%pi)ans=ans*i%pk;//计算冗余
    return ans*Mul(n/pi,pi,pk)%pk;//递归求最后一组
}

求出(n!) $ mod$ $ {pi}^{ki}$之后

只需要继续求出$ (m!)$ $ mod$ $ {pi}^{ki}$和$ ((n-m)!)$ $ mod$ $ {pi}^{ki}$，然后对后两项求 $ mod$ $ {pi}^{ki}$的逆元即可

显然剔除pi之后结果一定和$ {pi}^{ki}$互质，所以可以通过扩展欧几里德求逆元

最后加入被剔除的pi的贡献即可

求出pi在$ C(n,m)$中出现的次数k，然后把之前的答案乘上$ pi^k$即可

时间复杂度($ {pi}^{ki}$(求一组的时间)log(n))

注意到对于相同的模数p, 一组的答案相同

因此可以预处理前缀积，然后O(1)的返回每组的答案以及冗余答案

时间复杂度($ {pi}^{ki}$(求一组的时间)+log(n))

最后把所有答案CRT合并即可

 

总复杂度：                                                                                                                              

对于同一模数$ pi^{ki}$：O($ pi^{ki}$)预处理+log(n)求值

全程序复杂度只需要再乘上P的不同质因子个数即可

 

全代码：                                                                                                

 

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define rt register int
#define l putchar('\n')
#define ll long long
#define r read()
using namespace std;
inline ll read()
{
    register ll x = 0; char zf = 1; char ch;
    while (ch != '-' && !isdigit(ch)) ch = getchar();
    if (ch == '-') zf = -1, ch = getchar();
    while (isdigit(ch)) x = x * 10 + ch - '0', ch = getchar(); return x * zf;
}
void write(ll y)
{
    if (y < 0) putchar('-'), y = -y;
    if (y > 9) write(y / 10);
    putchar(y % 10 + '0');
}
inline void writeln(const ll y){write(y);putchar('\n');}
inline int min(const int x,const int y){return x<y?x:y;}
inline int max(const int x,const int y){return x>y?x:y;}
inline ll ksm(ll x,int y,int p)//快速幂 
{
    if(!y)return 1;int ew=1;
    while(y>1)
    {
        if(y&1)y--,ew=x*ew%p;
        else y>>=1,x=x*x%p;
    }
    return x*ew%p;
}
int i,j,k,m,n,x,y,z,cnt,p,P;
void exgcd(int a,int b,int &x,int &y)//扩展欧几里得 
{
    if(!b)x=1,y=0;
    else exgcd(b,a%b,y,x),y-=a/b*x;
}
int inv(int A,int p)//求A mod p意义下的逆元 
{
    if(!A)return 0;int x=0,y=0;
    exgcd(A,p,x,y);x=x%p+p;
    if (x>p)x-=p;return x;
}
ll qz[1000010];//预处理前缀积 
int Mul(const int n,const int pi,const int pk)//求n! % pi^ki pk=pi^ki
{
    if(n<=1)return 1;ll ans=1;
    if(n>=pk)
    {
        ans=qz[pk-1];
        ans=ksm(ans,n/pk,pk);
    }
    if(n%pk)ans=ans*qz[n%pk]%pk;
    return ans*Mul(n/pi,pi,pk)%pk;
}
int C(const int n,const int m,const int pi,const int pk)//求C(n,m)%pk  pk=pi^ki 
{
    qz[0]=qz[1]=1;
    for(rt i=2;i<pk;i++)
    {
        qz[i]=qz[i-1];
        if(i%pi)qz[i]=qz[i]*i%pk;
    }
    int a=Mul(n,pi,pk),b=inv(Mul(m,pi,pk),pk),c=inv(Mul(n-m,pi,pk),pk);
    ll ans=(ll)a*b%pk*c%pk;
    int k=0;//k记录pi出现次数 
    for(rt i=n/pi;i;i/=pi)k+=i;
    for(rt i=m/pi;i;i/=pi)k-=i;
    for(rt i=(n-m)/pi;i;i/=pi)k-=i;
    return ans*ksm(pi,k,pk)%pk;
}
ll ansc;
int main()
{
    for(rt t=read();t;t--)//t组询问，每组询问C(n,m)%p 
    {
        n=read();m=read();p=read();P=p;ansc=0;
        for(rt i=2;i<=p;i++)
        if(p%i==0)
        {
            int pi=i,pk=1;
            while(p%i==0)p/=i,pk*=i;
            (ansc+=(ll)C(n,m,pi,pk)*(P/pk)%P*inv(P/pk,pk))%=P;//CRT合并 
        }
        writeln(ansc%P);
    }
    return 0;
}