莫比乌斯反演学习笔记2

前置知识：莫比乌斯反演学习笔记1

本文介绍一些莫比乌斯反演的灵活应用，包括杜教筛等内容。

一、整除分块（数论分块）

补充一点基础知识...

具体来说就是：如果我们要求(sum limits_{i=1}^nlfloorfrac{n}{i} floor)，并不需要用朴素算法(O(n))去求，可以加快速度。

定理1：所有(1 leq i leq n)的(lfloorfrac{n}{i} floor)的结果最多只有(2sqrt n)种。

证明：对于(i in [1,sqrt n])，显然最多只有(sqrt n)种取值，而对于(i in [sqrt n+1,n])，(lfloorfrac{n}{i} floorleqsqrt n)，因此也最多有(sqrt n)种取值，两个部分相加，因此最多有(2sqrt n)种取值。

定理2：对于确定的(i)，能使得(lfloorfrac{n}{x} floor=lfloorfrac{n}{i} floor)的(x_{max}=Biglfloorcfrac{n}{lfloorfrac{n}{i} floor}Big floor)。

这里选取了两种证明方法。

证法1：记(f(i)=Biglfloorcfrac{n}{lfloorfrac{n}{i} floor}Big floor)，则因为(lfloor frac{n}{i} floor)单调不减，而且(f(i) ge i)，所以必有(lfloorfrac{n}{f(i)} floor leq lfloor frac{n}{i} floor)，又因为(f(i)leqcfrac{n}{lfloorfrac{n}{i} floor})，所以有：

(lfloorfrac{n}{f(i)} floor ge frac{n}{frac{n}{lfloorfrac{n}{i} floor}}=lfloorfrac{n}{i} floor)

同时证明大于等于和小于等于，因此(lfloorfrac{n}{f(i)} floor=lfloorfrac{n}{i} floor)。

另外，若(f(i))不是满足条件的最大值，则(lfloorfrac{n}{f(i)+1} floor=lfloorfrac{n}{i} floor)，则(f(i)=f(f(i)+1))，而又(f(i)ge i)，矛盾，因此(f(i))就是满足要求的最大值。

证法2：记(n=di+p=dx+qquad(0leq p<i,0leq q<x))，则(lfloorfrac{n}{x} floor=lfloorfrac{n}{i} floor=d)，由前式两等式相减得：

(q=p-d(x-i)qquadqquad q=p-dkquad(k=x-i))

要使(x)最大，也就是要让(k)尽量大，由(qge 0)可得(k)的最大值为(k_{max}=lfloorfrac{p}{d} floor)。

此时有：(x=i+d_{max}=i+lfloorfrac{p}{d} floor)，将(p,d)代换得：(x=i+Biglfloorfrac{n mod i}{lfloorfrac{n}{i} floor}Big floor)，将(n mod i)代换，并将(i)乘上去可得：

(x=i+Biglfloorfrac{n-ilfloorfrac{n}{i} floor}{lfloorfrac{n}{i} floor}Big floor=Biglfloorfrac{ilfloorfrac{n}{i} floor+n-ilfloorfrac{n}{i} floor}{lfloorfrac{n}{i} floor}Big floor=Biglfloorcfrac{n}{lfloorfrac{n}{i} floor}Big floor)

同样得到了结论。

因此我们只要通过上面对商的讨论，就可以在(O(sqrt n))时间复杂度内快速计算(sum limits_{i=1}^nlfloorfrac{n}{i} floor)。

for (int l=1,r;l<=n;l=r+1) {
	r=n/(n/l);
	ans+=(n/l)*(r-l+1);
}

来看一道题目。

洛谷2261 题目大意：给定(n,k leq 10^9)，求出(sum limits_{i=1}^nk mod i=k mod 1+...+k mod n)。

分析：直接(O(n))做是肯定不行的，这时我们想到利用余数的定义：(n mod i=n-ilfloorfrac{n}{i} floor)，于是我们就可以开心化简了：

(sum limits_{i=1}^n k mod i=kn-sumlimits_{i=1}^n ilfloorfrac{k}{i} floor)

而这个式子的后面一部分非常像整除分块的形式，但是还不完全一样。不过这也不难，我们只要将上面的代码中原本表示(l-r)之间数的个数的((r-l+1))替换成(l-r)的和(frac{(l+r)(r-l+1)}{2})即可。

有一点需要注意，此时求和上限和被除数不同，需要加上一句额外判断，当(Biglfloorcfrac{k}{lfloorfrac{k}{i} floor}Big floor)超过了(n)时也只能取到(n)。

