莫比乌斯反演题解汇总

Foreword

题目 + 题解

$(i, j) \Leftrightarrow \gcd (i, j)$

[SDOI2017] 数字表格

Description

$$
T \le 10 ^3, 1 \le n, m \le 10 ^6
$$

Solution

不妨设 $n \le m$

$$
\begin{aligned}
& \prod _{i = 1} ^n \prod _{j = 1} ^m f((i, j)) \\
=& \prod _{d = 1} ^n \prod _{i = 1} ^n \prod _{j = 1} ^m [(i, j) = d] f(d) \\
=& \prod _{d = 1} ^n f(d) ^{\sum _{i = 1} ^n \sum _{j = 1} ^m [(i, j) = d]} \\
\end{aligned}
$$

指数部分提出来 单独计算

$$
\begin{aligned}
& \sum _{i = 1} ^n \sum _{j = 1} ^m [(i, j) = d] \\
=& \sum _{i = 1} ^{\lfloor \frac n d \rfloor} \sum _{j = 1} ^{\lfloor \frac m d \rfloor} [(i, j) = 1] \\
=& \sum _{i = 1} ^{\lfloor \frac n d \rfloor} \sum _{j = 1} ^{\lfloor \frac m d \rfloor} \sum _{x | (i, j)} \mu (x) \\
=& \sum _{x = 1} ^{\lfloor \frac n d \rfloor} \mu (x) \lfloor \frac n {dx} \rfloor \lfloor \frac m {dx} \rfloor
\end{aligned}
$$

考虑先枚举 $dx$ 的乘积再枚举 $d$ ,可以证明这样得出的答案不变(简易理解就是可以遍历到所有可能的 $d, x$ 的取值集合,不重不漏)

这应该也算是一种套路了

令 $T = dx$ ,$x$ 可用 $\frac T d$ 代替
$$
\begin{aligned}
& \prod _{T = 1} ^n \prod _{d | T} f(d) ^{\mu (\frac T d) \lfloor \frac n T \rfloor \lfloor \frac m T \rfloor} \\
=& \prod _{T = 1} ^n \left( \prod _{d | T} f(d) ^{\mu (\frac T d) } \right) ^{ \lfloor \frac n T \rfloor \lfloor \frac m T \rfloor}
\end{aligned}
$$
预处理 斐波那契数列及其逆元 和 里面东西的前缀积,筛出 $\mu$ 之后直接枚举每个数对其倍数算贡献,复杂度 $O(n (\ln n + \log))$ ( $\log$ 是快速幂带的复杂度,这里 $\mu$ 只有三种取值,预处理逆元后就不需要快速幂了)

每次询问数论分块 + 快速幂,复杂度 $O(T (\sqrt n + \sqrt m) \log)$ ,其中 $T$ 是数据组数(注意逆元)

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
/**********************************************************
* Author : EndSaH
* Email : hjxhb1@gmail.com
* Created Time : 2019-05-31 21:11
* FileName : SDOI2017_shuzibiaoge.cpp
* Website : https://endsah.cf
* *******************************************************/

#include <cstdio>
#include <cctype>
#include <cmath>
#include <vector>
#include <bitset>

typedef long long LL;

#define debug(...) fprintf(stderr, __VA_ARGS__)
#define Debug(s) debug("The message in line %d, Function %s: %s\n", __LINE__, __FUNCTION__, s)
#define getchar() (ipos == iend and (iend = (ipos = _ibuf) + fread(_ibuf, 1, __bufsize, stdin), ipos == iend) ? EOF : *ipos++)
#define putchar(ch) (opos == oend ? fwrite(_obuf, 1, __bufsize, stdout), opos = _obuf : 0, *opos++ = (ch))
#define __bufsize (1 << 21 | 1)

char _ibuf[__bufsize], _obuf[__bufsize], _stk[20];
char *ipos = _ibuf, *iend = _ibuf, *opos = _obuf, *oend = _obuf + __bufsize, *stkpos = _stk;

struct END
{ ~END() { fwrite(_obuf, 1, opos - _obuf, stdout); } }
__;

inline int read()
{
register int x = 0;
register char ch;
while (!isdigit(ch = getchar()));
while (x = (x << 3) + (x << 1) + (ch & 15), isdigit(ch = getchar()));
return x;
}

template <typename _INT>
inline void write(_INT x)
{
while (*++stkpos = x % 10 ^ 48, x /= 10, x);
while (stkpos != _stk)
putchar(*stkpos--);
}

