数论函数练习记录（一）

参考资料：oi-wiki 莫比乌斯反演｜oi-wiki 筛法｜oi-wiki 欧拉函数｜blog｜blog

[CQOI2007]余数求和

题目链接

这道题主要是练习数论分块。

$$G(n,k)=\sum_{i=1}^nk\bmod i=\sum_{i=1}^nk-i\cdot\left\lfloor\frac{k}{i}\right\rfloor=nk-\sum_{i=1}^ni\cdot\left\lfloor\frac{k}{i}\right\rfloor$$

减号后面的值进行数论分块，每一块是常数 $\left\lfloor\frac{k}{l}\right\rfloor$ 乘以一个等差数列求和 $l+\cdots+r=\frac{(l+r)(r-l+1)}{2}$。

#include<algorithm>
#include<cstdio>

using namespace std;

long long n, k, ans;

int main(){
	scanf("%lld%lld", &n, &k);
	ans = 1ll * n * k;
	for(long long l = 1, r; l <= n; l = r + 1){
		if(k / l == 0)	r = n;
		else	r = min(n, k / (k / l));
		ans -= (k / l) * (l + r) * (r - l + 1) / 2;
	}
	printf("%lld\n", ans);
	return 0;
}

---

[POI2007]ZAP-Queries

题目链接

$$\begin{align} \sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=k]&=\sum_{i=1}^n\sum_{j=1}^m\left[\gcd\left(\frac{i}{k},\frac{j}{k}\right)=1\right]&&\gcd(i,j)=k\iff\gcd\left(\frac{i}{k},\frac{j}{k}\right)=1\\ &=\sum_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}[\gcd(i,j)=1]&&i\to ki,\;j\to kj\\ &=\sum_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}\epsilon\left(\gcd(i,j)\right)&&写作单位元函数\\ &=\sum_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}\sum_{d\mid\gcd(i,j)}\mu(d)&&\epsilon=\mu*1\;并展开\\ &=\sum_{d=1}^{\min\left(\left\lfloor\frac{n}{k}\right\rfloor,\left\lfloor\frac{m}{k}\right\rfloor\right)}\sum_{i=1}^{\left\lfloor\frac{n}{kd}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{kd}\right\rfloor}\mu(d)&&求和号换序并置\;i\to di,\,j\to dj\\ &=\sum_{d=1}^{\min\left(\left\lfloor\frac{n}{k}\right\rfloor,\left\lfloor\frac{m}{k}\right\rfloor\right)}\mu(d)\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor \end{align}$$

接下来进行数论分块，每一块是常数 $\left\lfloor\frac{n}{kl}\right\rfloor\left\lfloor\frac{m}{kl}\right\rfloor$ 乘上 $\mu(l)+\cdots+\mu(r)$，前缀和预处理即可。

这里 $r$ 的计算：$\left\lfloor\frac{n}{kl}\right\rfloor=\left\lfloor\frac{n}{kr}\right\rfloor\leqslant\frac{n}{kr}\implies r\leqslant\frac{\frac{n}{k}}{\left\lfloor\frac{n}{kl}\right\rfloor}\implies r=\left\lfloor\frac{\frac{n}{k}}{\left\lfloor\frac{n}{kl}\right\rfloor}\right\rfloor$.

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N = 50005;

int k;

int mu[N], pList[N], pID, muS[N];
bool notP[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	mu[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0)	pList[++pID] = i, mu[i] = -1;
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				mu[i * pList[j]] = 0;
				break;
			}
			else	mu[i * pList[j]] = -mu[i];
		}
	}
	for(int i = 1; i <= n; i++)	muS[i] = muS[i-1] + mu[i];
}

LL calc(int n, int m){
	LL res = 0;
	for(int l = 1, r; l <= min(n / k, m / k); l = r + 1){
		r = min((n / k) / ((n / k) / l), (m / k) / ((m / k) / l));
		res += 1ll * (muS[r] - muS[l-1]) * (n / k / l) * (m / k / l);
	}
	return res;
}

int main(){
	Euler(50000);
	int T; for(scanf("%d", &T); T; T--){
		int a, b; scanf("%d%d%d", &a, &b, &k);
		printf("%lld\n", calc(a, b));
	}
	return 0;
}

