可以先做这个题[SDOI2010]古代猪文
此算法和LUCAS定理没有半毛钱关系。
不保证P是质数。
$C_n^m=frac{n!}{m!(n-m)!}$
麻烦的是分母。
如果互质就有逆元了。
所以可以考虑把分子分母不互质的数单独提出来处理。
然鹅P太一般，直接处理要考虑的东西太多。
我们不妨令$p=p_1^{q_1}*p_2^{q_2}*...*p_k^{q_k}$
对每一个$p_i^{q_i}$分别求解(不妨叫这个数为$pk$)（这样会容易很多）
即求ai满足：$frac{n!}{m!(n-m)!} = a_ispace mod space pk$
然后可以$CRT$合并
（CRT可以合并的原因是，我们可以求出满足这些同余方程的通解。发现这些解mod lcm都是同一个数x。$C_n^m$一定是这些个解之一，不管是哪一个，$modspace lcm$即mod p都是x。我们也就求出了答案）
$frac{n!}{m!(n-m)!} = a_ispace mod space pk$
现在不互质的就是pi的倍数
首先我们可以把分子分母的所有$pi$质因子都提出来，然后上下次数相消。
对于$n!$中p的质因个数，就是不断除以p^i下取整。
剩下的都是和pk互质的了。存在逆元
以求$19!space modspace 3^2$为例
$19!=1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19$
提完质数3之后，变成：
$=(1*2*3*4*5*6)*3^2*(1*2*4*5*7*8*10*11*13*14*16*17*19)$
前面直接递归下去处理。后面的每一项，都可以%pk处理。
（不论分母还是分子位置，如果是分母位置，因为存在逆元，
10*inv =1 mod 9
1*inv = 1 mod 9
这两个inv显然是同一个inv.所以把所有项%pk没有问题
）
然后变成：
$=(1*2*3*4*5*6)*3^2*(1*2*4*5*7*8)^2*1)$
后面那个1是多出来的19
其实pk长度的循环节有n/pk个。直接算。剩下（例如这里的19），个数少于pk，直接算。
（所以，扩展LUCAS的重要适用条件是，$p_i^{q_i}$不能太大(1e5左右)）
递归算出来即可。
对于分母位置的两个阶乘，算出来结果之后，再处理inv
（这里可以先乘完之后再找inv，不用一边找inv一边乘。）
（注意inv处理要用exgcd，不保证质数，不能用费马）
(CRT可以不用保存结果，Mi=p/pk 可以一次到位）
#include<bits/stdc++.h> #define reg register int #define il inline #define int long long #define numb (ch^'0') using namespace std; typedef long long ll; il void rd(int &x){ char ch;x=0;bool fl=false; while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true); for(x=numb;isdigit(ch=getchar());x=x*10+numb); (fl==true)&&(x=-x); } namespace Miracle{ ll p; ll qm(ll x,ll y,ll pk){ x%=pk; ll ret=1; while(y){ if(y&1) ret=(ret*x)%pk; x=(x*x)%pk; y>>=1; } return ret; } ll calc(ll n,ll pi,ll pk){//计算阶乘部分 （质因子已经提前处理这里不予考虑） if(!n) return 1; ll res=1; for(reg i=2;i<pk;++i)//每个循环节 if(i%pi) res=(res*i)%pk; res=qm(res,n/pk,pk); for(reg i=2;i<=n%pk;++i) if(i%pi) res=(res*i)%pk; return res*calc(n/pi,pi,pk)%pk; } void exgcd(ll a,ll b,ll &x,ll &y){//exgcd if(!b){ x=1,y=0;return; } exgcd(b,a%b,y,x); y-=(a/b)*x; } ll inv(ll n,ll pk){//逆元 ll x,y;exgcd(n,pk,x,y); x=(x%pk+pk)%pk; return x; } ll C(ll n,ll m,ll pi,ll pk){//计算C(n,m)mod pi^k ll up=calc(n,pi,pk),d1=calc(m,pi,pk),d2=calc(n-m,pi,pk); ll k=0; for(reg i=n;i;i/=pi) k+=i/pi;//处理质因子个数 for(reg i=m;i;i/=pi) k-=i/pi; for(reg i=n-m;i;i/=pi) k-=i/pi; return up*inv(d1,pk)%pk*inv(d2,pk)%pk*qm(pi,k,pk)%pk; } ll CRT(ll b,ll mod){//CRT每步算出来了之后直接合并 return (b*inv(p/mod,mod)%p*(p/mod))%p; } ll EXLUCAS(ll n,ll m){//质因数分解+开始处理C ll ret=0; ll tmp=p; for(reg i=2;(ll)i*i<=tmp;++i){ if(tmp%i==0){ ll pi=i,pk=1; while(tmp%i==0) pk*=i,tmp/=i; (ret+=CRT(C(n,m,pi,pk),pk))%=p; } } if(tmp>1) (ret+=CRT(C(n,m,tmp,tmp),tmp))%=p; return ret; } int main(){ ll n,m; scanf("%lld%lld%lld",&n,&m,&p); printf("%lld",EXLUCAS(n,m)); return 0; } } signed main(){ Miracle::main(); return 0; } /* Author: *Miracle* Date: 2018/12/1 8:44:42 */
这个算法的核心思路是：
1.不互质的要提出来单独处理
2.直接处理P，不互质的太多了
3.分成质因子处理，CRT合并
4.对于阶乘上下提出不互质的部分（质因子），转化成互质存在逆元的情况
5.观察剩余部分，后面可以对pk取模。
发现一部分还是阶乘，递归处理。
另一部分发现有循环节，利用循环节加速处理。
剩下的边角考虑一下。
（5本质上就是对每个数提取pi质因子，剩下的再乘起来。不过用递归和循环节加速了一下）
还有一个无聊的题:
简单的组合数学，非要考你扩展LUCAS。。。。