首先
d
(
i
∗
j
)
=
∑
x
∣
i
∑
y
∣
j
[
g
c
d
(
x
,
y
)
=
1
]
d(i * j) = \sum_{x | i}\sum_{y|j}[gcd(x,y) = 1]
d(i∗j)=∑x∣i∑y∣j[gcd(x,y)=1],证明感性理解一下就好。然后整个式子变成:
∑
i
=
1
N
∑
j
=
1
M
∑
x
∣
i
∑
y
∣
j
[
g
c
d
(
x
,
y
)
=
1
]
\sum_{i = 1}^N\sum_{j = 1}^M\sum_{x | i}\sum_{y|j}[gcd(x,y) = 1]
i=1∑Nj=1∑Mx∣i∑y∣j∑[gcd(x,y)=1]
先改变一下枚举项,枚举因子 x,y:
∑
x
=
1
N
∑
y
=
1
M
∑
i
=
1
⌊
N
x
⌋
∑
j
=
1
⌊
M
y
⌋
[
g
c
d
(
x
,
y
)
=
1
]
\sum_{x = 1}^N\sum_{y = 1}^M\sum_{i = 1}^{\lfloor\frac{N}{x}\rfloor}\sum_{j = 1}^{\lfloor\frac{M}{y}\rfloor}[gcd(x,y) = 1]
x=1∑Ny=1∑Mi=1∑⌊xN⌋j=1∑⌊yM⌋[gcd(x,y)=1]
然后用
∑
d
∣
g
c
d
(
i
,
j
)
μ
(
d
)
\sum_{d | gcd(i,j)}\mu(d)
∑d∣gcd(i,j)μ(d)替换掉
[
g
c
d
(
i
,
j
)
=
1
]
[gcd(i,j) = 1]
[gcd(i,j)=1]:
∑
x
=
1
N
∑
y
=
1
M
∑
i
=
1
⌊
N
x
⌋
∑
j
=
1
⌊
M
y
⌋
∑
d
∣
g
c
d
(
x
,
y
)
μ
(
d
)
\sum_{x = 1}^N\sum_{y = 1}^M\sum_{i = 1}^{\lfloor\frac{N}{x}\rfloor}\sum_{j = 1}^{\lfloor\frac{M}{y}\rfloor}\sum_{d | gcd(x,y)}\mu(d)
x=1∑Ny=1∑Mi=1∑⌊xN⌋j=1∑⌊yM⌋d∣gcd(x,y)∑μ(d)
再次改变枚举项,枚举 d,将 d 提取到前面来:
∑
d
=
1
N
μ
(
d
)
∑
x
=
1
⌊
N
d
⌋
∑
y
=
1
⌊
M
d
⌋
∑
i
=
1
⌊
N
x
∗
d
⌋
∑
j
=
1
⌊
M
y
∗
d
⌋
\sum_{d = 1}^N\mu(d)\sum_{x = 1}^{\lfloor\frac{N}{d}\rfloor}\sum_{y = 1}^{\lfloor\frac{M}{d}\rfloor}\sum_{i = 1}^{\lfloor\frac{N}{x*d}\rfloor}\sum_{j = 1}^{\lfloor\frac{M}{y*d}\rfloor}
d=1∑Nμ(d)x=1∑⌊dN⌋y=1∑⌊dM⌋i=1∑⌊x∗dN⌋j=1∑⌊y∗dM⌋
把后面两个求和去掉:
∑
d
=
1
N
μ
(
d
)
∑
x
=
1
⌊
N
d
⌋
⌊
N
x
∗
d
⌋
∑
y
=
1
⌊
M
d
⌋
⌊
M
y
∗
d
⌋
\sum_{d = 1}^N\mu(d)\sum_{x = 1}^{\lfloor\frac{N}{d}\rfloor}{\lfloor\frac{N}{x*d}\rfloor}\sum_{y = 1}^{\lfloor\frac{M}{d}\rfloor}{\lfloor\frac{M}{y*d}\rfloor}
d=1∑Nμ(d)x=1∑⌊dN⌋⌊x∗dN⌋y=1∑⌊dM⌋⌊y∗dM⌋
到这一步就可以做了,后面两个可以数论分块,
O
(
n
∗
n
)
O(n * \sqrt n)
O(n∗n)预处理一下
∑
i
=
1
p
⌊
p
i
⌋
\sum_{i = 1}^{p}{\lfloor\frac{p}{i}\rfloor}
∑i=1p⌊ip⌋就可以O(1)算出后面一坨,再预处理一下莫比乌斯函数前缀和,就可以
O
(
n
)
O(\sqrt n)
O(n)处理单组,多组复杂度为:
O
(
T
∗
n
)
O(T * \sqrt n)
O(T∗n)
另外一种推法:我们要求的是
∑
i
=
1
N
∑
j
=
1
M
∑
x
∣
i
∑
y
∣
j
[
g
c
d
(
x
,
y
)
=
1
]
\sum_{i = 1}^N\sum_{j = 1}^M\sum_{x | i}\sum_{y|j}[gcd(x,y) = 1]
i=1∑Nj=1∑Mx∣i∑y∣j∑[gcd(x,y)=1]
改变枚举项,枚举 x,y:
∑
x
=
1
N
∑
y
=
1
M
∑
i
=
1
⌊
N
x
⌋
∑
j
=
1
⌊
M
y
⌋
[
g
c
d
(
x
,
y
)
=
1
]
\sum_{x = 1}^N\sum_{y = 1}^M\sum_{i = 1}^{\lfloor\frac{N}{x}\rfloor}\sum_{j = 1}^{\lfloor\frac{M}{y}\rfloor}[gcd(x,y) = 1]
x=1∑Ny=1∑Mi=1∑⌊xN⌋j=1∑⌊yM⌋[gcd(x,y)=1]
接下来是不同的地方:设
f
(
d
)
=
∑
x
=
1
N
∑
y
=
1
M
∑
i
=
1
⌊
N
x
⌋
∑
j
=
1
⌊
M
y
⌋
[
g
c
d
(
x
,
y
)
=
d
]
f(d) = \sum_{x = 1}^N\sum_{y = 1}^M\sum_{i = 1}^{\lfloor\frac{N}{x}\rfloor}\sum_{j = 1}^{\lfloor\frac{M}{y}\rfloor}[gcd(x,y) = d]
f(d)=x=1∑Ny=1∑Mi=1∑⌊xN⌋j=1∑⌊yM⌋[gcd(x,y)=d]
F
(
p
)
=
∑
p
∣
d
f
(
d
)
=
∑
x
=
1
N
∑
y
=
1
M
∑
i
=
1
⌊
N
x
⌋
∑
j
=
1
⌊
M
y
⌋
[
p
∣
g
c
d
(
x
,
y
)
]
F(p) = \sum_{p | d}f(d) = \sum_{x = 1}^N\sum_{y = 1}^M\sum_{i = 1}^{\lfloor\frac{N}{x}\rfloor}\sum_{j = 1}^{\lfloor\frac{M}{y}\rfloor}[p | gcd(x,y)]
F(p)=p∣d∑f(d)=x=1∑Ny=1∑Mi=1∑⌊xN⌋j=1∑⌊yM⌋[p∣gcd(x,y)]
把 p 提出来:枚举 x,y是p的倍数,i,j 是
p
∗
x
,
p
∗
y
p*x,p * y
p∗x,p∗y的倍数(因为i,j是 x,y的倍数,x,y现在是p的倍数,所以i,j是
p
∗
x
,
p
∗
y
p*x,p*y
p∗x,p∗y的倍数):
F
(
p
)
=
∑
p
∣
d
f
(
d
)
=
∑
x
=
1
⌊
N
p
⌋
∑
y
=
1
⌊
M
p
⌋
∑
i
=
1
⌊
N
x
∗
p
⌋
∑
j
=
1
⌊
M
y
∗
p
⌋
1
=
∑
x
=
1
⌊
N
p
⌋
∑
y
=
1
⌊
M
p
⌋
⌊
N
x
∗
p
⌋
∗
⌊
M
y
∗
p
⌋
F(p) = \sum_{p | d}f(d) = \sum_{x = 1}^{\lfloor\frac{N}{p}\rfloor}\sum_{y = 1}^{\lfloor\frac{M}{p}\rfloor}\sum_{i = 1}^{\lfloor\frac{N}{x * p}\rfloor}\sum_{j = 1}^{\lfloor\frac{M}{y * p}\rfloor} 1 = \sum_{x = 1}^{\lfloor\frac{N}{p}\rfloor}\sum_{y = 1}^{\lfloor\frac{M}{p}\rfloor}{\lfloor\frac{N}{x * p}\rfloor}*{\lfloor\frac{M}{y * p}\rfloor}
F(p)=p∣d∑f(d)=x=1∑⌊pN⌋y=1∑⌊pM⌋i=1∑⌊x∗pN⌋j=1∑⌊y∗pM⌋1=x=1∑⌊pN⌋y=1∑⌊pM⌋⌊x∗pN⌋∗⌊y∗pM⌋
反演一下:
f
(
d
)
=
∑
d
∣
p
μ
(
p
d
)
∗
F
(
p
)
f(d) = \sum_{d | p}\mu(\frac{p}{d}) * F(p)
f(d)=d∣p∑μ(dp)∗F(p)(这里p 最多枚举到 n)
我们要求的是:
f
(
1
)
=
∑
p
=
1
n
μ
(
p
)
∗
F
(
p
)
=
∑
p
=
1
n
μ
(
p
)
∗
∑
x
=
1
⌊
N
p
⌋
∑
y
=
1
⌊
M
p
⌋
⌊
N
x
∗
p
⌋
∗
⌊
M
y
∗
p
⌋
f(1) = \sum_{p = 1}^n\mu({p}) * F(p) =\sum_{p = 1}^n\mu({p}) * \sum_{x = 1}^{\lfloor\frac{N}{p}\rfloor}\sum_{y = 1}^{\lfloor\frac{M}{p}\rfloor}{\lfloor\frac{N}{x * p}\rfloor}*{\lfloor\frac{M}{y * p}\rfloor}
f(1)=p=1∑nμ(p)∗F(p)=p=1∑nμ(p)∗x=1∑⌊pN⌋y=1∑⌊pM⌋⌊x∗pN⌋∗⌊y∗pM⌋
这里已经上面一样,可以分块做了
代码:
#include<bits/stdc++.h>
using namespace std;
const int maxn = 1e7 + 10;
int t,n,m;
typedef long long ll;
bool ispri[maxn];
int pri[maxn],mu[maxn];
ll sum[maxn];
void sieve(int n) {
ispri[0] = ispri[1] = true;
pri[0] = 0;
mu[1] = 1;
for(int i = 2; i <= n; i++) {
if(!ispri[i]) pri[++pri[0]] = i,mu[i] = -1;
for(int j = 1; j <= pri[0] && i * pri[j] <= n; j++) {
ispri[i * pri[j]] = true;
if(i % pri[j] == 0) break;
mu[i * pri[j]] = -mu[i];
}
}
for(int i = 1; i <= n; i++)
mu[i] += mu[i - 1];
}
void prework(int n) {
int i,j;
for(i = 1; i <= n; i = j + 1) {
j = n / (n / i);
sum[n] += 1ll * (j - i + 1) * (n / i);
}
}
int main() {
scanf("%d",&t);
sieve(100000);
for(int i = 1; i <= 50000; i++)
prework(i);
while(t--) {
scanf("%d%d",&n,&m);
if(n > m) swap(n,m);
int i,j;
ll res = 0;
for(i = 1; i <= n; i = j + 1) {
j = min(n / (n / i),m / (m / i));
int p1 = n / i,p2 = m / i;
res += 1ll * (mu[j] - mu[i - 1]) * (sum[p1] * sum[p2]);
}
printf("%lld\n",res);
}
return 0;
}