#include <bits/stdc++.h>
using namespace std;
long long ans,n,k;
int main () {
	scanf("%lld%lld",&n,&k);
	ans=n*k;
	for (long long l=1,r;l<=min(n,k);l=r+1) {
		r=min(k/(k/l),n);
		ans-=(k/l)*(l+r)*(r-l+1)/2;
	}
	printf("%lld
",ans);
	return 0;
}

通过上面这道题目我们能看出来什么？其实整除分块并不只能求原本的标准形式。事实上，形如整除式后面乘上一个函数的值这样的求和式，都可以用整除分块的方法做出，即：

(sum limits_{i=1}^n lfloorfrac{n}{i} floor f(i)=sum lfloorfrac{n}{l} floor (s(r)-s(l-1)))

其中(s(k)=sum limits_{i=1}^kf(i))，即函数的前缀和，(l,r)与上面代码中的意义相同。其实上面的例子也是这样，只是(s(i))恰好是一个等差数列的形式，很容易算出。

二、gcd快速求和

由于这个部分比较好玩，而且与反演和卷积也有点关系，体现一些基本思想，就先讲一些了。

洛谷2303 一维gcd求和：求(sumlimits_{i=1}^Ngcd(N,i))。

分析：我们对不同的数据范围进行逐层递进的分析：

（1）(N leq 10^5)

这个非常简单，只要你会欧几里得算法，就可以直接计算(gcd(N,i))的值，然后将所有结果相加，暴力计算即可得到答案，时间复杂度(O(Nlog N))；

（2）(N leq 10^9)

主要的思维都体现在这里了......首先处理这种最大公因数问题的最常用方法之一，是枚举最大公因数的值(d)，也就是我们枚举每一个可能的(d)，然后去计数有多少个为(d)的结果，将(d)与计数量相乘即可得到结果为(d)的最大公因数和了，而(d=gcd(N,i))，则显然(d|N)，因此可以这样枚举：

(sum limits_{d|N} d sumlimits_{i=1}^N[gcd(N,i)=d])

这是第一步技巧，将gcd提取出来放在外面。

接下来还要继续化简，我们现在感觉这个方括号里的等式真的没什么用，那么遇到这种情况有两种方法，一种是用后面要讲的莫比乌斯反演，但这里还不必要；另一种就是“约分”，也就是将内侧求和的每一个d都约去，这样的话(N)就变成了(frac{N}{d})了：

(sum limits_{d|N} d sumlimits_{i=1}^frac{N}{d}[gcd(frac{N}{d},i)=1])

这时候你就应该意识到了：此时里面这个求和号其实就是欧拉函数！我们将它表示一下：

(sum limits_{d|N} d varphi(frac{N}{d}))

这样做(10^9)数据就足够了，我们枚举N的因子d，并算出N/d的欧拉函数值就可以了。

（3）(N leq 10^{12})

这时上面计算欧拉函数的时间代价就太高了，因此需要再次优化。

如果你掌握之前讲的卷积足够好的话，你就应该及时意识到，上面的这个式子是一个卷积式，其实是(id)与(varphi)的卷积。之前我们还学过：两个积性函数的卷积依然是积性函数，所以此时我们假设(f(N)=sum limits_{d|N} d varphi(frac{N}{d}))，此时(f)就是一个积性函数。那么我们将(n)分解素因数为(N=p_1^{c_1}...p_n^{c_n})，则：

(f(N)=f(p_1^{c_1})...f(p_n^{c_n}))

因此我们只需要算出对于一个素数的k次幂(p^k)，对应的函数值是多少即可。

直接用定义将(f(p^k))带入得：(f(p^k)=sum limits_{d|p^k} d varphi(frac{p^k}{d}))，而(p^k)的所有因数其实就是(p^0)一直到(p^k)，因此枚举d实质上就是枚举p的次数，于是原式可以化简为：

(f(p^k)=sum limits_{i=0}^k p^i varphi(frac{p^k}{p^i})=sum limits_{i=0}^k p^i varphi(p^{k-i}))

下面考验基本功了，有一个欧拉函数的结论：(varphi(p^x)=(p-1) imes p^{x-1}=p^x-p^{x-1})，这属于基础知识，就不证明了。

但是要注意，我们上面的式子中(x)是不为0的，在(f(p^k))的表达式中(k-i)却可能为0，因此(i=k)的情况单独讨论：

(f(p^k)=p^k varphi(1)+sum limits_{i=0}^{k-1} p^i varphi(p^{k-i})=p^k+sum limits_{i=0}^{k-1} p^i imes(p^{k-i}-p^{k-i-1}))

到这里思路已经很明确了，直接拆去括号利用乘法分配律计算：

(f(p^k)=p^k+sum limits_{i=0}^{k-1} p^k-p^{k-1}=(k+1)p^k-kp^{k-1}=(pk+p-k) imes p^{k-1})

最后这个式子非常简单，此时你已经可以直接对(N)分解素因数，对于每一个素因数用上面这行公式快速幂计算了，但是，也许用不到快速幂那多出来的(log)呢？

(f(N)=prodlimits_{i=1}^kf(p_i^{c_i})=prodlimits_{i=1}^k(p_ic_i+p_i-c_i) imesprodlimits_{i=1}^kp_i^{c_i-1})

碰巧后面这东西和N有关，所以进一步化简，最终得到很简单的一个式子：

(f(N)=N imesprodlimits_{i=1}^kfrac{p_ic_i+p_i-c_i}{p_i})

这样就可以直接分解做了，最终时间复杂度(O(sqrt{n}))。

#include <bits/stdc++.h>
#define ll long long
using namespace std;
ll n,cnt,ans;
ll g (ll p,ll k) {
    return (ans/p)*(p*k+p-k);
}
int main () {
    scanf("%lld",&n);
    ans=n;
    for (ll i=2;i*i<=n;i++) {
        cnt=0;
        while (n%i==0) {
            cnt++,n/=i;
        }
        if (cnt) {
    	    ans=g(i,cnt);
		}
    }
    if (n>1) {
    	ans=g(n,1);
	}
    printf("%lld
",ans);
    return 0;
}

一维gcd求和基本就到此为止了，可是下面才是真正的挑战...

二维gcd求和：求(sumlimits_{i=1}^Nsumlimits_{j=1}^Mgcd(i,j))。

我们将讨论这个问题的多种解法（当然暴力计算也算一种...），可以说是各有优劣，有些只适用于特殊的情况，有些可以处理多组数据等等。

（1）(N=M)，线性时间

这种情况还是比较简单的。我们只需要对这个式子稍作变形就能够直接转变为一维的形式，而我们发现gcd是内的两个数是可以交换的，所以我们只枚举(jleq i)的情况，就可以转化了：

(sumlimits_{i=1}^Nsumlimits_{j=1}^Ngcd(i,j)=2sumlimits_{i=1}^Nsumlimits_{j=1}^igcd(i,j)-sumlimits_{i=1}^Ngcd(i,i))

请注意这里因为(i=j)的情况这里被计算了两次，因此要减去一次。用上面的表示，若(f(N)=sumlimits_{i=1}^Ngcd(N,i))，则化简为：

(sumlimits_{i=1}^Nsumlimits_{j=1}^Ngcd(i,j)=2sumlimits_{i=1}^Nf(i)-frac{N imes(N+1)}{2})

我们上面已经发现，(f(N))是积性函数，因此直接做线性筛求前缀和即可。这种方法时间复杂度为(O(n))，优势是可以直接处理多组询问。

（2）无特殊限制，不妨设(N leq M)

这时上面的式子就没法推了，这时我们运用于一维gcd求和类似的方法，提取出gcd并枚举，同时利用类似“约分”的性质约去d：

(sumlimits_{i=1}^Nsumlimits_{j=1}^Mgcd(i,j)=sumlimits_{d=1}^Ndsumlimits_{i=1}^Nsumlimits_{j=1}^M[gcd(i,j)=d]=sumlimits_{d=1}^Ndsumlimits_{i=1}^{frac{N}{d}}sumlimits_{j=1}^{frac{M}{d}}[gcd(i,j)=1])

之前我们直接用欧拉函数解决后面一部分，这时就不行了。我们使用第二个技巧：莫比乌斯反演。

我们令(f(x)=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)=x])，与之对应，我们需要一个反演的辅助函数，而形如这种方括号里带等式的式子，我们通常可以变相等为整除构造，即：

