一个结论：（从二维扩展来的，三维也是对的，证明可以考虑质因数分解）
[d(ijk)=sum_{i'|i}sum_{j'|j}sum_{k'|k}[gcd(i',j')=1][gcd(i', k')=1][gcd(j', k')=1]
]
[sum_{i=1}^asum_{j=1}^bsum_{k=1}^csum_{i'|i}sum_{j'|j}sum_{k'|k}[gcd(i',j')=1][gcd(i', k')=1][gcd(j', k')=1]
]
[sum_{i'=1}^asum_{j'=1}^bsum_{k'=1}^c[gcd(i',j')=1][gcd(i', k')=1][gcd(j', k')=1]lfloorfrac{a}{i'}
floorlfloorfrac{b}{j'}
floorlfloorfrac{c}{k'}
floor
]
[sum_{i'=1}^asum_{j'=1}^blfloorfrac{a}{i'}
floorlfloorfrac{b}{j'}
floor[gcd(i',j')=1]sum_{k'=1}^c[gcd(i'j', k')=1]lfloorfrac{c}{k'}
floor
]
记：
[f(x)=sum_{k'=1}^c[gcd(x, k')=1]lfloorfrac{c}{k'}
floor
]
[f(x)=sum_{k'=1}^csum_{d|x,d|k'}mu(d)lfloorfrac{c}{k'}
floor
]
[f(x)=sum_{d|x}mu(d)sum_{d|k'}lfloorfrac{c}{k'}
floor
]
把先把(g(x)=sum_{x|k'}lfloorfrac{c}{k'} floor)预处理出来。
然后可以预处理出来(f(x)=sum_{d|x}mu(d)g(d),xin[1,ab])
则答案为：
[sum_{i'=1}^asum_{j'=1}^blfloorfrac{a}{i'}
floorlfloorfrac{b}{j'}
floor[gcd(i',j')=1] f(i'j')
]
#include <algorithm>
#include <cstdio>
using namespace std;
const int m = (1 << 30) - 1;
const int N = 4e6 + 5;
int p[N], mu[N], f[N], c;
bool tag[N];
void sieve(int n) {
mu[1] = 1;
for(int i = 2; i <= n; i ++) {
if(!tag[i]) { p[++ c] = i; mu[i] = -1; }
for(int j = 1; j <= c && i * p[j] <= n; j ++) {
tag[i * p[j]] = 1;
if(i % p[j] == 0) break ;
mu[i * p[j]] = - mu[i];
}
}
}
int main() {
int a, b, c;
scanf("%d%d%d", &a, &b, &c);
if(a > b) swap(a, b);
if(a > c) swap(a, c);
if(b > c) swap(b, c);
int ans = 0; sieve(a * b);
for(int i = 1; i <= c; i ++) {
int g = 0;
for(int j = i; j <= c; j += i) {
g += c / j;
}
for(int j = i; j <= a * b; j += i) {
(f[j] += mu[i] * g) &= m;
}
}
for(int i = 1; i <= a; i ++) {
for(int j = 1; j <= b; j ++) {
if(__gcd(i, j) == 1) {
(ans += (1ll * (a / i) * (b / j) * f[i * j]) & m) &= m;
}
}
}
printf("%d
", ans);
return 0;
}