\(\text{Solution}\)
容易发现答案为
\[\sum_{x = 1}^n\sum_{y = 1}^m[\gcd(x,y) = 1][\gcd(y,k) = 1]
\]
\[=\sum_{x = 1}^n\sum_{y = 1}^m\sum_{d|\gcd(x,y)}\mu(d)\sum_{g|\gcd(y,k)}\mu(g)
\]
考虑枚举\(\gcd(x,y)\)
\[=\sum_{i = 1}^{\min(n,m)}\mu(i)[i \perp k]\sum_{i | x}\sum_{j = 1}^{\lfloor \frac mi \rfloor}[j \perp k]
\]
设\(g(n) = \sum_{i = 1}^n [i \perp k]\)
\[=\sum_{i = 1}^{\min(n,m)}\mu(i)[i \perp k]\lfloor\dfrac ni\rfloor g(\lfloor \frac mi \rfloor)
\]
如何求\(g(n)\)？
对于互质有一个性质
当\(a\perp b\)时， \((a + c * b) \perp b\)(易证)。
所以\(g(n) = \lfloor\dfrac{n}{k}\rfloor * g(k) + g(n\mod k)\)
现在只需考虑如何求\(\sum_{i = 1}^{\min(n,m)}\mu(i)[i \perp k]\)
设\(f(n) = \sum_{i = 1}^n\mu(i)[i \perp k]\)
那么
\[f(n) = \sum_{i = 1}^n f(\lfloor \frac{n}{i} \rfloor)[i \perp k] - \sum_{i = 2}^n f(\lfloor \frac{n}{i} \rfloor)[i \perp k]
\]
发现
\[\sum_{i = 1}^n f(\lfloor \frac{n}{i} \rfloor)[i \perp k]
\]
\[=\sum_{i = 1}^n \sum_{j = 1}^{\lfloor \frac{n}{i} \rfloor}\mu(j)[i * j \perp k]
\]
设\(T = i * j\)，枚举\(T\)
\[=\sum_{T = 1}^n [T \perp k] \sum_{d|T}\mu(d)
\]
\[=\sum_{T = 1}^n [T \perp k][T = 1]
\]
\[=1
\]
最后柿子就是
\[f(n) = 1 - \sum_{i = 2}^n f(\lfloor \frac{n}{i} \rfloor)[i \perp k]
\]
就可以愉快地杜教筛了，时间复杂度\(O(n^{\frac{2}{3}})\)
\(\text{Code}\)
#include<cstdio>
#include<map>
#include<algorithm>
#define LL long long
using namespace std;
const int M = 2e6;
int n,m,k,vis[M + 5],p[M + 5],tot;
LL f[M + 5],cop[M + 5],mu[M + 5],g[M + 5];
map<int,LL> fm;
LL getG(int x) {return (LL)(x / k) * g[k] + g[x % k];}
int gcd(int x,int y){return y == 0 ? x : gcd(y,x % y);}
void init()
{
mu[1] = cop[1] = 1;
for (int i = 2; i <= M; i++)
{
if (!vis[i]) p[++tot] = i,mu[i] = -1,cop[i] = (k % i ? 1 : 0);
for (int j = 1; j <= tot && p[j] * i <= M; j++)
{
vis[p[j] * i] = 1,cop[p[j] * i] = cop[i] * cop[p[j]];
if (i % p[j] == 0) break;
mu[p[j] * i] = -mu[i];
}
}
for (int i = 1; i <= k; i++) g[i] = g[i - 1] + cop[i];
for (int i = 1; i <= M; i++) f[i] = f[i - 1] + cop[i] * mu[i];
}
LL gets(int x)
{
if (x <= M) return f[x];
if (fm[x]) return fm[x];
LL res = 1LL;
for (int l = 2,r; l <= x; l = r + 1)
r = x / (x / l),res -= gets(x / l) * (getG(r) - getG(l - 1));
return fm[x] = res;
}
int main()
{
scanf("%d%d%d",&n,&m,&k),init(); LL ans = 0;
for (int l = 1,r; l <= min(n,m); l = r + 1)
r = min(n / (n / l),m / (m / l)),ans += (LL)(gets(r) - gets(l - 1)) * (n / l) * getG(m / l);
printf("%lld\n",ans);
}