Mobius反演方法

<更新提示>

前置知识：建议有一点DirichletDirichletDirichlet卷积的基础，会线性筛求积性函数。本文可能会持续更新。

可以看我以前在博客园的博客。

<正文>

0. 符号及约定

1.1.1. [P][P][P]是指正则表达式，当PPP为truetruetrue时，[P]=1[P]=1[P]=1，当PPP为falsefalsefalse时，[P]=0[P]=0[P]=0。

2.2.2. 约数个数函数定义为τ(n)=∑d∣n1\tau(n)=\sum_{d|n}1τ(n)=∑d∣n1。

3.3.3. 约数和函数定义为σ(n)=∑d∣nd\sigma(n)=\sum_{d|n}dσ(n)=∑d∣nd。

4.4.4. 元函数定义为e(n)=[n=1]e(n)=[n=1]e(n)=[n=1]。

5.5.5. 恒等函数定义为I(n)=1I(n)=1I(n)=1。

6.6.6. 单位函数定义为ϵk(n)=nk\epsilon_k(n)=n^kϵk(n)=nk。

7.7.7. 欧拉函数定义为φ(n)=∑i=1n[gcd⁡(i,n)=1]\varphi(n)=\sum_{i=1}^n[\gcd(i,n)=1]φ(n)=∑i=1n[gcd(i,n)=1]。

