擴充套件盧卡斯定理
再看下面講解之前,請大家熟悉一下乘法逆元,中國剩餘定理,拓展gcd,費馬小定理 關於拓展盧卡斯,也就是盧卡斯的升級版,盧卡斯限定只能取餘一個素數,而拓展盧卡斯則沒有限定
若不是素數,將p分解質因數,將C(n,m)分別按照Lucas定理中的方法求對p的質因數的模,然後用中國剩餘定理合併。比如計算C(10,3)%14。C(10,3)=120,14有兩個質因數2和7,120%2=0,120%7=1,這樣用(2,0)(7,1)找到最小的正整數8即是答案,即C(10,3)%14=8。注意,這裡只適用於p分解完質因數後每個質因數只出現一次,例如12=2*2*3就不行,因為2出現了兩次。若p分解完質因數後,含有某個質因數出現多次,比如C(10,3)%98,其中98=2*7*7,此時就要把7*7看做一個數,即:120%2=0,120%49=22,用(2,0)(49,22)和中國剩餘定理得到答案22,即C(10,3)%98=22。此時,你又會有疑問,C(10,3)%49不也是模一個非素數嗎?此時不同的是這個非素數不是一般的非素數,而是某個素數的某次方。
如何計算C(n,m)%p^t(t>=2,p為素數)。 計算C(n,m)%p^t。我們知道,C(n,m)=n!/m!/(n-m)!,若我們可以計算出n!%p^t,我們就能計算出m!%p^t以及(n-m)!%p^t。我們不妨設x=n!%p^t,y=m!%p^t,z=(n-m)!%p^t,那麼答案就是x*reverse(y,p^t)*reverse(z,p^t)(reverse(a,b)計算a對b的乘法逆元)。那麼下面問題就轉化成如何計算n!%p^t。比如p=3,t=2,n=19, n!=1*2*3*4*5*6*7*8* ……*19 =[1*2*4*5*7*8*… 16*17*19](3*6*9*12*15*18) =[1*2*4*5*7*8*… *16*17*19]*3^6(1*2*3*4*5*6) 然後發現後面的是(n/p)!,於是遞迴即可。前半部分是以p^t為週期的[1*2*4*5*7*8]=[10*11*13*14*16*17](mod 9)。下面是孤立的19,可以知道孤立出來的長度不超過 p^t,於是暴力即可。那麼最後剩下的3^6啊這些數怎麼辦呢?我們只要計算出n!,m!,(n-m)!裡含有多少個p(不妨設a,b,c),那麼a-b-c就是C(n,m)中p的個數,直接算一下就行。 至此整個計算C(n,m)%p(p為任意數)的問題完美解決。
程式碼:
#include <iostream> #include <stdio.h> #include <algorithm> #include<cstring> using namespace std; typedef long long ll; ll quick_mod(ll a,ll b,ll m){ ll ans=1ll; while(b){ if(b&1) ans=ans*a%m; b>>=1; a=a*a%m; } return ans; } ll exgcd(ll a,ll b,ll &x,ll &y){ if(a%b==0){ x=0ll;y=1ll; return b; } ll v,tx,ty; v=exgcd(b,a%b,tx,ty); x=ty; y=tx-a/b*ty; return v; } ll inv(ll a,ll p){ if(!a) return 0ll; ll x,y; exgcd(a,p,x,y); x=(x%p+p)%p; return x; } ll Mul(ll n,ll pi,ll pk){ if(!n) return 1ll; ll ans=1ll; for(ll i=2;i<=pk;i++) if(i%pi) ans=ans*i%pk; ans=quick_mod(ans,n/pk,pk); for(ll i=2;i<=n%pk;i++){ if(i%pi) ans=ans*i%pk; } return ans*Mul(n/pi,pi,pk)%pk; } ll exlucas(ll n,ll m,ll p,ll pi,ll pk){ if(m>n) return 0ll; ll a=Mul(n,pi,pk),b=Mul(m,pi,pk),c=Mul(n-m,pi,pk); ll k=0ll,ans=0ll; for(ll i=n;i;i/=pi) k+=i/pi; for(ll i=m;i;i/=pi) k-=i/pi; for(ll i=n-m;i;i/=pi) k-=i/pi; ans=a*inv(b,pk)%pk*inv(c,pk)%pk*quick_mod(pi,k,pk)%pk; return ans*(p/pk)%p*inv(p/pk,pk)%p; //中國剩餘定理 a[i]*M*x 餘數*其他個個素數的乘積*x } int main() { ll n,m,p,ans=0; while(cin>>n>>m>>p){ for(ll x=p,i=2;i<=p;i++){ if(x%i==0){ ll pk=1ll; while(x%i==0) pk*=i,x/=i; ans=(ans+exlucas(n,m,p,i,pk))%p; } } cout<<ans<<endl; ans=0; } return 0; }