Lucas/exLucas 略解

Lucas 是一種用來求解 C n m m o d    p C_n^m \mod p Cnm​modp（ p p p 為質數）的算法。

利用 exLucas 可以求解 p p p 不為質數的情況。

Lucas

前置知識：

• 求乘法逆元
• 暴力求組合數

這種算法是基于 Lucas 定理的：

設 p p p 為素數，将正整數 n , m n,m n,m 表示為 p p p 進制如下：

n = ( n 0 n 1 … n k ‾ ) p , m = ( m 0 m 1 … m k ‾ ) p n=(\overline{n_0n_1\dots n_k})_p,m=(\overline{m_0m_1\dots m_k})_p n=(n0​n1​…nk​​)p​,m=(m0​m1​…mk​​)p​

其中允許首位為零，且規定當 m > n m>n m>n 時 C n m = 0 C_n^m=0 Cnm​=0，則：

C n m ≡ C n 0 m 0 ⋅ C n 1 m 1 ⋅ ⋯ ⋅ C n k m k ( m o d    p   ) C_n^m\equiv C_{n_0}^{m_0}\cdot C_{n_1}^{m_1}\cdot \dots \cdot C_{n_k}^{m_k} (\mod p\ ) Cnm​≡Cn0​m0​​⋅Cn1​m1​​⋅⋯⋅Cnk​mk​​(modp )

也可以換個遞歸的表達方式：

C n m ≡ C n / p m / p ⋅ C n   m o d   p m   m o d   p ( m o d p ) C_n^m\equiv C_{n/p}^{m/p}\cdot C_{n\bmod p}^{m\bmod p}\pmod p Cnm​≡Cn/pm/p​⋅Cnmodpmmodp​(modp)

于是對于 n , m < p n,m< p n,m<p 的情況，我們采用暴力求解；否則我們把它拆成較小的問題。

模闆題：洛谷 P3807 盧卡斯定理

#include<bits/stdc++.h>
using namespace std;
#define ll long long
ll getint(){
    ll ans=0,f=1;
    char c=getchar();
    while(c<'0'||c>'9'){
        if(c=='-')f=-1;
        c=getchar();
    }
    while(c>='0'&&c<='9'){
        ans=ans*10+c-'0';
        c=getchar();
    }
    return ans*f;
}
int inf=1926081700;
pair<ll,ll>exgcd(ll a,ll b){
    if(b==0){
        if(a==1)return make_pair(1,0);
        if(a==-1)return make_pair(-1,0);
        return make_pair(-inf,-inf);
    }
    pair<int,int>ans=exgcd(b,a%b);
    return make_pair(ans.second,ans.first-a/b*ans.second);
}
ll inv(ll x,ll mod){
    return (exgcd(x,mod).first%mod+mod)%mod;
}
ll qpow(ll x,ll y,ll p){
    ll ans=1;
    while(y){
        if(y&1)ans=ans*x%p;
        x=x*x%p;
        y>>=1;
    }
    return ans;
}
ll c(ll x,ll y,ll p){
    if(x<y)return 0;
    if(x==y)return 1;
    if(y>x-y)y=x-y;
    ll a=1,b=1;
    for(int i=0;i<y;i++)a=a*(x-i)%p,b=b*(y-i)%p;
    return a*inv(b,p);
}
ll luc(ll x,ll y,ll p){
    if(x<y)return 0;
    if(x==y)return 1;
    if(y>x-y)y=x-y;
    ll ans=1;
    while(x&&y){
        ans=ans*c(x%p,y%p,p)%p;
        x/=p,y/=p;
    }
    return ans;
}
int main(){
    int t=getint();
    while(t--){
        int x=getint(),y=getint(),z=getint();
        printf("%lld\n",luc(x+y,y,z));
    }
    return 0;
}
           

exLucas

前置知識：

• CRT（中國剩餘定理）
• 質因數分解（相信學到 Lucas 的沒人不會吧）
• 快速幂
• 求乘法逆元

當 p p p 不為質數時，我們采用拓展 Lucas（exLucas）

首先我們把 p p p 分解為 p 1 a 1 p 2 a 2 … p_1^{a_1}p_2^{a_2}\dots p1a1​​p2a2​​…。我們隻需要求出 C n m m o d    p i a i C_n^m\mod p_i^{a_i} Cnm​modpiai​​，用 CRT 就能求出答案。

這次我們用 C n m = n ! m ! ( n − m ) ! C_n^m=\frac{n!}{m!(n-m)!} Cnm​=m!(n−m)!n!​ 計算組合數。注意到這幾個階乘與 p p p 可能不互質，于是把所有 p k p^k pk 單獨提出來。

