自动搬运

查看原文

来自洛谷，原作者为

	ArrTue ac

搬运于2025-08-24 22:12:40，当前版本为作者最后更新于2021-11-19 12:32:32，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

update: 修正了错误代码。

题目传送门

前置知识：

欧拉函数，分块。

题意：

• 给出 $n$ 个数 $a_{1\dots n}$，$q$ 次询问，每次给出 $l,r$，求 $\sum\limits_{i=l}^r\sum\limits_{j=l}^r\gcd(a_i,a_j)$。

• 保证数列、询问均随机生成，$n,q,a_i\le10^5$。

分析：

$\begin{aligned} &\sum\limits_{i=l}^r\sum\limits_{j=l}^r\gcd(a_i,a_j)\\ =&\sum\limits_{i=l}^r\sum\limits_{j=l}^r\sum\limits_{d\mid a_i\land d\mid a_j}\varphi(d)&(1)\\ =&\sum\limits_{d=1}^{\max(a)}\varphi(d)\sum\limits_{i=l}^r\sum\limits_{j=l}^r[d\mid a_i\land d\mid a_j]\\ =&\sum\limits_{d=1}^{\max(a)}\varphi(d)(\sum\limits_{i=l}^r[d\mid a_i])^2 \end{aligned}$

其中 $(1)$ 用到了欧拉函数的性质 $n=\sum\limits_{d\mid n}\varphi(d)$（也可用莫反推出一个新函数 $f(d)$，实际上 $f$ 就是 $\varphi$）

官方解法中使用了莫队维护 $\sum\limits_{i=l}^r[d\mid a_i]$（不会莫队可以跳过这一段）：由于 $a_i$ 随机生成，故其因数个数期望较少（另一篇题解中认为是一个常数 $\alpha$，本题解认为是值域内的总因子个数 $v\ln v$ 除以值域 $v$，即 $\ln v$），移动莫队指针时可以直接枚举所有 $a_i$ 的因数 $d$ 并更新对应的 $\sum\limits_{i=l}^r[d\mid a_i]$。期望复杂度 $\mathcal O(n\sqrt n\ln v)$，惊人地可过。

但翻了 AC 代码，发现本题实际上有一种期望 $\mathcal O((n+q)\sqrt n)$ 的分块写法。

首先对值域分块，所有不大于 $V$ 的 $d$ 直接暴力求出 $[d\mid a_i]$ 的前缀和，询问时 $\mathcal O(V)$ 即可处理这一部分的数，总复杂度 $\mathcal O((n+q)V)$。

对于大于 $V$ 的 $d$，由于值比较大，所以随机情况下 $a$ 中 $d$ 的倍数较少。假设所有为 $d$ 的倍数的 $a$ 为 $a_{k_1},a_{k_2},\dots,a_{k_m}$。当询问左端点 $l\in(k_{x-1},k_x]$，右端点 $r\in[k_y,k_{y+1})$ 时，$\sum\limits_{i=l}^r[d\mid a_i]$ 即为 $y-x+1$；对答案贡献为 $\varphi(d)(y-x+1)^2$。对于每一个询问，显然不能循环每一个 $d$ 计算 $x,y$。但注意到因为 $m$ 较小，有序对 $(x,y)$ 只有 $\mathcal O(m^2)$ 种情况，可以从这里入手。

不妨将询问按右端点排序，对于 $y$ 相同的询问，$d$ 的贡献只与左端点有关。而左端点落在一个区间内时，贡献都是一样的。考虑差分，则 $d$ 只对 $y+1$ 个位置有影响。每次 $y$ 改变只需修改这几处的值即可。总修改次数为 $\mathcal O(m^2)$。

可以证明 $m^2$ 的期望为 $\mathcal O(\dfrac{n^2}{d^2})$（证明 link）。则对于所有 $d$，总修改次数期望为 $\mathcal O(\sum\limits_{d=V}^{\max(a)}\dfrac{n^2}{d^2})=\mathcal O(\dfrac{n^2}V)$。

对每个询问需要计算出差分数组的前缀和。考虑分块，使修改 $\mathcal O(1)$，查询 $\mathcal O(\sqrt n)$。总时间复杂度 $\mathcal O((n+q)V+\dfrac{n^2}V+q\sqrt n)$，取 $V=\sqrt n$，复杂度即为 $\mathcal O((n+q)\sqrt n)$。