template <typename _Tp>
inline bool Chkmax(_Tp& x, const _Tp& y)
{ return x < y ? x = y, true : false; }

template <typename _Tp>
inline bool Chkmin(_Tp& x, const _Tp& y)
{ return x > y ? x = y, true : false; }

const int maxN = 1e6 + 2;
const int mod = 1e9 + 7;

int T, n, m;
int mu[maxN], fib[maxN], fibinv[maxN], F[maxN];
std::bitset<maxN> vis;
std::vector<int> prim;

inline void Mod(int& x)
{ x >= mod ? x -= mod : 0; }

inline int QP(int a, LL n)
{
if (n == -1)
return QP(a, mod - 2);
if (a == 1)
return 1;
register int res = 1;
while (n)
{
if (n & 1)
res = (LL)res * a % mod;
a = (LL)a * a % mod;
n >>= 1;
}
return res;
}

inline void Pre(int n)
{
fib[1] = fibinv[1] = mu[1] = F[1] = F[0] = 1;
for (register int i = 2; i <= n; ++i)
{
F[i] = 1, Mod(fib[i] = fib[i - 1] + fib[i - 2]);
fibinv[i] = QP(fib[i], -1);
if (!vis[i])
prim.push_back(i), mu[i] = -1;
for (register auto j : prim)
{
if (i * j > n)
break;
vis.set(i * j);
if (i % j)
mu[i * j] = -mu[i];
else
break;
}
}
for (register int i = 1; i <= n; ++i)
for (register int j = 1; i * j <= n; ++j)
if (~mu[j])
{
if (mu[j])
F[i * j] = (LL)F[i * j] * fib[i] % mod;
}
else
F[i * j] = (LL)F[i * j] * fibinv[i] % mod;
for (register int i = 2; i <= n; ++i)
F[i] = (LL)F[i] * F[i - 1] % mod;
}

int main()
{
#ifndef ONLINE_JUDGE
freopen("SDOI2017_shuzibiaoge.in", "r", stdin);
freopen("SDOI2017_shuzibiaoge.out", "w", stdout);
#endif
Pre(maxN - 2);
T = read();
while (T--)
{
n = read(), m = read();
if (n > m)
std::swap(n, m);
register int ans = 1;
for (register int l = 1, r; l <= n; l = r + 1)
{
r = std::min(n / (n / l), m / (m / l));
ans = (LL)ans * QP( (LL)F[r] * QP(F[l - 1], -1) % mod, (LL)(n / l) * (m / l) ) % mod;
}
write(ans), putchar('\n');
}
return 0;
}

[SDOI2014] 数表

Description


$$
\sum _{i = 1} ^n \sum _{j = 1} ^m [\sigma _1 ((i, j)) \le a] \sigma _1 ((i, j)) \\
T \le 2 \times 10 ^4, 1 \le n, m \le 10 ^5, a \le 10 ^9
$$
每组数据有不同的 $a$

Solution

不妨设 $n \le m$

先考虑去掉 $\le a$ 的限制怎么做

可以知道 $\sigma _1 ((i, j)) = \sum \limits _{x | (i, j)} x$ ,于是不需要反演,可以直接这样推:

$$
\begin{aligned}
& \sum _{i = 1} ^n \sum _{j = 1} ^m \sigma _1 ((i, j)) \\
=& \sum _{i = 1} ^n \sum _{j = 1} ^m \sum _{x | (i, j)} x \\
=& \sum _{x = 1} ^n x \lfloor \frac n x \rfloor \lfloor \frac m x \rfloor
\end{aligned}
$$

然后就无路可走了

于是找另一条路走:

$$
\begin{aligned}
& \sum _{i = 1} ^n \sum _{j = 1} ^m \sigma _1 ((i, j)) \\
=& \sum _{d = 1} ^n \sigma _1 (d) \sum _{i = 1} ^n \sum _{j = 1} ^m [(i, j) = d] \\
=& \sum _{d = 1} ^n \sigma _1 (d) \sum _{x = 1} ^{\lfloor \frac n d \rfloor} \mu (x) \lfloor \frac n {dx} \rfloor \lfloor \frac m {dx} \rfloor
\end{aligned}
$$

(第二行最后那块实在是太常见了,就直接变了)

依然枚举 $T = dx$ 变为最终式:
$$
\sum _{T = 1} ^n \lfloor \frac n T \rfloor \lfloor \frac m T \rfloor \sum _{d | T} \sigma _1 (d) \mu (\frac T d)
$$
再回到题目的 $\le a$ 的限制