C n m = n ! m ! ( n − m ) ! = n ! p k 1 m ! p k 2 ⋅ ( n − m ) ! p k 3 ⋅ p k 1 − k 2 − k 3 C_n^m=\frac{n!}{m!(n-m)!}=\frac{\frac{n!}{p^{k_1}}}{\frac{m!}{p^{k_2}}\cdot\frac{(n-m)!}{p^{k_3}}}\cdot p^{k_1-k_2-k_3} Cnm​=m!(n−m)!n!​=pk2​m!​⋅pk3​(n−m)!​pk1​n!​​⋅pk1​−k2​−k3​

求 k k k 很容易，考慮 1 到 p 中有多少數被 p 1 , p 2 … p^1,p^2\dots p1,p2… 整除即可。

也就是說我們隻需要解決 n ! p k m o d    p a \frac{n!}{p^k}\mod p^a pkn!​modpa。

我們将 1 到 n n n 的數分成兩部分：是 p p p 的倍數的數與不是 p p p 的倍數的數。

p p p 的倍數的乘積可以遞歸求解：每個數除以 p p p，所求的變成 ⌊ n p ⌋ ! \lfloor\frac{n}{p}\rfloor! ⌊pn​⌋!

剩下的數字   m o d   p \bmod p modp 顯然是 p a p^a pa 個一循環的，我們把每個循環 %p 暴力求出來，然後使用快速幂。當然這些數字在循環之外可能有剩餘，暴力算就行了。

這樣求出了 n ! p k m o d    p a \frac{n!}{p^k}\mod p^a pkn!​modpa，就可以求 C n m m o d    p i a i C_n^m\mod p_i^{a_i} Cnm​modpiai​​，在拿中國剩餘定理算一下，答案就出來了。

模闆題：洛谷 P4720 【模闆】擴充盧卡斯

#include<bits/stdc++.h>
using namespace std;
#define ll long long
inline ll getint(){
    ll ans=0,f=1;
    char c=getchar();
    while(c<'0'||c>'9'){
        if(c=='-')f=-1;
        c=getchar();
    }
    while(c>='0'&&c<='9'){
        ans=ans*10+c-'0';
        c=getchar();
    }
    return ans*f;
}

int inf=1926081700;
inline pair<ll,ll>exgcd(ll a,ll b){
    if(b==0){
        if(a==1)return make_pair(1,0);
        if(a==-1)return make_pair(-1,0);
        return make_pair(-inf,-inf);
    }
    pair<ll,ll>ans=exgcd(b,a%b);
    return make_pair(ans.second,ans.first-a/b*ans.second);
}
inline ll inv(ll x,ll mod){
    return (exgcd(x,mod).first%mod+mod)%mod;
}

inline ll qmulti(ll a,ll b,ll mod){
    a%=mod,b%=mod;
    if(a<2e9&&b<2e9)return a*b%mod;
    ll ans=0;
    while(b){
        if(b&1)ans+=a,ans%=mod;
        a+=a,a%=mod;
        b>>=1;
    }
    return ans;
}

inline ll qpow(ll x,ll y,ll p){
    ll ans=1;
    while(y){
        if(y&1)ans=qmulti(ans,x,p);
        x=qmulti(x,x,p);
        y>>=1;
    }
    return ans;
}
inline ll qfac(ll x,ll p,ll pk){
    //(x!)/(p^t)%(p^k)
    if(!x)return 1;
    int ans=1;
    for(int i=1;i<pk;i++)if(i%p)ans=qmulti(ans,i,pk);
    ans=qpow(ans,x/pk,pk);
    for(int i=1;i<=x%pk;i++)if(i%p)ans=qmulti(ans,i,pk);
    return qmulti(ans,qfac(x/p,p,pk),pk);
}
inline ll countp(ll x,ll p){
    //(x!)至多有多少個 p
    int cnt=0;
    while(x)x/=p,cnt+=x;
    return cnt;
}
inline ll c(ll x,ll y,ll p,ll pk){
    if(x==y)return 1;
    if(x<y)return 0;
    return qmulti(qmulti(qmulti(qfac(x,p,pk),inv(qfac(y,p,pk),pk),pk),inv(qfac(x-y,p,pk),pk),pk),
           qpow(p,countp(x,p)-countp(y,p)-countp(x-y,p),pk),pk);
}
ll a[1000],b[1000],p[1000];
int main(){
    ll N=getint(),M=getint(),mod=getint();
    ll n=0;
    ll mul=mod,ans=0;
    for(ll i=2;i*i<=mod;i++){
        if(mod%i)continue;
        b[n]=1;p[n]=i;
        while(mod%i==0)mod/=i,b[n]*=i;
        //cout<<"> "<<i<<' '<<b[n]<<endl;
        a[n]=c(N,M,i,b[n]);
        //cout<<"! "<<a[n]<<endl;
        ++n;
    }
    if(mod!=1){
        b[n]=mod;p[n]=mod;
        a[n]=c(N,M,mod,mod);
        ++n;
    }
    //cout<<"."<<endl;
    for(int i=0;i<n;i++)ans+=qmulti(a[i],qmulti(mul/b[i],inv(mul/b[i],b[i]),mul),mul),ans%=mul;
    cout<<ans;
}