WA 3个点 50pts

其中第一组wa的是

919816379170448980 685276 36663

输出应该是0 但是我输出了4848，看了半天没看出个究竟。 excrt部分应该没有问题，我认为bi求错了，但是看不出来。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
/*
x 同余 bi mod ai
t*M同余bi-ans mod ai
M*t+ai*y=bi-ans
bi=f[n]*inv(f[m])*inv(f[n-m])* P^(g[n]-g[m]-g[n-m]);
ai=P^k;

inv

1/v 同余 inv %p
1 同余 inv*v%p
v*inv 同余 1 %p
v*inv+p*y=1;
v*x+p*y=1
显然v,p互质
*/
const LL N=1e6+10;
LL k,ai[N],bi[N],pi[N];

LL exgcd(LL a,LL b,LL &x,LL &y){
    if(!b){ x=1,y=0;return a; }
    LL gcd=exgcd(b,a%b,y,x);y-=a/b*x;
    return gcd;
}
LL inv(LL v,LL mod){
    LL  x,y;
    exgcd(v,mod,x,y);
    return (x%mod+mod)%mod;
}

LL mul(LL a,LL b,LL mod){
    LL ret=0;
    while(b){
        if(b&1) ret=(ret+a)%mod;
        b>>=1;a=(a+a)%mod;
    }
    return (ret+mod)%mod;
}
LL power(LL a,LL b,LL mod){
    LL ret=1;
    while(b){
        if(b&1) ret=ret*a%mod;
        a=a*a%mod;b>>=1; 
    }
    return (ret+mod)%mod;
}

LL calc(LL r,LL p,LL mod){
    LL ret=1;
    for(LL i=1;i<=r;i++)
        if(i%p) ret=ret*(i%mod)%mod;
    return (ret+mod)%mod;
}//计算1-r中%p不为0的数的乘积
inline LL f(LL n,LL p,LL pk){
    if(n==0) return 1;
    return (f(n/p,p,pk) * power( calc(pk,p,pk) , n/pk , pk ) % pk  * calc( n%pk,p,pk )%pk+pk )%pk ;
}//1-n中非p因子的乘积
inline LL g(LL n,LL p){
    if(n<p) return 0;
    return n/p+g(n/p,p);
}//计算p因子个数

LL excrt(){
    LL x,y,M=1,ans=0;
    for(LL i=1;i<=k;i++){
        LL a=M,b=ai[i],c=((bi[i]-ans)%ai[i]+ai[i])%ai[i];
        LL gcd=exgcd(a,b,x,y),ag=ai[i]/gcd;
        if(c%gcd) return -1;//扩欧有解条件
        x=mul(x,c/gcd,ag);//x*(c/gcd)%ag
        ans+=x*M;M*=ag;
        ans=(ans+M)%M;
    }
    return ans;
}
LL exLucas(LL n,LL m,LL p){
    k=0;
    for(LL i=2;i*i<=p;i++){
        if(p%i==0){
            ++k;pi[k]=i;ai[k]=1;
            while(p%i==0) p/=i,ai[k]*=i;
        }
    }
    //计算bi
    if(p!=1){
        ++k;pi[k]=p;ai[k]=p;
    }
    for(LL i=1;i<=k;i++){
        bi[i]= f(n,pi[i],ai[i]) * inv( f(m,pi[i],ai[i]) ,ai[i] ) %ai[i] * inv( f(n-m,pi[i],ai[i]) ,ai[i] ) %ai[i] 
                    *power( pi[i], g(n,pi[i]) - g(m,pi[i]) - g(n-m,pi[i]) , ai[k] )%ai[i];
        cout<<ai[i]<<" "<<bi[i]<<endl;
    }
    return excrt();
}

int main(){
    LL n,m,p;
    scanf("%lld%lld%lld",&n,&m,&p);
    printf("%lld\n",exLucas(n,m,p));
    return 0;
}
/*
p=3,t=2,n=19
19! =19*18*17*16*15*14*13*12*11*10*9*8*7*6*5*4*3*2*1
        =(19*17*16*14*13*11*10*8*7*5*4*2*1)*(18*15*12*9*6*3)
        =(19*17*16*14*13*11*10*8*7*5*4*2*1) * 3^6 * 6!
                                                                                            p^(n/p) *(n/p)!
        =19*(8*7*5*4*2*1)*(17*16*14*13*11*10)*…
        以p^t为周期，计算一个周期的+不满一周期的

mod P^t
=19*(8*7*5*4*2*1)^2*p^(n/p)*(n/p)
*/