本文的主要目的是记录莫比乌斯反演的常见的几个推导优化，数学不好就要多背诵。主要来源于OY的github，同时也使用了OY的算法模板。如果你恰好了解狄利克雷、莫比乌斯等相关概念，又不会自己推导，可以看一下。
示性函数
首先熟悉一下反演推导中常用的示性函数的写法：
[
c
o
n
d
i
t
i
o
n
]
=
{
1
c
o
n
d
i
t
i
o
n
=
t
r
u
e
0
c
o
n
d
i
t
i
o
n
=
f
a
l
s
e
[condition] =\left\{ \begin{array}{ll} 1 & condition = \boldsymbol{true} \\ 0 & condition = \boldsymbol{false} \end{array}\right.
[condition]={10condition=truecondition=false
因此，诸如 ∑ i = 1 N [ i ∣ N ] \sum_{i=1}^{N}[i|N] ∑i=1N[i∣N]就是求 N N N的所有约数的数量，而 ∑ i = 1 N i ⋅ [ i ∣ N ] \sum_{i=1}^{N}i\cdot[i|N] ∑i=1Ni⋅[i∣N]就是求约数和。
求和符号交换位置
然后要了解常见的求和符号交换位置。为了方便起见，本文假定求和均从1开始（除非交换需要），即下限固定，只有上限。则最普适的交换：当上限与求和项无关，可以交换。
∑ i = 1 N ∑ j = 1 M A i , j = ∑ j = 1 M ∑ i = 1 N A i , j \sum_{i=1}^{N}\sum_{j=1}^{M}A_{i,j}=\sum_{j=1}^{M}\sum_{i=1}^{N}A_{i,j} i=1∑Nj=1∑MAi,j=j=1∑Mi=1∑NAi,j
第二种比较常见的一般情况是：
∑ i = 1 N ∑ j = 1 i A i , j = ∑ j = 1 N ∑ i = j N A i , j \sum_{i=1}^{N}\sum_{j=1}^{i}A_{i,j}=\sum_{j=1}^{N}\sum_{i=j}^{N}A_{i,j} i=1∑Nj=1∑iAi,j=j=1∑Ni=j∑NAi,j
此公式可以借助矩阵形式进行演示，也可以使用示性函数结合第一种交换进行推导，具体如下：
∑ i = 1 N ∑ j = 1 i A i , j = ∑ i = 1 N ∑ j = 1 N A i , j ⋅ [ j ≤ i ] \sum_{i=1}^{N}\sum_{j=1}^{i}A_{i,j}=\sum_{i=1}^{N}\sum_{j=1}^{N}A_{i,j}\cdot[j\le{i}] i=1∑Nj=1∑iAi,j=i=1∑Nj=1∑NAi,j⋅[j≤i]
注意到右边的求和符号，上限已经与求和项无关，因此可以交换。有：
∑
i
=
1
N
∑
j
=
1
N
A
i
,
j
⋅
[
j
≤
i
]
=
∑
j
=
1
N
∑
i
=
1
N
A
i
,
j
⋅
[
j
≤
i
]
\sum_{i=1}^{N}\sum_{j=1}^{N}A_{i,j}\cdot[j\le{i}]=\sum_{j=1}^{N}\sum_{i=1}^{N}A_{i,j}\cdot[j\le{i}]
i=1∑Nj=1∑NAi,j⋅[j≤i]=j=1∑Ni=1∑NAi,j⋅[j≤i]
再把示性函数去掉，就可以得到上述等式。
常见优化推导
互质公式与互质枚举
[ g c d ( i , j ) = = 1 ] = ∑ d ∣ g c d ( i , j ) μ ( d ) \begin{equation} [gcd(i,j)==1]=\sum_{d|gcd(i,j)}\mu(d) \end{equation} [gcd(i,j)==1]=d∣gcd(i,j)∑μ(d)
这其实就是莫比乌斯函数的定义，因为 [ n = = 1 ] [n==1] [n==1]就是狄利克雷卷积的单位元，而莫比乌斯函数与全1函数互为逆元。
对于入门级别的反演，最重要的就是推导出包含gcd的式子，因为可以使用如下公式进行进一步推导。
∑ i = 1 N ∑ j = 1 M f ( g c d ( i , j ) ) g ( i , j ) = ∑ d = 1 N ∑ i = 1 N ∑ j = 1 M f ( d ) g ( i , j ) [ d = = g c d ( i , j ) ] = ∑ d = 1 N f ( d ) ( ∑ i = 1 N / d ∑ j = 1 M / j g ( i d , j d ) [ 1 = = g c d ( i , j ) ] ) \sum_{i=1}^{N}\sum_{j=1}^{M}f(gcd(i,j))g(i,j)=\sum_{d=1}^{N}\sum_{i=1}^{N}\sum_{j=1}^{M}f(d)g(i,j)[d==gcd(i,j)] \\=\sum_{d=1}^{N}f(d)\big(\sum_{i=1}^{N/d}\sum_{j=1}^{M/j}g(id,jd)[1==gcd(i,j)]\big) i=1∑Nj=1∑Mf(gcd(i,j))g(i,j)=d=1∑Ni=1∑Nj=1∑Mf(d)g(i,j)[d==gcd(i,j)]=d=1∑Nf(d)(i=1∑N/dj=1∑M/jg(id,jd)[1==gcd(i,j)])
再利用互质化简，可以进一步进行推导，直接记结论：
∑ i = 1 N ∑ j = 1 M f ( g c d ( i , j ) ) g ( i , j ) = ∑ d = 1 N ∑ i = 1 N / d ∑ j = 1 M / d g ( i d , j d ) ⋅ ∑ t ∣ d f ( t ) μ ( d t ) \begin{equation} \sum_{i=1}^{N}\sum_{j=1}^{M}f(gcd(i,j))g(i,j)=\sum_{d=1}^{N}\sum_{i=1}^{N/d}\sum_{j=1}^{M/d}g(id,jd)\cdot\sum_{t|d}f(t)\mu(\frac{d}{t}) \end{equation} i=1∑Nj=1∑Mf(gcd(i,j))g(i,j)=d=1∑Ni=1∑N/dj=1∑M/dg(id,jd)⋅t∣d∑f(t)μ(td)
也就是说一个二元函数求和最后可以转换成只在两个参数互质时才有贡献。这是最重要的一个公式。
先看一个没有用的推导，即又反过来求
∑
i
∑
j
[
1
=
=
g
c
d
(
i
,
j
)
]
\sum_{i}\sum_{j}[1==gcd(i,j)]
i∑j∑[1==gcd(i,j)]
令 f ( g c d ( i , j ) ) = [ 1 = = g c d ( i , j ) ] f(gcd(i,j))=[1==gcd(i,j)] f(gcd(i,j))=[1==gcd(i,j)]， g ( i , j ) ≡ 1 g(i,j)\equiv1 g(i,j)≡1，根据公式有：
∑ i ∑ j [ 1 = = g c d ( i , j ) ] = ∑ d = 1 N ∑ i = 1 N / d ∑ j = 1 N / d g ( i d , j d ) ⋅ ∑ t ∣ d f ( t ) μ ( d t ) \sum_{i}\sum_{j}[1==gcd(i,j)]=\sum_{d=1}^{N}\sum_{i=1}^{N/d}\sum_{j=1}^{N/d}g(id,jd)\cdot\sum_{t|d}f(t)\mu(\frac{d}{t}) i∑j∑[1==gcd(i,j)]=d=1∑Ni=1∑N/dj=1∑N/dg(id,jd)⋅t∣d∑f(t)μ(td)
显然只在
t
=
1
t=1
t=1时有贡献，又因为
g
g
g函数恒等于1，因此原式变为
∑
d
=
1
N
∑
i
=
1
N
/
d
∑
j
=
1
N
/
d
μ
(
d
)
\sum_{d=1}^{N}\sum_{i=1}^{N/d}\sum_{j=1}^{N/d}\mu(d)
d=1∑Ni=1∑N/dj=1∑N/dμ(d)
注意到
μ
(
d
)
\mu(d)
μ(d)与后两个求和变量无关，因此可以提取出来，有
∑
d
=
1
N
μ
(
d
)
∑
i
=
1
N
/
d
∑
j
=
1
N
/
d
1
\sum_{d=1}^{N}\mu(d)\sum_{i=1}^{N/d}\sum_{j=1}^{N/d}\mathit{1}
d=1∑Nμ(d)i=1∑N/dj=1∑N/d1
为了简便起见，这里一律省略向下取整符号。
约数公式与因数枚举
令
d
(
x
)
d(x)
d(x)表示
x
x
x的约数的数量，则
d
(
x
⋅
y
)
=
∑
i
∣
x
∑
i
∣
y
[
g
c
d
(
i
⋅
j
=
=
1
]
\begin{equation} d(x\cdot{y})=\sum_{i|x}\sum_{i|y}[gcd(i\cdot{j}==1] \end{equation}
d(x⋅y)=i∣x∑i∣y∑[gcd(i⋅j==1]
因此如果需要枚举因数，可以写成
∑
d
=
1
x
⋅
y
[
d
∣
x
⋅
y
]
=
∑
i
∣
x
∑
j
∣
y
[
g
c
d
(
i
,
j
)
=
=
1
]
\begin{equation} \sum_{d=1}^{x\cdot{y}}[d|x\cdot{y}]=\sum_{i|x}\sum_{j|y}[gcd(i,j)==1] \end{equation}
d=1∑x⋅y[d∣x⋅y]=i∣x∑j∣y∑[gcd(i,j)==1]
数组元素转化值域
对于数组 A A A，令 1 ≤ A i ≤ L 1\le{A_i}\le{L} 1≤Ai≤L，则
∑ i = 1 N ∑ j = 1 N f ( A i , A j ) = ∑ i = 1 L ∑ j = 1 L f ( i , j ) ⋅ c i ⋅ c j \begin{equation} \sum_{i=1}^{N}\sum_{j=1}^{N}f(A_i,A_j)=\sum_{i=1}^{L}\sum_{j=1}^{L}f(i,j)\cdot{c_i}\cdot{c_j} \end{equation} i=1∑Nj=1∑Nf(Ai,Aj)=i=1∑Lj=1∑Lf(i,j)⋅ci⋅cj
其中， c i c_i ci表示数值 i i i在数组中出现的次数。当满足 f ( i , j ) = f ( j , i ) f(i,j)=f(j,i) f(i,j)=f(j,i)时，一个常见的小变形：
∑ i = 1 N − 1 ∑ j = i + 1 N f ( A i , A j ) = 1 2 ( ∑ i = 1 N ∑ j = 1 N f ( A i , A j ) − ∑ i = 1 N f ( A i , A i ) ) \begin{equation} \sum_{i=1}^{N-1}\sum_{j=i+1}^{N}f(A_i,A_j)=\frac{1}{2}\big(\sum_{i=1}^{N}\sum_{j=1}^{N}f(A_i,A_j)-\sum_{i=1}^{N}f(A_i,A_i)\big) \end{equation} i=1∑N−1j=i+1∑Nf(Ai,Aj)=21(i=1∑Nj=1∑Nf(Ai,Aj)−i=1∑Nf(Ai,Ai))
洛谷P3327
求：
∑
i
=
1
N
∑
j
=
1
M
d
(
i
j
)
\sum_{i=1}^{N}\sum_{j=1}^{M}d(ij)
i=1∑Nj=1∑Md(ij)
首先利用约数公式，原式变为：
∑ i = 1 N ∑ j = 1 M ∑ u ∣ i ∑ v ∣ j [ g c d ( u ⋅ v ) = = 1 ] \sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{u|i}\sum_{v|j}[gcd(u\cdot{v})==1] i=1∑Nj=1∑Mu∣i∑v∣j∑[gcd(u⋅v)==1]
对后两个求和符号，可以利用第二种交换模式，即利用示性函数：
∑ i = 1 N ∑ j = 1 M ∑ u = 1 N ∑ v = 1 M [ g c d ( u ⋅ v ) = = 1 ] ⋅ [ u ∣ i ] ⋅ [ v ∣ j ] \sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{u=1}^{N}\sum_{v=1}^{M}[gcd(u\cdot{v})==1]\cdot{[u|i]}\cdot{[v|j]} i=1∑Nj=1∑Mu=1∑Nv=1∑M[gcd(u⋅v)==1]⋅[u∣i]⋅[v∣j]
注意到求和项与上限无关，因此可以把后两个求和符号与前两个求和符号交换位置，有
∑ u = 1 N ∑ v = 1 M ∑ i = 1 N ∑ j = 1 M [ g c d ( u ⋅ v ) = = 1 ] ⋅ [ u ∣ i ] ⋅ [ v ∣ j ] \sum_{u=1}^{N}\sum_{v=1}^{M}\sum_{i=1}^{N}\sum_{j=1}^{M}[gcd(u\cdot{v})==1]\cdot{[u|i]}\cdot{[v|j]} u=1∑Nv=1∑Mi=1∑Nj=1∑M[gcd(u⋅v)==1]⋅[u∣i]⋅[v∣j]
又因为 g c d ( u , v ) = = 1 gcd(u,v)==1 gcd(u,v)==1与 i , j i,j i,j均无关，因此可以提取出来，得到
∑ u = 1 N ∑ v = 1 M [ g c d ( u ⋅ v ) = = 1 ] ⋅ ∑ i = 1 N ∑ j = 1 M [ u ∣ i ] ⋅ [ v ∣ j ] \sum_{u=1}^{N}\sum_{v=1}^{M}[gcd(u\cdot{v})==1]\cdot\sum_{i=1}^{N}\sum_{j=1}^{M}{[u|i]}\cdot{[v|j]} u=1∑Nv=1∑M[gcd(u⋅v)==1]⋅i=1∑Nj=1∑M[u∣i]⋅[v∣j]
最后的求和其实就是问
[
1
,
N
]
[1,N]
[1,N]区间内有多少个
u
u
u的倍数，自然就是
N
u
\frac{N}{u}
uN，因此原式最终变为（顺便改变一下求和下标）
∑
i
=
1
N
∑
j
=
1
M
[
g
c
d
(
i
⋅
j
)
=
=
1
]
⋅
N
i
⋅
M
j
\sum_{i=1}^{N}\sum_{j=1}^{M}[gcd(i\cdot{j})==1]\cdot\frac{N}{i}\cdot\frac{M}{j}
i=1∑Nj=1∑M[gcd(i⋅j)==1]⋅iN⋅jM
利用互质枚举公式，令
f
(
x
)
=
[
x
=
=
1
]
f(x)=[x==1]
f(x)=[x==1]，
g
(
x
,
y
)
=
N
x
⋅
M
y
g(x,y)=\frac{N}{x}\cdot\frac{M}{y}
g(x,y)=xN⋅yM，原式化为：
∑
d
=
1
N
∑
i
=
1
N
/
d
∑
j
=
1
M
/
d
g
(
i
d
,
j
d
)
⋅
∑
t
∣
d
f
(
t
)
μ
(
d
t
)
\sum_{d=1}^{N}\sum_{i=1}^{N/d}\sum_{j=1}^{M/d}g(id,jd)\cdot\sum_{t|d}f(t)\mu(\frac{d}{t})
d=1∑Ni=1∑N/dj=1∑M/dg(id,jd)⋅t∣d∑f(t)μ(td)
注意到最后的和式只在 t = 1 t=1 t=1时有贡献，所以有所求为
∑ d = 1 N ∑ i = 1 N / d ∑ j = 1 M / d g ( i d , j d ) ⋅ μ ( d ) \sum_{d=1}^{N}\sum_{i=1}^{N/d}\sum_{j=1}^{M/d}g(id,jd)\cdot\mu(d) d=1∑Ni=1∑N/dj=1∑M/dg(id,jd)⋅μ(d)
又因为 μ ( d ) \mu(d) μ(d)与 i , j i,j i,j无关，可以提前，再代入 g g g函数有：
∑ d = 1 N μ ( d ) ∑ i = 1 N / d ∑ j = 1 M / d N i ⋅ d ⋅ M j ⋅ d \sum_{d=1}^{N}\mu(d)\sum_{i=1}^{N/d}\sum_{j=1}^{M/d}\frac{N}{i\cdot{d}}\cdot\frac{M}{j\cdot{d}} d=1∑Nμ(d)i=1∑N/dj=1∑M/di⋅dN⋅j⋅dM
最后一步比较简单，注意到
i
,
j
i,j
i,j互不影响，因此又可以将
i
i
i项提前，有
∑
d
=
1
N
μ
(
d
)
∑
i
=
1
N
/
d
N
i
⋅
d
∑
j
=
1
M
/
d
M
j
⋅
d
\sum_{d=1}^{N}\mu(d)\sum_{i=1}^{N/d}\frac{N}{i\cdot{d}}\sum_{j=1}^{M/d}\frac{M}{j\cdot{d}}
d=1∑Nμ(d)i=1∑N/di⋅dNj=1∑M/dj⋅dM
这就是最终结论，对于多样例情况可以使用整除分块做预处理，在
N
N
N\sqrt{N}
NN内解决。
以下是代码。这个代码用到了上面提到的old-yan的算法模板，因此比较长。实际上如果自行实现的话，只需要先筛一下求出莫比乌斯函数及其前缀和，然后对所有的 1 ≤ j ≤ 50000 1\le{j}\le{50000} 1≤j≤50000预处理出所有 ∑ i = 1 j j i \sum_{i=1}^{j}\frac{j}{i} ∑i=1jij。然后对每个case，用 N \sqrt{N} N的时间就能解决。
#include <bits/stdc++.h>
using namespace std;
#ifdef _MSC_VER
#define CPP_STANDARD _MSVC_LANG
#else
#define CPP_STANDARD __cplusplus
#endif
#if CPP_STANDARD >= 201402L
#define CONSTEXPR_IF_CPP14 constexpr
#else
#define CONSTEXPR_IF_CPP14
#endif
namespace std {
namespace my_bit_ops {
static constexpr int s32[32] = {31, 0, 1, 5, 2, 10, 6, 15, 3, 13, 11, 20, 7, 22, 16, 25, 30, 4, 9, 14, 12, 19, 21, 24, 29, 8, 18, 23, 28, 17, 27, 26};
static constexpr int s64[64] = {63, 0, 1, 6, 2, 12, 7, 18, 3, 24, 13, 27, 8, 33, 19, 39, 4, 16, 25, 37, 14, 45, 28, 47, 9, 30, 34, 53, 20, 49, 40, 56, 62, 5, 11, 17, 23, 26, 32, 38, 15, 36, 44, 46, 29, 52, 48, 55, 61, 10, 22, 31, 35, 43, 51, 54, 60, 21, 42, 50, 59, 41, 58, 57};
static constexpr int s256[256] = {8, 7, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
static constexpr int ctz32(uint32_t x) noexcept { return s32[(x & -x) * 0b00001000110010100111010110111110 >> 27]; }
static constexpr int ctz64(uint64_t x) noexcept { return s64[(x & -x) * 0b0000010000110001010001110010010110011010011110101011101101111110 >> 58]; }
static constexpr int clz8(uint32_t x) noexcept { return s256[x]; }
static constexpr int clz16(uint32_t x) noexcept { return x >> 8 ? clz8(x >> 8) : 8 + clz8(x); }
static constexpr int clz32(uint32_t x) noexcept { return x >> 16 ? clz16(x >> 16) : 16 + clz16(x); }
static constexpr int clz64(uint64_t x) noexcept { return x >> 32 ? clz32(x >> 32) : 32 + clz32(x); }
static constexpr uint32_t f1_32(uint32_t x) noexcept { return x - ((x >> 1) & 0x55555555); }
static constexpr uint32_t f2_32(uint32_t x) noexcept { return (x & 0x33333333) + ((x >> 2) & 0x33333333); }
static constexpr int f3_32(uint32_t x) noexcept { return (x + (x >> 4) & 0x0f0f0f0f) * 0x01010101 >> 24; }
static constexpr int popcnt32(uint32_t x) noexcept { return f3_32(f2_32(f1_32(x))); }
static constexpr uint64_t f1_64(uint64_t x) noexcept { return x - ((x >> 1) & 0x5555555555555555); }
static constexpr uint64_t f2_64(uint64_t x) noexcept { return (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333); }
static constexpr int f3_64(uint64_t x) noexcept { return (x + (x >> 4) & 0x0f0f0f0f0f0f0f0f) * 0x0101010101010101 >> 56; }
static constexpr int popcnt64(uint64_t x) noexcept { return f3_64(f2_64(f1_64(x))); }
constexpr uint32_t bit_ceil_msvc_32(uint32_t x) noexcept { return x ? uint32_t(1) << (32 - clz32(x - 1)) : 1; }
constexpr uint64_t bit_ceil_msvc_64(uint64_t x) noexcept { return x ? uint64_t(1) << (64 - clz64(x - 1)) : 1; }
template <typename Tp>
constexpr Tp bit_ceil_msvc(Tp x) noexcept { return sizeof(Tp) == 4 ? bit_ceil_msvc_32(x) : bit_ceil_msvc_64(x); }
constexpr uint32_t bit_floor_msvc_32(uint32_t x) noexcept { return x ? uint32_t(1) << (31 - clz32(x)) : 0; }
constexpr uint64_t bit_floor_msvc_64(uint64_t x) noexcept { return x ? uint64_t(1) << (63 - clz64(x)) : 0; }
template <typename Tp>
constexpr Tp bit_floor_msvc(Tp x) noexcept { return sizeof(Tp) == 4 ? bit_floor_msvc_32(x) : bit_floor_msvc_64(x); }
constexpr uint32_t bit_width_msvc_32(uint32_t x) noexcept { return 32 - clz32(x); }
constexpr uint64_t bit_width_msvc_64(uint64_t x) noexcept { return 64 - clz64(x); }
template <typename Tp>
constexpr Tp bit_width_msvc(Tp x) noexcept { return sizeof(Tp) == 4 ? bit_width_msvc_32(x) : bit_width_msvc_64(x); }
constexpr int countl_zero_msvc_32(uint32_t x) noexcept { return x ? clz32(x) : 32; }
constexpr int countl_zero_msvc_64(uint64_t x) noexcept { return x ? clz64(x) : 64; }
template <typename Tp>
constexpr int countl_zero_msvc(Tp x) noexcept { return sizeof(Tp) == 4 ? countl_zero_msvc_32(x) : countl_zero_msvc_64(x); }
constexpr int countr_zero_msvc_32(uint32_t x) noexcept { return x ? ctz32(x) : 32; }
constexpr int countr_zero_msvc_64(uint64_t x) noexcept { return x ? ctz64(x) : 64; }
template <typename Tp>
constexpr int countr_zero_msvc(Tp x) noexcept { return sizeof(Tp) == 4 ? countr_zero_msvc_32(x) : countr_zero_msvc_64(x); }
template <typename Tp>
constexpr int popcount_msvc(Tp x) noexcept { return sizeof(Tp) == 4 ? popcnt32(x) : popcnt64(x); }
#ifdef _MSC_VER
template <typename Tp>
constexpr Tp bit_ceil(Tp x) noexcept { return bit_ceil_msvc(x); }
template <typename Tp>
constexpr Tp bit_floor(Tp x) noexcept { return bit_floor_msvc(x); }
template <typename Tp>
constexpr Tp bit_width(Tp x) noexcept { return bit_width_msvc(x); }
template <typename Tp>
constexpr int countl_zero(Tp x) noexcept { return countl_zero_msvc(x); }
template <typename Tp>
constexpr int countr_zero(Tp x) noexcept { return countr_zero_msvc(x); }
template <typename Tp>
constexpr int popcount(Tp x) noexcept { return popcount_msvc(x); }
#else
constexpr uint32_t bit_ceil_32(uint32_t x) noexcept { return x > 1 ? uint32_t(1) << (32 - __builtin_clz(x - 1)) : 1; }
constexpr uint64_t bit_ceil_64(uint64_t x) noexcept { return x > 1 ? uint64_t(1) << (64 - __builtin_clzll(x - 1)) : 1; }
template <typename Tp>
constexpr Tp bit_ceil(Tp x) noexcept { return sizeof(Tp) == 4 ? bit_ceil_32(x) : bit_ceil_64(x); }
constexpr uint32_t bit_floor_32(uint32_t x) noexcept { return x ? uint32_t(1) << (31 - __builtin_clz(x)) : 0; }
constexpr uint64_t bit_floor_64(uint64_t x) noexcept { return x ? uint64_t(1) << (63 - __builtin_clzll(x)) : 0; }
template <typename Tp>
constexpr Tp bit_floor(Tp x) noexcept { return sizeof(Tp) == 4 ? bit_floor_32(x) : bit_floor_64(x); }
constexpr uint32_t bit_width_32(uint32_t x) noexcept { return x ? 32 - __builtin_clz(x) : 0; }
constexpr uint64_t bit_width_64(uint64_t x) noexcept { return x ? 64 - __builtin_clzll(x) : 0; }
template <typename Tp>
constexpr Tp bit_width(Tp x) noexcept { return sizeof(Tp) == 4 ? bit_width_32(x) : bit_width_64(x); }
constexpr int countl_zero_32(uint32_t x) noexcept { return x ? __builtin_clz(x) : 32; }
constexpr int countl_zero_64(uint64_t x) noexcept { return x ? __builtin_clzll(x) : 64; }
template <typename Tp>
constexpr int countl_zero(Tp x) noexcept { return sizeof(Tp) == 4 ? countl_zero_32(x) : countl_zero_64(x); }
constexpr int countr_zero_32(uint32_t x) noexcept { return x ? __builtin_ctz(x) : 32; }
constexpr int countr_zero_64(uint64_t x) noexcept { return x ? __builtin_ctzll(x) : 64; }
template <typename Tp>
constexpr int countr_zero(Tp x) noexcept { return sizeof(Tp) == 4 ? countr_zero_32(x) : countr_zero_64(x); }
template <typename Tp>
constexpr int popcount(Tp x) noexcept { return sizeof(Tp) == 4 ? __builtin_popcount(x) : __builtin_popcountll(x); }
#endif
}
#ifndef __cpp_lib_bitops
using my_bit_ops::bit_ceil;
using my_bit_ops::bit_floor;
using my_bit_ops::bit_width;
using my_bit_ops::countl_zero;
using my_bit_ops::countr_zero;
using my_bit_ops::popcount;
#endif
}
namespace OY {
namespace MOBIUS {
using size_type = uint32_t;
struct CompoundPlus {
template <typename Tp>
void operator()(Tp &x, const Tp &y) const { x += y; }
};
struct CompoundMinus {
template <typename Tp>
void operator()(Tp &x, const Tp &y) const { x -= y; }
};
template <size_type MAX_RANGE, bool Flag>
struct OptionalArray {
size_type m_data[MAX_RANGE + 1];
template <typename Mapping>
void prepare_presum(size_type range, Mapping mapping) {
m_data[0] = 0;
for (size_type i = 1; i <= range; i++) m_data[i] = m_data[i - 1] + mapping(i);
}
void set(size_type i, size_type val) { m_data[i] = val; }
size_type query(size_type left, size_type right) const { return m_data[right] - m_data[left - 1]; }
size_type operator[](size_type i) const { return m_data[i]; }
};
template <size_type MAX_RANGE>
struct OptionalArray<MAX_RANGE, false> {};
constexpr size_type get_estimated_ln(size_type x) {
return x <= 7 ? 1
: x <= 32 ? 2
: x <= 119 ? 3
: x <= 359 ? 4
: x <= 1133 ? 5
: x <= 3093 ? 6
: x <= 8471 ? 7
: x <= 24299 ? 8
: x <= 64719 ? 9
: x <= 175196 ? 10
: x <= 481451 ? 11
: x <= 1304718 ? 12
: x <= 3524653 ? 13
: x <= 9560099 ? 14
: x <= 25874783 ? 15
: x <= 70119984 ? 16
: x <= 189969353 ? 17
: x <= 514278262 ? 18
: x <= 1394199299 ? 19
: 20;
}
constexpr size_type get_estimated_Pi(size_type x) { return x / get_estimated_ln(x); }
template <size_type MAX_RANGE, bool RangeQuery = false>
struct Table {
static constexpr size_type max_pi = get_estimated_Pi(MAX_RANGE);
char m_val[MAX_RANGE + 1];
std::bitset<MAX_RANGE + 1> m_notp;
size_type m_primes[max_pi + 1], m_pcnt;
OptionalArray<MAX_RANGE, RangeQuery> m_mobius_presum;
Table(size_type range = MAX_RANGE) { resize(range); }
void resize(size_type range) {
if (!range) return;
m_val[1] = 1;
m_notp.set(1);
m_pcnt = 0;
m_primes[m_pcnt++] = 2;
for (size_type i = 3; i <= range; i += 2) {
if (!m_notp[i]) m_primes[m_pcnt++] = i, m_val[i] = -1;
for (size_type j = 1; j != m_pcnt && i * m_primes[j] <= range; j++) {
m_notp.set(i * m_primes[j]);
if (i % m_primes[j] == 0) break;
m_val[i * m_primes[j]] = -m_val[i];
}
}
if constexpr (RangeQuery) m_mobius_presum.prepare_presum(range, [&](size_type i) { return query_mobius(i); });
m_primes[m_pcnt] = std::numeric_limits<size_type>::max() / 2;
}
int query_mobius(size_type x) const { return (x & 3) ? ((x & 1) ? m_val[x] : -m_val[x >> 1]) : 0; }
int query_mobius(size_type left, size_type right) const {
static_assert(RangeQuery, "RangeQuery Must Be True");
return m_mobius_presum.query(left, right);
}
bool is_prime(size_type x) const { return (x & 1) ? !m_notp[x] : x == 2; }
size_type query_kth_prime(size_type k) const { return m_primes[k]; }
};
template <size_type MAX_RANGE>
struct Multiplicative {
static constexpr size_type max_pi = get_estimated_Pi(MAX_RANGE);
struct node {
size_type m_val, m_cnt, m_low, m_prod;
};
node m_minp[MAX_RANGE / 2];
size_type m_primes[max_pi];
template <typename Tp, typename CalcPrime, typename CalcPrimePow>
void solve(size_type range, Tp *array, CalcPrime &&calc_prime, CalcPrimePow &&calc_prime_pow) {
if (range >= 1) array[1] = Tp(1);
if (range >= 2) array[2] = calc_prime(2);
for (size_type x = 4, c = 2; x <= range; x <<= 1, c++) array[x] = calc_prime_pow(2, c, x >> 1);
for (size_type i = 3, odd_pcnt = 0; i <= range; i += 2) {
if (!m_minp[i / 2].m_val) {
m_primes[odd_pcnt++] = i;
m_minp[i / 2] = {i, 1, 1, i};
array[i] = calc_prime(i);
} else
array[i] = m_minp[i / 2].m_prod == i ? calc_prime_pow(m_minp[i / 2].m_val, m_minp[i / 2].m_cnt, m_minp[i / 2].m_low) : array[m_minp[i / 2].m_prod] * array[m_minp[i / 2].m_low];
for (size_type j = 0; j != odd_pcnt && i * m_primes[j] <= range; j++) {
size_type p = m_primes[j];
m_minp[i * p / 2].m_val = p;
if (i % p == 0) {
m_minp[i * p / 2].m_cnt = m_minp[i / 2].m_cnt + 1;
m_minp[i * p / 2].m_low = m_minp[i / 2].m_prod == i ? i : m_minp[i / 2].m_low;
m_minp[i * p / 2].m_prod = m_minp[i / 2].m_prod * p;
break;
}
m_minp[i * p / 2].m_cnt = 1;
m_minp[i * p / 2].m_low = i;
m_minp[i * p / 2].m_prod = p;
}
if (i + 1 <= range) {
size_type ctz = std::countr_zero(i + 1);
if (i + 1 != 1 << ctz) array[i + 1] = array[1 << ctz] * array[i + 1 >> ctz];
}
}
}
template <typename Tp, typename CalcPrime, typename CalcPrimePow>
std::vector<Tp> solve(size_type range, CalcPrime &&calc_prime, CalcPrimePow &&calc_prime_pow) {
std::vector<Tp> res(range + 1);
solve(range, res.data(), calc_prime, calc_prime_pow);
return res;
}
};
template <typename Tp, typename FindKthPrime, typename Operation = CompoundPlus>
inline void partial_sum_Dirichlet_divisor(size_type range, Tp *array, FindKthPrime &&find_kth_prime, Operation &&op = Operation()) {
for (size_type i = 0;; i++) {
const size_type p = find_kth_prime(i);
if (p > range) break;
for (size_type j = 1, k = p; k <= range; j++, k += p) op(array[k], array[j]);
}
}
template <typename Tp, typename FindKthPrime, typename Operation = CompoundMinus>
inline void adjacent_difference_Dirichlet_divisor(size_type range, Tp *array, FindKthPrime &&find_kth_prime, Operation &&op = Operation()) {
for (size_type i = 0;; i++) {
const size_type p = find_kth_prime(i);
if (p > range) break;
for (size_type j = range / p, k = j * p; k; j--, k -= p) op(array[k], array[j]);
}
}
template <typename Tp, typename FindKthPrime, typename Operation = CompoundPlus>
inline void partial_sum_Dirichlet_multiple(size_type range, Tp *array, FindKthPrime &&find_kth_prime, Operation &&op = Operation()) {
for (size_type i = 0;; i++) {
const size_type p = find_kth_prime(i);
if (p > range) break;
for (size_type j = range / p, k = j * p; k; j--, k -= p) op(array[j], array[k]);
}
}
template <typename Tp, typename FindKthPrime, typename Operation = CompoundMinus>
inline void adjacent_difference_Dirichlet_multiple(size_type range, Tp *array, FindKthPrime &&find_kth_prime, Operation &&op = Operation()) {
for (size_type i = 0;; i++) {
const size_type p = find_kth_prime(i);
if (p > range) break;
for (size_type j = 1, k = p; k <= range; j++, k += p) op(array[j], array[k]);
}
}
}
}
namespace OY {
template <typename Tp, typename = typename std::enable_if<std::is_unsigned<Tp>::value>::type>
struct SqrtDecomposition {
struct iterator {
struct _range {
Tp m_quot, m_left, m_right;
Tp quot() const { return m_quot; }
Tp left() const { return m_left; }
Tp right() const { return m_right; }
template <typename Ostream>
friend Ostream &operator<<(Ostream &out, const _range &x) { return out << '(' << x.m_quot << ": " << x.m_left << '~' << x.m_right << ')'; }
} m_range;
Tp m_val;
iterator &operator--() {
if (m_range.m_quot == 1)
m_range = {0, m_val + 1, Tp(-1)};
else {
m_range.m_left = m_range.m_right + 1;
m_range.m_quot = m_val / m_range.m_left;
if (m_range.m_quot >= m_range.m_left)
m_range.m_right = m_range.m_left;
else
m_range.m_right = m_val / m_range.m_quot;
}
return *this;
}
iterator &operator++() {
m_range.m_right = m_range.m_left - 1;
if (m_range.m_quot < m_range.m_right)
m_range.m_left = m_val / (++m_range.m_quot + 1) + 1;
else if (m_range.m_right)
m_range.m_quot = m_val / --m_range.m_left;
else
m_range.m_quot = -1, m_range.m_left = 0;
return *this;
}
const _range &operator*() const { return m_range; }
const _range *operator->() const { return &m_range; }
bool operator!=(const iterator &rhs) const { return m_range.m_quot != rhs.m_range.m_quot; }
bool operator==(const iterator &rhs) const { return m_range.m_quot == rhs.m_range.m_quot; }
static iterator null(Tp val) { return {{Tp(-1)}, val}; }
};
struct reverse_iterator : iterator {
iterator base() const { return *this; }
reverse_iterator &operator--() {
iterator::operator++();
return *this;
}
reverse_iterator &operator++() {
iterator::operator--();
return *this;
}
static reverse_iterator null(Tp val) {
reverse_iterator res;
res.m_range = {0, val + 1, Tp(-1)}, res.m_val = val;
return res;
}
};
Tp m_val;
SqrtDecomposition(Tp val) : m_val(val) {}
Tp size() const {
Tp r;
if constexpr (sizeof(Tp) == 8) {
#ifdef _MSC_VER
r = std::sqrt(val);
r -= r * r > val;
#else
r = std::sqrt((long double)m_val);
#endif
} else
r = std::sqrt((double)m_val);
return r * 2 - (m_val * 2 + 1 < (r + 1) * r * 2);
}
iterator lower_bound(Tp quot) const {
if (!quot) return rend().base();
if (quot > m_val) return end();
Tp right = m_val / quot;
if (quot < right) {
Tp left = m_val / (quot + 1) + 1;
return iterator{{quot, left, right}, m_val};
} else
return iterator{{m_val / right, right, right}, m_val};
}
iterator upper_bound(Tp quot) const { return lower_bound(quot + 1); }
iterator find_by_divisor(Tp divisor) const {
if (!divisor) return end();
if (divisor > m_val) return rend().base();
Tp quot = m_val / divisor;
if (divisor <= quot) return iterator{{quot, divisor, divisor}, m_val};
Tp right = m_val / quot;
Tp left = m_val / (quot + 1) + 1;
return iterator{{quot, left, right}, m_val};
}
iterator begin() const { return iterator{{1, m_val / 2 + 1, m_val}, m_val}; }
iterator end() const { return iterator::null(m_val); }
reverse_iterator rbegin() const {
reverse_iterator res;
res.m_range = {m_val, 1, 1}, res.m_val = m_val;
return res;
}
reverse_iterator rend() const { return reverse_iterator::null(m_val); }
};
template <typename Ostream, typename Tp>
Ostream &operator<<(Ostream &out, const SqrtDecomposition<Tp> &x) {
out << '{';
for (auto i : x) {
if (i.m_quot != 1) out << ", ";
out << i;
}
return out << '}';
}
}
using llt = long long;
int const SZ = 50000 + 1;
vector<llt> Vec(SZ, 0);
OY::MOBIUS::Table<SZ, true> Mobius;
void init(){
for(int i=1;i<SZ;++i){
llt ans = 0;
for(const auto & p : OY::SqrtDecomposition<uint32_t>(i)){
ans += (p.right() - p.left() + 1LL) * p.quot();
}
Vec[i] = ans;
}
return;
}
llt N, M;
void work(){
cin >> N >> M;
llt ans = 0;
for(int r,l=1,n=min(N,M);l<=n;l=r+1){
r = min(N/(N/l), M/(M/l));
ans += Mobius.query_mobius(l, r) * Vec[N / l] * Vec[M / l];
}
cout << ans << endl;
}
int main(){
#ifndef ONLINE_JUDGE
freopen("z.txt", "r", stdin);
#endif
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
init();
int nofkase = 1;
cin >> nofkase;
while(nofkase--) work();
return 0;
}
洛谷P3455
求：
∑
i
=
1
N
∑
j
=
1
M
[
g
c
d
(
i
,
j
)
=
=
K
]
\sum_{i=1}^{N}\sum_{j=1}^{M}[gcd(i,j)==K]
i=1∑Nj=1∑M[gcd(i,j)==K]
其中，
N
,
M
,
K
N,M,K
N,M,K是固定值。原式简单转化一下，变为：
∑
i
=
1
N
K
∑
j
=
1
M
K
[
g
c
d
(
i
,
j
)
=
=
1
]
\sum_{i=1}^{\frac{N}{K}}\sum_{j=1}^{\frac{M}{K}}[gcd(i,j)==1]
i=1∑KNj=1∑KM[gcd(i,j)==1]
令 f ( x ) = [ x = = 1 ] , g ( x ) ≡ 1 f(x)=[x==1],g(x)\equiv{1} f(x)=[x==1],g(x)≡1，套用互质枚举的公式得到：
∑ d = 1 N ∑ i = 1 N K d ∑ j = 1 M K d μ ( d ) \sum_{d=1}^{N}\sum_{i=1}^{\frac{N}{Kd}}\sum_{j=1}^{\frac{M}{Kd}}\mu(d) d=1∑Ni=1∑KdNj=1∑KdMμ(d)
很显然三个求和下标互不相关，可以分开，原式变为：
∑ d = 1 N μ ( d ) ∑ i = 1 N K d 1 ∑ j = 1 M K d 1 \sum_{d=1}^{N}\mu(d)\sum_{i=1}^{\frac{N}{Kd}}\mathit{1}\sum_{j=1}^{\frac{M}{Kd}}\mathit{1} d=1∑Nμ(d)i=1∑KdN1j=1∑KdM1
因为是多case，所以预处理莫比乌斯函数的和，然后整除分块求每个case即可.
#include <bits/stdc++.h>
using namespace std;
#ifdef _MSC_VER
#define CPP_STANDARD _MSVC_LANG
#else
#define CPP_STANDARD __cplusplus
#endif
#if CPP_STANDARD >= 201402L
#define CONSTEXPR_IF_CPP14 constexpr
#else
#define CONSTEXPR_IF_CPP14
#endif
namespace std {
namespace my_bit_ops {
static constexpr int s32[32] = {31, 0, 1, 5, 2, 10, 6, 15, 3, 13, 11, 20, 7, 22, 16, 25, 30, 4, 9, 14, 12, 19, 21, 24, 29, 8, 18, 23, 28, 17, 27, 26};
static constexpr int s64[64] = {63, 0, 1, 6, 2, 12, 7, 18, 3, 24, 13, 27, 8, 33, 19, 39, 4, 16, 25, 37, 14, 45, 28, 47, 9, 30, 34, 53, 20, 49, 40, 56, 62, 5, 11, 17, 23, 26, 32, 38, 15, 36, 44, 46, 29, 52, 48, 55, 61, 10, 22, 31, 35, 43, 51, 54, 60, 21, 42, 50, 59, 41, 58, 57};
static constexpr int s256[256] = {8, 7, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
static constexpr int ctz32(uint32_t x) noexcept { return s32[(x & -x) * 0b00001000110010100111010110111110 >> 27]; }
static constexpr int ctz64(uint64_t x) noexcept { return s64[(x & -x) * 0b0000010000110001010001110010010110011010011110101011101101111110 >> 58]; }
static constexpr int clz8(uint32_t x) noexcept { return s256[x]; }
static constexpr int clz16(uint32_t x) noexcept { return x >> 8 ? clz8(x >> 8) : 8 + clz8(x); }
static constexpr int clz32(uint32_t x) noexcept { return x >> 16 ? clz16(x >> 16) : 16 + clz16(x); }
static constexpr int clz64(uint64_t x) noexcept { return x >> 32 ? clz32(x >> 32) : 32 + clz32(x); }
static constexpr uint32_t f1_32(uint32_t x) noexcept { return x - ((x >> 1) & 0x55555555); }
static constexpr uint32_t f2_32(uint32_t x) noexcept { return (x & 0x33333333) + ((x >> 2) & 0x33333333); }
static constexpr int f3_32(uint32_t x) noexcept { return (x + (x >> 4) & 0x0f0f0f0f) * 0x01010101 >> 24; }
static constexpr int popcnt32(uint32_t x) noexcept { return f3_32(f2_32(f1_32(x))); }
static constexpr uint64_t f1_64(uint64_t x) noexcept { return x - ((x >> 1) & 0x5555555555555555); }
static constexpr uint64_t f2_64(uint64_t x) noexcept { return (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333); }
static constexpr int f3_64(uint64_t x) noexcept { return (x + (x >> 4) & 0x0f0f0f0f0f0f0f0f) * 0x0101010101010101 >> 56; }
static constexpr int popcnt64(uint64_t x) noexcept { return f3_64(f2_64(f1_64(x))); }
constexpr uint32_t bit_ceil_msvc_32(uint32_t x) noexcept { return x ? uint32_t(1) << (32 - clz32(x - 1)) : 1; }
constexpr uint64_t bit_ceil_msvc_64(uint64_t x) noexcept { return x ? uint64_t(1) << (64 - clz64(x - 1)) : 1; }
template <typename Tp>
constexpr Tp bit_ceil_msvc(Tp x) noexcept { return sizeof(Tp) == 4 ? bit_ceil_msvc_32(x) : bit_ceil_msvc_64(x); }
constexpr uint32_t bit_floor_msvc_32(uint32_t x) noexcept { return x ? uint32_t(1) << (31 - clz32(x)) : 0; }
constexpr uint64_t bit_floor_msvc_64(uint64_t x) noexcept { return x ? uint64_t(1) << (63 - clz64(x)) : 0; }
template <typename Tp>
constexpr Tp bit_floor_msvc(Tp x) noexcept { return sizeof(Tp) == 4 ? bit_floor_msvc_32(x) : bit_floor_msvc_64(x); }
constexpr uint32_t bit_width_msvc_32(uint32_t x) noexcept { return 32 - clz32(x); }
constexpr uint64_t bit_width_msvc_64(uint64_t x) noexcept { return 64 - clz64(x); }
template <typename Tp>
constexpr Tp bit_width_msvc(Tp x) noexcept { return sizeof(Tp) == 4 ? bit_width_msvc_32(x) : bit_width_msvc_64(x); }
constexpr int countl_zero_msvc_32(uint32_t x) noexcept { return x ? clz32(x) : 32; }
constexpr int countl_zero_msvc_64(uint64_t x) noexcept { return x ? clz64(x) : 64; }
template <typename Tp>
constexpr int countl_zero_msvc(Tp x) noexcept { return sizeof(Tp) == 4 ? countl_zero_msvc_32(x) : countl_zero_msvc_64(x); }
constexpr int countr_zero_msvc_32(uint32_t x) noexcept { return x ? ctz32(x) : 32; }
constexpr int countr_zero_msvc_64(uint64_t x) noexcept { return x ? ctz64(x) : 64; }
template <typename Tp>
constexpr int countr_zero_msvc(Tp x) noexcept { return sizeof(Tp) == 4 ? countr_zero_msvc_32(x) : countr_zero_msvc_64(x); }
template <typename Tp>
constexpr int popcount_msvc(Tp x) noexcept { return sizeof(Tp) == 4 ? popcnt32(x) : popcnt64(x); }
#ifdef _MSC_VER
template <typename Tp>
constexpr Tp bit_ceil(Tp x) noexcept { return bit_ceil_msvc(x); }
template <typename Tp>
constexpr Tp bit_floor(Tp x) noexcept { return bit_floor_msvc(x); }
template <typename Tp>
constexpr Tp bit_width(Tp x) noexcept { return bit_width_msvc(x); }
template <typename Tp>
constexpr int countl_zero(Tp x) noexcept { return countl_zero_msvc(x); }
template <typename Tp>
constexpr int countr_zero(Tp x) noexcept { return countr_zero_msvc(x); }
template <typename Tp>
constexpr int popcount(Tp x) noexcept { return popcount_msvc(x); }
#else
constexpr uint32_t bit_ceil_32(uint32_t x) noexcept { return x > 1 ? uint32_t(1) << (32 - __builtin_clz(x - 1)) : 1; }
constexpr uint64_t bit_ceil_64(uint64_t x) noexcept { return x > 1 ? uint64_t(1) << (64 - __builtin_clzll(x - 1)) : 1; }
template <typename Tp>
constexpr Tp bit_ceil(Tp x) noexcept { return sizeof(Tp) == 4 ? bit_ceil_32(x) : bit_ceil_64(x); }
constexpr uint32_t bit_floor_32(uint32_t x) noexcept { return x ? uint32_t(1) << (31 - __builtin_clz(x)) : 0; }
constexpr uint64_t bit_floor_64(uint64_t x) noexcept { return x ? uint64_t(1) << (63 - __builtin_clzll(x)) : 0; }
template <typename Tp>
constexpr Tp bit_floor(Tp x) noexcept { return sizeof(Tp) == 4 ? bit_floor_32(x) : bit_floor_64(x); }
constexpr uint32_t bit_width_32(uint32_t x) noexcept { return x ? 32 - __builtin_clz(x) : 0; }
constexpr uint64_t bit_width_64(uint64_t x) noexcept { return x ? 64 - __builtin_clzll(x) : 0; }
template <typename Tp>
constexpr Tp bit_width(Tp x) noexcept { return sizeof(Tp) == 4 ? bit_width_32(x) : bit_width_64(x); }
constexpr int countl_zero_32(uint32_t x) noexcept { return x ? __builtin_clz(x) : 32; }
constexpr int countl_zero_64(uint64_t x) noexcept { return x ? __builtin_clzll(x) : 64; }
template <typename Tp>
constexpr int countl_zero(Tp x) noexcept { return sizeof(Tp) == 4 ? countl_zero_32(x) : countl_zero_64(x); }
constexpr int countr_zero_32(uint32_t x) noexcept { return x ? __builtin_ctz(x) : 32; }
constexpr int countr_zero_64(uint64_t x) noexcept { return x ? __builtin_ctzll(x) : 64; }
template <typename Tp>
constexpr int countr_zero(Tp x) noexcept { return sizeof(Tp) == 4 ? countr_zero_32(x) : countr_zero_64(x); }
template <typename Tp>
constexpr int popcount(Tp x) noexcept { return sizeof(Tp) == 4 ? __builtin_popcount(x) : __builtin_popcountll(x); }
#endif
}
#ifndef __cpp_lib_bitops
using my_bit_ops::bit_ceil;
using my_bit_ops::bit_floor;
using my_bit_ops::bit_width;
using my_bit_ops::countl_zero;
using my_bit_ops::countr_zero;
using my_bit_ops::popcount;
#endif
}
namespace OY {
namespace MOBIUS {
using size_type = uint32_t;
struct CompoundPlus {
template <typename Tp>
void operator()(Tp &x, const Tp &y) const { x += y; }
};
struct CompoundMinus {
template <typename Tp>
void operator()(Tp &x, const Tp &y) const { x -= y; }
};
template <size_type MAX_RANGE, bool Flag>
struct OptionalArray {
size_type m_data[MAX_RANGE + 1];
template <typename Mapping>
void prepare_presum(size_type range, Mapping mapping) {
m_data[0] = 0;
for (size_type i = 1; i <= range; i++) m_data[i] = m_data[i - 1] + mapping(i);
}
void set(size_type i, size_type val) { m_data[i] = val; }
size_type query(size_type left, size_type right) const { return m_data[right] - m_data[left - 1]; }
size_type operator[](size_type i) const { return m_data[i]; }
};
template <size_type MAX_RANGE>
struct OptionalArray<MAX_RANGE, false> {};
constexpr size_type get_estimated_ln(size_type x) {
return x <= 7 ? 1
: x <= 32 ? 2
: x <= 119 ? 3
: x <= 359 ? 4
: x <= 1133 ? 5
: x <= 3093 ? 6
: x <= 8471 ? 7
: x <= 24299 ? 8
: x <= 64719 ? 9
: x <= 175196 ? 10
: x <= 481451 ? 11
: x <= 1304718 ? 12
: x <= 3524653 ? 13
: x <= 9560099 ? 14
: x <= 25874783 ? 15
: x <= 70119984 ? 16
: x <= 189969353 ? 17
: x <= 514278262 ? 18
: x <= 1394199299 ? 19
: 20;
}
constexpr size_type get_estimated_Pi(size_type x) { return x / get_estimated_ln(x); }
template <size_type MAX_RANGE, bool RangeQuery = false>
struct Table {
static constexpr size_type max_pi = get_estimated_Pi(MAX_RANGE);
char m_val[MAX_RANGE + 1];
std::bitset<MAX_RANGE + 1> m_notp;
size_type m_primes[max_pi + 1], m_pcnt;
OptionalArray<MAX_RANGE, RangeQuery> m_mobius_presum;
Table(size_type range = MAX_RANGE) { resize(range); }
void resize(size_type range) {
if (!range) return;
m_val[1] = 1;
m_notp.set(1);
m_pcnt = 0;
m_primes[m_pcnt++] = 2;
for (size_type i = 3; i <= range; i += 2) {
if (!m_notp[i]) m_primes[m_pcnt++] = i, m_val[i] = -1;
for (size_type j = 1; j != m_pcnt && i * m_primes[j] <= range; j++) {
m_notp.set(i * m_primes[j]);
if (i % m_primes[j] == 0) break;
m_val[i * m_primes[j]] = -m_val[i];
}
}
if constexpr (RangeQuery) m_mobius_presum.prepare_presum(range, [&](size_type i) { return query_mobius(i); });
m_primes[m_pcnt] = std::numeric_limits<size_type>::max() / 2;
}
int query_mobius(size_type x) const { return (x & 3) ? ((x & 1) ? m_val[x] : -m_val[x >> 1]) : 0; }
int query_mobius(size_type left, size_type right) const {
static_assert(RangeQuery, "RangeQuery Must Be True");
return m_mobius_presum.query(left, right);
}
bool is_prime(size_type x) const { return (x & 1) ? !m_notp[x] : x == 2; }
size_type query_kth_prime(size_type k) const { return m_primes[k]; }
};
template <size_type MAX_RANGE>
struct Multiplicative {
static constexpr size_type max_pi = get_estimated_Pi(MAX_RANGE);
struct node {
size_type m_val, m_cnt, m_low, m_prod;
};
node m_minp[MAX_RANGE / 2];
size_type m_primes[max_pi];
template <typename Tp, typename CalcPrime, typename CalcPrimePow>
void solve(size_type range, Tp *array, CalcPrime &&calc_prime, CalcPrimePow &&calc_prime_pow) {
if (range >= 1) array[1] = Tp(1);
if (range >= 2) array[2] = calc_prime(2);
for (size_type x = 4, c = 2; x <= range; x <<= 1, c++) array[x] = calc_prime_pow(2, c, x >> 1);
for (size_type i = 3, odd_pcnt = 0; i <= range; i += 2) {
if (!m_minp[i / 2].m_val) {
m_primes[odd_pcnt++] = i;
m_minp[i / 2] = {i, 1, 1, i};
array[i] = calc_prime(i);
} else
array[i] = m_minp[i / 2].m_prod == i ? calc_prime_pow(m_minp[i / 2].m_val, m_minp[i / 2].m_cnt, m_minp[i / 2].m_low) : array[m_minp[i / 2].m_prod] * array[m_minp[i / 2].m_low];
for (size_type j = 0; j != odd_pcnt && i * m_primes[j] <= range; j++) {
size_type p = m_primes[j];
m_minp[i * p / 2].m_val = p;
if (i % p == 0) {
m_minp[i * p / 2].m_cnt = m_minp[i / 2].m_cnt + 1;
m_minp[i * p / 2].m_low = m_minp[i / 2].m_prod == i ? i : m_minp[i / 2].m_low;
m_minp[i * p / 2].m_prod = m_minp[i / 2].m_prod * p;
break;
}
m_minp[i * p / 2].m_cnt = 1;
m_minp[i * p / 2].m_low = i;
m_minp[i * p / 2].m_prod = p;
}
if (i + 1 <= range) {
size_type ctz = std::countr_zero(i + 1);
if (i + 1 != 1 << ctz) array[i + 1] = array[1 << ctz] * array[i + 1 >> ctz];
}
}
}
template <typename Tp, typename CalcPrime, typename CalcPrimePow>
std::vector<Tp> solve(size_type range, CalcPrime &&calc_prime, CalcPrimePow &&calc_prime_pow) {
std::vector<Tp> res(range + 1);
solve(range, res.data(), calc_prime, calc_prime_pow);
return res;
}
};
template <typename Tp, typename FindKthPrime, typename Operation = CompoundPlus>
inline void partial_sum_Dirichlet_divisor(size_type range, Tp *array, FindKthPrime &&find_kth_prime, Operation &&op = Operation()) {
for (size_type i = 0;; i++) {
const size_type p = find_kth_prime(i);
if (p > range) break;
for (size_type j = 1, k = p; k <= range; j++, k += p) op(array[k], array[j]);
}
}
template <typename Tp, typename FindKthPrime, typename Operation = CompoundMinus>
inline void adjacent_difference_Dirichlet_divisor(size_type range, Tp *array, FindKthPrime &&find_kth_prime, Operation &&op = Operation()) {
for (size_type i = 0;; i++) {
const size_type p = find_kth_prime(i);
if (p > range) break;
for (size_type j = range / p, k = j * p; k; j--, k -= p) op(array[k], array[j]);
}
}
template <typename Tp, typename FindKthPrime, typename Operation = CompoundPlus>
inline void partial_sum_Dirichlet_multiple(size_type range, Tp *array, FindKthPrime &&find_kth_prime, Operation &&op = Operation()) {
for (size_type i = 0;; i++) {
const size_type p = find_kth_prime(i);
if (p > range) break;
for (size_type j = range / p, k = j * p; k; j--, k -= p) op(array[j], array[k]);
}
}
template <typename Tp, typename FindKthPrime, typename Operation = CompoundMinus>
inline void adjacent_difference_Dirichlet_multiple(size_type range, Tp *array, FindKthPrime &&find_kth_prime, Operation &&op = Operation()) {
for (size_type i = 0;; i++) {
const size_type p = find_kth_prime(i);
if (p > range) break;
for (size_type j = 1, k = p; k <= range; j++, k += p) op(array[j], array[k]);
}
}
}
}
using llt = long long;
int const SZ = 50001;
OY::MOBIUS::Table<SZ, true> Mobius;
llt N, M, K;
void work(){
cin >> N >> M >> K;
if(N > M) swap(N, M);
llt ans = 0;
llt nk = N / K, mk = M / K;
for(int j,i=1;i<=N;i=j+1){
llt x = N / i, y = M / i;
j = min(N/x, M/y);
x = nk / i, y = mk / i;
ans += Mobius.query_mobius(i, j) * x * y;
}
cout << ans << "\n";
}
int main(){
#ifndef ONLINE_JUDGE
freopen("z.txt", "r", stdin);
#endif
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
int nofkase = 1;
cin >> nofkase;
while(nofkase--) work();
return 0;
}
洛谷2257
求：
∑
i
=
1
N
∑
j
=
1
M
[
g
c
d
(
i
,
j
)
i
s
p
r
i
m
e
]
\sum_{i=1}^{N}\sum_{j=1}^{M}[gcd(i,j)\ is\ prime]
i=1∑Nj=1∑M[gcd(i,j) is prime]
同样套用互质枚举，原式变为
∑ d = 1 N ∑ i = 1 N / d ∑ j = 1 M / d ∑ t ∣ d [ t i s p r i m e ] μ ( d t ) \sum_{d=1}^{N}\sum_{i=1}^{N/d}\sum_{j=1}^{M/d}\sum_{t|d}[t\ is\ prime]\mu(\frac{d}{t}) d=1∑Ni=1∑N/dj=1∑M/dt∣d∑[t is prime]μ(td)
将与 i , j i,j i,j无关的提前，且 i , j i,j i,j也可以分离，得到：
∑ d = 1 N ∑ t ∣ d [ t i s p r i m e ] μ ( d t ) ∑ i = 1 N / d 1 ∑ j = 1 M / d 1 \sum_{d=1}^{N}\sum_{t|d}[t\ is\ prime]\mu(\frac{d}{t})\sum_{i=1}^{N/d}\mathbb{1}\sum_{j=1}^{M/d}\mathbb{1} d=1∑Nt∣d∑[t is prime]μ(td)i=1∑N/d1j=1∑M/d1
注意到第二个求和符号专门为质数准备，因此可以使用筛法预处理。对于多个case的情况下，先处理和再使用整除分块即可。
#include <bits/stdc++.h>
using namespace std;
using llt = long long;
using vi = vector<int>;
using si = set<int>;
using vll = vector<llt>;
int const SZ = 10000001;
/// Fx[i] = SIGMA{mu(i/p)}
vector<llt> Fx;
vector<llt> Sum;
vector<bool> isComp; // isComp[i]表示i是否为合数
vector<llt> Primes; // Primes[i]是第i个质数
vector<int> Mobius; // 莫比乌斯函数
void init(){
Primes.reserve(SZ);
isComp.assign(SZ, false);
Mobius.assign(SZ, 0);
Mobius[1] = 1;
long long tmp;
for(int i=2;i<SZ;++i){
if(!isComp[i]){
Primes.push_back(i);
Mobius[i] = -1;
}
for(int j=0,n=Primes.size();j<n&&(tmp=i*(long long)Primes[j])<SZ;++j){
isComp[tmp] = true;
if(0 == i % Primes[j]){
Mobius[tmp] = 0;
break;
}else{
Mobius[tmp] = -Mobius[i];
}
}
}
Fx.assign(SZ, 0);
for(auto p : Primes){
for(auto i=p;i<SZ;i+=p){
Fx[i] += Mobius[i / p];
}
}
Sum.assign(SZ, 0);
for(int i=1;i<SZ;++i){
Sum[i] = Sum[i - 1] + Fx[i];
}
return;
}
int N, M;
void work(){
cin >> N >> M;
if(N > M) swap(N, M);
llt ans = 0;
for(int j,i=1;i<=N;i=j+1){
auto x = N / i, y = M / i;
j = min(N / x, M / y);
ans += (Sum[j] - Sum[i - 1]) * x * y;
}
cout << ans << "\n";
return;
}
int main(){
#ifndef ONLINE_JUDGE
freopen("z.txt", "r", stdin);
#endif
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
init();
int nofkase = 1;
cin >> nofkase;
while(nofkase--) work();
return 0;
}
洛谷P1829
求：
∑ i = 1 N ∑ j = 1 M l c m ( i , j ) \sum_{i=1}^{N}\sum_{j=1}^{M}lcm(i,j) i=1∑Nj=1∑Mlcm(i,j)
将原式变为：
∑ i = 1 N ∑ j = 1 M i ⋅ j g c d ( i , j ) \sum_{i=1}^{N}\sum_{j=1}^{M}\frac{i\cdot{j}}{gcd(i,j)} i=1∑Nj=1∑Mgcd(i,j)i⋅j
利用互质枚举公式，令
f
(
x
)
=
1
x
f(x)=\frac{1}{x}
f(x)=x1，
g
(
i
,
j
)
=
i
⋅
j
g(i,j)=i\cdot{j}
g(i,j)=i⋅j，套公式有
∑
d
=
1
N
∑
i
=
1
N
/
d
∑
j
=
1
M
/
d
i
d
⋅
j
d
⋅
∑
t
∣
d
1
t
⋅
μ
(
d
t
)
\sum_{d=1}^{N}\sum_{i=1}^{N/d}\sum_{j=1}^{M/d}{id\cdot{jd}}\cdot\sum_{t|d}\frac{1}{t}\cdot{\mu(\frac{d}{t})}
d=1∑Ni=1∑N/dj=1∑M/did⋅jd⋅t∣d∑t1⋅μ(td)
令
x
=
d
t
x=\frac{d}{t}
x=td，最后一个求和符号可以改写为
∑
x
∣
d
x
d
⋅
μ
(
x
)
\sum_{x|d}\frac{x}{d}\cdot{\mu(x)}
∑x∣ddx⋅μ(x)。注意到
d
d
d只与第一个求和符号的下标有关，因此可以提到最前。同时最后一个求和符号连带其和式均可以提到
i
,
j
i,j
i,j之前；
i
i
i本身也可以提到
j
j
j之前。因此原式变为：
∑
d
=
1
N
(
d
⋅
∑
x
∣
d
x
μ
(
x
)
⋅
∑
i
i
⋅
∑
j
j
)
\sum_{d=1}^{N}\big(d\cdot{\sum_{x|d}x\mu(x)\cdot{\sum_{i}i\cdot{\sum_{j}j}}}\big)
d=1∑N(d⋅x∣d∑xμ(x)⋅i∑i⋅j∑j)
注意到后两个求和符号已经非常简单了，因此直接省略上下限，简单记录一下。再考虑第二个求和符合，令
F
(
x
)
=
∑
i
∣
x
i
⋅
μ
(
i
)
F(x)=\sum_{i|x}i\cdot{\mu(i)}
F(x)=i∣x∑i⋅μ(i)
可以证明其是积性函数，有
F
(
1
)
=
1
F(1)=1
F(1)=1
F
(
p
)
=
F
(
p
k
)
=
1
−
p
F(p)=F(p^k)=1-p
F(p)=F(pk)=1−p
F
(
p
1
⋅
p
2
)
=
F
(
p
1
)
⋅
F
(
p
2
)
F(p_1\cdot{p_2})=F(p_1)\cdot{F(p_2)}
F(p1⋅p2)=F(p1)⋅F(p2)
因此 F ( x ) F(x) F(x)的值可以筛出来，对于单case的题目就已经可以解了。如果是多个case的情况，可以使用整除分块使得每个case在根号时间内解出。
#include <bits/stdc++.h>
using namespace std;
namespace OY {
template <uint64_t P, bool IsPrime, typename = typename std::enable_if<(P > 1 && P < uint64_t(1) << 63)>::type>
struct StaticModInt64 {
using mint = StaticModInt64<P, IsPrime>;
using mod_type = uint64_t;
mod_type m_val;
static mod_type _mul(mod_type a, mod_type b) {
#ifdef _MSC_VER
mod_type high, low, res;
low = _umul128(a, b, &high);
_udiv128(high, low, mod(), &res);
return res;
#else
int64_t res = a * b - mod_type((long double)(a)*b / mod()) * mod();
if (res < 0)
res += mod();
else if (res >= mod())
res -= mod();
return res;
#endif
}
StaticModInt64() = default;
template <typename Tp, typename std::enable_if<std::is_signed<Tp>::value>::type * = nullptr>
StaticModInt64(Tp val) : m_val{} {
auto x = val % int64_t(mod());
if (x < 0) x += mod();
m_val = x;
}
template <typename Tp, typename std::enable_if<std::is_unsigned<Tp>::value>::type * = nullptr>
StaticModInt64(Tp val) : m_val(val % mod()) {}
static mint raw(mod_type val) {
mint res;
res.m_val = val;
return res;
}
static constexpr mod_type mod() { return P; }
mod_type val() const { return m_val; }
mint pow(uint64_t n) const {
mod_type res = 1, b = m_val;
while (n) {
if (n & 1) res = _mul(res, b);
b = _mul(b, b), n >>= 1;
}
return raw(res);
}
mint inv() const {
if (IsPrime)
return inv_Fermat();
else
return inv_exgcd();
}
mint inv_exgcd() const {
mod_type x = mod(), y = m_val, m0 = 0, m1 = 1;
while (y) {
mod_type z = x / y;
x -= y * z, m0 -= m1 * z, std::swap(x, y), std::swap(m0, m1);
}
if (m0 >= mod()) m0 += mod() / x;
return raw(m0);
}
mint inv_Fermat() const { return pow(mod() - 2); }
mint &operator++() {
if (++m_val == mod()) m_val = 0;
return *this;
}
mint &operator--() {
if (!m_val) m_val = mod();
m_val--;
return *this;
}
mint operator++(int) {
mint old(*this);
++*this;
return old;
}
mint operator--(int) {
mint old(*this);
--*this;
return old;
}
mint &operator+=(const mint &rhs) {
m_val += rhs.m_val;
if (m_val >= mod()) m_val -= mod();
return *this;
}
mint &operator-=(const mint &rhs) {
m_val += mod() - rhs.m_val;
if (m_val >= mod()) m_val -= mod();
return *this;
}
mint &operator*=(const mint &rhs) {
m_val = _mul(m_val, rhs.m_val);
return *this;
}
mint &operator/=(const mint &rhs) { return *this *= rhs.inv(); }
mint operator+() const { return *this; }
mint operator-() const { return raw(m_val ? mod() - m_val : 0); }
bool operator==(const mint &rhs) const { return m_val == rhs.m_val; }
bool operator!=(const mint &rhs) const { return m_val != rhs.m_val; }
bool operator<(const mint &rhs) const { return m_val < rhs.m_val; }
bool operator>(const mint &rhs) const { return m_val > rhs.m_val; }
bool operator<=(const mint &rhs) const { return m_val <= rhs.m_val; }
bool operator>=(const mint &rhs) const { return m_val <= rhs.m_val; }
template <typename Tp>
explicit operator Tp() const { return Tp(m_val); }
friend mint operator+(const mint &a, const mint &b) { return mint(a) += b; }
friend mint operator-(const mint &a, const mint &b) { return mint(a) -= b; }
friend mint operator*(const mint &a, const mint &b) { return mint(a) *= b; }
friend mint operator/(const mint &a, const mint &b) { return mint(a) /= b; }
};
template <typename Istream, uint64_t P, bool IsPrime>
Istream &operator>>(Istream &is, StaticModInt64<P, IsPrime> &x) { return is >> x.m_val; }
template <typename Ostream, uint64_t P, bool IsPrime>
Ostream &operator<<(Ostream &os, const StaticModInt64<P, IsPrime> &x) { return os << x.val(); }
using mint4611686018427387847 = StaticModInt64<4611686018427387847, true>;
using mint9223372036854775783 = StaticModInt64<9223372036854775783, true>;
}
using llt = long long;
using mint = OY::StaticModInt64<20101009, true>;
int const SZ = 10000001;
struct Sieve{ // 线性筛法
using llt = int;
vector<bool> isComp; // isComp[i]表示i是否为合数
vector<llt> primes; // primes[i]是第i个质数
vector<mint> Fx; // F(x) = SIGMA{ixmu(i),i|x}
Sieve(int SZ = 35000){
primes.reserve(SZ);
isComp.assign(SZ, false);
Fx.assign(SZ, 0);
Fx[1] = 1;
long long tmp;
for(int i=2;i<SZ;++i){
if(!isComp[i]){
primes.push_back(i);
Fx[i] = mint(1 - i);
}
for(int j=0,n=primes.size();j<n&&(tmp=i*(long long)primes[j])<SZ;++j){
isComp[tmp] = true;
if(0 == i % primes[j]){
Fx[tmp] = Fx[i];
break;
}else{
Fx[tmp] = Fx[i] * Fx[primes[j]];
}
}
}
}
}Sie(SZ);
vector<mint> Sum;
void init(){
Sum.assign(SZ, 0);
for(int i=1;i<SZ;++i){
Sum[i] = Sum[i - 1] + mint(i) * Sie.Fx[i];
}
return;
}
llt N, M;
void work(){
cin >> N >> M;
if(N > M) swap(N, M);
mint ans = 0;
for(int j,i=1;i<=N;i=j+1){
llt x = N / i, y = M / i;
j = min(N/x, M/y);
x = x * (x + 1) / 2;
y = y * (y + 1) / 2;
ans += (Sum[j] - Sum[i - 1]) * x * y;
}
cout << ans << endl;
}
int main(){
#ifndef ONLINE_JUDGE
freopen("z.txt", "r", stdin);
#endif
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
init();
int nofkase = 1;
// // cin >> nofkase;
while(nofkase--) work();
return 0;
}
AGC038C
给定数组
A
N
A_N
AN，求
∑
i
=
1
N
−
1
∑
j
=
i
+
1
N
l
c
m
(
A
i
,
A
j
)
\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}lcm(A_i,A_j)
i=1∑N−1j=i+1∑Nlcm(Ai,Aj)
首先将数组元素求和转为值域求和，变为：
1
2
(
∑
i
=
1
L
∑
j
=
1
L
l
c
m
(
i
,
j
)
⋅
c
+
i
⋅
c
+
j
−
∑
i
=
1
N
A
i
)
\frac{1}{2}\big(\sum_{i=1}^{L}\sum_{j=1}^{L}lcm(i,j)\cdot{c+i}\cdot{c+j}-\sum_{i=1}^{N}A_i\big)
21(i=1∑Lj=1∑Llcm(i,j)⋅c+i⋅c+j−i=1∑NAi)
令 f ( x ) = 1 x f(x)=\frac{1}{x} f(x)=x1， g ( i , j ) = i ⋅ j ⋅ c i ⋅ c j g(i,j)=i\cdot{j}\cdot{c_i}\cdot{c_j} g(i,j)=i⋅j⋅ci⋅cj，仍然使用互质枚举，最后转化为：
∑ d = 1 L d ∑ t ∣ d t ⋅ μ ( t ) ∑ i = 1 L / d i ⋅ c i d ∑ j = 1 L / d j ⋅ c j d \sum_{d=1}^{L}d\sum_{t|d}t\cdot{\mu(t)}\sum_{i=1}^{L/d}i\cdot{c_{id}}\sum_{j=1}^{L/d}j\cdot{c_{jd}} d=1∑Ldt∣d∑t⋅μ(t)i=1∑L/di⋅cidj=1∑L/dj⋅cjd
第二个求和符号仍然用筛法，最后两个求和符号其实是一样的，暴力方法即可。
#include <bits/stdc++.h>
using namespace std;
#if __cpp_constexpr >= 201304L
#define CONSTEXPR14 constexpr
#else
#define CONSTEXPR14
#endif
namespace OY {
template <uint32_t P, bool IsPrime, typename = typename std::enable_if<(P > 1 && P < uint32_t(1) << 31)>::type>
struct StaticModInt32 {
using mint = StaticModInt32<P, IsPrime>;
using mod_type = uint32_t;
mod_type m_val;
static constexpr mod_type _reduce_norm(int32_t x) { return x < 0 ? x + mod() : x; }
static constexpr mod_type _mul(mod_type a, mod_type b) { return uint64_t(a) * b % mod(); }
constexpr StaticModInt32() = default;
template <typename Tp, typename std::enable_if<std::is_signed<Tp>::value>::type * = nullptr>
constexpr StaticModInt32(Tp val) : m_val(_reduce_norm(val % int32_t(mod()))) {}
template <typename Tp, typename std::enable_if<std::is_unsigned<Tp>::value>::type * = nullptr>
constexpr StaticModInt32(Tp val) : m_val(val % mod()) {}
static CONSTEXPR14 mint raw(mod_type val) {
mint res{};
res.m_val = val;
return res;
}
static constexpr mod_type mod() { return P; }
constexpr mod_type val() const { return m_val; }
CONSTEXPR14 mint pow(uint64_t n) const {
mod_type res = 1, b = m_val;
while (n) {
if (n & 1) res = _mul(res, b);
b = _mul(b, b), n >>= 1;
}
return raw(res);
}
CONSTEXPR14 mint inv() const {
if constexpr (IsPrime)
return inv_Fermat();
else
return inv_exgcd();
}
CONSTEXPR14 mint inv_exgcd() const {
mod_type x = mod(), y = m_val, m0 = 0, m1 = 1;
while (y) {
mod_type z = x / y;
x -= y * z, m0 -= m1 * z, std::swap(x, y), std::swap(m0, m1);
}
if (m0 >= mod()) m0 += mod() / x;
return raw(m0);
}
constexpr mint inv_Fermat() const { return pow(mod() - 2); }
CONSTEXPR14 mint &operator++() {
if (++m_val == mod()) m_val = 0;
return *this;
}
CONSTEXPR14 mint &operator--() {
if (!m_val) m_val = mod();
m_val--;
return *this;
}
CONSTEXPR14 mint operator++(int) {
mint old(*this);
++*this;
return old;
}
CONSTEXPR14 mint operator--(int) {
mint old(*this);
--*this;
return old;
}
CONSTEXPR14 mint &operator+=(const mint &rhs) {
m_val += rhs.m_val;
if (m_val >= mod()) m_val -= mod();
return *this;
}
CONSTEXPR14 mint &operator-=(const mint &rhs) {
m_val += mod() - rhs.m_val;
if (m_val >= mod()) m_val -= mod();
return *this;
}
CONSTEXPR14 mint &operator*=(const mint &rhs) {
m_val = _mul(m_val, rhs.m_val);
return *this;
}
CONSTEXPR14 mint &operator/=(const mint &rhs) { return *this *= rhs.inv(); }
constexpr mint operator+() const { return *this; }
constexpr mint operator-() const { return raw(m_val ? mod() - m_val : 0); }
constexpr bool operator==(const mint &rhs) const { return m_val == rhs.m_val; }
constexpr bool operator!=(const mint &rhs) const { return m_val != rhs.m_val; }
constexpr bool operator<(const mint &rhs) const { return m_val < rhs.m_val; }
constexpr bool operator>(const mint &rhs) const { return m_val > rhs.m_val; }
constexpr bool operator<=(const mint &rhs) const { return m_val <= rhs.m_val; }
constexpr bool operator>=(const mint &rhs) const { return m_val <= rhs.m_val; }
template <typename Tp>
constexpr explicit operator Tp() const { return Tp(m_val); }
friend CONSTEXPR14 mint operator+(const mint &a, const mint &b) { return mint(a) += b; }
friend CONSTEXPR14 mint operator-(const mint &a, const mint &b) { return mint(a) -= b; }
friend CONSTEXPR14 mint operator*(const mint &a, const mint &b) { return mint(a) *= b; }
friend CONSTEXPR14 mint operator/(const mint &a, const mint &b) { return mint(a) /= b; }
};
template <typename Istream, uint32_t P, bool IsPrime>
Istream &operator>>(Istream &is, StaticModInt32<P, IsPrime> &x) { return is >> x.m_val; }
template <typename Ostream, uint32_t P, bool IsPrime>
Ostream &operator<<(Ostream &os, const StaticModInt32<P, IsPrime> &x) { return os << x.val(); }
using mint998244353 = StaticModInt32<998244353, true>;
using mint1000000007 = StaticModInt32<1000000007, true>;
}
using llt = long long;
using vi = vector<int>;
using vll = vector<llt>;
using mint = OY::mint998244353;
int const SZ = 1E6 + 1;
vi Primes;
vi isComp;
/// Fx = SIGMA{t*mu(t), t|x}
vector<mint> Fx;
void init(){
Primes.reserve(SZ);
isComp.assign(SZ, false);
Fx.assign(SZ, {});
Fx[1] = 1;
long long tmp;
for(int i=2;i<SZ;++i){
if(!isComp[i]){
Primes.push_back(i);
Fx[i] = mint(1 - i);
}
for(int j=0,n=Primes.size();j<n&&(tmp=i*(long long)Primes[j])<SZ;++j){
isComp[tmp] = true;
if(0 == i % Primes[j]){
Fx[tmp] = Fx[i];
break;
}else{
Fx[tmp] = Fx[i] * Fx[Primes[j]];
}
}
}
}
int N;
vi A;
vi Cnt;
void work(){
cin >> N;
A.assign(N, {});
Cnt.assign(SZ, 0);
for(auto & i : A){
cin >> i;
Cnt[i] += 1;
}
mint ans = mint(-accumulate(A.begin(), A.end(), 0LL));
int L = *max_element(A.begin(), A.end());
for(int d=1;d<=L;++d){
auto m = L / d;
auto t = mint(0);
for(int i=1;i<=m&&i*d<=L;i+=1){
t += mint(i) * Cnt[i*d];
}
ans += mint(d) * Fx[d] * t * t;
}
ans /= 2;
cout << ans << endl;
return;
}
int main(){
#ifndef ONLINE_JUDGE
freopen("z.txt", "r", stdin);
#endif
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
init();
int nofkase = 1;
// cin >> nofkase;
while(nofkase--) work();
return 0;
}
洛谷3768
给定正整数
N
N
N和质数
P
P
P，求
(
∑
i
=
1
N
∑
j
=
1
N
i
⋅
j
⋅
g
c
d
(
i
,
j
)
)
m
o
d
p
\big(\sum_{i=1}^{N}\sum_{j=1}^{N}i\cdot{j}\cdot{gcd(i,j)}\big)\mod{p}
(i=1∑Nj=1∑Ni⋅j⋅gcd(i,j))modp
显然利用互质枚举，令 f ( x ) = x f(x)=x f(x)=x， g ( i , j ) = i ⋅ j g(i,j)=i\cdot{j} g(i,j)=i⋅j，并且交换求和符号，有（省略上限）
∑ d d ∑ t ∣ d t ⋅ μ ( t ) ∑ i i ∑ j j \sum_{d}d\sum_{t|d}t\cdot{\mu(t)}\sum_{i}i\sum_{j}j d∑dt∣d∑t⋅μ(t)i∑ij∑j
注意到本题 N N N到 5 × 1 0 8 5\times{10^8} 5×108，因此还需要杜教筛。就此停下，以后再说。
扫一扫
851
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
