「狗屁不会」exlucas

Posted 2020-11-24 hzoi-kx

填坑
考试的时候不会自闭了，考完板子打不对自闭了。
lucas用来处理一类问题 (C_n^{m}mod P) 其中P为质数。
可是如果P变成GK非质数了咋办？？
就要用exlucas了。
首先，如果设P=(prod {p_i}^{c_i}),那么如果能分别求出组合数模质数幂的答案，就可以用CRT合并了。
现在问题转化为求(C_n^m mod p^k)。
有基本式子(C_n^m=frac{n!}{m!(n-m)!})，然后唯一的问题就是如果n!,(n-m)!模(p^k)意义下无逆元怎么办。
无逆元是因为不互质，所以考虑把p都提出来，也就是(C_n^m=frac{frac{n!}{p^x}}{frac{m!}{p^y}frac{(n-m)!}{p^z}}p^{x-y-z})也就是求参考lucas的形式，可以递归实现求(frac{n!}{p^x})。
具体来说，对于n先把含p因子的提出来。
(n!=p^{left lfloor frac{n}{p} ight floor}left lfloor frac{n}{p} ight floor!prodlimits_{i=1&&iperp p}^n i)
质因子只会来自于第一项与第二项。第一项直接搞掉，第二项递归下去，第三项可以预处理。
具体看代码吧。

#include<cstdio>
#include<iostream>
#define LL long long
#define int LL
#define Pair pair<int,int>
#define fr first
#define sc second
#define mp make_pair
using namespace std;
const int N=35;
int pk[N],pc[N],tt,p,k,kx,tk[N],js[N][1000005];
inline int qpow(int a,int b,int mod=k,int c=1) {
    for(;b;b>>=1,a=1ll*a*a%mod) if(b&1) c=1ll*c*a%mod;
    return c;
}
Pair calc(int x) {
    if(!x) return make_pair(1,0);Pair tmp=calc(x/p);
    return make_pair(tmp.fr*qpow(js[kx][k-1],x/k)%k*js[kx][x%k]%k,tmp.sc+x/p);
}
inline int C(int n,int m) {
    if(m>n) return 0; Pair x=calc(n),y=calc(m),z=calc(n-m);
    return x.fr*qpow(y.fr,k-k/p-1)%k*qpow(z.fr,k-k/p-1)%k*qpow(p,x.sc-y.sc-z.sc)%k;
}
inline int CRT(int ans=0,int M=1) {
    for(int i=1;i<=tt;++i) M=M*pc[i];
    for(int i=1;i<=tt;++i) k=pc[i],ans=(ans+qpow(M/k,pc[i]-pc[i]/pk[i]-1)*tk[i]%k*(M/k))%M;
    return ans;
}
signed main() {
    LL n,m;int P; scanf("%lld%lld%lld",&n,&m,&P);int x=P;
    for(int i=2;i*i<=x;++i) {
        if(x%i==0) {
            ++tt;pk[tt]=pc[tt]=i;x/=i;
            while(x%i==0) x/=i,pc[tt]*=i;
        }
    }
    if(x^1) ++tt,pk[tt]=pc[tt]=x;
    for(kx=1;kx<=tt;++kx) {
        js[kx][0]=1;p=pk[kx],k=pc[kx];
        for(int j=1;j<k;++j) {
            if(j%p==0) js[kx][j]=js[kx][j-1];
            else js[kx][j]=1ll*js[kx][j-1]*j%k;
        }
        tk[kx]=C(n,m);
    }
    printf("%lld
",CRT());
    return 0;
}