【数论】HDU6390 GuGuFishtion

分析：

//详情参见problem B
很明显，如果$(x,y)=1$(“（）”表示最大公因数)
那么$\varphi(x*y)=\varphi(x)*\varphi(y)$

那么如果$(x,y)>1$
将$x、y$唯一分解：
$x=p_{x1}^{k_{x1}}\times p_{x2}^{k_{x2}}\times ……\times p_{xm}^{k_{xm}}$
$y=p_{y1}^{k_{y1}}\times p_{y2}^{k_{y2}}\times ……\times p_{yn}^{k_{yn}}$

然后发现，若x存在某个质因子$p^a$，且y存在同样的一个质因子$p^b$，那么它们的乘积中就有$p^{a+b}$这个质因子。
根据$\varphi$函数的定义：$\varphi(n)=(p_1^{k_1}-p_1^{k_1-1})*(p_2^{k_2}-p_2^{k_2-1})*(p_3^{k_3}-p_3^{k_3-1})……$
那么对于$\varphi(x),p$的贡献为$p^a-p^{a-1}$
对于$\varphi(y),p$的贡献为$p^b-p^{b-1}$
对于$\varphi(x*y),p$的贡献为$p^{a+b}-p^{a+b-1}$
发现：$(p^a-p^{a-1})\times (p^b-p^{b-1})=p^{a+b}-2\times p^{a+b-1}+p^{a+b-2}=(p^{a+b}-p^{a+b-1})\times \frac {p-1} {p}$

所以，定义$g=(x,y)$
$g=p_{g1}^{k_{g1}}\times p_{g2}^{k_{g2}}\times ……\times p_{gs}^{k_{gs}}$
那么$\frac {\varphi(x*y)}{\varphi(x)*\varphi(y)}=\frac {p_{g1}} {p_{g1}-1}\times \frac {p_{g2}} {p_{g2}-1}\times \frac {p_{g3}} {p_{g3}-1}……$

所以只需要枚举一下gcd=g，然后求所有gcd为g的方案数即可。莫比乌斯反演一下就出来了。。（其实用不着反演，直接用莫比乌斯函数的性质上去容斥就可以了）

#include<cstdio>
#include<cstring>
#include<algorithm>
#define SF scanf
#define PF printf
#define MAXN 1000000
using namespace std;
typedef long long ll;
int tot;
int primes[MAXN+10],isprime[MAXN+10],minp[MAXN+10],mu[MAXN+10];
ll fspx[MAXN+10];
ll p,n,m;
ll fsp(ll x,int y){
    ll res=1;
    while(y){
        if(y&1)
            res=res*x%p;
        x=x*x%p;
        y>>=1;    
    }
    return res;
}
void prepare(){
    mu[1]=1;
    for(int i=2;i<=MAXN;i++){
        if(isprime[i]==0){
            primes[++tot]=i;
            minp[i]=i;
            mu[i]=-1;
        }
        for(int j=1;primes[j]*i<=MAXN;j++){
            isprime[i*primes[j]]=1;
            if(i%primes[j]==0){
                minp[i*primes[j]]=primes[j];
                mu[i*primes[j]]=0;
                break;
            }
            minp[i*primes[j]]=primes[j];
            mu[i*primes[j]]=-mu[i];
        }
    }
    for(int i=1;i<=MAXN;i++)
        mu[i]=mu[i-1]+mu[i];
}
int t;
ll ans;
ll work(ll x,ll y){
    if(x>y)
        swap(x,y);
    ll last=0,res=0;
    for(ll i=1;i<=x;i=last+1){
        last=min(x/(x/i),y/(y/i));
        res+=1ll*(x/i)*(y/i)%p*(mu[last]-mu[i-1])%p;
        res%=p;
    }
    return (res+p)%p;
}
ll solve(ll x){
    ll x1=x;
    ll s=1;
    ll s1=1;
    while(x!=1ll){
        if(s1%minp[x]!=0){
            (s*=1ll*minp[x]*fspx[minp[x]]%p)%=p;
            s1*=minp[x];
        }
        x/=minp[x];
    }
    ll sumx=work(n/x1,m/x1);
    ll res=sumx*s%p;
    return res;
}
int main(){
    prepare();
    SF("%d",&t);
    while(t--){
        ans=0;
        SF("%I64d%I64d%I64d",&n,&m,&p);
        for(int i=1;i<=tot&&primes[i]<=min(n,m);i++){
            fspx[primes[i]]=fsp(primes[i]-1,p-2);
        }
        for(int i=2;i<=min(n,m);i++){
            //PF("{%d %lld}\n",i,solve(i));
            (ans+=solve(i))%=p;
        }
        (ans+=work(n,m))%=p;
        ans=(ans+p)%p;
        PF("%I64d\n",ans);
    }
}