bzoj3481 DZY Loves Math III

当x一定时，方程相当于一个模线性方程组，当且仅当gcd(x,P)|Q时有gcd(x,P)个解。于是我们枚举d=gcd(x,P)，答案是

=∑d[d|P][d|Q]d∑x[gcd(x,Pd)=1]∑d|gcd(P,Q)φ(Pd)d

容易发现这个式子具有积性，因此我们只需要对每个质因子单独考虑。
记P=pqx，Q=pq′x，对于p的答案是

=∑i=0q′φ(pq−i)pi(q′+1)(p−1)(pq−1)

注意当i=q时展开欧拉函数并不会产生(p−1)一项，需要特判。
质因数分解用到rho，注意Q=0的情况。

#include<cstdio>#include<algorithm>using namespace std;#define LL long longconst int mod=1000000007;LL p[1000];int q[1000],q1[1000],tot;LL Rand(LL n){    return (((LL)rand()<<31)|rand())%(n-1)+1;}LL inc(LL x,LL y,LL p){    x+=y;    return x>=p?x-p:x;}LL dec(LL x,LL y,LL p){    x-=y;    return x<0?x+p:x;}LL mul(LL x,LL y,LL p){    LL ret=0;    for (;y;y>>=1,x=inc(x,x,p))        if (y&1) ret=inc(ret,x,p);    return ret;}LL pow(LL x,LL y,LL p){    LL ret=1;    for (;y;y>>=1,x=mul(x,x,p))        if (y&1) ret=mul(ret,x,p);    return ret;}LL gcd(LL x,LL y){    return y?gcd(y,x%y):x;}void add(LL x,int flag){    if (flag)    {        for (int i=1;i<=tot;i++)            if (p[i]==x)            {                q[i]++;                return;            }        p[++tot]=x;        q[tot]=1;    }    else    {        for (int i=1;i<=tot;i++)            if (p[i]==x)            {                if (q1[i]<q[i]) q1[i]++;                return;            }    }}int check(LL n){    LL x=n-1,y,z;    int k=0;    while (!(x&1)) x>>=1,k++;    for (int i=1;i<=8;i++)    {        y=pow(Rand(n),x,n);        for (int j=1;j<=k;j++)        {            z=mul(y,y,n);            if (z==1&&y!=1&&y!=n-1) return 0;            y=z;        }        if (y!=1) return 0;    }    return 1;}void solve(LL n,int flag){    if (n==1) return;    if (check(n))    {        add(n,flag);        return;    }    LL x,y,d,c=Rand(n);    x=y=Rand(n);    for (int i=1,k=1;;i++)    {        x=inc(mul(x,x,n),c,n);        d=gcd(dec(x,y,n),n);        if (d>1)        {            solve(d,flag);            solve(n/d,flag);            return;        }        if (i==k) k<<=1,y=x;    }}int main(){    //freopen("c.in","r",stdin);    srand(2333);    int n;    LL x,ans=1;    scanf("%d",&n);    for (int i=1;i<=n;i++)    {        scanf("%lld",&x);        solve(x,1);    }    for (int i=1;i<=n;i++)    {        scanf("%lld",&x);        if (!x)        {            for (int j=1;j<=tot;j++) q1[j]=q[j];            break;        }        else solve(x,0);    }    for (int i=1;i<=tot;i++)        ans=mul(ans,mul(pow(p[i]%mod,q[i]-1,mod),inc(mul((p[i]-1)%mod,q1[i]+1,mod),(q[i]==q1[i]),mod),mod),mod);    printf("%lld\n",ans);}