惊喜的发现这个式子让这个题目变得可做了

将后面这一块单独提出来看:令 $g(T) = \sum \limits _{d | T} \sigma _1 (d) \mu (\frac T d)$

我们一开始枚举的 $d$ 的意义就是 $i, j$ 的 $\gcd$ ,中间就算经过了变换,它的意义依然没有改变;也就是说,$\sigma _1 ((i, j)) \le a$ 也对应着上面 $g(T)$ 中所枚举的 $d$ 需要满足 $\sigma _1 (d) \le a$

于是我们离线将 $a$ 排序,同时将 $10 ^5$ 内的数线性筛出其约数和并按约数和排序

每次询问将满足约数和小于 $a$ 的数对其倍数算贡献,也就是动态的更新 $g(T)$ ,考虑到整除分块时需要查询前缀和,而修改只需要单点加,于是便用树状数组动态维护 $g(T)$

对 $2 ^{31}$ 取模,只需要算的时候用 unsigned int ,将答案再与 $2 ^{31} - 1$ 相与即可(注意 (1 << 31) - 1 此处会爆 int ,因为 $2 ^{31}$ 超出了 int 范围,故用 0x7FFFFFFF 为宜)

复杂度的话,视作 $n, m$ 同阶,一开始排序和线性筛即为 $O(T \log T + n \log n)$ ,因为对约数和和询问排了序;处理询问时,修改处枚举倍数和树状数组有两个 $\log$ (其中有一个是 $\ln$ ,粗略算成 $\log$ 算了),回答用到整除分块和树状数组,综上,复杂度约为 $O(T \sqrt n \log n + n \log ^2 n + T \log T)$

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
/**********************************************************
* Author : EndSaH
* Email : hjxhb1@gmail.com
* Created Time : 2019-06-04 19:13
* FileName : SDOI2014_number.cpp
* Website : https://endsah.cf
* *******************************************************/

#include <cstdio>
#include <cctype>
#include <cmath>
#include <bitset>
#include <vector>
#include <algorithm>

typedef std::pair<int, int> pii;
typedef unsigned int uint;

#define debug(...) fprintf(stderr, __VA_ARGS__)
#define Debug(s) debug("The message in line %d, Function %s: %s\n", __LINE__, __FUNCTION__, s)
#define getchar() (ipos == iend and (iend = (ipos = _ibuf) + fread(_ibuf, 1, __bufsize, stdin), ipos == iend) ? EOF : *ipos++)
#define putchar(ch) (opos == oend ? fwrite(_obuf, 1, __bufsize, stdout), opos = _obuf : 0, *opos++ = (ch))
#define __bufsize (1 << 21 | 1)

char _ibuf[__bufsize], _obuf[__bufsize], _stk[20];
char *ipos = _ibuf, *iend = _ibuf, *opos = _obuf, *oend = _obuf + __bufsize, *stkpos = _stk;

struct END
{ ~END() { fwrite(_obuf, 1, opos - _obuf, stdout); } }
__;

inline int read()
{
register int x = 0;
register char ch;
while (!isdigit(ch = getchar()));
while (x = (x << 3) + (x << 1) + (ch & 15), isdigit(ch = getchar()));
return x;
}

template <typename _INT>
inline void write(_INT x)
{
while (*++stkpos = x % 10 ^ 48, x /= 10, x);
while (stkpos != _stk)
putchar(*stkpos--);
}

template <typename _Tp>
inline bool Chkmax(_Tp& x, const _Tp& y)
{ return x < y ? x = y, true : false; }

template <typename _Tp>
inline bool Chkmin(_Tp& x, const _Tp& y)
{ return x > y ? x = y, true : false; }

const int maxN = 1e5;
const int maxQ = 2e4;

int T;
int ans[maxN + 2], low[maxN + 2], mu[maxN + 2];
uint BIT[maxN + 2];
std::bitset<maxN + 2> vis;
std::vector<int> prim;
pii sigma[maxN + 2];

struct Ask
{
int id, n, m, a;

Ask() { }

Ask(int id, int n, int m, int a) : id(id), n(n), m(m), a(a) { }

bool operator<(const Ask& x) const
{ return a < x.a; }
} ask[maxQ + 2];

inline void Add(int x, int addval)
{
while (x <= maxN)
BIT[x] += addval, x += x & -x;
}

inline uint Query(int x)
{
register uint res = 0;
while (x)
res += BIT[x], x -= x & -x;
return res;
}