思路：

1. 求欧拉函数。

2. 对值域分块，分别处理。

具体说明值域分块后大于 $V$ 的 $d$ 如何处理。

开 $\max(a)$ 个 vector，存下每个 $d$ 的 $k_i$。考虑右端点所在区间由 $[k_{y_0},k_{y_0+1})$ 到 $[k_{y_0+1},k_{y_0+2})$，差分数组如何变化：$\forall x\in[1,y_0+1]$，左端点在 $(k_{x-1},k_x]$ 的贡献由 $\varphi(d)(y_0-x+1)^2$ 变为 $\varphi(d)(y_0-x+2)^2$，增加 $\varphi(d)(2(y_0-x+1)+1)$。

在差分数组上，体现为位置 1 增加 $\varphi(d)(2y_0+1)$，位置 $k_x+1\ (x\in[1,y_0])$ 减少 $2\varphi(d)$，位置 $k_{y_0}+1$ 减少 $\varphi(d)$。

分块求和，记得处理对应位置块和。

#include <bits/stdc++.h>
#define rep(i, l, r) for(int i=l, _=r; i<=_; ++i)
typedef long long ll;
using namespace std;
inline int read() {
	int res=0; bool k=1; char ch;
	while(!isdigit(ch=getchar())) k=ch^45;
	do res=(res<<1)+(res<<3)+(ch&15); while(isdigit(ch=getchar()));
	return k? res: -res;
}
const int mN=1e5+9, mP=1e4, mL=320;
int n, q, L, V, a[mN], sum[mN][mL];
ll ans[mN];

int pri[mP], cntp, phi[mN];
inline void init(int n) {
	phi[1]=1;
	for(int i=2, j, u; i<=n; ++i) {
		if(!phi[i]) pri[++cntp]=i, phi[i]=i-1;
		for(j=1; (u=pri[j]*i)<=n && i/pri[j]*pri[j]^i; ++j) phi[u]=phi[i]*(pri[j]-1);
		if(u<=n) phi[u]=phi[i]*pri[j];
	}
}

struct Que {int x, y, id;} qn[mN];
bool operator < (Que q1, Que q2) {return q1.y<q2.y;}

ll b[mN], sumb[mL];	//b 为差分数组，sumb 为分块的差分数组
vector<int> num[mN];	//用于存 ki
inline void sol(int x, int p) {
	if(x<=V) return ++sum[p][x], void();	//小于 V，直接前缀和
	int y0=(int) num[x].size(), t=phi[x];
	b[1]+=(ll) (y0<<1|1)*t, sumb[0]+=(ll) (y0<<1|1)*t;
	b[p+1]-=t, sumb[(p+1)/L]-=t;
	t<<=1;
	for(int i=y0-1; i>=0; --i) b[num[x][i]+1]-=t, sumb[(num[x][i]+1)/L]-=t;
	num[x].push_back(p);
}

inline ll sqr_(int x) {return (ll) x*x;}

int main() {
	init(1e5), n=read(), q=read(), L=sqrt(n), V=sqrt(1e5);
	if(L<2) L=2;
	rep(i, 1, n) a[i]=read();
	rep(i, 1, q) qn[i].x=read(), qn[i].y=read(), qn[i].id=i;
	sort(qn+1, qn+q+1);

	rep(i, 1, q) {
		rep(r, qn[i-1].y+1, qn[i].y) {	//扩展 y
			memcpy(sum[r], sum[r-1], sizeof sum[r]);
			for(int j=1; j*j<=a[r]; ++j) if(a[r]%j==0) {
				sol(j, r);
				if(j*j^a[r]) sol(a[r]/j, r);
			}
		}
		int x=qn[i].x, y=qn[i].y, id=qn[i].id;
		rep(j, 0, x/L-1) ans[id]+=sumb[j];	//大于 V，值域分块
		rep(j, x/L*L, x) ans[id]+=b[j];
		rep(v, 1, V) ans[id]+=phi[v]*sqr_(sum[y][v]-sum[x-1][v]);	//小于 V，直接前缀和差分
	}
	rep(i, 1, q) printf("%lld\n", ans[i]);
	return 0;
}