(g(x)=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[x|gcd(i,j)]=sumlimits_{x|d}f(d))

而由(x|gcd(i,j))得(x|i,x|j)，因此其实(g(x)=lfloorfrac{n}{x} floorlfloorfrac{m}{x} floor)。因此利用反演，我们得到：

(f(x)=sumlimits_{x|d}mu(frac{d}{x})lfloorfrac{n}{d} floorlfloorfrac{m}{d} floor)

原题中，我们只要求(f(1))，那我们来推一推：(f(1)=sumlimits_{i=1}^nmu(i)lfloorfrac{n}{i} floorlfloorfrac{m}{i} floor)。那么我们回代原式：

(sumlimits_{i=1}^Nsumlimits_{j=1}^Mgcd(i,j)=sumlimits_{d=1}^Ndsumlimits_{i=1}^{frac{N}{d}}mu(i)iglfloorfrac{lfloorfrac{N}{d} floor}{i}ig flooriglfloorfrac{lfloorfrac{M}{d} floor}{i}ig floor)

推到这里算是第一个层次，接下来看到这个式子就至少会有一个初步的想法了。

这个式子的内层求和其实是整除分块的形式，如果没看出来的话建议再看一下第一部分...这里是一个二元的整除分块，但是形式是一样的，时间复杂度变为了两元的根号和而已，依然存在若干对两个被除数来说商都一样的区间。那么，里面我们就可以进行整除分块了。那么外面怎么办呢？事实上，我们发现对于里面算出的分块结果，只与(lfloorfrac{N}{d} floor)以及(lfloorfrac{M}{d} floor)有关，从函数的角度出发，这还是一个整除分块，即：(sumlimits_{i=1}^nf(lfloorfrac{n}{i} floor))，这里还是增为二元，本质不变，因此外层是另一个整除分块，所以我们至少得到了一种用两次整除分块得到结果的方法，时间复杂度其实是两次分块的积：(O((sqrt n+sqrt m)^2)=O(n + m))。能做到线性。

但是这还不够好，我们需要低于线性复杂度的做法。更好的算法将在杜教筛中给出。

三、快速求和举例

