https://www.luogu.com.cn/problem/P1587
首先思考我们要求的是什么？
{ x k l y } = { x y } \large \begin{Bmatrix}\frac{xk^l}{y}\end{Bmatrix} =\begin{Bmatrix}\frac{x}{y}\end{Bmatrix} {yxkl}={yx}
x
k
l
−
y
[
x
k
l
y
]
=
x
−
y
[
x
y
]
\large xk^l-y \begin{bmatrix}\frac{xk^l}{y} \end{bmatrix}=x-y \begin{bmatrix}\frac{x}{y} \end{bmatrix}
xkl−y[yxkl]=x−y[yx]
k
l
≡
1
(
m
o
d
y
)
k^l \equiv 1 (\mod y)
kl≡1(mody)
可以得到我们实际上要求的是
∑ i = 1 n ∑ j = 1 m [ i ⊥ j ] [ j ⊥ k ] \sum_{i=1}^n\sum_{j=1}^m [i\perp j] [j \perp k] i=1∑nj=1∑m[i⊥j][j⊥k]
按照套路
∑ d = 1 m i n ( n , m ) μ ( d ) ∑ i = 1 n / d ∑ j = 1 m / d [ j d ⊥ k ] \large \sum\limits_{d=1}^{min(n, m)} \mu(d) \sum\limits_{i=1}^{n/d} \sum\limits_{j=1}^{m/d} [jd \perp k] d=1∑min(n,m)μ(d)i=1∑n/dj=1∑m/d[jd⊥k]
= ∑ d = 1 m i n ( n , m ) μ ( d ) ( n / d ) ∑ j = 1 m / d [ j ⊥ k ] [ d ⊥ k ] \large = \sum\limits_{d=1}^{min(n, m)} \mu(d) (n/d) \sum\limits_{j=1}^{m/d} [j \perp k][d \perp k] =d=1∑min(n,m)μ(d)(n/d)j=1∑m/d[j⊥k][d⊥k]
设 g ( n ) = ∑ i = 1 n [ i ⊥ k ] g(n)=\sum\limits_{i=1}^n [i\perp k] g(n)=i=1∑n[i⊥k]
容易发现
g
(
n
)
=
(
n
/
k
)
ϕ
(
k
)
+
g
(
n
m
o
d
k
)
g(n)=(n/k) \phi(k)+g(n \mod k)
g(n)=(n/k)ϕ(k)+g(nmodk)
上式子
=
∑
d
=
1
m
i
n
(
n
,
m
)
(
n
/
d
)
g
(
m
/
d
)
μ
(
d
)
[
d
⊥
k
]
\large = \sum\limits_{d=1}^{min(n, m)} (n/d) g(m/d) \mu(d) [d \perp k]
=d=1∑min(n,m)(n/d)g(m/d)μ(d)[d⊥k]
前两项显然可以整除分块，现在要计算的是后两项的区间和
记 f ( n ) = ∑ i = 1 n μ ( i ) [ i ⊥ k ] f(n)=\sum\limits_{i=1}^n \mu(i)[i \perp k] f(n)=i=1∑nμ(i)[i⊥k]
区间和就可以表示为 f ( r ) − f ( l − 1 ) f(r)-f(l-1) f(r)−f(l−1)
考虑怎么快速计算 f f f
发现是两个积性函数卷在一起,还是一个积性函数
我们考虑让它卷上 [ i ⊥ k ] [i \perp k] [i⊥k]
得到杜教筛的那条式子
f ( n ) = ∑ i = 1 n [ i ⊥ k ] f ( n / i ) − ∑ i = 2 n [ i ⊥ k ] f ( n / i ) \large f(n)=\sum\limits_{i=1}^n[i \perp k] f(n/i)-\sum\limits_{i=2}^n[i\perp k] f(n/i) f(n)=i=1∑n[i⊥k]f(n/i)−i=2∑n[i⊥k]f(n/i)
考虑前面那一部分怎么算？
∑
i
=
1
n
[
i
⊥
k
]
∑
j
=
1
n
/
i
μ
(
j
)
[
j
⊥
k
]
\sum\limits_{i=1}^n[i \perp k] \sum\limits_{j=1}^{n/i}\mu(j)[j\perp k]
i=1∑n[i⊥k]j=1∑n/iμ(j)[j⊥k]
=
∑
i
=
1
n
∑
j
=
1
n
/
i
μ
(
j
)
[
i
j
⊥
k
]
=\sum\limits_{i=1}^n\sum\limits_{j=1}^{n/i}\mu(j)[ij\perp k]
=i=1∑nj=1∑n/iμ(j)[ij⊥k]
考虑枚举 i j ij ij
∑
i
=
1
n
∑
d
∣
i
μ
(
d
)
[
i
⊥
k
]
\sum\limits_{i=1}^n \sum\limits_{d|i}\mu(d)[i\perp k]
i=1∑nd∣i∑μ(d)[i⊥k]
∑
i
=
1
n
[
i
⊥
k
]
∑
d
∣
i
μ
(
d
)
=
1
\sum\limits_{i=1}^n[i\perp k] \sum\limits_{d|i} \mu(d)=1
i=1∑n[i⊥k]d∣i∑μ(d)=1
于是乎就得到了
f ( n ) = 1 − ∑ i = 2 n [ i ⊥ k ] f ( n / i ) \large f(n)=1-\sum\limits_{i=2}^n[i\perp k] f(n/i) f(n)=1−i=2∑n[i⊥k]f(n/i)
而
f
(
n
/
i
)
f(n/i)
f(n/i)前面那个在整除分块的时候就是
g
g
g
大力杜教筛即可
code:
#include<bits/stdc++.h>
#define ll long long
#define N 5000050
using namespace std;
int gcd(int x, int y) {
return y? gcd(y, x % y) : x;
}
int g[N], f[N], mu[N], n, m, k, prime[N], sz, vis[N];
ll G(int n) {
return 1ll * (n / k) * g[k] + g[n % k];
}
const int lim = 5000000;
void init(int n) {
for(int i = 1; i <= k; i ++) g[i] = g[i - 1] + (gcd(i, k) == 1);
f[1] = mu[1] = 1;
for(int i = 2; i <= n; i ++) {
if(!vis[i]) {
prime[++ sz] = i;
mu[i] = -1;
}
for(int j = 1; j <= sz && i * prime[j] <= n; j ++) {
vis[i * prime[j]] = 1;
if(i % prime[j] == 0) break;
mu[i * prime[j]] = - mu[i];
}
f[i] = f[i - 1] + mu[i] * (G(i) - G(i - 1));
}
}
map<int, ll> sf;
ll F(int n) {
if(n <= lim) return f[n];
if(sf[n]) return sf[n];
ll ret = 1;
for(int l = 2, r = 0; l <= n; l = r + 1) {
r = n / (n / l);
ret -= F(n / l) * (G(r) - G(l - 1));
} return sf[n] = ret;
}
int main() {
scanf("%d%d%d", &n, &m, &k);
init(lim);
ll ans = 0;
for(int l = 1, r = 0; l <= min(n, m); l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans += 1ll * (n / l) * G(m / l) * (F(r) - F(l - 1));
}
printf("%lld", ans);
return 0;
}