---

[HAOI2011]Problem b

题目链接

通过二维前缀和就可以把问题转换成上一道题。

不能乱开 long long，否则常数太大会 TLE。

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N = 50005;

int k;

int mu[N], pList[N], pID, muS[N];
bool notP[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	mu[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0)	pList[++pID] = i, mu[i] = -1;
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				mu[i * pList[j]] = 0;
				break;
			}
			else	mu[i * pList[j]] = -mu[i];
		}
	}
	for(int i = 1; i <= n; i++)	muS[i] = muS[i-1] + mu[i];
}

LL calc(int n, int m){
	LL res = 0;
	for(int l = 1, r; l <= min(n / k, m / k); l = r + 1){
		r = min((n / k) / ((n / k) / l), (m / k) / ((m / k) / l));
		res += 1ll * (muS[r] - muS[l-1]) * (n / k / l) * (m / k / l);
	}
	return res;
}

int main(){
	Euler(50000);
	int T; for(scanf("%d", &T); T; T--){
		int a, b, c, d; scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
		printf("%lld\n", calc(b, d) - calc(a-1, d) - calc(b, c-1) + calc(a-1, c-1));
	}
	return 0;
}

---

[SPOJ5971]LCMSUM

题目链接

$$\begin{align} \sum_{i=1}^n\text{lcm(i,n)}&=\sum_{i=1}^n\frac{i\cdot n}{\gcd(i,n)}=n+\sum_{i=1}^{n-1}\frac{i\cdot n}{\gcd(i,n)}=n+\sum_{i=1}^{n-1}\frac{(n-i)\cdot n}{\gcd(n-i,n)}&&根据对称性做变换\;i\to n-i\\ &=n+\frac{1}{2}\sum_{i=1}^{n-1}\frac{n^2}{\gcd(i,n)}=\frac{n}{2}+\frac{1}{2}\sum_{i=1}^n\frac{n^2}{\gcd(i,n)}&&根据\;\gcd(i,n)=\gcd(n-i,n)\;进行合并\\ &=\frac{n}{2}+\frac{n^2}{2}\sum_{d\mid n}\frac{\varphi\left(\frac{n}{d}\right)}{d}&&枚举\;\gcd(i,n)\;的值，乘上\;\gcd(i,n)=d\;的个数，\\&&&即\;\gcd(i/d,n/d)=1\;的个数，即\;\varphi(n/d)\\ &=\frac{n}{2}+\frac{n^2}{2}\sum_{d\mid n}\frac{d\cdot\varphi(d)}{n}=\frac{n}{2}\left(\sum_{d\mid n}d\cdot\varphi(d)+1\right)&&根据对称性做变换\;d\to\frac{n}{d} \end{align}$$

令 $g(n)=\sum\limits_{d\mid n}d\cdot\varphi(d)$，如果我们可以预处理出 $g(n)$，那么我们就可以 $O(1)$ 回答了。我们这样预处理：枚举 $i$，将所有 $i\mid n$ 的 $g(n)$ 加上贡献 $i\cdot\varphi(i)$，这样预处理的复杂度是调和级数，其阶为 $O(n\lg n)$. 【这种预处理方式很重要！】

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N = 1000005;

int T, n;

int phi[N], pList[N], pID;
bool notP[N];
LL g[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	phi[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0)	pList[++pID] = i, phi[i] = i - 1;
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				phi[i * pList[j]] = phi[i] * pList[j];
				break;
			}
			else	phi[i * pList[j]] = phi[i] * (pList[j] - 1);
		}
	}
}

int main(){
	Euler(1000000);
	for(int i = 1; i <= 1000000; i++)
		for(int j = 1; i * j <= 1000000; j++)
			g[i*j] += 1ll * i * phi[i];
	for(scanf("%d", &T); T; T--){
		scanf("%d", &n);
		printf("%lld\n", (g[n] + 1) * n / 2);
	}
	return 0;
}

---

[国家集训队]Crash的数字表格

题目链接