inline void Shuffle(int n)
{
prim.reserve(n / log(n));
sigma[1] = pii(1, 1), mu[1] = low[1] = 1;
for (register int i = 2; i <= n; ++i)
{
if (!vis[i])
{
prim.push_back(i);
mu[i] = -1, low[i] = i, sigma[i] = pii(i + 1, i);
}
for (register auto j : prim)
{
if (i * j > n)
break;
vis.set(i * j);
if (i % j)
low[i * j] = j, mu[i * j] = -mu[i], sigma[i * j] = pii(sigma[i].first * sigma[j].first, i * j);
else
{
low[i * j] = low[i] * j;
if (low[i] == i)
sigma[i * j] = pii(sigma[i].first + i * j, i * j);
else
sigma[i * j] = pii(sigma[i / low[i]].first * sigma[j * low[i]].first, i * j);
break;
}
}
}
std::sort(sigma + 1, sigma + n + 1);
}

int main()
{
#ifndef ONLINE_JUDGE
freopen("SDOI2014_number.in", "r", stdin);
freopen("SDOI2014_number.out", "w", stdout);
#endif
Shuffle(maxN);
T = read();
for (register int i = 1, n, m, a; i <= T; ++i)
{
n = read(), m = read(), a = read();
if (n > m)
std::swap(n, m);
ask[i] = Ask(i, n, m, a);
}
std::sort(ask + 1, ask + T + 1);
register pii* j = sigma + 1;
for (register Ask* p = ask + 1; p - ask <= T; ++p)
{
while (j - sigma <= maxN and j->first <= p->a)
{
for (register int k = 1; k * j->second <= maxN; ++k)
Add(k * j->second, j->first * mu[k]);
++j;
}
uint cans = 0;
for (register int l = 1, r; l <= p->n; l = r + 1)
{
r = std::min(p->n / (p->n / l), p->m / (p->m / l));
cans += (uint)(p->n / l) * (p->m / l) * (Query(r) - Query(l - 1));
}
ans[p->id] = cans & 0x7FFFFFFF;
}
for (register int i = 1; i <= T; ++i)
write(ans[i]), putchar('\n');
return 0;
}

简单的数学题

Description

给定 $n, p$,求

$$
\sum _{i = 1} ^n \sum _{j = 1} ^n i \cdot j \cdot (i, j)
$$

$$
1 \le n \le 10 ^{10}, 5 \times 10 ^8 \le p \le 1.1 \times 10 ^9, p \in \mathbb{P}
$$

Solution

依然先提一个 $(i, j)$ 出来:

$$
\sum _{d = 1} ^n d \sum _{i = 1} ^n \sum _{j = 1} ^n i \cdot j \cdot [(i, j) = d]
$$

$i, j$ 改为枚举 $d$ 的倍数,将 $id \cdot jd$ 中的两个 $d$ 提出去:

$$
\sum _{d = 1} ^n d ^3 \sum _{i = 1} ^{\lfloor \frac n d \rfloor} \sum _{j = 1} ^{\lfloor \frac n d \rfloor} i \cdot j \cdot [(i, j) = 1]
$$

对后面那一块反演一下:

$$
\begin{aligned}
& \sum _{d = 1} ^n d ^3 \sum _{x = 1} ^{\lfloor \frac n d \rfloor} \mu (x) \cdot x ^2 \sum _{i = 1} ^{\lfloor \frac n {dx} \rfloor} i \cdot \sum _{j = 1} ^{\lfloor \frac n {dx} \rfloor} j \\
=& \sum _{d = 1} ^n d ^3 \sum _{x = 1} ^{\lfloor \frac n d \rfloor} \mu (x) x ^2 S ^2(\lfloor \frac n {dx} \rfloor)
\end{aligned}
$$

注意:提 $x$ 出来的时候也得带上 $x ^2$。刚刚忘记带然后完全错了…

其中 $S(i) = \frac {i (i + 1)} 2$。

再枚举 $T = dx$:

$$
\begin{aligned}
& \sum _{T = 1} ^n S ^2 (\lfloor \frac n T \rfloor) T ^2 \sum _{d | T} d \mu (\frac T d) \\
=& \sum _{T = 1} ^n S ^2 (\lfloor \frac n T \rfloor) T ^2 \varphi (T)
\end{aligned}
$$

神奇的是后面就是 $\varphi$。那说明可能走了弯路。

