莫比乌斯反演——hdu6390推公式

/*
首先要把原始化简成 k/phi[k] 的格式，然后把有关k的sigma提出来，
后面就是求gcd(i,j)==k的莫比乌斯反演
这里要用整除分块加下速
*/
#include<bits/stdc++.h>
using namespace std;
#define ll long long 
#define maxn 2000005
ll n,m,mod;

bool vis[maxn+10];
int prime[maxn+10],mm,phi[maxn+10],mu[maxn+10],sum[maxn+10];
void primes(){
    phi[1]=1;mu[1]=1;
    for(int i=2;i<maxn;i++){
        if(!vis[i]){mu[i]=-1;prime[++mm]=i;phi[i]=i-1;}
        for(int j=1;j<=mm;j++){
            if(i*prime[j]>=maxn)break;
            vis[i*prime[j]]=1;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j],mu[i*prime[j]]=0; 
                break;
            }
            else phi[i*prime[j]]=phi[i]*(prime[j]-1),mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=maxn;i++)
        sum[i]=sum[i-1]+mu[i];
} 
//处理i/phi[i] 
ll x[maxn],inv[maxn];
void init(){
    inv[0]=inv[1]=1;
    for(ll i=2;i<=n;i++)
        inv[i]=(mod-mod/i)*inv[mod%i]%mod;
    for(ll i=1;i<=n;i++)
        x[i]=i*inv[phi[i]]%mod;
}
inline ll calc(ll i){//求莫比乌斯公式值 
    ll d1=n/i,d2=m/i,res=0;
    for(ll l=1,r;l<=d1;l=r+1){
        r=min(d1/(d1/l),d2/(d2/l));
        res=(res+(d1/l)*(d2/l)%mod*(sum[r]-sum[l-1])%mod+mod)%mod;
    }
    
    return res; 
}

int main(){
    int t;
    scanf("%d",&t);
    primes();
    while(t--){
        scanf("%lld%lld%lld",&n,&m,&mod);
        init();
        ll ans=0;
        if(n>m)swap(n,m);
        for(int i=1;i<=n;i++)
            ans=(ans+x[i]*calc(i)%mod+mod)%mod;
        printf("%lld
",ans);
    }
return 0;
}