8.8.8. 设n=∏picin=\prod p_i^{c_i}n=∏pici，则莫比乌斯函数定义为μ(n)={1n=1(−1)k∀ci=10∃ci>1\mu(n)=\begin{cases}1&n=1\\(-1)^k& \forall c_i=1\\0&\exist c_i>1\end{cases}μ(n)=⎩⎪⎨⎪⎧1(−1)k0n=1∀ci=1∃ci>1。

9.9.9. 对于数论函数fff，若满足gcd⁡(a,b)=1\gcd(a,b)=1gcd(a,b)=1时，有f(ab)=f(a)f(b)f(ab)=f(a)f(b)f(ab)=f(a)f(b)，则称函数fff为积性函数。

10.10.10. 对于两个函数f,gf,gf,g，定义他们的dirichletdirichletdirichlet卷积为：(f×g)(n)=∑d∣nf(d)g(nd)(f\times g)(n)=\sum_{d|n}f(d)g(\frac{n}{d})(f×g)(n)=∑d∣nf(d)g(dn)，函数f,gf,gf,g不必要是积性函数。

1. 欧拉筛法

前置知识里已经提到过，以下给出本文可能会使用的线性筛代码：

int flag[N],Prime[N],cnt;
inline void EularSieve(void)
{for (int i=2;i<=Lim;i++){if ( !flag[i] ) Prime[++cnt] = i;for (int j=1;j<=cnt&&i*Prime[j]<=Lim;j++){flag[ i*Prime[j] ] = true;if ( i % Prime[j] == 0 ) break;}}
}

这是最朴素的欧拉筛法，不过当我们要筛一些比较复杂的积性函数的时候，线性筛方程的推导可能会非常复杂。

因此我们考虑引入一种更好的线性筛法，目的是能够方便的筛出积性函数的值。

首先，我们注意到积性函数的性质：当gcd⁡(a,b)=1\gcd(a,b)=1gcd(a,b)=1时，有f(ab)=f(a)f(b)f(ab)=f(a)f(b)f(ab)=f(a)f(b)。那么对于一个任意的正整数nnn，我们可以用算术基本定理进行分解：

n=∏i=1kpicin=\prod_{i=1}^kp_i^{c_i}n=i=1∏kpici

此时，任意的i,j∈[1,k]i,j\in[1,k]i,j∈[1,k]都满足gcd⁡(pici,pjcj)=1\gcd(p_i^{c_i},p_j^{c_j})=1gcd(pici,pjcj)=1。那么我们就可以得到：

f(n)=f(∏i=1kpici)=∏i=1kf(pici)f(n)=f\left(\prod_{i=1}^kp_i^{c_i}\right)=\prod_{i=1}^kf(p_i^{c_i})f(n)=f(i=1∏kpici)=i=1∏kf(pici)

那么我们就有一个想法：利用定义求出积性函数fff在所有素数幂出的取值，然后直接递推出所有函数值。

具体地说，可以这样递推：

f[1] = 1;
for (int i=2;i<=Lim;i++)if ( i == Pow[i] ) f[i] = calc(p[i],e[i]);else f[i] = f[i/Pow[i]] * f[Pow[i]];

其中p[i]p[i]p[i]代表数iii的最小素因子，e[i]e[i]e[i]代表iii的分解式中p[i]p[i]p[i]这个素因子的指数，calccalccalc就是求素数幂处取值的函数，而Pow[i]=p[i]e[i]Pow[i]=p[i]^{e[i]}Pow[i]=p[i]e[i]，，这样是不是就符合我们的要求了呢？

那么现在我们的问题就是如何求出p,e,Powp,e,Powp,e,Pow这三个数组，幸运的是，他们都可以在线性筛的过程中求。

根据欧拉筛法每次用一个数的最小素因子筛去这个合数，就可以更新这三个数组的值了。

代码如下：

int p[N],e[N],Pow[N],Prime[N],cnt;
inline void EularSieve(void)
{Pow[1] = p[1] = 1 , e[1] = 0;for (int i=2;i<=Lim;i++){if ( p[i] == 0 )p[i] = Pow[i] = Prime[++cnt] = i , e[i] = 1;// 判定一个新素数for (int j=1;j<=cnt&&i*Prime[j]<=Lim;j++){int Next = i * Prime[j];if ( p[i] == Prime[j] ){e[Next] = e[i] + 1;Pow[Next] = Pow[i] * Prime[j];// 有相同的素因子break;}else e[Next] = 1 , Pow[Next] = Prime[j];// 没有相同的素因子}}f[1] = 1;for (int i=2;i<=Lim;i++)if ( i == Pow[i] ) f[i] = calc(p[i],e[i]);else f[i] = f[i/Pow[i]] * f[Pow[i]];
}

2. DirichletDirichletDirichlet卷积

2.1 DirichletDirichletDirichlet卷积的性质

1.1.1. 两个积性函数fff和ggg的dirichletdirichletdirichlet卷积仍为积性函数。

证明：

设有两个积性函数fff和ggg，则它们的dirichletdirichletdirichlet卷积为：

h=f×g=∑d∣nf(d)g(nd)h=f\times g=\sum_{d|n}f(d)g(\frac{n}{d})h=f×g=d∣n∑f(d)g(dn)

对于函数hhh则可以得到：

h(x)h(y)=(∑d1∣xf(d1)g(xd1))(∑d2∣yf(d2)g(yd2))=∑d1∣x,d2∣yf(d1d2)g(xyd1d2)=∑d∣xyf(d)g(xyd)=h(xy)h(x)h(y)=\left (\sum_{d_1|x}f(d_1)g(\frac{x}{d_1})\right)\left(\sum_{d_2|y}f(d_2)g(\frac{y}{d_2})\right) \\ \ \\=\sum_{d_1|x,d_2|y}f(d_1d_2)g\left(\frac{xy}{d_1d_2}\right)\\=\sum_{d|xy}f(d)g\left(\frac{xy}{d}\right)=h(xy)h(x)h(y)=⎝⎛d1∣x∑f(d1)g(d1x)⎠⎞⎝⎛d2∣y∑f(d2)g(d2y)⎠⎞ =d1∣x,d2∣y∑f(d1d2)g(d1d2xy)=d∣xy∑f(d)g(dxy)=h(xy)

故函数hhh为积性函数。

2.2.2. dirichletdirichletdirichlet卷积满足交换律。

证明：

设有数论函数fff和ggg，则有：

f×g=∑d∣nf(n)g(nd)=∑d∣ng(n)f(nd)=g×ff\times g=\sum_{d|n}f(n)g(\frac{n}{d})\\=\sum_{d|n}g(n)f(\frac{n}{d})=g\times ff×g=d∣n∑f(n)g(dn)=d∣n∑g(n)f(dn)=g×f

3.3.3. dirichletdirichletdirichlet卷积满足结合律。

证明：

设有数论函数fff，ggg和hhh，则有：

(f×g)×h=f×g×h=g×h×f=(g×h)×f=f×(g×h)(f\times g)\times h=f\times g \times h\\=g\times h \times f=(g\times h)\times f\\=f\times (g\times h)(f×g)×h=f×g×h=g×h×f=(g×h)×f=f×(g×h)

4.4.4. dirichletdirichletdirichlet卷积满足分配律。

证明：

设有数论函数fff，ggg和hhh，则有：

(g+h)×f=∑d∣n(g(d)+h(d))f(nd)=∑d∣ng(d)f(nd)+∑d∣nh(d)f(nd)=g×f+h×f(g+h)\times f=\sum_{d|n}(g(d)+h(d))f(\frac{n}{d}) \\=\sum_{d|n}g(d)f(\frac{n}{d})+\sum_{d|n}h(d)f(\frac{n}{d}) \\=g\times f+h\times f(g+h)×f=d∣n∑(g(d)+h(d))f(dn)=d∣n∑g(d)f(dn)+d∣n∑h(d)f(dn)=g×f+h×f

2.2 常见的DirichletDirichletDirichlet卷积

1.1.1. f×e=ff\times e=ff×e=f

2.2.2. ϵ=φ×I\epsilon=\varphi \times Iϵ=φ×I

3.3.3. τ=I×I\tau=I\times Iτ=I×I

4.4.4. σ=ϵ×I\sigma=\epsilon\times Iσ=ϵ×I

5.5.5. e=μ×Ie=\mu\times Ie=μ×I

6.6.6. φ=ϵ×μ\varphi=\epsilon\times \muφ=ϵ×μ

7.7.7. σ=τ×φ\sigma=\tau \times \varphiσ=τ×φ

这些卷积的证明多数在『简单积性函数和dirichlet卷积』一文中有，此处篇幅受限，顾略去证明。

3. 下取整函数的性质

1.1.1. 对于i∈[x,⌊k⌊kx⌋⌋]i\in \left [x, \left \lfloor \frac{k}{ \lfloor \frac{k}x{} \rfloor } \right \rfloor \right ]i∈[x,⌊⌊xk⌋k⌋]，⌊ki⌋\lfloor \frac{k}{i} \rfloor⌊ik⌋的值都相等。

证明：

设f(x)=⌊k⌊kx⌋⌋f(x)= \left \lfloor \frac{k}{ \lfloor \frac{k}{x} \rfloor } \right \rfloorf(x)=⌊⌊xk⌋k⌋，显然有f(x)≥⌊k(kx)⌋=xf(x)\geq \left \lfloor \frac{k}{ ( \frac{k}{x} ) } \right \rfloor=xf(x)≥⌊(xk)k⌋=x，则可得⌊kf(x)⌋≤⌊kx⌋\left \lfloor \frac{k}{f(x)} \right \rfloor\leq \left \lfloor \frac{k}{x} \right \rfloor⌊f(x)k⌋≤⌊xk⌋。

从另一方面考虑，则有⌊kf(x)⌋≥⌊kk⌊k/x⌋⌋=⌊kx⌋\left \lfloor \frac{k}{f(x)} \right \rfloor\geq\left \lfloor \frac{k}{ \frac{k}{\left \lfloor k/x \right \rfloor } } \right \rfloor=\lfloor \frac{k}{x} \rfloor⌊f(x)k⌋≥⌊⌊k/x⌋kk⌋=⌊xk⌋，则可得：⌊kf(x)⌋=⌊kx⌋\left \lfloor \frac{k}{f(x)} \right \rfloor=\left \lfloor \frac{k}{x} \right \rfloor⌊f(x)k⌋=⌊xk⌋。

2.2.2. ⌊nk⌋\left\lfloor\frac{n}{k}\right\rfloor⌊kn⌋最多只有2n2\sqrt n2n种取值。

证明：

当k≤nk\leq\sqrt nk≤n时，至多只有n\sqrt nn个kkk，所以对应的取值只有n\sqrt nn种。

当k>nk>\sqrt nk>n时，有1≤⌊nk⌋≤n1\leq\left\lfloor\frac{n}{k}\right\rfloor\leq\sqrt n1≤⌊kn⌋≤n，所以对应的取值也只有n\sqrt nn种。

那么当我们要对某个有关下取整函数的值求和的时候，就可以使用如下的整除分块算法，时间复杂度O(n)O(\sqrt n)O(n)。

inline int calc(void)
{int res = 0;for (int l=1,r;l<=n;l=r+1){r = n/l ? min( n/(n/l) , n ) : n;res += f(n/l);}return res;
}

对于二元的情况，其实也是一样的，我们考虑间断点合并，就可以证明其时间复杂度为O(n+m)O(\sqrt n + \sqrt m)O(n+m)。

inline int calclim(int n,int k,int l) { return k/l ? min( k/(k/l) , n ) : n; }
inline int calc(void)
{int res = 0;for (int l=1,r;l<=min(n,m);l=r+1){r = calclim( min(n,m) , n , l );r = min( r , calclim( min(n,m) , m , r ) );// 此时，区间[l,r]中所有的 [n/l] 都相等,所有的 [m/l] 也都相等。 res += f( n/l , m/l );}return res;
}

3.3.3. 若mmm是正整数，则有⌊⌊x⌋m⌋=⌊xm⌋\left\lfloor\frac{\lfloor x \rfloor}{m}\right\rfloor=\left\lfloor\frac{x}{m}\right\rfloor⌊m⌊x⌋⌋=⌊mx⌋。

证明：

设x=km+r(r<m)x=km+r(r<m)x=km+r(r<m)，则⌊⌊x⌋m⌋=⌊k+⌊r⌋m⌋=k⌊xm⌋=⌊k+rm⌋=k\left\lfloor\frac{\lfloor x \rfloor}{m}\right\rfloor=\left\lfloor k+ \frac{\lfloor r \rfloor}{m}\right\rfloor=k\\ \ \\ \left\lfloor\frac{x}{m}\right\rfloor=\left\lfloor k+ \frac{r}{m}\right\rfloor=k⌊m⌊x⌋⌋=⌊k+m⌊r⌋⌋=k ⌊mx⌋=⌊k+mr⌋=k

推论： ⌊⌊kn⌋m⌋=⌊knm⌋\left\lfloor\frac{\lfloor \frac{k}{n} \rfloor}{m}\right\rfloor=\left\lfloor\frac{k}{nm}\right\rfloor⌊m⌊nk⌋⌋=⌊nmk⌋

4. 反演原理

反演形式111：

∑i=1n∑j=1m[gcd⁡(i,j)=1]⟺∑i=1n∑j=1m∑d∣gcd⁡(i,j)μ(d)\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=1]\Longleftrightarrow\sum_{i=1}^n\sum_{j=1}^m\sum_{d|\gcd(i,j)}\mu(d)i=1∑nj=1∑m[gcd(i,j)=1]⟺i=1∑nj=1∑md∣gcd(i,j)∑μ(d)

证明：由e=μ×Ie=\mu \times Ie=μ×I可得。

反演形式222（因数形式莫比乌斯定理）：

F(n)=∑d∣nf(d)⟺f(n)=∑d∣nμ(d)F(nd)F(n)=\sum_{d|n}f(d)\Longleftrightarrow f(n)=\sum_{d|n}\mu(d)F\left(\frac{n}{d}\right)F(n)=d∣n∑f(d)⟺f(n)=d∣n∑μ(d)F(dn)

证明1：dirichletdirichletdirichlet卷积

已知F=I×fF=I \times fF=I×f，则有μ×F=μ×I×F\mu\times F=\mu\times I\times Fμ×F=μ×I×F，故μ×F=e×f=f\mu\times F=e\times f=fμ×F=e×f=f。

证明2：代数反演

∑d∣nμ(d)F(nd)=∑d∣nμ(d)∑d′∣ndf(d′)=∑d′∣nf(d′)∑d∣nd′μ(d)\sum_{d|n}\mu(d)F\left(\frac{n}{d}\right)=\sum_{d|n}\mu(d)\sum _{d'|\frac{n}{d}}f\left(d'\right)\\ \ \\ =\sum_{d'|n}f(d')\sum _{d|\frac{n}{d'}}\mu\left(d\right)d∣n∑μ(d)F(dn)=d∣n∑μ(d)d′∣dn∑f(d′) =d′∣n∑f(d′)d∣d′n∑μ(d)

由μ×I=e\mu\times I=eμ×I=e可知∑d∣nd′μ(d)\sum_{d|\frac{n}{d'}}\mu(d)∑d∣d′nμ(d)非000当且仅当nd′=1\frac{n}{d'}=1d′n=1，所以d′=nd'=nd′=n，原式即为f(n)f(n)f(n)。

反演形式333（倍数形式莫比乌斯定理）：

F(n)=∑n∣df(d)⟺f(n)=∑n∣dμ(dn)F(d)F(n)=\sum_{n|d}f(d)\Longleftrightarrow f(n)=\sum_{n|d}\mu\left(\frac{d}{n}\right)F(d)F(n)=n∣d∑f(d)⟺f(n)=n∣d∑μ(nd)F(d)

证明：

好像没有找到dirichletdirichletdirichlet卷积的证明，就先写一个代数证明吧。

设dn=k\frac{d}{n}=knd=k，则可以得到：

∑n∣dμ(dn)F(d)=∑k=1∞μ(k)F(nk)=∑k=1∞μ(k)∑nk∣tf(t)=∑n∣tf(t)∑k∣tnμ(k)\sum_{n|d}\mu\left(\frac{d}{n}\right)F(d)=\sum_{k=1}^{\infty}\mu(k)F(nk)\\ \ \\ =\sum_{k=1}^{\infty}\mu(k)\sum_{nk|t}f(t)=\sum_{n|t}f(t)\sum_{k|\frac{t}{n}}\mu(k)n∣d∑μ(nd)F(d)=k=1∑∞μ(k)F(nk) =k=1∑∞μ(k)nk∣t∑f(t)=n∣t∑f(t)k∣nt∑μ(k)

由μ×I=e\mu\times I=eμ×I=e可知∑k∣tnμ(k)\sum_{k|\frac{t}{n}}\mu(k)∑k∣ntμ(k)非000当且仅当tn=1\frac{t}{n}=1nt=1，所以t=nt=nt=n，原式即为f(n)f(n)f(n)。

5. 例题

5.1 前言

好了，终于到了例题环节，这一部分对于真正学会莫比乌斯反演才是最重要的。这一节将会通过例题着重讲解反演原理的运用。对于反演原理一节，形式1,31,31,3最为重要，即使其本质基本相同，也希望读者能够时刻牢记。

5.2 BZOJ2301 Problem b

有TTT组询问，求

∑i=ab∑j=cd[gcd⁡(i,j)=k]\sum_{i=a}^b\sum_{j=c}^d[\gcd(i,j)=k]i=a∑bj=c∑d[gcd(i,j)=k]

a,b,c,d,T,k≤50000a,b,c,d,T,k\leq 50000a,b,c,d,T,k≤50000。

首先，我们不妨将这个式子看成一个数表，第iii行第jjj列的格子符合要求gcd⁡(i,j)=k\gcd(i,j)=kgcd(i,j)=k，就有111的价值。于是原式求的就是数表的子矩阵(a,c)−(b,d)(a,c)-(b,d)(a,c)−(b,d)的价值和。

这样我们就可以用二维前缀和的思想容斥

矩阵的价值和。即求：

∑i=1n∑j=1m[gcd⁡(i,j)=k]\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=k]i=1∑nj=1∑m[gcd(i,j)=k]

想必这个函数没法直接求吧，我们可以考虑运用倍数形式莫比乌斯定理，设：

f(k)=∑i=1n∑j=1m[gcd⁡(i,j)=k],F(x)=∑x∣kf(k)f(k)=\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=k],F(x)=\sum_{x|k}f(k)f(k)=i=1∑nj=1∑m[gcd(i,j)=k],F(x)=x∣k∑f(k)

思考一下F(x)F(x)F(x)是什么？函数fff在xxx所有倍数处的取值之和。也就是说一对(i,j)(i,j)(i,j)的最大公约数只要是xxx的倍数，就会对F(x)F(x)F(x)有一个111的贡献。

我们知道x∣gcd⁡(i,j)x|\gcd(i,j)x∣gcd(i,j)等价于i,ji,ji,j都是xxx的倍数，既然i∈[1,n],j∈[1,m]i\in[1,n],j\in[1,m]i∈[1,n],j∈[1,m]，那么函数FFF的值就很好求了。

F(x)=∑x∣kf(k)=⌊nx⌋⌊mx⌋F(x)=\sum_{x|k}f(k)=\left\lfloor\frac{n}{x}\right\rfloor\left\lfloor\frac{m}{x}\right\rfloorF(x)=x∣k∑f(k)=⌊xn⌋⌊xm⌋

我们知道了函数FFF的求法，就可以大胆地套用倍数形式莫比乌斯定理了：

f(k)=∑k∣xμ(xk)F(x)=∑k∣xμ(xk)⌊nx⌋⌊mx⌋f(k)=\sum_{k|x}\mu\left(\frac{x}{k}\right)F(x)=\sum_{k|x}\mu\left(\frac{x}{k}\right)\left\lfloor\frac{n}{x}\right\rfloor\left\lfloor\frac{m}{x}\right\rfloorf(k)=k∣x∑μ(kx)F(x)=k∣x∑μ(kx)⌊xn⌋⌊xm⌋

我们可以枚举kkk的倍数O(n)O(n)O(n)计算答案了，但似乎这样还会超时。

回想起整除分块的优化了吗？我们不妨设t=xkt=\frac{x}{k}t=kx，那么就有：

f(k)=∑t=1min(n,m)kμ(t)⌊ntk⌋⌊mtk⌋=∑t=1min(n,m)kμ(t)⌊⌊n/k⌋t⌋⌊⌊m/k⌋t⌋f(k)=\sum_{t=1}^{\frac{min(n,m)}{k}}\mu(t)\left\lfloor\frac{n}{tk}\right\rfloor\left\lfloor\frac{m}{tk}\right\rfloor=\sum_{t=1}^{\frac{min(n,m)}{k}}\mu(t)\left\lfloor\frac{\lfloor n/k\rfloor}{t}\right\rfloor\left\lfloor\frac{\lfloor m/k \rfloor}{t}\right\rfloorf(k)=t=1∑kmin(n,m)μ(t)⌊tkn⌋⌊tkm⌋=t=1∑kmin(n,m)μ(t)⌊t⌊n/k⌋⌋⌊t⌊m/k⌋⌋

注意，上面的推导还用了下取整函数的性质333。

好了，我们可以线性筛预处理莫比乌斯函数，并求出他的前缀和，然后运用整除分块算法，O(n+m)O(\sqrt n+\sqrt m)O(n+m)计算每一组询问的答案。

参考代码如下：

#include <bits/stdc++.h>
using namespace std;
const int N = 50020;
int a,b,c,d,k,sum[N];
int cnt,Prime[N],p[N],e[N],Pow[N],mu[N];
inline void EularSieve(void)
{p[1] = Pow[1] = 1 , e[1] = 0;for (int i=2;i<=N-20;i++){if ( p[i] == 0 )p[i] = Pow[i] = Prime[++cnt] = i , e[i] = 1;for (int j=1;j<=cnt&&i*Prime[j]<=N-20;j++){int Next = i * Prime[j];p[Next] = Prime[j];if ( p[i] == Prime[j] ){e[Next] = e[i] + 1;Pow[Next] = Pow[i] * Prime[j];break;}else e[Next] = 1 , Pow[Next] = Prime[j];}}mu[1] = 1;for (int i=2;i<=N-20;i++)if ( Pow[i] == i ) mu[i] = e[i] >= 2 ? 0 : -1;else mu[i] = mu[i/Pow[i]] * mu[Pow[i]];// 欧拉筛法+递推求莫比乌斯函数
}
inline void init(void)
{for (int i=1;i<=N-20;i++)sum[i] = sum[i-1] + mu[i];// 计算一下莫比乌斯函数的前缀和
}
inline int calclim(int n,int k,int l) { return k/l ? min( k/(k/l) , n ) : n; }
inline long long Mobius(int n,int m,int k)
{long long res = 0; n /= k , m /= k;for (int l=1,r;l<=min(n,m);l=r+1){r = calclim( min(n,m) , n , l );r = min( r , calclim( min(n,m) , m , l ) );res += 1LL * ( sum[r] - sum[l-1] ) * (n/l) * (m/l);// 整除分块，直接计算答案}return res;
}
int main(void)
{EularSieve() , init();int T; scanf("%d",&T);while ( T --> 0 ){scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);long long ans = Mobius( b , d , k );ans -= Mobius( b , c-1 , k );ans -= Mobius( a-1 , d , k );ans += Mobius( a-1 , c-1 , k );// 容斥一下，求原问题的答案printf("%lld\n",ans);}return 0; // 拜拜程序
}

5.3 BZOJ2820 YY的GCD

有TTT组询问，求

∑i=1n∑j=1m[gcd⁡(i,j)∈Prime]\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)\in Prime]i=1∑nj=1∑m[gcd(i,j)∈Prime]

n,m≤107,T≤10000n,m\leq 10^7,T\leq 10000n,m≤107,T≤10000。

有了上一题的经验，相信大家可以轻松推导到这一步：

ans=∑p∈Prime∑t=1min(n,m)pμ(t)⌊ntp⌋⌊mtp⌋ans=\sum_{p\in Prime}\sum_{t=1}^{\frac{min(n,m)}{p}}\mu(t)\left\lfloor\frac{n}{tp}\right\rfloor\left\lfloor\frac{m}{tp}\right\rfloorans=p∈Prime∑t=1∑pmin(n,m)μ(t)⌊tpn⌋⌊tpm⌋

直接枚举质数ppp求和，由素数定理可知时间复杂度为O(nln⁡n(n+m))O(\frac{n}{\ln n}(\sqrt n+\sqrt m))O(lnnn(n+m))，仍然会超时。

由于ppp也是需要枚举的未知量，此时我们不妨设tp=Ttp=Ttp=T，则原式可以化为：

∑T=1min(n,m)⌊nT⌋⌊mT⌋∑p∈Prime∣Tμ(Tp)\sum_{T=1}^{min(n,m)}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum_{p\in Prime|T}\mu\left(\frac{T}{p}\right)T=1∑min(n,m)⌊Tn⌋⌊Tm⌋p∈Prime∣T∑μ(pT)

再设f(T)=∑p∈Prime∣Tμ(Tp)f(T)=\sum_{p\in Prime|T}\mu\left(\frac{T}{p}\right)f(T)=∑p∈Prime∣Tμ(pT)，那么就有：

ans=∑T=1min(n,m)⌊nT⌋⌊mT⌋f(T)ans=\sum_{T=1}^{min(n,m)}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor f(T)ans=T=1∑min(n,m)⌊Tn⌋⌊Tm⌋f(T)

借鉴上一题的思路，如果我们能够预处理出函数fff的值，并求出前缀和，就也可以整除分块回答询问了。

不幸的是，函数fff和质数有关，不是积性函数，不能用我们之前的方法递推。

但是我们发现nnn以内素数大约只有nln⁡n\frac{n}{\ln n}lnnn个，直接枚举每一个质数去更新它的倍数，均摊每个质数也只会更新ln⁡n\ln nlnn次，所以暴力求值的时间复杂度也不会超时，大约是O(n)O(n)O(n)的。

这个函数的递推其实有严格线性的方法，不过鉴于要推导更复杂的公式，在这道题中完全没有必要，故本文不再赘述。

参考代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 10000020;
int cnt,Prime[N],p[N],e[N],Pow[N],mu[N],f[N],a,b;
inline void EularSieve(void)
{p[1] = Pow[1] = 1 , e[1] = 0;for (int i=2;i<=N-20;i++){if ( p[i] == 0 )p[i] = Pow[i] = Prime[++cnt] = i , e[i] = 1;for (int j=1;j<=cnt&&i*Prime[j]<=N-20;j++){int Next = i * Prime[j];p[Next] = Prime[j];if ( p[i] == Prime[j] ){e[Next] = e[i] + 1;Pow[Next] = Pow[i] * Prime[j];break;}else e[Next] = 1 , Pow[Next] = Prime[j];}}mu[1] = 1;for (int i=2;i<=N-20;i++)if ( Pow[i] == i ) mu[i] = e[i] >= 2 ? 0 : -1;else mu[i] = mu[i/Pow[i]] * mu[Pow[i]];for (int i=1;i<=cnt;i++)for (int j=1;j*Prime[i]<=N-20;j++)f[j*Prime[i]] += mu[j];// 枚举质数暴力更新for (int i=1;i<=N-20;i++) f[i] += f[i-1];
}
inline int calclim(int n,int k,int l) { return k/l ? min( k/(k/l) , n ) : n; }
inline ll Mobius(void)
{ll res = 0;for (int l=1,r;l<=min(a,b);l=r+1){r = calclim( min(a,b) , a , l );r = min( r , calclim( min(a,b) , b , l ) );res += 1LL * ( f[r] - f[l-1] ) * (a/l) * (b/l);}return res;
}
int main(void)
{EularSieve();int T; scanf("%d",&T);while ( T --> 0 )scanf("%d%d",&a,&b),printf("%lld\n",Mobius());return 0;
}

<后记>