給出一個表示式

求解

首先，設p為a，b共同的質因數

phi(n) = p^α = (p-1)*p^(α-1)

phi(ab) = (p-1)*p^(α1+α2-1)

phi(a) = (p-1)*p^(α1-1)

phi(b) = (p-1)*p^(α2-1)

phi(ab)/(phi(a)*phi(b)) = p-1/p;

這樣原式子就能轉化為

若p只為a或b的質因子，可以發現約分之後都為1

所以上式成立。

所以問題就轉化為

進一步轉化為

可以通過莫比烏斯反演來求解gcd(a,b)=k的個數，最後乘上逆元即為答案

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <cmath>

using namespace std;
typedef long long ll;
const int maxn = 1e6+10;
bool check[maxn];
int prime[maxn];
int mu[maxn];
int phi[maxn];
ll inv[maxn];
int p;


void init()
{
    for(int i = 2; i < maxn; i++)
        phi[i] = 0;
    phi[1] = 1;
    for(int i = 2; i < maxn; i++)
    {
        if(!phi[i])
        {
            for(int j = i; j < maxn; j+=i)
            {
                if(!phi[j]) phi[j] = j;
                phi[j] = phi[j]/i*(i-1);
            }
        }
    }
}

void Moblus()
{
    memset(check,false,sizeof(check));
    mu[1] = 1;
    int tot = 0;
    for(int i = 2; i < maxn; i++)
    {
        if(!check[i])
        {
            prime[tot++] = i;
            mu[i] = -1;
        }
        for(int j = 0; j < tot; j++)
        {
            if(i*prime[j]>maxn)
                break;
            check[i*prime[j]] = true;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]] = 0;
                break;
            }
            else
            {
                mu[i*prime[j]] = -mu[i];
            }
        }
    }
}

ll cal(int b,int d,int k)
{
    ll ans1 = 0;
    b = b/k;
    d = d/k;
    if(b>d)
        swap(b,d);
    for(int i = 1; i <= b; i++)
    {
        ans1 += (ll)mu[i]*(b/i)*(d/i);
    }

    return ans1;
}

void init2()
{
    inv[1] = 1;
    for(int i = 2; i < maxn; i++)
    {
        inv[i] = inv[p%i]*(ll)(p-p/i)%p;
    }
}

int main()
{
    init();
    Moblus();
    int T;
    scanf("%d",&T);
    while(T--)
    {
        int n,m;
        scanf("%d%d%d",&n,&m,&p);
        ll ans = 0;
        init2();
        for(int k = 1; k <= min(n,m); k++)
        {
            ll temp = cal(n,m,k)%p;
            temp = temp*k%p*inv[phi[k]];
            ans = (ans+temp)%p;
        }
        printf("%lld\n",ans);
    }
    return 0;
}