事实上由于 $id = \varphi \ast I$,所以一开始的 $(i, j)$ 乘积项可以直接变为 $\sum _{x | (i, j)} \varphi(x)$。后续一样的套路,就不写了。

$I \ast \varphi = Id$,其中 $I(n) = 1, Id(n) = n$,根据这个式子杜教筛即可

然而又错了。现在要筛的函数是 $f(i) = i ^2 \varphi(i)$。

我是什么垃圾 竟然会觉得对 T ^2 phi(T) 分别求前缀和就能乘起来

然后发现这个玩意不会构造 $g$ 了,于是去看题解。

其实仔细想想,狄利克雷卷积是 $f(d) g(\frac T d)$ 的形式。

如果构造 $g(i) = i ^2$,那么刚好平方项会被约为 $T ^2$,并且 $\sum _{d | T} \varphi (d) = i$,于是得到 $f \ast g = h, h(i) = i ^3$。杜教筛即可。

设 $g(x) = x ^3, f(x) = \sum _{i | x} i ^3 \mu (\frac x i) = \mu \ast g$。因为需要杜教筛,并且看到了 $\mu$,那么显然左右两边同卷上 $I(x) = 1$,得到 $f \ast I = g$。$g, I$ 的前缀和都可以快速算出,杜教筛即可

上面这个是刚刚没带 $x ^2$ 的时候写的。带了 $x ^2$ 就没这么麻烦了。

另附:

$$
\sum _{i = 1} ^n i ^2 = \frac {n (n + 1) (2n + 1)} 6 \\
\sum _{i = 1} ^n i ^3 = \frac {n ^2 (n + 1) ^2} 4
$$

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
#include <iostream>
#include <vector>
#include <bitset>
#include <algorithm>

using namespace std;
using LL = long long;

const int maxN = 1e7;
const int sqr = 1e5 + 5;

int p, ans, inv2, inv4, inv6;
LL n;
int sphi[maxN + 2];
int mem[sqr];
vector<int> prim;
bitset<maxN + 2> vis;

inline void Mod(int& x)
{ x >= p ? x -= p : x < 0 ? x += p : 0; }

int QPow(int bas, int ind)
{
int res = 1;
while (ind)
{
if (ind & 1)
res = (LL)res * bas % p;
bas = (LL)bas * bas % p;
ind >>= 1;
}
return res;
}

inline int S(LL x)
{
x %= p;
return x * (x + 1) % p * inv2 % p;
}

inline int Q(LL x)
{
x %= p;
return x * (x + 1) % p * (x * 2 + 1) % p * inv6 % p;
}

inline int P(LL x)
{
x %= p;
return (x + 1) * (x + 1) % p * x % p * x % p * inv4 % p;
}

int F(LL x)
{
if (x <= maxN)
return sphi[x];
int& cur = mem[n / x];
if (~cur)
return cur;
cur = P(x);
for (LL l = 2, r; l <= x; l = r + 1)
{
r = x / (x / l);
Mod(cur -= LL(Q(r) - Q(l - 1)) * F(x / l) % p);
}
return cur;
}

int main()
{
#ifndef ONLINE_JUDGE
freopen("wib.in", "r", stdin);
freopen("wib.out", "w", stdout);
#endif
cin >> p >> n;
inv2 = (p + 1) >> 1, inv4 = QPow(4, p - 2), inv6 = QPow(6, p - 2);
// at first, sphi[i] = phi[i]
sphi[1] = 1;
for (int i = 2; i <= maxN; ++i)
{
if (!vis[i])
{
prim.push_back(i);
sphi[i] = i - 1;
}
for (int j : prim)
{
int x = i * j;
if (x > maxN)
break;
vis.set(x);
if (i % j)
sphi[x] = sphi[i] * (j - 1);
else
{
sphi[x] = sphi[i] * j;
break;
}
}
}
for (int i = 1; i <= maxN; ++i)
Mod(sphi[i] = sphi[i - 1] + (LL)i * i % p * sphi[i] % p);
fill(mem, mem + sqr, -1);
for (LL l = 1, r; l <= n; l = r + 1)
{
r = n / (n / l);
Mod(ans += (LL)S(n / l) * S(n / l) % p * (F(r) - F(l - 1)) % p);
}
cout << ans << endl;
return 0;
}

2019 南京网络赛 E K sum

Description

定义函数

$$
f _n (k) = \sum _{l _1 = 1} ^n \sum _{l _2 = 1} ^n \cdots \sum _{l _k = 1} ^n \gcd(l _1, l _2, \cdots, l _k) ^2
$$