$$\begin{align} \sum_{i=1}^n\sum_{j=1}^m\text{lcm}(i,j)&=\sum_{i=1}^n\sum_{j=1}^m\frac{i\cdot j}{\gcd(i,j)}\\ &=\sum_{d=1}^{\min(n,m)}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]d\cdot i\cdot j&&考虑\;d=\gcd(i,j)\;的贡献，并置\;i\to di,\,j\to dj\\ &=\sum_{d=1}^{\min(n,m)}d\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\epsilon(\gcd(i,j))\cdot i\cdot j&&写作单位元函数\\ &=\sum_{d=1}^{\min(n,m)}d\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\sum_{r\mid\gcd(i,j)}\mu(r)\cdot i\cdot j&&\epsilon=1*\mu\;并展开\\ &=\sum_{d=1}^{\min(n,m)}d\sum_{r=1}^{\min\left(\left\lfloor\frac{n}{d}\right\rfloor,\left\lfloor\frac{m}{d}\right\rfloor\right)}\sum_{i=1}^{\left\lfloor\frac{n}{dr}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{dr}\right\rfloor}r^2\mu(r)\cdot i\cdot j&&经典的求和号换序并置\;i\to ri,\,j\to rj\\ &=\sum_{d=1}^{\min(n,m)}d\sum_{r=1}^{\min\left(\left\lfloor\frac{n}{d}\right\rfloor,\left\lfloor\frac{m}{d}\right\rfloor\right)}r^2\mu(r)\sum_{i=1}^{\left\lfloor\frac{n}{dr}\right\rfloor}i\sum_{j=1}^{\left\lfloor\frac{m}{dr}\right\rfloor}j\\ &=\sum_{d=1}^{\min(n,m)}d\sum_{r=1}^{\min\left(\left\lfloor\frac{n}{d}\right\rfloor,\left\lfloor\frac{m}{d}\right\rfloor\right)}r^2\mu(r)\cdot\frac{\left(\left\lfloor\frac{n}{dr}\right\rfloor+1\right)\left\lfloor\frac{n}{dr}\right\rfloor}{2}\cdot\frac{\left(\left\lfloor\frac{m}{dr}\right\rfloor+1\right)\left\lfloor\frac{m}{dr}\right\rfloor}{2}&&后两项是独立的，可求和 \end{align}$$

我们设

$$S(n,m)=\sum_{r=1}^{\min(n,m)}r^2\mu(r)\cdot\frac{\left(\left\lfloor\frac{n}{r}\right\rfloor+1\right)\left\lfloor\frac{n}{r}\right\rfloor}{2}\cdot\frac{\left(\left\lfloor\frac{m}{r}\right\rfloor+1\right)\left\lfloor\frac{m}{r}\right\rfloor}{2}$$

它是可以通过数论分块计算的，每一块是常数 $\frac{\left(\left\lfloor\frac{n}{l}\right\rfloor+1\right)\left\lfloor\frac{n}{l}\right\rfloor}{2}\cdot\frac{\left(\left\lfloor\frac{m}{l}\right\rfloor+1\right)\left\lfloor\frac{m}{l}\right\rfloor}{2}$ 乘上 $l^2\mu(l)+\cdots+r^2\mu(r)$，前缀和预处理即可。

那么原式可写作：

$$\sum_{d=1}^{\min(n,m)}d\cdot S\left(\left\lfloor\frac{n}{d}\right\rfloor,\left\lfloor\frac{m}{d}\right\rfloor\right)$$

而它又是一个可以数论分块计算的，每一块是常数 $S\left(\left\lfloor\frac{n}{l}\right\rfloor,\left\lfloor\frac{m}{l}\right\rfloor\right)$ 乘上 $l+\cdots+r=\frac{(l+r)(r-l+1)}{2}$.

所以这道题是数论分块的嵌套，复杂度为 $O(\sqrt n\cdot\sqrt n)=O(n)$.

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N = 10000005;
const LL MOD = 20101009;

int mu[N], pList[N], pID;
bool notP[N];
LL sum[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	mu[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0)	pList[++pID] = i, mu[i] = -1;
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				mu[i * pList[j]] = 0;
				break;
			}
			else	mu[i * pList[j]] = -mu[i];
		}
	}
	for(int i = 1; i <= n; i++)
		sum[i] = (sum[i-1] + 1ll * i * i % MOD * mu[i] % MOD) % MOD;
}

