bzoj3481 DZY Loves Math III

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

$$\begin{align} & \sum_{d} [d|P][d|Q]d\sum_{x}[\gcd(x,\frac Pd)=1] \\ = & \sum_{d|\gcd(P,Q)} \varphi(\frac P d) d \end{align}$$
容易发现这个式子具有积性，因此我们只需要对每个质因子单独考虑。
记 $P=p^qx$， $Q=p^{q'}x$，对于 $p$的答案是
$$\begin{align} & \sum_{i=0}^{q'}\varphi(p^{q-i})p^i \\ = & (q'+1)(p-1)(p^{q-1}) \end{align}$$
注意当 $i=q$时展开欧拉函数并不会产生 $(p-1)$一项，需要特判。
质因数分解用到rho，注意 $Q=0$的情况。

#include<cstdio>
#include<algorithm>
using namespace std;
#define LL long long
const 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);
}