考慮對於每一個x有多少個合法解。得到ax+by=c形式的方程。如果gcd(x,y)|c，則a在0~y-1範圍內的解的個數為gcd(x,y)。也就是說現在所要求的是Σ[gcd(x,P)|Q]*gcd(x,P)。

　　對這個式子套路地枚舉gcd，可以得到Σdφ(P/d) (d|gcd(P,Q))。質因子間相互獨立，考慮每個質因子的貢獻再累乘。如果d取完了P的某項質因子，那麽該質因子的貢獻為p_i^q_i，否則為(p_i-1)p_i^q_i-1。於是rho分解完質因數就可以算了。

　　註意特判Q=0。

#include<iostream> 
#include<cstdio>
#include
 
<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
#define P 1000000007
#define ll long long
ll read()
{
    ll x=0,f=1;char c=getchar();
    while (c<‘0‘||c>‘9‘) {if (c==‘-‘) f=-1;c=getchar();}
    while (c>=‘0‘&&c<=‘9‘) x=(x<<1
 
)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
int n,cntp=0,cntq=0,ans=1;
ll p[11],q[11],a[500],b[500],c[500],d[500];
ll gcd(ll n,ll m){return m==0?n:gcd(m,n%m);}
ll ksc(ll a,ll b,ll p)
{
    ll t=a*b-(ll)((long double)a*b/p+0.5)*p;
    return (t<0)?t+p:t; 
}
ll ksm(int a,ll k,ll p)
{
    if (k==0) return
 
 1;
    ll tmp=ksm(a,k>>1,p);tmp=ksc(tmp,tmp,p);
    if (k&1) return ksc(tmp,a,p);else return tmp;
}
bool check(int k,ll n)
{
    if (k>=n) return 1;
    ll p=n-1;
    ll t=ksm(k,p,n);
    if (t==n-1) return 1;
    if (t!=1) return 0;
    while (!(p&1))
    {
        p>>=1;
        ll t=ksm(k,p,n);
        if (t==n-1) return 1;
        if (t!=1) return 0;    
    }
    return 1;
}
bool Miller_Rabin(ll n)
{
    if (n==1) return 0;
    for (int i=2;i*i<=min(n,100ll);i++)
    if (n%i==0) return n==i;
    if (n<=100) return 1;
    else return check(2,n)&&check(3,n)&&check(7,n)&&check(61,n)&&check(24251,n)&&n!=46856248255981;
}
ll f(ll x,ll n,int c){return (ksc(x,x,n)+c)%n;}
void Pollard_Rho(ll n,ll *a,int &cnt)
{
    if (n==1) return;
    if (Miller_Rabin(n)) {a[++cnt]=n;return;}
    if (n<=100) for (int i=2;i<=n;i++) if (n%i==0&&Miller_Rabin(n/i)) {a[++cnt]=n/i;Pollard_Rho(i,a,cnt);return;}
    while (1)
    {
        int c=rand()%(n-1)+1;
        ll x=(rand()%n+c)%n,y=x;
        do
        {
            ll z=gcd(abs(x-y),n);
            if (z>1&&z<n) {Pollard_Rho(n/z,a,cnt),Pollard_Rho(z,a,cnt);return;}
        }while ((x=f(x,n,c))!=(y=f(f(y,n,c),n,c)));
    }
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("bzoj3481.in","r",stdin);
    freopen("bzoj3481.out","w",stdout);
    const char LL[]="%I64d\n";
#else
    const char LL[]="%lld\n";
#endif
    n=read();srand(20020509);
    cntp=0;for (int i=1;i<=n;i++) p[i]=read(),Pollard_Rho(p[i],a,cntp);
    cntq=0;for (int i=1;i<=n;i++) q[i]=read(),Pollard_Rho(q[i],b,cntq);
    sort(a+1,a+cntp+1);sort(b+1,b+cntq+1);
    for (int i=1;i<=cntp;i++)
    {
        int t=i;
        while (a[t+1]==a[i]) t++;
        c[i]=t-i+1;i=t;
    }
    for (int i=1;i<=cntq;i++)
    {
        int t=i;
        while (b[t+1]==b[i]) t++;
        d[i]=t-i+1;i=t;
    }
    for (int i=1;i<=cntp;i++)
    if (c[i]&&!c[i-1])
        for (int j=i-1;j&&!c[j];j--) c[j]=c[j+1],c[j+1]=0;
    cntp=unique(a+1,a+cntp+1)-a-1;
    for (int i=1;i<=cntq;i++)
    if (d[i]&&!d[i-1])
        for (int j=i-1;j&&!d[j];j--) d[j]=d[j+1],d[j+1]=0;
    cntq=unique(b+1,b+cntq+1)-b-1;
    for (int i=1;i<=cntp;i++) a[i]%=P;
    for (int i=1;i<=cntq;i++) b[i]%=P;
    if (b[1]==0)
    {
        cntq=cntp;
        for (int i=1;i<=cntp;i++) b[i]=a[i],d[i]=c[i];
    }
    for (int j=1;j<=cntp;j++)
    {
        int x=0;
        for (int i=1;i<=cntq;i++)
        if (b[i]==a[j]) {x=d[i];break;}
        if (x<c[j]) ans=1ll*ans*ksm(a[j],c[j]-1,P)%P*(x+1)%P*(a[j]-1)%P;
        else ans=1ll*ans*ksm(a[j],c[j]-1,P)%P*(1ll*c[j]*(a[j]-1)%P+a[j])%P;
    }
    cout<<ans<<endl; 
    return 0;
}

BZOJ3481 DZY Loves Math III（數論+Pollard_Rho）