擴展盧卡斯

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

這是蒟蒻對擴展盧卡斯的一些看法
若有錯誤歡迎指出，不勝感激

普通盧卡斯                                                                                                                      

快速求出$ 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;
}