有时我们遇见的不是gcd的求和，而是其余一些函数的求和式，通常的思路有转化成gcd求和或者根据公式展开，然后交换枚举顺序，或者进行莫比乌斯反演等。

  1. 洛谷3327 计算(sumlimits_{i=1}^nsumlimits_{j=1}^md(ij))，其中d为因数个数。

    这里采用的是直接将d的求和转换成gcd求和，首先证明引理：(d(ij)=sumlimits_{x|i}sumlimits_{y|j}[gcd(x,y)=1])。我比较喜欢上面链接中一个题解的证明方法：构造一一对应。

    对于ij的一个因数k，考虑k的某个素因数p，假设k中p的次数为a，i中p的次数为b，j中p的次数为c，则：

    （1）对于(aleq b)，取(x=p^a,y=1)，显然(gcd(x,y)=1)

    （2）对于(bleq aleq c)，取(x=1,y=p^{b-a})，显然(gcd(x,y)=1)

    而上面两种情况都对应了一个(x|i,y|j)的数对((x,y))，对于k的不同素因数分别独立计算出来的(gcd(x,y)=1)当然也是成立的。因此我们建立了等式两边的一一对应，即证明了两边相等。（不过这不是我们的重点）

    于是我们要求的式子变成了：

    (sumlimits_{i=1}^nsumlimits_{j=1}^msumlimits_{x|i}sumlimits_{y|j}[gcd(x,y)=1])

    这时看上去没有什么好办法了，那么在这种因数枚举在内层的求和式，我们就可以尝试转换一下枚举顺序：

    (sumlimits_{x=1}^nsumlimits_{y=1}^mlfloorfrac{n}{x} floorlfloorfrac{m}{y} floor[gcd(x,y)=1])

    明显求和号减少了，于是我们看到方括号里等号的形式，想到要进行莫比乌斯反演化简，下面我们进行反演。前面讲过，遇见方括号里为等式的，转换成整除通常就可以用一个简单的式子表示出来，令(f(k)=sumlimits_{x=1}^nsumlimits_{y=1}^mlfloorfrac{n}{x} floorlfloorfrac{m}{y} floor[gcd(x,y)=k])，则我们可以这样定义：

    (g(k)=sumlimits_{x=1}^nsumlimits_{y=1}^mlfloorfrac{n}{x} floorlfloorfrac{m}{y} floor[k|gcd(x,y)]=sumlimits_{k|d}f(d))

    下面重点就是计算出(g(k))的值了，看上去不是很容易，但是联想到我们在做gcd求和的时候可以约去枚举的因数d，在这里其实也可以把k约掉（有点神奇的一步）：

    (g(k)=sumlimits_{x=1}^{lfloorfrac{n}{k} floor}sumlimits_{y=1}^{lfloorfrac{m}{k} floor}lfloorfrac{n}{xk} floorlfloorfrac{m}{yk} floor)

    记(p=lfloorfrac{n}{k} floor,q=lfloorfrac{m}{k} floor)，则可化为：

    (g(k)=sumlimits_{x=1}^plfloorfrac{p}{x} floorsumlimits_{y=1}^qlfloorfrac{q}{y} floor)

    到这里已经非常明确了：整除分块。我们可以预处理所有这样的求和式，这样即可(O(1))计算g的值。如果只有一组数据，也可以遇到了再算。回到题目，题中要求的实际上就是(f(1))，所以我们带入即可：

    (f(1)=sumlimits_{i=1}^nmu(i)(sumlimits_{x=1}^{lfloorfrac{n}{i} floor}lfloorfrac{lfloorfrac{n}{i} floor}{x} floorsumlimits_{y=1}^{lfloorfrac{m}{i} floor}lfloorfrac{lfloorfrac{m}{i} floor}{y} floor))

    看似是一个(O(n))的枚举算法，但是实际上后面这堆明显是个两个整除分块，所以按照上面说的预处理的商求和式的值就可以做到(O(sqrt n))的复杂度。但是这种算法需要初始化，因此只适用于多组数据，而不适用n真的很大的情况。最终时间复杂度(O(nsqrt n+Tsqrt n))，T为数据组数。（懒得线性筛莫比乌斯函数了，直接暴力计算，反正不会影响最终复杂度）

    #include <bits/stdc++.h>
#define ll long long
using namespace std;
const int MAXN=50010;
ll k,n,m,p[MAXN],s[MAXN],ans;
ll calcp (int x) {
	ll ans=0;
	for (int l=1,r;l<=x;l=r+1) {
		r=(x/(x/l));
		ans+=(x/l)*(r-l+1);
	}
	return ans;
}
ll calcm (int x) {
	ll mu=1;
	for (int i=2;i*i<=x;i++) {
		if (x%i==0) {
			mu*=-1,x/=i;
			if (x%i==0) {
				return 0;
			}
		}
	}
	if (x>1) {
		mu*=-1;
	}
	return mu;
}
int main () {
	for (int i=1;i<=50000;i++) {
		p[i]=calcp(i);
		s[i]=s[i-1]+calcm(i);
	}
	scanf("%lld",&k);
	for (int ii=1;ii<=k;ii++) {
		scanf("%lld%lld",&n,&m);
		ans=0;
		for (int l=1,r;l<=min(n,m);l=r+1) {
			r=min(n/(n/l),m/(m/l));
			ans+=(s[r]-s[l-1])*p[n/l]*p[m/l];
		}
		printf("%lld
",ans);
	}
	return 0;
}

  2. 洛谷2257 计算(sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)=p])，(p)为任意素数。

    首先我们按照惯用套路，枚举gcd的值并约去p：

    (sumlimits_{pin prim}sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)=p]=sumlimits_{pin prim}sumlimits_{i=1}^{lfloorfrac{n}{p} floor}sumlimits_{j=1}^{lfloorfrac{m}{p} floor}[gcd(i,j)=1])

    于是最后部分直接莫比乌斯反演（参见gcd求和部分）：

    (sumlimits_{pin prim}sumlimits_{d=1}^{lfloorfrac{n}{p} floor}mu(d)lfloorfrac{n}{pd} floorlfloorfrac{m}{pd} floor)

    于是呢？直接开始枚举素数吗？然而(nleq 10^7)会让你T爆，所以需要继续化简。比较神奇的一步是：枚举pd=t，然后直接转换整个结构。这里用到这个技巧主要是因为注意到最后的整除部分只和pd有关，因此转换后可以转化为一个整除分块形式，这也是常用技巧：

    (sumlimits_{t=1}^nlfloorfrac{n}{t} floorlfloorfrac{m}{t} floorsumlimits_{p|t,pin prim}mu(frac{t}{p}))

    于是这道题到这里就告一段落了，从整体来看，这是一个典型的整除分块，内层求和看，这是一个线性筛的题目，我们可以直接筛出内层求和式子。

    （代码未经卡常，O2过）

    #include <bits/stdc++.h>
