题意：

有$G_u(a,b)=\frac{\varphi(ab)}{\varphi(a)\varphi(b)}$。给定m,n,p，求$(\sum_{a=1}^m\sum_{b=1}^nG_u(a,b))\ mod\ p$。其中$1≤m,n≤1000000$，$max(m,n)<p≤1000000007$，p为质数。

思路：

可以推导出$\frac{\varphi(ab)}{\varphi(a)\varphi(b)}=\frac{gcd(a,b)}{\varphi(gcd(a,b))}$。

举个例子形象一点，

因为有欧拉公式$φ(n) = n *(1-\frac{1}{x})* (1-\frac{1}{y})$，其中x,y为n的质因子。

$$\frac{\varphi(72)}{\varphi(6)\varphi(12)}=\frac{72*(1-\frac{1}{2})*(1-\frac{1}{3})}{6*(1-\frac{1}{2})*(1-\frac{1}{3})*12*(1-\frac{1}{2})*(1-\frac{1}{3})}=\frac{1}{(1-\frac{1}{2})*(1-\frac{1}{3})}$$

化简后分母只剩下$gcd$的质因子的部分，上下同乘$gcd$，得

$$\frac{6}{(6-3)*(1-\frac{1}{3})}=\frac{6}{6-2-3+1}$$

这里分子就是$gcd(6,12)$，分母则是$gcd(6,12)$中减去2的倍数的个数，减去3的倍数的个数，加上6的倍数的个数，这个容斥就是求出了比$gcd(6,12)$小的数中与它互质的数的个数，即$\varphi(gcd(6,12))$。

所以问题就转化为了求$(\sum_{a=1}^m\sum_{b=1}^n\frac{gcd(a,b)}{\varphi(gcd(a,b))})\ mod\ p$。这里用到了容斥。

枚举$gcd(a,b)$。设f[i]:gcd为i的倍数的(a,b)对数，g[i]:gcd为i的(a,b)对数。

以a=6,b=8为例，

$$f[1]=6*8=48,f[2]=3*4=12,f[3]=2*2=4,$$ $$f[4]=1*2=2,f[5]=1*1=1,f[6]=1*1=1$$。 $$g[6]=f[6]=1,g[5]=f[5]=1,g[4]=f[4]=2,g[3]=f[3]-g[6]=4-1=3,$$ $$g[2]=f[2]-g[4]-g[6]=12-2-1=9,$$

$g[1]=f[1]-g[2]-g[3]-g[4]-g[5]-g[6]=48-9-3-2-1-1=32​$。

所以就可以求出g[i]了。

代码：

#define MAXN 1000000
typedef long long ll;
ll eu[MAXN+10];
ll pri[MAXN+10];
bool vis[MAXN+10];
ll f[MAXN+10];//gcd为i的倍数的对数
ll g[MAXN+10];//gcd为i的对数
ll inv[MAXN+10];
void init()
{
    eu[1]=1;
    ll cnt=0;
    for(ll i=2;i<=MAXN;i++)
    {
        if(!vis[i]){pri[cnt++]=i;eu[i]=i-1;}
        for(int j=0;j<=cnt&&pri[j]*i<=MAXN;j++)
        {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0)
            {
                eu[i*pri[j]]=eu[i]*pri[j];
                break;
            }
            eu[i*pri[j]]=eu[i]*(pri[j]-1);
        }
    }
}
long long PowerMod(ll a,ll b,ll c)
{
    long long ans=1;
    long long k;
    k=a;
    while(b>0)
    {
        if(b%2==1)
        ans=(ans*k)%c;
        b=b/2;
        k=(k*k)%c;
    }
    return ans;
}
int main()
{
    int t;
    ll m,n,p;
    init();
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lld%lld%lld",&m,&n,&p);
        ll ans=0;
        for(ll i=1;i<=min(m,n);i++)
            f[i]=(m/i)*(n/i);
        ll tot=min(m,n);
        inv[1]=1;
        for(ll i=2;i<=tot;i++)
           inv[i]=(p-(p/i))*inv[p%i]%p;
        for(ll i=tot;i>=1;i--)
        {
            g[i]=f[i];
            for(ll j=2;i*j<=tot;j++)
                g[i]=(g[i]-g[i*j]+p)%p;
            ll temp=g[i];
            temp=((temp*i)%p*inv[eu[i]])%p;
            ans=(ans+temp)%p;
        }
        printf("%lld\n",ans);
    }
    return 0;
}