现给定 $n, k$,需要求出 $\sum _{i = 2} ^k f _n (i)$,答案对 $10 ^9 + 7$ 取模。

$T$ 组数据。

$$
1 \le T \le 10, 1 \le n \le 10 ^9, 2 \le k \le 10 ^{10 ^5}, \text{3000ms}
$$

Solution

一样的先提出 $\gcd$:

$$
f _n (k) = \sum _{d = 1} ^n d ^2 \sum _{l _1 = 1} ^n \sum _{l _2 = 1} ^n \cdots \sum _{l _k = 1} ^n [\gcd(l _1, l _2, \cdots, l _k) = d]
$$

一开始发现自己做到这里就不会了… 后面那一块怎么处理啊…菜的真实

实际上这就是只会板子不理解原理的后果。$[n = 1] = \sum _{d | n} \mu (d)$,就算有 $k$ 个 $\gcd$ 一样可以变为 $\sum _{d | \gcd(l _1, l _2, \cdots, l _k)} \mu (d)$。可以和平常一样的变化出来:

$$
f _n (k) = \sum _{d = 1} ^n d ^2 \sum _{x = 1} ^{\lfloor \frac n d \rfloor} \mu (x) \lfloor \frac n {dx} \rfloor ^k
$$

再枚举 $T = dx$:

$$
f _n (k) = \sum _{T = 1} ^n \lfloor \frac n T \rfloor ^k \sum _{d | T} d ^2 \mu (\frac T d)
$$

正常来讲到这里就可以直接杜教筛了。但是这个坑题还要求前缀和… 还得另想办法。

由欧拉定理可以得到 $f _n(k) = f _n(k \bmod (10 ^9 + 6))$,尝试利用这个性质做一下

然后做不了。考虑在外层枚举 $i$ 从 2 到 $k$:

$$
\begin{aligned}
ans &= \sum _{i = 2} ^k \sum _{T = 1} ^n \lfloor \frac n T \rfloor ^i \sum _{d | T} d ^2 \mu (\frac T d) \\
&= \sum _{T = 1} ^n \left( \sum _{i = 2} ^k \lfloor \frac n T \rfloor ^i \right) \sum _{d | T} d ^2 \mu (\frac T d)
\end{aligned}
$$

想到了不会算等比数列求和菜的真实

令 $P(x) = \sum _{i = 2} ^k x ^i, f(x) = \sum _{d | x} d ^2 \mu (\frac T d)$。$f$ 可以与 $I(x) = 1$ 卷积并杜教筛,$P$ 可以等比数列求和得到 $(x - 1) P(x) = x ^{k + 1} - x ^2$,特判 $x$ 为 1 的情况并用欧拉定理降幂即可快速求和。因为 $\lfloor \frac n T \rfloor$ 取值种数不超过 $2 \sqrt n$,直接快速幂可以接受。

另外线性筛 $f$ 的时候可以证明的是 $f(px) = f(x) \cdot p ^2$,比较简单就略去证明过程了。利用这个线性筛,或者用 low 存下最小质因子的幂次也行。代码里面用的是 low

不会分析杜教筛的题的复杂度。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
#include <cctype>
#include <iostream>
#include <vector>
#include <bitset>
#include <unordered_map>

using namespace std;
using LL = long long;

const int maxN = 1e6;
const int mod = 1e9 + 7;

int n, T, k, tmpk, ans, inv6;
int sf[maxN + 2], low[maxN + 2];
bitset<maxN + 2> vis;
vector<int> prim;
unordered_map<int, int> f;

inline void Mod(int& x)
{ x >= mod ? x -= mod : x < 0 ? x += mod : 0; }

int QPow(int bas, int ind)
{
int res = 1;
while (ind)
{
if (ind & 1)
res = (LL)res * bas % mod;
bas = (LL)bas * bas % mod;
ind >>= 1;
}
return res;
}

inline int P(int x)
{
if (x == 1)
return tmpk;
int res = LL(QPow(x, k + 1) - QPow(x, 2)) * QPow(x - 1, mod - 2) % mod;
Mod(res);
return res;
}

int Sf(int x)
{
if (x <= maxN)
return sf[x];
if (f.count(x))
return f[x];
int& cur = f[x];
cur = 1ll * x * (x + 1) % mod * (2 * x + 1) % mod * inv6 % mod;
for (int l = 2, r; l <= x; l = r + 1)
{
r = x / (x / l);
Mod(cur -= LL(r - l + 1) * Sf(x / l) % mod);
}
return cur;
}

