扩展卢卡斯
-------------------------------------------------------------------------------------------
这是蒟蒻对扩展卢卡斯的一些见解
如有错误欢迎指出，不胜感激
普通卢卡斯
快速求出$ 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; }
