考慮枚舉x，分情況討論

當Q=0時，無論x取什麽都有對應的y合法，個數為gcd(x,P)個

否則對於x有合法的y當且僅當gcd(x,P)|Q，並且數量為gcd(x,P)，這個仔細思考還是不難理解的

ans=sigema(1~P-1)x [gcd(x,P)|Q]*gcd(x,P)

=sigema(d|gcd(P,Q))d d*(sigema(1~P/d)x [gcd(x,P/d)==1]) (如果Q等於0把d的限制去掉就完事了)

後面明顯是歐拉函數，=sigema(d|gcd(P,Q))d d*phi(P/d)

然後是一個很類似狄利克雷卷積的形式，當積性函數好了

那麽拆分出來，把每個質因子分開搞，設A為gcd(P,Q)在全局第k個質因子的出現次數，B為P在全局第k個質因子出現的次數

=product(1~plen)k [sigema(0~A)i pk^i*phi(pk^(B-i))]

=product(1~plen)k [sigema(0~A)i pk^i* (pk^(B-i)*(pk-1)/pk)]

=product(1~plen)k [sigema(0~A)i pk^B-1*(pk-1)]

=product(1~plen)k [sigema(0~A)i pk^B-1*(pk-1)]

=product(1~plen)k (A+1)*pk^B-1*(pk-1)

這樣就可以快速算完了，但是有一點需要註意，當i==B時，phi(1)不能表示成1*(1-1)/1的形式，要特判一下A==B的情況

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long LL;
const LL mod=1e9+7;
LL gcd(LL a,LL b)
{
    if(a<0)a=-a;
    if(a==0)return b;
    return gcd(b%a,a);
}
LL quick_mul(LL A,LL p,LL mo)
{
    LL ret
 
=0;
    while(p!=0)
    {
        if(p%2==1)ret=(ret+A)%mo;
        A=(A+A)%mo;p/=2;
    }
    return ret;
}
LL quick_pow(LL A,LL p,LL mo)
{
    LL ret=1;
    while(p!=0)
    {
        if(p%2==1)ret=quick_mul(ret,A,mo);
        A=quick_mul(A,A,mo);p/=2;
    }
    return ret;
}

//-----------------------------------def----------------------------------------------

bool check(LL a,LL t,LL mi,LL n)
{
    LL d=quick_pow(a,t,n);
    if(d==1||d==n-1)return true;
    while(mi--)
    {
        d=quick_mul(d,d,n);
        if(d==n-1)return true;
    }
    return false;
}
int prr[10]={2,3,5,7,11,13,17,19,23,29};
bool miller_rabin(LL n)
{
    if(n<=1)return false;
    for(int i=0;i<=9;i++)
    {
        if(n==prr[i])return true;
        if(n%prr[i]==0)return false;
    }
    
    LL t=n-1;int mi=0;
    while(t%2==0)t/=2,mi++;
    for(int i=1;i<=20;i++)
    {
        LL a=rand()%(n-1)+1;
        if(!check(a,t,mi,n))return false;
    }
    return true;
}
//~~~~~~~~~~~~~~~miller_rabin~~~~~~~~~~~~~~

LL pollard_rho(LL n)
{
    LL a=rand()%n,b=a,c=rand()%n;
    LL i=1,k=2;
    while(1)
    {
        a=(quick_mul(a,a,n)+c)%n;
        LL d=gcd(b-a,n);
        if(d!=1&&d!=n)return d;
        if(a==b)return n;
        
        i++;
        if(i==k)b=a,k*=2;
    }
}
//~~~~~~~~~~~~~~~pollard_rho~~~~~~~~~~~~~~~

int pr;LL prime[110];
void getprime(LL n)
{
    if(n==1)return ;
    if(miller_rabin(n)){prime[++pr]=n;return ;}
    
    LL d=n;
    while(d>=n)d=pollard_rho(d);
    
    getprime(d);
    getprime(n/d);
}
LL plen,p[1100],tt[1100];LL c[2][1100];
void merge()
{
    int i=1,j=1,k=0;
    while(i<=pr&&j<=plen)
    {
             if(prime[i]<p[j])tt[++k]=prime[i++];
        else if(prime[i]>p[j])tt[++k]=p[j++];
        else tt[++k]=p[j],i++,j++;
    }
    while(i<=pr)tt[++k]=prime[i++];
    while(j<=plen)tt[++k]=p[j++];
    memcpy(p,tt,sizeof(p));plen=k;
}
int n;LL s[2][110];
void getci(int w)
{
    for(int i=1;i<=n;i++)
    {
        LL u=s[w][i];
        for(int j=1;j<=plen;j++)
            while(u%p[j]==0)u/=p[j],c[w][j]++;
    }
    if(w==1)for(int j=1;j<=plen;j++)c[w][j]=min(c[0][j],c[1][j]);
}

//------------------------------------------divi--------------------------------------------------

int main()
{
    freopen("4.in","r",stdin);
    freopen("a.out","w",stdout);
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        scanf("%lld",&s[0][i]);
        pr=0,getprime(s[0][i]);
        sort(prime+1,prime+pr+1);
        pr=unique(prime+1,prime+pr+1)-prime-1;
        merge();
    }
    bool bk=false;
    for(int i=1;i<=n;i++)
    {
        scanf("%lld",&s[1][i]); if(s[1][i]==0){bk=true;break;}
        pr=0,getprime(s[1][i]);
        sort(prime+1,prime+pr+1);
        pr=unique(prime+1,prime+pr+1)-prime-1;
        merge();
    }
    if(bk==true)
    {
        getci(0);
        memcpy(c[1],c[0],sizeof(c[1]));
    }
    else getci(0),getci(1);
    
    LL ans=1;
    for(int i=1;i<=plen;i++)
    {
        LL A=c[1][i],B=c[0][i];
        if(A<B)ans=ans*(quick_pow(p[i],B-1,mod)*(p[i]-1)%mod*(A+1)%mod)%mod;
        else ans=ans*( (quick_pow(p[i],B-1,mod)*(p[i]-1)%mod*A%mod+quick_pow(p[i],A,mod))%mod )%mod;
    }
    printf("%lld\n",ans);
        
    return 0;
}

bzoj3481: DZY Loves Math III