LL S(int n, int m){
	LL res = 0;
	if(n > m)	swap(n, m);
	for(int l = 1, r; l <= n; l = r + 1){
		r = min(n / (n / l), m / (m / l));
		res += 1ll * (1 + n / l) * (n / l) % MOD * (1 + m / l) % MOD * (m / l) % MOD * ((sum[r] - sum[l-1]) % MOD + MOD) % MOD * 15075757 % MOD;
		res %= MOD;
	}
	return res;
}

int n, m;
LL ans;

int main(){
	Euler(10000000);
	scanf("%d%d", &n, &m);
	if(n > m)	swap(n, m);
	for(int dl = 1, dr; dl <= n; dl = dr + 1){
		dr = min(n / (n / dl), m / (m / dl));
		ans += (1ll * (dl + dr) * (dr - dl + 1) % MOD * 10050505 % MOD) * S(n / dl, m / dl) % MOD;
		ans %= MOD;
	}
	printf("%lld\n", ans);
	return 0;
}

其实这道题还能继续优化：先鸽一下（雾

---

[luoguP2398]GCD SUM

题目链接

解法一：和上一道题差不多。

$$\begin{align} \sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)&=\sum_{d=1}^n\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}d[\gcd(i,j)=1]&&考虑\;d=\gcd(i,j)\;的贡献，并置\;i\to di,\,j\to dj\\ &=\sum_{d=1}^nd\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{r\mid \gcd(i,j)}\mu(r)&&写作单位元函数并展开\\ &=\sum_{d=1}^nd\sum_{r=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{i=1}^{\left\lfloor\frac{n}{dr}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{dr}\right\rfloor}\mu(r)&&求和号换序并置\;i\to ri,\,j\to rj\\ &=\sum_{d=1}^nd\sum_{r=1}^{\left\lfloor\frac{n}{d}\right\rfloor}{\left\lfloor\frac{n}{dr}\right\rfloor}^2\mu(r) \end{align}$$

数论分块嵌套。

复杂度：$O(n)$

解法二：

$$\begin{align} \sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)&=\sum_{i=1}^n\sum_{j=1}^n\varphi(\gcd(i,j))*1(\gcd(i,j))&&由\;\text{id}=\varphi*1\\ &=\sum_{i=1}^n\sum_{j=1}^n\sum_{d\mid\gcd(i,j)}\varphi(d)&&把卷积打开\\ &=\sum_{d=1}^n\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\varphi(d)&&求和号换序并置\;i\to di,\,j\to dj\\ &=\sum_{d=1}^n{\left\lfloor\frac{n}{d}\right\rfloor}^2\varphi(d) \end{align}$$

数论分块。

复杂度：$O(\sqrt n)$ （不算线性筛）

注：其实解法一继续推导可得解法二的表达式，不过推导的下一步不是很显然（要卷回去）：

$$\begin{align} \sum_{d=1}^nd\sum_{r=1}^{\left\lfloor\frac{n}{d}\right\rfloor}{\left\lfloor\frac{n}{dr}\right\rfloor}^2\mu(r)&=\sum_{d=1}^nd\sum_{r=1}^{n}[d\mid r]{\left\lfloor\frac{n}{r}\right\rfloor}^2\mu\left(\frac{r}{d}\right)&&置\;r\to\frac{r}{d}\\ &=\sum_{r=1}^n{\left\lfloor\frac{n}{r}\right\rfloor}^2\sum_{d\mid r}\text{id}(d)\cdot\mu\left(\frac{r}{d}\right)&&求和号换序\\ &=\sum_{r=1}^n{\left\lfloor\frac{n}{r}\right\rfloor}^2\varphi(r)&&由\;\text{id}*\mu=\varphi \end{align}$$

这个博客讲的也是这个，并给这种方法起了个名字叫做欧拉反演，哈哈。

解法一：

#include<bits/stdc++.h>

using namespace std;

const int N = 100005;