#define ll long long
using namespace std;
const int MAXN=10000010;
int t,tot,mu[MAXN],q[MAXN],p[MAXN];
ll f[MAXN],s[MAXN],n,m;
void init () {
	mu[1]=1;
	for (int i=2;i<=MAXN-10;i++) {
		if (!q[i]) {
			p[++tot]=i;
			mu[i]=-1;
		}
		for (int j=1;j<=tot&&i*p[j]<=MAXN-10;j++) {
			q[i*p[j]]=1;
			if (i%p[j]==0) {
				mu[i*p[j]]=0;
				break;
			}
			mu[i*p[j]]=-mu[i];
		}
	}
	for (int i=1;i<=tot;i++) {
		for (int j=p[i];j<=MAXN-10;j+=p[i]) {
			f[j]+=mu[j/p[i]];
		}
	}
	for (int i=1;i<=MAXN-10;i++) {
		s[i]=s[i-1]+f[i];
	}
}
int main () {
	init();
	scanf("%d",&t);
	for (int ii=1;ii<=t;ii++) {
		scanf("%lld%lld",&n,&m);
		ll ans=0;
		for (int l=1,r;l<=min(n,m);l=r+1) {
			r=min(n/(n/l),m/(m/l));
			ans+=(n/l)*(m/l)*(s[r]-s[l-1]);
		}
		printf("%lld
",ans);
	}
	return 0;
}

  3. 洛谷1829 计算(sumlimits_{i=1}^nsumlimits_{j=1}^mlcm(i,j))。

    这道题有点毒瘤，来看看如何化简。

    首先lcm直接做肯定不好做，要转化成gcd的问题方便求解，转换也很简单，我们知道(ij=gcd(i,j)lcm(i,j))，因此自然而然的将(lcm(i,j))用(frac{ij}{gcd(i,j)})代换即可。

    然后是惯用的套路，我们将gcd进行枚举，注意这次gcd在分母上，所以枚举出来也在分母：

    (sumlimits_{d=1}^nfrac{1}{d}sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)=d]ij)

    这时来看里面的这个式子，我们将d约去：

    (sumlimits_{d=1}^nfrac{1}{d}sumlimits_{i=1}^{lfloorfrac{n}{d} floor}sumlimits_{j=1}^{lfloorfrac{m}{d} floor}[gcd(i,j)=1]id imes jd)

    注意里面的ij要乘上d，保留原本的结果。做到这里，我们发现里面的两个d是可以直接拿出来的，于是分母上的d就直接被约分了：

    (sumlimits_{d=1}^ndsumlimits_{i=1}^{lfloorfrac{n}{d} floor}sumlimits_{j=1}^{lfloorfrac{m}{d} floor}[gcd(i,j)=1]ij)

    我们设(f(n,m)=sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)=1]ij)，着手对它进行化简。对于这样的式子，我们已经有了很多经验，可以直接把([gcd(i,j)=1])转换成莫比乌斯函数的求和，具体反演步骤请看上文gcd求和。

    (f(n,m)=sumlimits_{d=1}^nmu(d)sumlimits_{d|i}sumlimits_{d|j}ij)

    做到这里遇到了一个小瓶颈，这种时候通常要么是转换枚举顺序，要么就要转换枚举对象，这里我们选择后者。设(p=frac{i}{d},q=frac{j}{d})。枚举p和q，直接将ij代换就得到了以下的式子：

    (f(n,m)=sumlimits_{d=1}^nmu(d)sumlimits_{p=1}^{lfloorfrac{n}{d} floor}sumlimits_{q=1}^{lfloorfrac{m}{d} floor}pd imes qd)

    这个形式很好，我们直接将(d^2)放到外层求和，内层的式子也变得十分简单了：

    (f(n,m)=sumlimits_{d=1}^nd^2mu(d)sumlimits_{p=1}^{lfloorfrac{n}{d} floor}sumlimits_{q=1}^{lfloorfrac{m}{d} floor}pq)

    这时里面的式子直接用等差数列求和公式代换掉：

    (f(n,m)=sumlimits_{d=1}^nd^2 imes mu(d) imesfrac{lfloorfrac{n}{d} floor imes(lfloorfrac{n}{d} floor+1)}{2} imesfrac{lfloorfrac{m}{d} floor imes(lfloorfrac{m}{d} floor+1)}{2})

    式子当中的(f(n,m))可以用整除分块的方式做，(d^2mu(d))的前缀和是完全可以线性筛的。而我们再来看看整个求和式子的值，将(f(n,m))放进式子里其实就是：

    (sumlimits_{d=1}^nd imes f(lfloorfrac{n}{d} floor,lfloorfrac{m}{d} floor))

    这是一个很好的形式，因为外层还是一个整除分块...于是我们用分块套分块的方式就完成了这道题目。时间复杂度为(O((sqrt n+sqrt m)^2)=O(n+m))

    #include <bits/stdc++.h>