int main()
{
#ifndef ONLINE_JUDGE
freopen("wib.in", "r", stdin);
freopen("wib.out", "w", stdout);
#endif
ios::sync_with_stdio(false);
cin >> T;
inv6 = QPow(6, mod - 2);
sf[1] = 1;
for (int i = 2; i <= maxN; ++i)
{
if (!vis[i])
{
low[i] = i;
Mod(sf[i] = (LL)i * i % mod - 1);
prim.push_back(i);
}
for (int j : prim)
{
int x = i * j;
if (x > maxN)
break;
vis.set(x);
if (i % j)
low[x] = j, sf[x] = (LL)sf[i] * sf[j] % mod;
else
{
low[x] = low[i] * j;
if (low[i] == i)
sf[x] = (LL)sf[i] * j % mod * j % mod;
else
sf[x] = (LL)sf[i / low[i]] * sf[j * low[i]] % mod;
break;
}
}
}
//for (int i = 1; i <= 12; ++i)
// cout << sf[i] << endl;
for (int i = 1; i <= maxN; ++i)
Mod(sf[i] += sf[i - 1]);
while (T--)
{
cin >> n;
k = tmpk = ans = 0;
char ch;
while (!isdigit(ch = cin.get()));
while (k = (k * 10ll + (ch & 15)) % (mod - 1),
tmpk = (tmpk * 10ll + (ch & 15)) % mod,
isdigit(ch = cin.get()));
Mod(--tmpk);
for (int l = 1, r; l <= n; l = r + 1)
{
r = n / (n / l);
Mod(ans += (LL)P(n / l) * (Sf(r) - Sf(l - 1)) % mod);
}
cout << ans << endl;
}
return 0;
}

[PE530]GCD of Divisors

Description

$$
f(n) = \sum \limits _{d | n} (d, \frac n d)
$$

(省略了 $\gcd$)

$$
F(k) = \sum _{n = 1} ^k f(n), k = 10 ^{15}
$$

Solution

首先这个 $10 ^{15}$ 的范围应该不是任何一个筛能够承受的…得考虑直接从 $F$ 开始展开才行。于是先得到最初的式子:

$$
\sum _{n = 1} ^k \sum _{d | n} (d, \frac n d)
$$

将有贡献的 $(d, \frac n d)$ 提出来,设为 $g$:

$$
\sum _{g = 1} ^k g \sum _{n = 1} ^k \sum _{d | n} [(d, \frac n d) = g]
$$

$(d, \frac n d) = g$ 说明 $n$ 应当是 $g ^2$ 的倍数,于是改变 $n$ 的枚举意义与范围,变为枚举 $g ^2$ 的倍数:

$$
\sum _{g = 1} ^k g \sum _{n = 1} ^{\lfloor \frac k {g ^2} \rfloor} \sum _{d | n} [(d, \frac n d) = 1]
$$

对后面应用莫比乌斯反演:

$$
\sum _{g = 1} ^k g \sum _{n = 1} ^{\lfloor \frac k {g ^2} \rfloor} \sum _{d | n} \sum _{x | (d, \frac n d)} \mu (x)
$$

将 $x$ 提到前面:

$$
\sum _{g = 1} ^k \sum _{x = 1} ^{\lfloor \frac k {g ^2} \rfloor} \mu(x) g \sum _{n = 1} ^{\lfloor \frac k {g ^2} \rfloor} \sum _{d | n} [x | d, x | \frac n d]
$$

Update:$x$ 的枚举上界 $k$ 错了,应当是 $\lfloor \frac k {g ^2} \rfloor$,已修正。之前的变成伪证了…

同理,此时 $n$ 应当是 $x ^2$ 的倍数,再次更改枚举意义:

$$
\sum _{g = 1} ^k \sum _{x = 1} ^{\lfloor \frac k {g ^2} \rfloor} \mu(x) g \sum _{n = 1} ^{\lfloor \frac k {(gx) ^2} \rfloor} \sum _{d | n} 1
$$

这貌似不太好化,应当考虑如何扩界。主要是看到 $x$ 的上界 $\lfloor \frac k {g ^2} \rfloor$ 这里,如果能够扩成 $\lfloor \frac k g \rfloor$,那么便可以构造卷积使其变为 $\varphi$。

如果对为什么扩成 $\lfloor \frac k g \rfloor$ 感到疑惑,请继续往下看。

