dls真厉害，快速求$\sum_{a=1}^n \sum_{b=1}^m gcd(a,b) $的个数，我想的方法是根据上节课dls讲的方法，要容过来容过去，这次不用了。
则$f[d]=（n/d)\times (m/d)$。
而$g[d]=f[d]-\sum_{d|x} g[x]$，从大到小枚举因子就可以了。
#include <bits/stdc++.h> using namespace std; const int maxn = 1e6 + 50; typedef long long ll; bool flag[maxn]; //标记数组 ll phi[maxn]; //欧拉函数值 int prime[maxn]; //同时得到素数筛 int cnt = 0; void Get_phi(int n) { cnt = 0; memset(flag,true,sizeof(flag)); phi[1] = 1; for(int i=2;i<=n;i++) { if(flag[i]) //素数 { prime[cnt++] = i; phi[i] = i-1; //素数的欧拉函数值是i-1 } for(int j=0;j<cnt;j++) { if(i*prime[j]>n) { break; } flag[i*prime[j]] = false;//素数的倍数不是素数 if(i%prime[j]==0) //i%mod prime = 0,那么phi(i*p) = p*phi(i) { phi[i*prime[j]] = prime[j]*phi[i]; break; } else phi[i*prime[j]] = (prime[j]-1)*phi[i];//i mod prime != 0, 那么 phi(i * prime) == phi(i) * (prime-1) } } } int inv[maxn]; ll f[maxn]; int main() { int T; scanf("%d", &T); Get_phi(1e6 + 10); while(T--) { int n, m, mod; scanf("%d %d %d", &n, &m, &mod); if(n > m) swap(n, m); inv[0] = inv[1] = 1; for(int i = 2; i <= n; i++) { inv[i] = (ll)(mod - mod / i) * inv[mod % i] % (ll)mod; } ll ans = 0; for(int i = n; i >= 1; i--) ///枚举 d | (a, b) { f[i] = (ll)n / i * (m / i); for(int j = i + i; j <= n; j += i) ///枚举 能够整除d的 { f[i] = f[i] - f[j]; } ans = (ans + f[i] * i % mod * inv[phi[i]]) % mod; //printf("%d %lld phi = %lld %lld\n", i, f[i], phi[i], ans); } printf("%lld\n", ans); } return 0; }