EOJ Monthly 2019.11
∑
i
=
1
n
∑
a
1
=
1
i
∑
a
2
=
1
i
∑
a
3
=
1
i
⋯
∑
a
k
−
1
i
∑
a
k
i
[
g
c
d
(
a
1
,
a
2
,
a
3
,
…
,
a
k
−
1
,
a
k
,
i
)
=
=
1
]
=
∑
i
=
1
n
∑
d
∣
i
μ
(
d
)
⌊
i
d
⌋
k
=
∑
d
=
1
n
μ
(
d
)
∑
d
∣
i
⌊
i
d
⌋
k
=
∑
d
=
1
n
μ
(
d
)
∑
t
=
1
n
d
t
k
\sum_{i = 1} ^{n} \sum_{a_1 = 1} ^{i} \sum_{a_2 = 1} ^{i} \sum_{a_3 = 1} ^{i} \dots \sum_{a_{k - 1}} ^{i} \sum_{a_k} ^{i} [gcd(a_1, a_2, a_3, \dots, a_{k - 1}, a_{k}, i) == 1]\\ = \sum_{i = 1} ^{n} \sum_{d \mid i} \mu(d) \lfloor \frac{i}{d} \rfloor ^ k\\ = \sum_{d = 1} ^{n} \mu(d) \sum_{d \mid i} \lfloor \frac{i}{d}\rfloor ^ k\\ = \sum_{d = 1} ^{n} \mu(d) \sum_{t = 1} ^{\frac{n}{d}} t ^ k\\
i=1∑na1=1∑ia2=1∑ia3=1∑i⋯ak−1∑iak∑i[gcd(a1,a2,a3,…,ak−1,ak,i)==1]=i=1∑nd∣i∑μ(d)⌊di⌋k=d=1∑nμ(d)d∣i∑⌊di⌋k=d=1∑nμ(d)t=1∑dntk
然后杜教筛筛出
∑
i
=
1
n
μ
(
i
)
\sum\limits_{i = 1} ^{n} \mu(i)
i=1∑nμ(i)的前缀和，用拉格朗日插值得到
∑
i
=
1
n
i
k
\sum\limits_{i = 1} ^{n} i ^ k
i=1∑nik这个式子，再加上数论分块即可完美解决。
代码
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>
using namespace std
;
typedef long long ll
;
const int inf
= 0x3f3f3f3f;
const double eps
= 1e-7;
const int N
= 1e6 + 10, mod
= 998244353;
ll fac
[N
], pre
[N
], suc
[N
], inv
[N
], prime
[N
], sum
[N
], mu
[N
], n
, k
, cnt
;
bool st
[N
];
ll
quick_pow(ll a
, int n
) {
ll ans
= 1;
while(n
) {
if(n
& 1) ans
= ans
* a
% mod
;
a
= a
* a
% mod
;
n
>>= 1;
}
return ans
;
}
void init() {
sum
[1] = mu
[1] = 1;
for(int i
= 2; i
< N
; i
++) {
if(!st
[i
]) {
prime
[cnt
++] = i
;
mu
[i
] = -1;
sum
[i
] = quick_pow(i
, k
);
}
for(int j
= 0; j
< cnt
&& i
* prime
[j
] < N
; j
++) {
st
[i
* prime
[j
]] = 1;
sum
[i
* prime
[j
]] = 1ll * sum
[i
] * sum
[prime
[j
]] % mod
;
if(i
% prime
[j
] == 0) break;
mu
[i
* prime
[j
]] = -mu
[i
];
}
}
fac
[0] = inv
[0] = 1;
for(int i
= 1; i
< N
; i
++) {
mu
[i
] = (mu
[i
] + mu
[i
- 1]) % mod
;
sum
[i
] = (sum
[i
] + sum
[i
- 1]) % mod
;
fac
[i
] = 1ll * fac
[i
- 1] * i
% mod
;
}
inv
[N
- 1] = quick_pow(fac
[N
- 1], mod
- 2);
for(int i
= N
- 2; i
>= 1; i
--) {
inv
[i
] = 1ll * inv
[i
+ 1] * (i
+ 1) % mod
;
}
}
ll
solve(ll n
) {
ll ans
= 0;
pre
[0] = suc
[k
+ 3] = 1;
for(int i
= 1; i
<= k
+ 2; i
++) pre
[i
] = 1ll * pre
[i
- 1] * (n
- i
) % mod
;
for(int i
= k
+ 2; i
>= 1; i
--) suc
[i
] = 1ll * suc
[i
+ 1] * (n
- i
) % mod
;
for(int i
= 1; i
<= k
+ 2; i
++) {
ll a
= 1ll * pre
[i
- 1] * suc
[i
+ 1] % mod
, b
= 1ll * inv
[i
- 1] * inv
[k
+ 2 - i
] % mod
;
if((k
+ 2 - i
) & 1) b
*= -1;
ans
= ((ans
+ 1ll * sum
[i
] * a
% mod
* b
% mod
) % mod
+ mod
) % mod
;
}
return ans
;
}
unordered_map
<ll
, ll
> ans_s
;
ll
S(ll n
) {
if(n
< N
) return mu
[n
];
if(ans_s
.count(n
)) return ans_s
[n
];
ll ans
= 1;
for(ll l
= 2, r
; l
<= n
; l
= r
+ 1) {
r
= n
/ (n
/ l
);
ans
= ((ans
- 1ll *(r
- l
+ 1) * S(n
/ l
)) % mod
+ mod
) % mod
;
}
return ans_s
[n
] = ans
;
}
int main() {
scanf("%lld %lld", &n
, &k
);
init();
ll ans
= 0;
for(ll l
= 1, r
; l
<= n
; l
= r
+ 1) {
r
= n
/ (n
/ l
);
ans
= (ans
+ 1ll * ((S(r
) - S(l
- 1)) % mod
+ mod
) % mod
* solve(n
/ l
) % mod
) % mod
;
}
printf("%lld\n", ans
);
return 0;
}
