[BZOJ4338][BJOI2015]糖果(扩展Lucas)

先求出式子$P_{C_{K+m-1}^{m}}^{n}$，然后对于排列直接$O(n)$求解，对于组合用扩展Lucas求解。

但这题数据并没有保证任何一个模数的质因子的$p^k$在可线性处理的范围内，于是并不会标准解法，只会面向数据编程。

数据中保证了如果某个质因子p的次数不为1，则它的$p^k$一定在可线性处理的范围内，于是只要特判次数为1的质数即可。

次数为1就可以直接求逆元$O(m)$处理了，于是问题解决，虽然随便出组数据就能卡掉。

#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (ll i=(l); i<=(r); i++)
typedef long long ll;
using namespace std;
 
const int N=100010;
ll n,m,K,P,tot,num[N],pi[N],pk[N],sm[N];
 
ll ksm(ll a,ll b,ll p){
    ll res=1;
    for (; b; a=a*a%p,b>>=1)
        if (b & 1) res=res*a%p;
    return res;
}
 
void exgcd(ll a,ll b,ll &x,ll &y){
    if (!b){ x=1; y=0; return; }
        else exgcd(b,a%b,y,x),y-=a/b*x;
}
 
ll inv(ll a,ll p){ ll x,y; exgcd(a,p,x,y); return (x%p+p)%p; }
 
ll F(ll n,ll pi,ll pk){ return n ? ksm(sm[pk],n/pk,pk)*sm[n%pk]%pk*F(n/pi,pi,pk)%pk : 1; }
 
ll exLucas(ll n,ll m,ll pi,ll pk){
    if (n<m) return 0;
    sm[0]=sm[1]=1; rep(i,2,pk) sm[i]=sm[i-1]*((i%pi)?i:1)%pk;
    ll a=F(n,pi,pk),b=F(m,pi,pk),c=F(n-m,pi,pk),k=0;
    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;
    return a*inv(b*c%pk,pk)%pk*ksm(pi,k,pk)%pk;
}
 
ll CRT(ll r[],ll m[]){
    ll res=0;
    rep(i,1,tot) res=(res+P/m[i]*inv(P/m[i],m[i])%P*r[i])%P;
    return res;
}
 
void Frac(ll n){
    tot=0;
    for (ll i=2; i*i<=n; i++) if (n%i==0){
        pi[++tot]=i; pk[tot]=i; n/=i;
        while (n%i==0) n/=i,pk[tot]=pk[tot]*i;
    }
    if (n>1) pi[++tot]=n,pk[tot]=n;
}
 
ll Pe(ll n,ll m,ll P){ ll res=1; rep(i,1,m) res=res*(n-i+P+1)%P; return res; }
 
int main(){
    scanf("%lld%lld%lld%lld",&n,&m,&K,&P);
    Frac(P); ll res=0;
    rep(i,1,tot){
        ll s=1;
        if (pi[i]==pk[i] && pi[i]>m)
            rep(j,1,m) s=s*(K+m-j)%pi[i]*ksm(j,pi[i]-2,pi[i])%pi[i];
        else s=exLucas(K+m-1,m,pi[i],pk[i]);
        s=Pe(s,n,pi[i]); res=(res+P/pi[i]*inv(P/pi[i],pi[i])%P*s)%P;
    }
    printf("%lld
",res);
    return 0;
}