bzoj4891: [Tjoi2017]龙舟

求$frac{b_1b_2b_3...b_m}{a_1a_2a_3...a_m}\%M$

M<=1e18,m<=100000,数据组数<=50

用pollard-rho分解M的质因数，提取出$b_i,a_i$与M不互质的部分处理一下

#include<cstdio>
#include<algorithm>
typedef long long i64;
typedef long double ld;
char buf[65536],*ptr=buf+65536;
int G(){
    if(buf+65536==ptr)fread(ptr=buf,1,65536,stdin);
    return *ptr++;
}
i64 _(){
    i64 x=0;
    if(ptr-buf<65000){
        while(*ptr<48)++ptr;
        while(*ptr>47)x=x*10+(*ptr++)-48;
    }else{
        int c=G();
        while(c<48)c=G();
        while(c>47)x=x*10+c-48,c=G();
    }
    return x;
}
int n,m,k;
i64 b[100007],a[22][100007];
int fp,ts[100];
i64 fs[100];
i64 mul(i64 a,i64 b,i64 c){
    if(c<=2000000000ll)return a*b%c;
    i64 r=a*b-i64(ld(a)/c*b)*c;
    if(r>=c||r<=-c)r%=c;
    return r>=0?r:r+c;
}
i64 pw(i64 a,i64 n,i64 P){
    i64 v=1;
    for(;n;n>>=1,a=mul(a,a,P))if(n&1)v=mul(v,a,P);
    return v;
}
i64 gcd(i64 a,i64 b){
    if(a<0)a=-a;
    for(i64 c;b;c=a,a=b,b=c%b);
    return a;
}
bool mr(i64 n){
    i64 z=n-1;
    int t=0;
    while(~z&1)z>>=1,++t;
    for(int i=0;i<15;++i){
        i64 a=rand()%(n-1)+1;
        i64 x=pw(a,z,n);
        for(int j=0;j<t;++j){
            i64 y=mul(x,x,n);
            if(y==1&&x!=1&&x!=n-1)return 0;
            x=y;
        }
        if(x!=1)return 0;
    }
    return 1;
}
i64 get(i64 x,int c){
    int i=1,j=2;
    i64 a=(rand()^i64(rand())<<31)%(x-1)+1,b=a;
    while(1){
        a=mul(a,a,x);
        if((a+=c)>=x)a%=x;
        i64 p=gcd(a-b,x);
        if(p!=1)return p;
        if((++i)==j)j<<=1,b=a;
    }
}
void calc(i64 n){
    if(n==1)return;
    if(mr(n)){
        fs[fp++]=n;
        return;
    }
    for(int c=12347;;++c){
        i64 a=get(n,c);
        if(a!=n){
            i64 b=gcd(a,n/a);
            calc(a/b);
            calc(n/a/b);
            calc(b);
            return;
        }
    }
}
i64 cal(i64*b,i64*a,i64 M){
    fp=0;
    calc(M);
    std::sort(fs,fs+fp);
    fp=std::unique(fs,fs+fp)-fs;
    i64 phi_M=M;
    for(int i=0;i<fp;++i)ts[i]=0,phi_M=phi_M/fs[i]*(fs[i]-1);
    i64 B=1,A=1,x;
    for(int i=1;i<=m;++i){
        x=b[i];
        if(fs[0]==2)for(;~x&1;x>>=1,++ts[0]);
        for(int j=(fs[0]==2);j<fp;++j){
            for(i64 p=fs[j],y=x/p;y*p==x;x=y,y=y/p,++ts[j]);
        }
        B=mul(B,x,M);
    }
    for(int i=1;i<=m;++i){
        x=a[i];
        if(fs[0]==2)for(;~x&1;x>>=1,--ts[0]);
        for(int j=(fs[0]==2);j<fp;++j){
            for(i64 p=fs[j],y=x/p;y*p==x;x=y,y=y/p,--ts[j]);
        }
        A=mul(A,x,M);
    }
    B=mul(B,pw(A,phi_M-1,M),M);
    for(int i=0;i<fp;++i){
        if(ts[i]<0)return -1;
        B=mul(B,pw(fs[i],ts[i],M),M);
    }
    return B;
}
int main(){
    srand(77151);
    n=_();m=_();k=_();
    for(int i=1;i<=m;++i)b[i]=_();
    for(int t=1;t<=n;++t)
    for(int i=1;i<=m;++i)a[t][i]=_();
    for(int i=0;i<k;++i){
        int x=_();
        i64 M=_();
        printf("%lld
",cal(b,a[x],M));
    }
    return 0;
}