int mu[N], pList[N], pID;
bool notP[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	mu[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0){
			pList[++pID] = i;
			mu[i] = -1;
		}
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				mu[i * pList[j]] = 0;
				break;
			}
			else	mu[i * pList[j]] = -mu[i];
		}
	}
}

int n, muS[N];

int main(){
	Euler(100000);
	for(int i = 1; i <= 100000; i++)	muS[i] = muS[i-1] + mu[i];
	scanf("%d", &n);
	long long ans = 0;
	for(int dl = 1, dr; dl <= n; dl = dr + 1){
		dr = n / (n / dl);
		long long res = 0;
		for(int l = 1, r; l <= n / dl; l = r + 1){
			r = (n / dl) / (n / dl / l);
			res += 1ll * (muS[r] - muS[l-1]) * (n / dl / l) * (n / dl / l);
		}
		ans += 1ll * (dl + dr) * (dr - dl + 1) / 2 * res;
	}
	printf("%lld\n", ans);
	return 0;
}

解法二：

#include<bits/stdc++.h>

using namespace std;

const int N = 100005;

int phi[N], pList[N], pID;
bool notP[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	phi[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0){
			pList[++pID] = i;
			phi[i] = i - 1;
		}
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				phi[i * pList[j]] = phi[i] * pList[j];
				break;
			}
			else	phi[i * pList[j]] = phi[i] * (pList[j] - 1);
		}
	}
}

int n;
long long phiS[N];

int main(){
	Euler(100000);
	for(int i = 1; i <= 100000; i++)	phiS[i] = phiS[i-1] + phi[i];
	scanf("%d", &n);
	long long ans = 0;
	for(int l = 1, r; l <= n; l = r + 1){
		r = n / (n / l);
		ans += 1ll * (phiS[r] - phiS[l-1]) * (n / l) * (n / l);
	}
	printf("%lld\n", ans);
	return 0;
}

---

[NOI2010]能量采集

题目链接

题目要求的就是

$$\sum_{i=1}^n\sum_{j=1}^m\left(2\gcd(i,j)-1\right)=2\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)-nm$$

和上一题一样了。

#include<bits/stdc++.h>

using namespace std;

const int N = 100005;

int phi[N], pList[N], pID;
bool notP[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	phi[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0){
			pList[++pID] = i;
			phi[i] = i - 1;
		}
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				phi[i * pList[j]] = phi[i] * pList[j];
				break;
			}
			else	phi[i * pList[j]] = phi[i] * (pList[j] - 1);
		}
	}
}

int n, m;
long long phiS[N];

int main(){
	Euler(100000);
	for(int i = 1; i <= 100000; i++)	phiS[i] = phiS[i-1] + phi[i];
	scanf("%d%d", &n, &m); if(n > m) swap(n, m);
	long long ans = 0;
	for(int l = 1, r; l <= n; l = r + 1){
		r = min(n / (n / l), m / (m / l));
		ans += 1ll * (phiS[r] - phiS[l-1]) * (n / l) * (m / l);
	}
	printf("%lld\n", 2 * ans - 1ll * n * m);
	return 0;
}

---

[SDOI2015]约数个数和

题目链接

引理：

$$d(i\cdot j)=\sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)=1]$$

证明：约数个数公式：$d(n)=d\left(\prod\limits_{i=1}^k{p_i}^{r_i}\right)=\prod\limits_{i=1}^k(r_i+1)$。考虑 $ij$ 的一个质因子 $p$，设 $i=i’\cdot{p}^{r_i}$ ，$j=j’\cdot{p}^{r_j}$，则 $p$ 对答案有 $r_i+r_j+1$ 的贡献。又考虑右边的式子，由于 $\gcd(x,y)=1$ 的约束，$p$ 只能从 $i$ 或 $j$ 其中一边取出，或者不取出，贡献仍是 $r_i+r_j+1$。证毕。

于是：