考虑这个 $x$ 的界:因为后面的 $n$ 枚举到的是 $\lfloor \frac k {(gx) ^2} \rfloor$,相当于 $\lfloor \frac {\lfloor \frac k {g ^2} \rfloor} {x ^2} \rfloor$,所以 $x$ 的上界应当是 $\lfloor \frac {\sqrt k} g \rfloor$。也就是说,$x$ 可以任意扩到 $\lfloor \frac {\sqrt k} g \rfloor$ 以上的界。

将界修改为 $\lfloor \frac k g \rfloor$,可知现在 $x$ 的意义应当是枚举 $g$ 的 $x$ 倍。于是令 $T = gx$,仅枚举 $T, g$,则 $x = \frac T g$,得到

$$
\sum _{T = 1} ^{\sqrt k} \sum _{g | T} g \mu(\frac T g) \sum _{n = 1} ^{\lfloor \frac k {T ^2} \rfloor} \sum _{d | n} 1
$$

(相当于是:原本枚举所有数的倍数,现在枚举所有数的因数,易知这两个枚举顺序等价)

因为 $T$ 在 $\sqrt k$ 以上时,后面那一块值为 0,所以可以将 $T$ 的范围缩减到 $\sqrt k$。

将 $\mu \ast Id$(这里的 $\ast$ 指狄利克雷卷积)变为 $\varphi$,并将后面的 $d$ 提到 $n$ 前面:

(虽然 $\sum \limits _{d | n} 1$ 是 $\sigma _0 (n)$,但是这个范围怕是筛不动…)

$$
\sum _{T = 1} ^{\sqrt k} \varphi (T) \sum _{d = 1} ^{\lfloor \frac k {T ^2} \rfloor} \lfloor \frac {\lfloor \frac k {T ^2} \rfloor} d \rfloor
$$

$T$ 线性筛,后面的对 $\lfloor \frac k {T ^2} \rfloor$ 整除分块即可解决问题。

复杂度:

$$
\begin{aligned}
& \sum _{T = 1} ^{\sqrt k} \sqrt {\frac k {T ^2}} \\
=& \sqrt k \sum _{T = 1} ^{\sqrt k} O(\frac 1 T) \\
=& \sqrt k \ln \sqrt k
\end{aligned}
$$

开 O2 大概要跑 30s。$F(10 ^{15}) = 207366437157977206$。

Update:实际上令 $h(n) = \sum _{d = 1} ^n \lfloor \frac n d \rfloor$,原式变为 $\sum _{T = 1} ^{\sqrt k} \varphi (T) h(\lfloor \frac k {T ^2} \rfloor)$,可以对 $\lfloor \frac k {T ^2} \rfloor$ 整除分块(并不知道对于平方的整除分块复杂度怎么证),跑起来大概要 22s。

具体分块方法:对于当前端点 $l$,先得到对于 $l ^2$ 而言的右端点 $r’$,那么取值一样的区间应当是 $[l, \lfloor \sqrt {r’} \rfloor]$。orz jambow

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#include <cmath>
#include <iostream>
#include <vector>
#include <bitset>

using namespace std;
using LL = long long;

const int N = 3.2e7;
const int maxN = N + 2;

LL k, ans;
LL phi[maxN];
bitset<maxN> vis;
vector<int> prim;

void Euler_Sieve()
{
phi[1] = 1;
for (int i = 2; i <= N; ++i)
{
if (!vis[i])
phi[i] = i - 1, prim.push_back(i);
for (int j : prim)
{
int x = i * j;
if (x > N)
break;
vis.set(x);
if (i % j)
phi[x] = phi[i] * (j - 1);
else
{
phi[x] = phi[i] * j;
break;
}
}
}
for (int i = 2; i <= N; ++i)
phi[i] += phi[i - 1];
}

inline LL F(LL n)
{
LL res = 0;
for (LL l = 1, r; l <= n; l = r + 1)
{
r = n / (n / l);
res += (r - l + 1) * (n / l);
}
return res;
}

int main()
{
freopen("wib.in", "r", stdin);
freopen("wib.out", "w", stdout);
ios::sync_with_stdio(false);
cin >> k;
Euler_Sieve();
int t = sqrt((long double)k);
for (int l = 1, r; l <= t; l = r + 1)
{
r = sqrt((long double)(k / (k / ((LL)l * l))));
LL p = k / ((LL)l * l);
ans += F(p) * (phi[r] - phi[l - 1]);
}
cout << ans << endl;
return 0;
}

To be continued