#define ll long long
using namespace std;
const int MAXN=10000010,P=20101009;
ll n,m,ans,tot,s[MAXN];
int mu[MAXN],q[MAXN],p[MAXN];
void init () {
	mu[1]=1;
	for (int i=2;i<=MAXN-10;i++) {
		if (!q[i]) {
			p[++tot]=i,mu[i]=-1;
		}
		for (int j=1;j<=tot&&i*p[j]<=MAXN-10;j++) {
			if (i%p[j]==0) {
				mu[i*p[j]]=0,q[i*p[j]]=1;
				break;
			}
			mu[i*p[j]]=mu[i]*(-1),q[i*p[j]]=1;
		}
	}
	for (int i=1;i<=MAXN-10;i++) {
		s[i]=(s[i-1]+1LL*mu[i]*i*i)%P;
	}
	return;
}
ll gt (ll n,ll m) {
	return ((((n*n+n)/2)%P)*(((m*m+m)/2)%P))%P;
}
ll calc (ll n,ll m) {
	ll ans=0;
	for (ll l=1,r;l<=min(n,m);l=r+1) {
		r=min(n/(n/l),m/(m/l));
		ans=(ans+((gt(n/l,m/l)%P)*((s[r]-s[l-1]+P)%P))%P)%P;
	}
	return ans;
}
int main () {
	init();
	scanf("%d%d",&n,&m);
	for (ll l=1,r;l<=min(n,m);l=r+1) {
		r=min(n/(n/l),m/(m/l));
		ans=(ans+((((l+r)*(r-l+1))/2)%P*calc(n/l,m/l)%P)%P)%P;
	}
	printf("%d
",ans);
	return 0;
}

  4. 洛谷3768 计算(sumlimits_{i=1}^nsumlimits_{j=1}^mijgcd(i,j))。

    本题和上题基本相同，只是一开始d出现在分子的位置，因此我们直接写出最后的计算式：

    (sumlimits_{d=1}^nd^3 imes f(lfloorfrac{n}{d} floor,lfloorfrac{m}{d} floor))

    接下来是小学奥数部分，你也许学过(sumlimits_{i=1}^ni^3=(sumlimits_{i=1}^ni)^2)。即使没学过，也很容易用数学归纳法证明：

    n=1时显然成立，n>1时假定n-1的情况成立，则等式左侧的增量为(n^3)，右侧的增量为：

    (n^2+2n(sumlimits_{i=1}^{n-1}i)=n^2+n^2(n-1)=n^3)

    两侧增量相等，因此等式成立。

    于是我们就可以解决(d^3)的前缀和了。代码和lcm求和非常相近，就不发了。（当然上一道题和这道题数据范围完全不同，本题做到这里只能解决部分数据，满分需要用到杜教筛处理(f(i,j))解决）

总结：

通过几道题目的练习，我们掌握了一些处理快速求和式的基本方法：

第一步：一般都是转换成我们较熟悉gcd求和的问题，这里就要求掌握一些转换的公式法则；

第二步：遇见了gcd求和式，我们一般把gcd提取出来进行枚举，然后将内层求和转换成方括号形式；

第三步：如果需要，在内层求和约去d，将n和m分别除掉d，这样可以使方括号内变成等于1的形式；

第四步：讨论内层方括号的求值，通常是整除分块类型，有等差数列/线性筛/继续化简几种情况；

第五步：带回整个式子进行计算，通常也是线性筛+整除分块类型。

另外，要熟练掌握这个式子，这可以帮助你快速化简方括号内的值：

(sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)=x]=sumlimits_{x|d}mu(frac{d}{x})lfloorfrac{n}{d} floorlfloorfrac{m}{d} floor)

通常在x=1的时候有特例：

(sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)=1]=sumlimits_{d=1}^nmu(i)lfloorfrac{n}{d} floorlfloorfrac{m}{d} floor)

另外也要熟练掌握整除分块以及和莫比乌斯函数有关的能够线性筛出前缀和的函数，也就是要掌握下面这种形式的式子如何在线性初始化后由根号的时间复杂度内求出：

(Ans=sumlimits_{d=1}^ng(d)f(lfloorfrac{n}{d} floor,lfloorfrac{m}{d} floor))

其中g(d)是与莫比乌斯函数有关的某个可线性筛函数...

最后，如果我们要求和的式子有左边界，也就是形如(aleq ileq c,bleq jleq d)，那么我们简单容斥一下，其实就是一个二维前缀和：

(sumlimits_{i=a}^csumlimits_{j=b}^df(i,j)=sumlimits_{i=1}^csumlimits_{j=1}^df(i,j)+sumlimits_{i=1}^{a-1}sumlimits_{j=1}^{b-1}f(i,j)-sumlimits_{i=1}^{a-1}sumlimits_{j=1}^df(i,j)-sumlimits_{i=1}^csumlimits_{j=1}^{b-1}f(i,j))

四、杜教筛

杜教筛是一类用来算积性函数前缀和的工具，可以做到低于线性的时间复杂度。通常作为线性筛的扩展。也就是说，我们现在要求：

(s(i)=sumlimits_{i=1}^nf(i), nleq10^9)

这样一来线筛就被卡死了......但是我们可以用杜教筛来求这个东西。首先在此之前，我们依然需要用线性筛预处理出一部分的$ s(i)$值，例如我们这里做到(10^6)，这是可以用线性筛解决的。

接下来，就是考验脑洞的时候了。我们找来一个函数(g)，使得(g)的前缀和与(f*g)的前缀和都是容易求的（其实杜教筛的对象函数通常比较固定，所以记一些选取辅助函数的模板也是可以的），设(h=f*g)，然后我们来看看如何用f和g来表示h的前缀和：

(sumlimits_{i=1}^nh(i)=sumlimits_{i=1}^nsumlimits_{d|i}g(d)f(frac{i}{d})=sumlimits_{d=1}^ng(d)sumlimits_{d|i}f(frac{i}{d}))

