题解
我们要求的其实是这个东西= =
(sum_{i = 1}^{n}sum_{j = 1}^{n}[(i,j) == 1][(j,k) == 1])
然后变一下形
(sum_{j = 1}^{n}[(j,k) == 1]sum_{i = 1}^{n}[(i,j) == 1])
(sum_{j = 1}^{n}[(j,k) == 1]sum_{i = 1}^{n}sum_{d|i,j}mu(d))
(sum_{j = 1}^{n}[(j,k) == 1]sum_{d | j} mu(d) lfloor frac{n}{d}
floor)
(sum_{d = 1}^{n} mu(d) lfloor frac{n}{d}
floor sum_{y = 1}^{lfloor frac{m}{d}
floor}[(yd,k) == 1])
(sum_{d = 1}^{n} mu(d) lfloor frac{n}{d}
floor sum_{y = 1}^{lfloor frac{m}{d}
floor}[(d,k) == 1][(y,k) == 1])
(sum_{d = 1}^{n} [(d,k) == 1] mu(d) lfloor frac{n}{d}
floor sum_{y = 1}^{lfloor frac{m}{d}
floor}[(y,k) == 1])
设(f(n) = sum_{i = 1}^{n}[(i,k) == 1])
我们可以预处理(f(1))到(f(k))
那么就有
(f(n) = lfloorfrac{n}{k}
floor f(k) + f(n mod k))
因为((a,b) == 1) 等价于((a mod b,b) == 1)
我们现在可以(O(n))解决这个式子了，但是还不够
我们可以用数论分块处理(lfloor frac{n}{d}
floor) 和 (lfloor frac{m}{d}
floor)
我们尝试算
(sum_{d = 1}^{n} [(d,k) == 1] mu(d))
设(g(n,k) = sum_{d = 1}^{n} [(d,k) == 1] mu(d))
我们选择一个(k)的质因子(p)，然后把(k)表示成(p^{c}q)的形式
我们从和(q)互质的数里除掉和(p)互质的数
(g(n,k) = sum_{d = 1}^{n} [(d,q) == 1] mu(d) - sum_{y = 1}^{n} [(yp,q) == 1] mu(yp))
由于(p)和(q)互质，所以我们只需要保证([(y,p) == 1][(y,q) == 1])
(g(n,k) = sum_{d = 1}^{n} [(d,q) == 1] mu(d) - mu(p) sum_{y = 1}^{n} [(y,p) == 1][(y,q) == 1] mu(y))
(g(n,k) = sum_{d = 1}^{n} [(d,q) == 1] mu(d) + sum_{y = 1}^{lfloor frac{n}{d}
floor} [(y,k) == 1]mu(y))
(g(n,k) = g(n,q) + g(lfloor frac{n}{p}
floor,k))
边界是(n <= 1)返回(n)，(k = 1)返回一个莫比乌斯函数的前缀和，可以杜教筛
代码
#include <bits/stdc++.h>
//#define ivorysi
#define enter putchar('
')
#define space putchar(' ')
#define fi first
#define se second
#define pb push_back
#define mp make_pair
#define eps 1e-8
#define mo 974711
#define MAXN 1000000
#define pii pair<int,int>
using namespace std;
typedef long long int64;
typedef double db;
template<class T>
void read(T &res) {
res = 0;char c = getchar();T f = 1;
while(c < '0' || c > '9') {
if(c == '-') f = -1;
c = getchar();
}
while(c >= '0' && c <= '9') {
res = res * 10 + c - '0';
c = getchar();
}
res *= f;
}
template<class T>
void out(T x) {
if(x < 0) {putchar('-');x = -x;}
if(x >= 10) {
out(x / 10);
}
putchar('0' + x % 10);
}
int N,M,K;
int mu[MAXN + 5],prime[MAXN + 5],tot,Mu[MAXN + 5];
bool nonprime[MAXN + 5];
int f[2005];
struct HASH {
struct node {
int64 x,v;
int next;
}E[MAXN * 2];
int head[mo + 5],sumE;
HASH() {
memset(head,0,sizeof(head));sumE = 0;
}
void add(int u,int64 x,int64 v) {
E[++sumE].x = x;E[sumE].v = v;E[sumE].next = head[u];
head[u] = sumE;
}
void Insert(int64 x,int64 v) {
add(x % mo,x,v);
}
int64 Query(int64 x){
for(int i = head[x % mo] ; i ; i = E[i].next) {
if(E[i].x == x) return E[i].v;
}
return -1;
}
}H[2];
int gcd(int a,int b) {
return b == 0 ? a : gcd(b,a % b);
}
int64 calcF(int x) {
return 1LL * (x / K) * f[K] + f[x % K];
}
int64 S(int x) {
if(x <= MAXN) return Mu[x];
int64 c = H[0].Query(x);
if(c != -1) return c;
int64 res = 0;
for(int i = 2 ; i <= x ; ++i) {
int r = x / (x / i);
res = res + 1LL * (r - i + 1) * S(x / i);
i = r;
}
res = 1 - res;
H[0].Insert(x,res);
return res;
}
int64 G(int n,int k) {
if(k == 1) return S(n);
else if(n <= 1) return n;
int64 c = H[1].Query(1LL * (n - 1) * K + k);
if(c != -1) return c;
for(int i = 1 ; i <= tot ; ++i) {
if(k % prime[i] == 0) {
int x = k;
while(x % prime[i] == 0) x /= prime[i];
int64 res = G(n,x) + G(n / prime[i],x * prime[i]);
H[1].Insert(1LL * (n - 1) * K + k,res);
return res;
}
}
}
void Solve() {
int t = min(N,M);
int64 res = 0;
for(int i = 1 ; i <= t ; ++i) {
int r = min(N / (N / i),M / (M / i));
int64 s = calcF(M / i) * (N / i);
res = res + s * (G(r,K) - G(i - 1,K));
i = r;
}
out(res);enter;
}
int main() {
#ifdef ivorysi
freopen("f1.in","r",stdin);
#endif
read(N);read(M);read(K);
mu[1] = 1;Mu[1] = 1;
for(int i = 2 ; i <= MAXN ; ++i) {
if(!nonprime[i]) {
prime[++tot] = i;
mu[i] = -1;
}
for(int j = 1 ; j <= tot ; ++j) {
if(prime[j] > MAXN / i) break;
nonprime[i * prime[j]] = 1;
if(i % prime[j] == 0) break;
else mu[i * prime[j]] = -mu[i];
}
Mu[i] = Mu[i - 1] + mu[i];
}
for(int i = 1 ; i <= K ; ++i) {
f[i] = f[i - 1] + (gcd(i,K) == 1);
}
Solve();
return 0;
}