從一個例題:【HDU 3037】 Saving Beans 來開始Lucas定理的應(yīng)用刁岸。
題目大意為:松鼠要從n棵樹上摘一共m個豆子,結(jié)果的方案數(shù)對素數(shù)p(不大于1e5)取模鹃唯,求解搔耕。
思路:
可以理解為m個豆子分為n份茫因,求分的方法個數(shù)。
由插板法來對m個數(shù)進行劃分牍疏,由于可能某棵樹沒有摘豆子蠢笋,可以理解為:x1+x2+x3+……+xn=m的解的個數(shù),即為C(m+n-1,n-1)鳞陨。(將m顆豆子加上n-1個板子的位置昨寞,得到的序列再從中取n-1個板子的位置)=C(m+n-1,m)。
由于m的值取0~m厦滤,那么就得sum=C(n-1,0)+C(n,1)+C(n+1,2)+C(n+2,3)+……+C(m+n-1,m)援岩。
利用公式C(n,r)=C(n-1馁害,r)+C(n-1窄俏,r-1)=C(n-1,r)+C(n-2碘菜,r-1)+C(n-3凹蜈,r-2)……
sum=C(n+m,m)忍啸。
也就是說仰坦,接下來的算法變成了C(n+m,m)%p计雌。
然后就是Lucas定理的運用:
Lucas(m悄晃,n,p)=C(m%p凿滤,n%p妈橄,p)?Lucas(m/p,n/p翁脆,p)眷蚓。
Lucas(x,0反番,p)=1沙热。
這里可以采用的方法是遞歸求解叉钥。
簡單的理解就是:
以求解n! % p 為例,把n分段篙贸,每p個一段投队,每一段求得結(jié)果是一樣的。但是需要單獨處理每一段的末尾p,2p,...,把p提取出來爵川,會發(fā)現(xiàn)剩下的數(shù)正好又是(n/p)! 敷鸦,相當于劃歸了一個子問題,這樣遞歸求解即可雁芙。
這個是單獨處理n轧膘!的情況,當然C(n,m)就是n兔甘!/(m! *(n-m);寻),每一個階乘都用上面的方法處理的話洞焙,就是Lucas定理了.
Lucas最大的數(shù)據(jù)處理能力是p在10^5左右蟆淀。
而C(a,b) =a! / ( b! ? (a-b)! ) mod p
其實就是求 ( a! / (a-b)!) ? ( b! )^(p-2) mod p
上面這一步變換是根據(jù)費馬小定理:假如p是質(zhì)數(shù),且a,p互質(zhì)澡匪,那么a的(p-1)次方除以p的余數(shù)恒為1,
那么a和a^(p-2)互為乘法逆元熔任,則(b / a) = (b * a^(p-2) ) mod p)
b!與b!(p-2)互為乘法逆元甸鸟,即b!?b!(p-2)=1惦费,那么,
//快速冪a^b % k
ll PowerMod(ll a, ll b, ll k) {
ll tmp = a, ret = 1;
while (b) {
if (b & 1) ret = ret * tmp % k;
tmp = tmp * tmp % k;
b >>= 1;
}
return ret;
}
//求C(n, m)%p p最大為10^5 n, m可以很大抢韭!
ll Lucas(ll n, ll m, ll p) {
ll ret = 1;
while (n && m) {
ll nn = n%p, mm = m%p;
if (nn < mm) return 0;
//fac[nn]為預(yù)處理的 fac[n] = n!%p
ret = ret*fac[nn]*PowerMod(fac[mm]*fac[nn-mm]%p, p-2, p)%p;
n /= p;
m /= p;
}
return ret;
}
//C(n, m) % p
Lucas(n, m, p);
用下面的Lucas定理程序?qū)崿F(xiàn)就能得出結(jié)果薪贫,實現(xiàn)過程中要注意乘法時的強制轉(zhuǎn)換
#include <iostream>
#include <cstdio>
typedef long long lld;
lld N,M,P;
int Pow(lld a,lld n,lld p)
{
lld x = a;
lld res = 1;
while(n)
{
if(n & 1)
{
res = ((lld)res * (lld)x) % p;
}
n >>= 1;
x = ((lld)x*(lld)x) % p;
}
return res;
}
int Cm(lld n,lld m,lld p)
{
lld a = 1,b = 1;
if(m > n) return 0;
//實現(xiàn)(a!/(a-b)!) * (b!)^(p-2)) mod p,由于n比較大,所以刻恭,此處不知道有什么好的優(yōu)化
while(m)
{
a = (a * n) % p;
b = (b * m) % p;
m--;
n--;
}
return ((lld)a * (lld)Pow(b,p-2,p))%p;
}
int Lucas(lld n,lld m,lld p)
{
if(m==0)
return 1;
return((lld)Cm(n%p,m%p,p)*(lld)Lucas(n/p,m/p,p))%p;
}
int main()
{
int t;
scanf("%d",&t);
while(t--)
{
scanf("%I64d%I64d%I64d",&N,&M,&P);
printf("%d\n",Lucas(N+M,M,P));
}
return 0;
}