$$\begin{align} \sum_{i=1}^n\sum_{j=1}^md(i\cdot j)&=\sum_{i=1}^n\sum_{j=1}^m\sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)=1]&&由引理\\ &=\sum_{x=1}^n\sum_{y=1}^m\sum_{i=1}^{\left\lfloor\frac{n}{x}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{y}\right\rfloor}\epsilon(\gcd(x,y))&&求和号换序并置\;i\to xi,\,j\to yj\\ &=\sum_{x=1}^n\sum_{y=1}^m{\left\lfloor\frac{n}{x}\right\rfloor}{\left\lfloor\frac{m}{y}\right\rfloor}\sum_{d\mid\gcd(x,y)}\mu(d)&&\epsilon=\mu*1\;并展开，i\,和\,j\,的和式计算出来\\ &=\sum_{d=1}^{\min(n,m)}\sum_{x=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{y=1}^{\left\lfloor\frac{m}{d}\right\rfloor}{\left\lfloor\frac{n}{dx}\right\rfloor}{\left\lfloor\frac{m}{dy}\right\rfloor}\mu(d)&&求和号换序并置\;x\to dx,\,y\to dy\\ &=\sum_{d=1}^{\min(n,m)}\mu(d)\sum_{x=1}^{\left\lfloor\frac{n}{d}\right\rfloor}{\left\lfloor\frac{n}{dx}\right\rfloor}\sum_{y=1}^{\left\lfloor\frac{m}{d}\right\rfloor}{\left\lfloor\frac{m}{dy}\right\rfloor}\\ &=\sum_{d=1}^{\min(n,m)}\mu(d)\cdot S\left(\left\lfloor\frac{n}{d}\right\rfloor\right)\cdot S\left(\left\lfloor\frac{m}{d}\right\rfloor\right)&&设\,S(i)=\sum_{r=1}^i\left\lfloor\frac{i}{r}\right\rfloor \end{align}$$

$S(1…n)$ 可以通过数论分块 $O(n\sqrt n)$ 地预处理出来，随后上式可通过数论分块 $O(\sqrt n)$ 地回答。

#include<bits/stdc++.h>

using namespace std;

const int N = 50005;

int mu[N], pList[N], pID;
bool notP[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	mu[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0){
			pList[++pID] = i;
			mu[i] = -1;
		}
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				mu[i * pList[j]] = 0;
				break;
			}
			else	mu[i * pList[j]] = -mu[i];
		}
	}
}

int T, n, m, muS[N];
long long S[N];

int main(){
	Euler(50000);
	for(int i = 1; i <= 50000; i++)	muS[i] = muS[i-1] + mu[i];
	for(int i = 1; i <= 50000; i++){
		for(int l = 1, r; l <= i; l = r + 1){
			r = i / (i / l);
			S[i] += (r - l + 1) * (i / l);
		}
	}
	for(scanf("%d", &T); T; T--){
		scanf("%d%d", &n, &m);
		if(n > m)	swap(n, m);
		long long ans = 0;
		for(int dl = 1, dr; dl <= n; dl = dr + 1){
			dr = min(n / (n / dl), m / (m / dl));
			ans += 1ll * (muS[dr] - muS[dl-1]) * S[n / dl] * S[m / dl];
		}
		printf("%lld\n", ans);
	}
	return 0;
}

---

[luoguP2257]YY的GCD

题目链接

$$\begin{align} \sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=\text{prime}]&=\sum_{p=1}^{\min(n,m)}\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=p]\quad(p\in\text{prime})&&枚举质数\;p\\ &=\sum_{p=1}^{\min(n,m)}\sum_{d=1}^{\min\left(\left\lfloor\frac{n}{p}\right\rfloor,\left\lfloor\frac{m}{p}\right\rfloor\right)}\mu(d)\left\lfloor\frac{n}{pd}\right\rfloor\left\lfloor\frac{m}{pd}\right\rfloor\quad(p\in\text{prime})&&和之前做过的题目一样的推导\\ &=\sum_{T=1}^{\min(n,m)}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum_{p\mid T}\mu\left(\frac{T}{p}\right)\quad(p\in\text{prime})&&置\;T=pd\\ \end{align}$$

最后一步的代换非常重要！详细逻辑如下：

$$\sum\limits_{p=1}^{\min(n,m)}\sum\limits_{d=1}^{\min\left(\left\lfloor\frac{n}{p}\right\rfloor,\left\lfloor\frac{m}{p}\right\rfloor\right)}f(p,d)=\sum\limits_{1\leqslant p\cdot d\leqslant\min(n,m)}f(p,d)\overset{T=p\cdot d}{===}\sum\limits_{T=1}^{\min(n,m)}\sum\limits_{p\mid T}f\left(p,\frac{T}{p}\right)$$