第二步直接转换了枚举顺序。接下来我们将后面一部分化简，设(p=frac{i}{d})，则：

(sumlimits_{i=1}^nh(i)=sumlimits_{d=1}^ng(d)sumlimits_{p=1}^{lfloorfrac{n}{d} floor}f(p)=sumlimits_{d=1}^ng(d)s(lfloorfrac{n}{d} floor))

看到这里你可能觉得没什么用...但是实际上我们要算的(s(n))就包含在了右边这个式子里面，我们只要将右边式子第一项分离出来：

(sumlimits_{i=1}^nh(i)=g(1)s(n)+sumlimits_{d=2}^ng(d)s(lfloorfrac{n}{d} floor))

于是(s(n))就在那里了，我们只要稍微移个项就可以得到(s(n))的式子了，这就是杜教筛的套路式：

(g(1)s(n)=sumlimits_{i=1}^nh(i)-sumlimits_{d=2}^ng(d)s(lfloorfrac{n}{d} floor))

前一部分使我们说好的可以快速算前缀和，后一部分实则是个递归，但是整除分块的形式，可以帮助我们快速计算出它的值。至于等式左边......只能说我们一般选取的g的(g(1))值都为1吧......

根据各种我不会的神奇的数学知识可以得到计算(s(n))的一般复杂度为(O(n^{frac{3}{4}}))，但是由于我们可以用线性筛筛出前k的(s(i))，因此复杂度变成了(O(k+frac{n}{sqrt k}))，取(k=n^{frac{2}{3}})时可以得到最优时间复杂度(O(n^{frac{2}{3}}))。

模板题：洛谷4213：杜教筛(Sum)

  1. 求(s(n)=sumlimits_{i=1}^nmu(i))

    根据套路式，我们需要找到一个合适的函数使得它与(mu)的卷积的前缀和很容易求，那么在将卷积的时候我们已经提到过一组逆元的关系：(mu*I=varepsilon)，因此我们只要令g=I即可，然后套进套路式里，(g(i)=1)作为因数就不写了，而(varepsilon)的前缀和无论如何都是1...因此作为被减数的一部分直接就是1了：

    (s(n)=1-sumlimits_{d=2}^ns(lfloorfrac{n}{d} floor))

    于是算法步骤就很简单，先线性筛出(10^6)以下的(mu)值，然后用杜教筛递归计算上面的式子即可。

  2. 求(s(n)=sumlimits_{i=1}^nvarphi(i))

    还记得我们在上一次讲到了非常有用的一个等式：(mu*id=varphi)，因此我们对此进行反演，即可得(varphi*I=id)，这是少数利用带莫比乌斯的式子转换成带I式子的反演之一。于是思路很清晰了，我们取g=I，带进套路式中，只要求出id的前缀和即可，而id的前缀和实际上就是等差数列。

    (s(n)=frac{n(n+1)}{2}-sumlimits_{d=2}^ns(lfloorfrac{n}{d} floor))

    先线性筛出(10^6)以下的(varphi)值，然后用杜教筛递归计算上面的式子即可。

    （由于题目中每个点有十组数据，因此线性筛的范围可以适当增大，下面代码中开到了(2 imes 10^6)）

    一般的杜教筛选用的是map来做记忆化，但是我对某些stl容器有心理障碍，所以开了一个数组，如果有一些其他的应用可能还是得用map才行。

