直接反演一下:
∑
i
=
1
n
∑
i
=
1
n
f
(
g
c
d
(
i
,
j
)
)
k
\sum_{i = 1}^n\sum_{i = 1}^nf(gcd(i,j))^k
i=1∑ni=1∑nf(gcd(i,j))k
=
∑
d
=
1
n
∑
i
=
1
⌊
n
d
⌋
∑
i
=
1
⌊
n
d
⌋
f
(
d
)
k
∗
[
g
c
d
(
i
,
j
)
=
1
]
=\sum_{d = 1}^n\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}f(d)^k * [gcd(i,j) = 1]
=d=1∑ni=1∑⌊dn⌋i=1∑⌊dn⌋f(d)k∗[gcd(i,j)=1]
=
∑
d
=
1
n
f
(
d
)
k
∑
p
=
1
⌊
n
d
⌋
μ
(
p
)
∗
⌊
n
p
∗
d
⌋
2
=\sum_{d = 1}^nf(d)^k\sum_{p = 1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)*{\lfloor\frac{n}{p*d}\rfloor}^2
=d=1∑nf(d)kp=1∑⌊dn⌋μ(p)∗⌊p∗dn⌋2
=
∑
T
=
1
n
⌊
n
T
⌋
2
∑
d
∣
T
f
(
d
)
k
μ
(
T
d
)
=\sum_{T = 1}^n{\lfloor\frac{n}{T}\rfloor}^2\sum_{d | T}f(d)^k\mu(\frac{T}{d})
=T=1∑n⌊Tn⌋2d∣T∑f(d)kμ(dT)这个式子可以分块,需要预处理右半边前缀和。
令
g
(
T
)
=
∑
d
∣
T
f
(
d
)
k
μ
(
T
d
)
g(T) = \sum_{d | T}f(d)^k\mu(\frac{T}{d})
g(T)=∑d∣Tf(d)kμ(dT),则
g
(
T
)
=
f
k
∗
μ
g(T) = f^k * \mu
g(T)=fk∗μ,
∗
*
∗代表卷积
按杜教筛的套路,将
g
g
g函数卷上
I
I
I 函数进行求和,令
h
=
g
∗
I
=
f
k
h = g * I =f^k
h=g∗I=fk
∑
i
=
1
n
h
(
i
)
\sum_{i = 1}^nh(i)
i=1∑nh(i)
=
∑
i
=
1
n
∑
d
∣
i
I
(
d
)
∗
g
(
i
d
)
= \sum_{i = 1}^n\sum_{d | i}I(d) * g(\frac{i}{d})
=i=1∑nd∣i∑I(d)∗g(di)
=
∑
d
=
1
n
I
(
d
)
∑
i
=
1
⌊
n
d
⌋
g
(
i
)
=\sum_{d = 1}^nI(d)\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}g(i)
=d=1∑nI(d)i=1∑⌊dn⌋g(i)
=
∑
d
=
1
n
I
(
d
)
S
(
⌊
n
d
⌋
)
=
∑
d
=
1
n
S
(
⌊
n
d
⌋
)
=\sum_{d = 1}^nI(d)S({\lfloor\frac{n}{d}\rfloor})=\sum_{d = 1}^nS({\lfloor\frac{n}{d}\rfloor})
=d=1∑nI(d)S(⌊dn⌋)=d=1∑nS(⌊dn⌋)提出第一项
I
(
1
)
S
(
n
)
=
S
(
n
)
I(1)S(n) = S(n)
I(1)S(n)=S(n) 并移项:
S
(
n
)
=
∑
i
=
1
n
h
(
i
)
−
∑
i
=
2
n
S
(
⌊
n
i
⌋
)
S(n) = \sum_{i = 1}^nh(i) - \sum_{i = 2}^nS(\lfloor\frac{n}{i}\rfloor)
S(n)=i=1∑nh(i)−i=2∑nS(⌊in⌋)
S
(
n
)
=
∑
i
=
1
n
f
(
i
)
k
−
∑
i
=
2
n
S
(
⌊
n
i
⌋
)
S(n) = \sum_{i = 1}^nf(i)^k - \sum_{i = 2}^nS(\lfloor\frac{n}{i}\rfloor)
S(n)=i=1∑nf(i)k−i=2∑nS(⌊in⌋)
f
(
i
)
k
f(i)^k
f(i)k可以用min_25筛筛出来,筛法和Uoj 188相同,但是这题
f
(
i
)
f(i)
f(i)质数部分贡献不为0,先筛出合数部分再加上质数部分即可。由于递归版没有记忆化,杜教筛里每次都要计算会T,因此需要用递推版,并且由于模数是
2
32
2^{32}
232,可以用 unsigned int 类型自然溢出来取代取模运算减少常数。
代码:
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 2e5 + 10;
const ll mod = 1ll << 32;
bool ispri[maxn];
ll pri[maxn], num, tot;
ll pw[maxn];
ll g[maxn], id1[maxn], id2[maxn], w[maxn], sqr;
ll sum[maxn];
ll n, k;
ll mp[maxn];
ll fpow(ll a, ll b) {
ll r = 1;
while (b) {
if (b & 1)
r = r * a % mod;
a = a * a % mod;
b >>= 1;
}
return r;
}
void sieve(int n) {
ispri[0] = ispri[1] = true;
num = 0;
for (int i = 2; i <= n; i++) {
if (!ispri[i])
pri[++num] = i;
for (int j = 1; j <= num && i * pri[j] <= n; j++) {
ispri[i * pri[j]] = true;
if (i % pri[j] == 0)
break;
}
}
}
/*ll S(ll x, ll y) { //递归版,T飞
if (pri[y] > x)
return 0;
ll z = x <= sqr ? id1[x] : id2[n / x];
ll res = 0;
for (ll i = y + 1; i <= num && pri[i] * pri[i] <= x; i++) {
for (ll e = 1, pe = pri[i]; pe * pri[i] <= x; e++, pe *= pri[i]) {
ll t = x / pe;
ll k = t <= sqr ? id1[t] : id2[n / t];
res += (S(t, i) + (g[k] - i + 1) % mod * pw[i] % mod) % mod;
res %= mod;
}
}
return res;
}*/
ll getsum(ll x) { //杜教筛
ll t = x <= sqr ? id1[x] : id2[n / x];
if (mp[t] != -1)
return mp[t];
ll res = (sum[t] + g[t]) % mod;
for (ll i = 2, j; i <= x; i = j + 1) {
j = x / (x / i);
res -= getsum(x / i) * (j - i + 1) % mod;
if (res < 0)
res += mod;
}
return mp[t] = res;
}
int main() {
sieve(maxn - 10);
scanf("%lld%lld", &n, &k);
memset(mp, -1, sizeof mp);
for (int i = 1; i <= num; i++) {
pw[i] = fpow(pri[i], k);
}
sqr = sqrt(n);
for (ll i = 1, j; i <= n; i = j + 1) { //min_25筛离散化
j = n / (n / i);
w[++tot] = n / i;
g[tot] = (w[tot] - 1) % mod;
if (n / i <= sqr)
id1[n / i] = tot;
else
id2[j] = tot;
}
for (int i = 1; i <= num; i++) { //min_25筛预处理 g
for (int j = 1; j <= tot && pri[i] * pri[i] <= w[j]; j++) {
ll t = w[j] / pri[i];
ll k = t <= sqr ? id1[t] : id2[n / t];
g[j] -= (g[k] - i + 1) % mod;
if (g[j] < 0)
g[j] += mod;
}
}
for (int i = num; i >= 1; i--) { //min_25筛 递推版
for (int j = 1; j <= tot && pri[i] * pri[i] <= w[j]; j++) {
for (ll e = 1, pe = pri[i]; pe * pri[i] <= w[j]; pe = pe * pri[i]) {
ll t = w[j] / pe;
ll k = t <= sqr ? id1[t] : id2[n / t];
sum[j] += (sum[k] + (g[k] - i + 1) % mod * pw[i] % mod) % mod;
sum[j] %= mod;
}
}
}
ll ans = 0;
for (ll i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
ll p = n / i % mod;
ans += (getsum(j) - getsum(i - 1) + mod) % mod * p * p % mod;
ans %= mod;
}
printf("%lld\n", ans);
return 0;
}