然后令 $g(T)=\sum\limits_{p\mid T}{\mu\left(\frac{T}{p}\right)}\quad(p\in\text{prime})$，这是可以预处理的，即枚举质数 $p$，将 $p\mid n$ 的 $g(n)$ 加上贡献 $\mu\left(\frac{n}{p}\right)$，复杂度是 $O(n\cdot H(n))\sim O(n\lg n)$

随后数论分块，每一块是常数 $\left\lfloor\frac{n}{l}\right\rfloor\left\lfloor\frac{m}{l}\right\rfloor$ 乘上 $g(l)+\cdots+g(r)$，前缀和预处理即可。

#include<bits/stdc++.h>

using namespace std;

const int N = 10000005;

int mu[N], pList[N], pID;
bool notP[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	mu[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0){
			pList[++pID] = i;
			mu[i] = -1;
		}
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				mu[i * pList[j]] = 0;
				break;
			}
			else	mu[i * pList[j]] = -mu[i];
		}
	}
}

long long g[N];

int main(){
	Euler(10000000);
	for(int i = 1; i <= pID; i++)
		for(int j = pList[i]; j <= 10000000; j += pList[i])
			g[j] += mu[j / pList[i]];
	for(int i = 1; i <= 10000000; i++)	g[i] += g[i-1];
	int T; for(scanf("%d", &T); T; T--){
		int n, m; scanf("%d%d", &n, &m); if(n > m) swap(n, m);
		long long ans = 0;
		for(int l = 1, r; l <= n; l = r + 1){
			r = min(n / (n / l), m / (m / l));
			ans += 1ll * (g[r] - g[l-1]) * (n / l) * (m / l);
		}
		printf("%lld\n", ans);
	}
	return 0;
}

---

[luoguP2586]GCD

题目链接

求的东西和上一道题基本一样，但是没有多组数据，时限加紧。

注意上一道题为了 $O(\sqrt n)$ 地回答而进行了 $O(n\lg n)$ 的预处理，这个预处理在这道题里会超时。但是这道题只需要 $O(n)$ 回答就行了。

重新推一波式子：

$$\begin{align} \sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=\text{prime}]&=\sum_{p=1}^{n}\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=p]\quad(p\in\text{prime})&&枚举质数\;p\\ &=\sum_{p=1}^n\sum_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{p}\right\rfloor}[\gcd(i,j)=1]\quad(p\in\text{prime})&&置\;i\to pi,\,j\to pj\\ &=\sum_{p=1}^n\left(2\sum_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor}\sum_{j=1}^{i}[\gcd(i,j)=1]-1\right)\quad(p\in\text{prime})&&由 \gcd 的对称性\\ &=\sum_{p=1}^n\left(2\sum_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor}\varphi(i)-1\right)\quad(p\in\text{prime})&&由\;\varphi(n)\;的定义式 \end{align}$$

预处理 $\varphi(n)$ 的前缀和，就可以 $O(n)$ 解决了。

上面的推导中最妙的一点是用上了 $\varphi(n)=\sum\limits_{i=1}^n[\gcd(i,n)=1]$ 这个定义式。

#include<bits/stdc++.h>

using namespace std;

const int N = 10000005;

int phi[N], pList[N], pID;
bool notP[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	phi[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0){
			pList[++pID] = i;
			phi[i] = i - 1;
		}
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				phi[i * pList[j]] = phi[i] * pList[j];
				break;
			}
			else	phi[i * pList[j]] = phi[i] * (pList[j] - 1);
		}
	}
}

long long phiS[N];

int main(){
	Euler(10000000);
	for(int i = 1; i <= 10000000; i++)	phiS[i] = phiS[i-1] + phi[i];
	int n; scanf("%d", &n);
	long long ans = 0;
	for(int i = 1; pList[i] <= n; i++)
		ans += 2 * phiS[n / pList[i]] - 1;
	printf("%lld\n", ans);
	return 0;
}