#include <bits/stdc++.h>
#define ll long long 
using namespace std;
const int MAXN=2000010,DSQRT=50010;
ll sphi[MAXN],mem[DSQRT],mep[DSQRT];
int t,n,tot,mu[MAXN],phi[MAXN],q[MAXN],smu[MAXN],p[MAXN];
void init () {
	mu[1]=phi[1]=smu[1]=sphi[1]=1;
	for (int i=2;i<=MAXN-10;i++) {
		if (!q[i]) {
			p[++tot]=i,phi[i]=i-1,mu[i]=-1;
		}
		for (int j=1;j<=tot&&i*p[j]<=MAXN-10;j++) {
			if (i%p[j]==0) {
				q[i*p[j]]=1,phi[i*p[j]]=phi[i]*p[j],mu[i*p[j]]=0;
				break;
			}
			q[i*p[j]]=1,phi[i*p[j]]=phi[i]*(p[j]-1),mu[i*p[j]]=-mu[i];
		}
		smu[i]=smu[i-1]+mu[i],sphi[i]=sphi[i-1]+phi[i];
	}
	return;
}
ll sp (int x) {
	if (x<=MAXN-10) {
		return sphi[x];
	}
	if (mep[n/x]) {
		return mep[n/x];
	}
	ll ans=x*(ll)(x+1)/2;
	for (int l=2,r;l<=x;l=r+1) {
		r=x/(x/l);
		ans-=sp(x/l)*(r-l+1);
	}
	return mep[n/x]=ans;
}
ll sm (int x) {
	if (x<=MAXN-10) {
		return smu[x];
	}
	if (mem[n/x]) {
		return mem[n/x];
	}
	ll ans=1;
	for (int l=2,r;l<=x;l=r+1) {
		r=x/(x/l);
		ans-=sm(x/l)*(r-l+1);
	}
	return mem[n/x]=ans;
}
int main () {
	init();
	scanf("%d",&t);
	for (int i=1;i<=t;i++) {
		memset(mem,0,sizeof(mem));
		memset(mep,0,sizeof(mep));
		scanf("%d",&n);
		printf("%lld %lld
",sp(n),sm(n));
	}
	return 0;
}

回过来看这题：(sumlimits_{i=1}^nsumlimits_{j=1}^mijgcd(i,j))，我们已经化简到了：

(sumlimits_{x=1}^nx^3sumlimits_{d=1}^{lfloorfrac{n}{x} floor}d^2 imes mu(d) imesfrac{lfloorfrac{n}{xd} floor imes(lfloorfrac{n}{xd} floor+1)}{2} imesfrac{lfloorfrac{m}{xd} floor imes(lfloorfrac{m}{xd} floor+1)}{2})

里面这个分数太烦，我们记(g(i,j)=frac{i imes(i+1)}{2} imesfrac{j imes(j+1)}{2})，然后我们发现里面这东西和(xd)有关，所以我们可以设(t=xd)，转换一下枚举顺序，枚举t的值，这时我们发现，外层的(x^3)可分两个给(d^2)，正好凑成(t^2)，剩下一个x放在内层：

(sumlimits_{t=1}^nt^2g(lfloorfrac{n}{t} floor,lfloorfrac{m}{t} floor)sumlimits_{x|t}x imes mu(frac{t}{x}))

这里又出现了一个很基础的式子：(sumlimits_{x|t}x imes mu(frac{t}{x}))，这东西是一个卷积式，而且非常喜闻乐见：(mu*id)，那不多说了，直接化成(varphi(t))即可。于是x就消失了：

(sumlimits_{t=1}^nvarphi(t) imes t^2 imes g(lfloorfrac{n}{t} floor,lfloorfrac{m}{t} floor))

最后就这么个东西，我们要做到(nleq 10^9)，就必须解决(x^2varphi(x))的前缀和，而且要使用杜教筛解决。下面我们的主要目标就是这个东西：

(sumlimits_{i=1}^nvarphi(i) imes i^2)

带个(i^2)始终很不方便，我们想要消掉它，于是下面介绍杜教筛中的一个技巧：如果要求和的函数出现了一个单独的数，我们可以配一个(frac{n}{i})之类的式子，去把它消掉。我们先看看，随意找一个(g)，得到的卷积形式h是：

(h(n)=sumlimits_{d|n}varphi(d) imes d^2 imes g(frac{n}{d}))

想让d消失掉？那么我们只需要让分母上出现(d^2)即可，于是我们可以设(g(x)=x^2)，这样(h(n))就可以进行化简了，只剩下一个关于(varphi)的式子：

(h(n)=sumlimits_{d|n}varphi(d) imes d^2 imes frac{n^2}{d^2}=n^2sumlimits_{d|n}varphi(d))

好了，要的就是这个效果，我们看到了(sumlimits_{d|n}varphi(d))，在筛欧拉函数的时候我们已经讲过了，(varphi*I=id)，那么这个求和式实际上就是n：

(h(n)=n^2 imes n=n^3)

因此原本的求和就可以用下面这个式子表示了：

(s(n)=sumlimits_{i=1}^ni^3-sumlimits_{i=2}^ni^2s(lfloorfrac{n}{i} floor))

而我们下面就是要求(1^3+...+n^3)以及(1^2+...+n^2)（不写sigma更接地气，毕竟小学题目），前面那个我们已经证明过了，等于((1+...+n)^2)，后面这个其实是(frac{n(n+1)(2n+1)}{6})，我们来证明一下：

(n=1)时显然成立；

(n>1)时若对于n-1成立，则等式左侧的增量为(n^2)，右侧增量为：

(frac{n(n+1)(2n+1)}{6}-frac{(n-1)n(2n-1)}{6}=frac{n}{6}((n+1)(2n+1)-(n-1)(2n-1)))
(=frac{n}{6}(6n)=n^2)

两侧增量相等，因此等式对于n也成立。

于是本题至此可以完全解决，时间复杂度为(O(n^frac{2}{3}))。

最后，我们看看之前还没完全解决的一个问题：

(sumlimits_{i=1}^Nsumlimits_{j=1}^Mgcd(i,j))

我们之前通过反演已经得到了这个结果：

(sumlimits_{d=1}^Ndsumlimits_{i=1}^{frac{N}{d}}mu(i)iglfloorfrac{lfloorfrac{N}{d} floor}{i}ig flooriglfloorfrac{lfloorfrac{M}{d} floor}{i}ig floor)

里面这个分母看起来不好，于是我们转换枚举顺序，设(k=i imes d)：

(sumlimits_{k=1}^Nlfloorfrac{N}{k} floorlfloorfrac{M}{k} floorsumlimits_{d|k}d imes mu(frac{k}{d}))

内层求和式又是一个非常喜闻乐见的式子，看到卷积式，又有id和(mu)，于是毫不犹豫地直接化简：(varphi=id*mu)：

(sumlimits_{k=1}^Nlfloorfrac{N}{k} floorlfloorfrac{M}{k} floorvarphi(k))

这个式子就非常简单了，首先一看是个整除分块的形式，可以处理欧拉函数的前缀和，然后快速求出答案。那么怎么求欧拉函数的前缀和呢？如果你还说线性筛的话，那么你这一段可能白学了。当然我们可以用杜教筛解决这个问题，于是我们就可以先筛欧拉函数前缀和，然后分块做了。

由此我们可见，通常我们要将待处理的式子化成较基本的函数的求和，这样才方便用杜教筛进行处理。

【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/ix35/p/11972989.html