自动搬运

查看原文

来自洛谷，原作者为

	wolfind **

搬运于2025-08-24 22:30:25，当前版本为作者最后更新于2021-03-17 13:26:56，作者可能在搬运后再次修改，您可在原文处查看最新版

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

以下是正文

---

$n\le 8$

直接枚举 $p$ 然后计算方差即可。

写起来意外的麻烦

复杂度 $O(n!)$。

#include<bits/stdc++.h>
#define N 10
#define MOD 1000000007
using namespace std;
void read(int &x)
{
	x=0;
	char ch=getchar();
	while(ch<'0'||ch>'9')ch=getchar();
	while(ch>='0'&&ch<='9')
	{
		x=x*10+ch-'0';
		ch=getchar();
	}
}
int n,invn,ans;
int a[N],b[N];
int p[N],v[N];
int ksm(int x,int k)
{
	int res=1;
	while(k)
	{
		if(k&1)res=1ll*res*x%MOD;
		x=1ll*x*x%MOD;
		k>>=1; 
	}
	return res;
}
int gcd(int x,int y)
{
	if(x==0)return y;
	return gcd(y%x,x);
}
void dfs(int x)
{
	if(x>n)
	{
		int pj=0;
		for(int i=1;i<=n;i++)
		{
			int res=gcd(a[i],b[p[i]]);
			pj=(pj+res)%MOD;
		}
		pj=1ll*pj*invn%MOD;
		for(int i=1;i<=n;i++)
		{
			int res=gcd(a[i],b[p[i]]);
			ans=(ans+1ll*(res-pj)*(res-pj)%MOD)%MOD;
		}
		return;
	}
	for(int i=1;i<=n;i++)
		if(!v[i])
		{
			p[x]=i;
			v[i]=1;
			dfs(x+1);
			p[x]=0;
			v[i]=0;
		}
}
int main()
{
	scanf("%d",&n);
	invn=ksm(n,MOD-2);
	for(int i=1;i<=n;i++)
		read(a[i]);
	for(int i=1;i<=n;i++)
		read(b[i]);
	dfs(1);
	ans=1ll*ans*invn%MOD;
	printf("%d",ans);
}

---

$n\le 100$

如果你想要直接递推，比如说 $dp_i$ 表示只考虑前 $i$ 位的时候的答案，那么由于原题是个随机的排列，你就得知道前 $i$ 位到底用了什么位置，必须要 $2^n$ 来记录，肯定是做不了的。

所以我们把方差拆掉：$\begin{aligned}\dfrac{1}{n}\sum\limits_{i=1}^{n}(x_i-\dfrac{S}{n})^2=\dfrac{1}{n}(\sum\limits_{i=1}^{n}x_i^2-\dfrac{S^2}{n})\end{aligned}$

那么这题就变成了求两个东西：

1.数列中的数的平方的和 $\begin{aligned}sum1=\sum\limits_{i=1}^{n}x_i^2\end{aligned}$。

2.数列中的数的和的平方 $sum2=S^2$。

$sum1$ 是好弄的，只要数列中的每个数独立考虑即可。

问题就是 $sum2$，没法独立考虑，因为这题的性质不允许我们 $dp$，所以就直接把它拆掉！

具体地， $\begin{aligned}S^2=\sum\limits_{i=1}^{n}x_i^2+\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[i\ne j]x_ix_j\end{aligned}$

那么我们记录两个值：

$\begin{aligned}res1=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\gcd(a_i,b_j)^2\end{aligned}$

$\begin{aligned}res2=\sum\limits_{i1=1}^{n}\sum\limits_{i2=1}^{n}\sum\limits_{j1=1}^{n}\sum\limits_{j2=1}^{n}[i1 \ne i2][j1 \ne j2]\gcd(a_{i1},b_{j1})\times \gcd(a_{i2},b_{j2})\end{aligned}$

就会有 $sum1=res1\times (n-1)!$，$sum2=res1\times (n-1)!+res2\times (n-2)!$

最后答案就是 $ans=\dfrac{1}{n}(sum1-\dfrac{sum2}{n})$。

这样直接做就是 $O(n^4)$ 的，可以拿 $25$ 分。

下面这个代码真的是直接明了，就不写注释了。

#include<bits/stdc++.h>
#define N 110
#define MOD 1000000007
using namespace std;
void read(int &x)
{
	x=0;
	char ch=getchar();
	while(ch<'0'||ch>'9')ch=getchar();
	while(ch>='0'&&ch<='9')
	{
		x=x*10+ch-'0';
		ch=getchar();
	}
}
int n,invn,sum1,sum2,ans;
int a[N],b[N],mul[N];
int c[N][N];
int ksm(int x,int k)
{
	int res=1;
	while(k)
	{
		if(k&1)res=1ll*res*x%MOD;
		x=1ll*x*x%MOD;
		k>>=1;
	}
	return res;
}
int gcd(int x,int y)
{
	if(x==0)return y;
	return gcd(y%x,x);
}
int main()
{
	scanf("%d",&n);
	invn=ksm(n,MOD-2);
	for(int i=1;i<=n;i++)
		read(a[i]);
	for(int i=1;i<=n;i++)
		read(b[i]);
	mul[0]=1;
	for(int i=1;i<=n;i++)
		mul[i]=1ll*mul[i-1]*i%MOD;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			c[i][j]=gcd(a[i],b[j]);
			sum1=(sum1+1ll*c[i][j]*c[i][j])%MOD;
		}
	}
	sum1=1ll*sum1*mul[n-1]%MOD;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			for(int ii=1;ii<=n;ii++)
			{
				if(ii==i)continue;
				for(int jj=1;jj<=n;jj++)
				{
					if(jj==j)continue;
					int res1=c[i][j],res2=c[ii][jj];
					sum2=(sum2+1ll*res1*res2)%MOD;
				}
			}
		}
	}
	sum2=1ll*sum2*mul[n-2]%MOD;
	sum2=(sum1+sum2)%MOD;
	ans=(sum1-1ll*invn*sum2%MOD+MOD)%MOD;
	ans=1ll*ans*invn%MOD;
	printf("%d",ans);
}

---

$n\le 10^6$

（O(n^2)都是些奇怪的东西，就跳过了。）

然后就是想办法加快求 $res1,res2$ 的速度了。

首先是 $res1$，非常经典的 $\gcd$ 容斥：（知道的可以跳过了）

我们记录 $f_x$ 表示 $\gcd(a,b)=x$ 的 $(a,b)$ 有多少组。

也就是：$\begin{aligned}f_x=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[\gcd(a_i,b_j)=x]\end{aligned}$。

不过想要 $\gcd(a,b)=x$ 的话，至少得有 $x|a,x|b$ 吧。

也就是： $\begin{aligned}f_x=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[x|a_i][x|b_j]\end{aligned}$

但是这样的话如果 $\gcd(a,b)=kx$ 也会被算进去，所以要减去 $f_{2x},f_{3x},f_{4x}...$ 就把错的去掉了。

最终会有：$\begin{aligned}f_x=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[x|a_i][x|b_j]-\sum\limits_{i=2}^{\left\lfloor\frac{n}{x}\right\rfloor}f_{ix}\end{aligned}$

直接预处理 $\begin{aligned}ca_x=\sum\limits_{i=1}^{n}[x|a_i]\end{aligned}$，$\begin{aligned}cb_x=\sum\limits_{i=1}^{n}[x|b_i]\end{aligned}$，就可以从大往小求出 $f_i$ 了。

最后 $\begin{aligned}res1=\sum\limits_{i=1}^{n}f_i\times i^2\end{aligned}$

复杂度为 $O(\left\lfloor\dfrac{n}{1}\right\rfloor+\left\lfloor\dfrac{n}{2}\right\rfloor+...\left\lfloor\dfrac{n}{n}\right\rfloor)=O(n\log n)$。

---

然后就是本题的难点所在，$res2$ 该怎么求？

仔细观察原式：

$\begin{aligned}res2=\sum\limits_{i1=1}^{n}\sum\limits_{i2=1}^{n}\sum\limits_{j1=1}^{n}\sum\limits_{j2=1}^{n}[i1 \ne i2][j1 \ne j2]\gcd(a_{i1},b_{j1})\times \gcd(a_{i2},b_{j2})\end{aligned}$

其中要满足两个条件 $[i1 \ne i2][j1 \ne j2]$。

但是在计数的时候保证选的两个不同的数是难的，而选两个相同的数或者是任选两个数都是好做的。

所以这就启发我们去容斥！

具体来说我们先求出无限制下的答案，然后减去满足 $[i1=i2]$（注意这里是相等）而 $(j1,j2)$ 无限制的答案，再减去 $(i1,i2)$ 无限制而满足 $[j1=j2]$ 的答案，最后还要加回来满足 $[i1=i2]$ 且 $[j1=j2]$ 的答案。

有点繁琐，不过意会一下应该能明白。

我们把上面四个情况顺次记作 $tot1,tot2,tot3,tot4$，下面我们一个个求。

---

$tot1$

先把式子写出来：

$\begin{aligned}tot1=\sum\limits_{i1=1}^{n}\sum\limits_{i2=1}^{n}\sum\limits_{j1=1}^{n}\sum\limits_{j2=1}^{n}\gcd(a_{i1},b_{j1})\times \gcd(a_{i2},b_{j2})\end{aligned}$

观察一下发现：

$\begin{aligned}tot1=(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\gcd(a_{i},b_{j}))^2\end{aligned}$

想到之前求 $res1$ 时用的 $f_i$ 了吗？借用一下即可。

$\begin{aligned}tot1=(\sum\limits_{i=1}^{n}f_i\times i)^2\end{aligned}$

---

$tot2,tot3$

因为它俩一模一样，所以就讲讲 $tot2$ 咋求。

$\begin{aligned}tot2=\sum\limits_{i=1}^{n}\sum\limits_{j1=1}^{n}\sum\limits_{j2=1}^{n}\gcd(a_{i},b_{j1})\times \gcd(a_{i},b_{j2})\end{aligned}$

（因为 $i1=i2$，所以就直接改成 $i$ 了）

再观察一下可以发现 $j1,j2$ 能用一个平方弄掉：

$\begin{aligned}tot2=\sum\limits_{i=1}^{n}(\sum\limits_{j=1}^{n}\gcd(a_{i},b_{j}))^2\end{aligned}$

然后就用到我们的欧拉反演了！

当然，反演的式子很简单：$\begin{aligned}n=\sum\limits_{d|n}\varphi(d)\end{aligned}$。

咋用呢，其实就是把gcd换掉！

$\begin{aligned}\gcd(a_i,b_j)=\sum\limits_{d|\gcd(a_i,b_j)}\varphi(d)\end{aligned}$

代入即有：

$\begin{aligned} tot2&= \sum\limits_{i=1}^{n}(\sum\limits_{j=1}^{n}\sum\limits_{d|\gcd(a_i,b_j)}\varphi(d))^2\\ &=\sum\limits_{i=1}^{n}(\sum\limits_{j=1}^{n}\sum\limits_{d|a_i,d|b_j}\varphi(d))^2\\ &=\sum\limits_{i=1}^{n}(\sum\limits_{d|a_i}\sum\limits_{j=1}^{n}[d|b_j]\varphi(d))^2\\ \end{aligned}$

我们记录一个辅助数组 $\begin{aligned}gb_x=\sum\limits_{i=1}^{n}[x|b_i]\varphi(x)\end{aligned}$。

那么上式就是：$\begin{aligned}tot2= \sum\limits_{i=1}^{n}(\sum\limits_{d|a_i}gb_d)^2\end{aligned}$

这样就可以 $O(n\log n)$ 啦。

如果你不用欧拉反演而是继续用之前的容斥的话，复杂度就是 $O(nlog^2n)$ 的，不过仔细想想把这个容斥的过程预处理出来，其实就是欧拉函数！

还是挺简单的。

---

$tot4$

$\begin{aligned}tot4=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\gcd(a_{i},b_{j})^2=res1\end{aligned}$

没了。

---

然后根据 $res2=tot1-tot2-tot3+tot4$，这题就做完了。

复杂度 $O(n\log n)$。

写代码的时候还要注意一个细节，就是所有的 $a$ 我们是要存进桶中的，这样不同的数才会去计算一遍，否则同一个数计算多次的话复杂度是错的。

代码长度也就暴力的两倍

#include<bits/stdc++.h>
#define N 1000010
#define MOD 1000000007
using namespace std;
void read(int &x)
{
	x=0;
	char ch=getchar();
	while(ch<'0'||ch>'9')ch=getchar();
	while(ch>='0'&&ch<='9')
	{
		x=x*10+ch-'0';
		ch=getchar();
	}
}
int n,m,invn,sum1,sum2,ans;
int a[N],b[N],mul[N];
int ca[N],cb[N],f[N],phi[N],ga[N],gb[N];
int ksm(int x,int k)
{
	int res=1;
	while(k)
	{
		if(k&1)res=1ll*res*x%MOD;
		x=1ll*x*x%MOD;
		k>>=1; 
	}
	return res;
}
int gcd(int x,int y)
{
	if(x==0)return y;
	return gcd(y%x,x);
}
int main()
{
	scanf("%d",&n);
	invn=ksm(n,MOD-2);
	for(int i=1;i<=n;i++)
	{
		int x;
		read(x);
		++a[x];//要存进桶里 
		m=max(m,x);
	}
	for(int i=1;i<=n;i++)
	{
		int x;
		read(x);
		++b[x];
		m=max(m,x);
	}
	mul[0]=1;
	for(int i=1;i<=n;i++)
		mul[i]=1ll*mul[i-1]*i%MOD;//前缀积 
	for(int i=m;i>=1;i--)
		for(int j=1;j*i<=m;j++)
			ca[i]+=a[j*i];//把多少个a是i的倍数算进ca[i]中 
	for(int i=m;i>=1;i--)
		for(int j=1;j*i<=m;j++)
			cb[i]+=b[j*i];
	for(int i=m;i>=1;i--)
	{
		f[i]=1ll*ca[i]*cb[i]%MOD;
		for(int j=2;j*i<=m;j++)
			(f[i]+=MOD-f[j*i])%=MOD;//容斥,求出有多少gcd(a,b)=i,存进f[i] 
		sum1=(sum1+1ll*f[i]*i%MOD*i)%MOD;//其实是res1
		sum2=(sum2+1ll*f[i]*i)%MOD;//tot1
	}
	sum2=1ll*sum2*sum2%MOD;
	
	for(int i=1;i<=m;i++)
	{
		phi[i]=(phi[i]+i)%MOD;
		//这也是phi的一种求法，这样写主要是为了体现出其容斥的内涵 
		int tmpa=1ll*phi[i]*ca[i]%MOD,tmpb=1ll*phi[i]*cb[i]%MOD;
		(gb[i]+=tmpb)%=MOD;
		(ga[i]+=tmpa)%=MOD;
		sum2=(sum2-1ll*gb[i]*gb[i]%MOD*a[i]%MOD+MOD)%MOD;//tot2
		sum2=(sum2-1ll*ga[i]*ga[i]%MOD*b[i]%MOD+MOD)%MOD;//tot3
		for(int j=2;j*i<=m;j++)
			(phi[j*i]+=MOD-phi[i])%=MOD,(gb[j*i]+=tmpb)%=MOD,(ga[j*i]+=tmpa)%=MOD;
	}
	
	sum2=(sum2+sum1)%MOD;//tot4
	
	sum1=1ll*sum1*mul[n-1]%MOD;
	sum2=1ll*sum2*mul[n-2]%MOD;
	sum2=(sum1+sum2)%MOD;
	ans=(sum1-1ll*invn*sum2%MOD+MOD)%MOD;
	ans=1ll*ans*invn%MOD;
	printf("%d",ans);
}

---

$n\le 5\times 10^6$

没有被放进去的数据范围。

这咋整啊！卡常？

其实是高维前缀和优化。

复杂度上看，高维前缀和的复杂度是 $\sum\limits_{p}\left\lfloor\dfrac{n}{p}\right\rfloor$ （$p$ 是质数）$=O(n\log \log n)$，卡5000000还是很轻松的。

那么高维前缀和优化到底是啥呢？

其实就是我们想要求一个东西 $\begin{aligned}f_i=\sum\limits_{j=1}^{ij\le n}g_{ij}\end{aligned}$，其中 $g$ 是已知的。

那么我们可以先把 $j=2$ 的求出来，再加上 $j=3$ 的，再加上 $j=5$ 的......

而在加的过程中我们可以把之前算过的 $f$ 拿来用，跟前缀和一样。而且为了不算重，枚举的 $j$ 都是质数。

而这可以看作是一个多维的空间，其中每一维都是由一个质数 $p$ 的倍数构成的，因此原理其实就是高维前缀和。

代码也很短，如下：

for(int i=1;i<=n;i++)f[i]=g[i];
	
for(int i=1;i<=p[0];i++)
	for(int j=n/p[i];j>=1;j--)
			f[j]+=f[j*p[i]];

那么我们原来的求 $ca[i]$ 的代码就可以改掉了，但是有些是无法改的！

我们举例一下原来所有要 $O(n\log n)$ 求的东西现在可以做到的复杂度：

1.$phi[i]$： 直接线性筛。$O(n)$

2.$ca[i],cb[i]$：高维前缀和。$O(n\log\log n)$

3.$ga[i],gb[i]$：高维前缀和。$O(n\log\log n)$

4.$f[i]$：因为要由自己容斥出自己，所以无法高维前缀和优化。$O(n\log n)$

好家伙，经典 $\gcd$ 容斥部分的求解成为了复杂度上限。

不过我们仍然有办法！

我们想，反正都是容斥，这里也用欧拉反演不好吗？

比如说：

$\begin{aligned} tot1&=(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\gcd(a_{i},b_{j}))^2\\ &=(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{d|a_i,d|b_j}\varphi(d))^2\\ &=(\sum\limits_{d=1}^{n}ca_dcb_d\varphi(d))^2 \end{aligned}$

要用 $f$ 求的还有一个 $res1$：

$\begin{aligned}res1=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\gcd(a_i,b_j)^2\end{aligned}$

啊这怎么欧拉反演。。。

我们参考欧拉反演的式子 $n=\sum\limits_{d|n}\varphi(d)$，也设计一个反演！

具体的，我们找到一个函数 $\phi(x)$，使得 $n^2=\sum\limits_{d|n}\phi(d)$ 就能像上面那样反演了。

这里根据狄利克雷卷积有 $id_2=\phi\ \ast\ I$ ，$\phi=id_2\ \ast\ \mu$ 所以 $\phi$ 是个积性函数，我们只要知道 $\phi(p^k)$ 是多少即可。

我们可以轻易发现 $\phi(p^k)=(p^k)^2-\phi(1)-\phi(p)-...-\phi(p^{k-1})=(p^k)^2-(p^{k-1})^2$

然后就没了。

复杂度 $O(n\log\log n)$。

全村最快的代码

在洛谷上跑 $n=5000000$ 只用 $500ms$ （不计读入）。

#include<bits/stdc++.h>
#define N 5000010
#define MOD 1000000007
using namespace std;
void read(int &x)
{
	x=0;
	char ch=getchar();
	while(ch<'0'||ch>'9')ch=getchar();
	while(ch>='0'&&ch<='9')
	{
		x=x*10+ch-'0';
		ch=getchar();
	}
}
int f[N],phi[N],v[N],p[N];
int n,m,invn,sum1,sum2,ans;
int a[N],b[N],mul[N];
int ca[N],cb[N],ga[N],gb[N];
int ksm(int x,int k)
{
	int res=1;
	while(k)
	{
		if(k&1)res=1ll*res*x%MOD;
		x=1ll*x*x%MOD;
		k>>=1; 
	}
	return res;
}
int gcd(int x,int y)
{
	if(x==0)return y;
	return gcd(y%x,x);
}
int main()
{
	read(n);
	invn=ksm(n,MOD-2);
	for(int i=1;i<=n;i++)
	{
		int x;
		read(x);
		++a[x];
		++ca[x];//这边要多一步，等价于ca[x]=a[x] 
		m=max(m,x);
	}
	for(int i=1;i<=n;i++)
	{
		int x;
		read(x);
		++b[x];
		++cb[x];
		m=max(m,x);
	}
	mul[0]=1;
	for(int i=1;i<=n;i++)
		mul[i]=1ll*mul[i-1]*i%MOD;
	
	phi[1]=1;f[1]=1;
	for(int i=2;i<=m;i++)
	{
		if(!v[i])p[++p[0]]=i,phi[i]=i-1,f[i]=1ll*i*i%MOD-1;//f[i]就是我们自己求的积性函数！ 
		for(int j=1;j<=p[0]&&p[j]*i<=m;j++)
		{
			v[i*p[j]]=1;
			if(i%p[j]==0)
			{
				phi[i*p[j]]=phi[i]*p[j];
				f[i*p[j]]=1ll*f[i]*p[j]%MOD*p[j]%MOD;
				break;
			}
			else phi[i*p[j]]=phi[i]*phi[p[j]],f[i*p[j]]=1ll*f[i]*f[p[j]]%MOD;
		}
	}
	
	for(int i=1;i<=p[0];i++)
		for(int j=m/p[i];j>=1;j--)
			ca[j]+=ca[j*p[i]],cb[j]+=cb[j*p[i]];//高维前缀和！ 
	
	for(int i=m;i>=1;i--)
	{
		int tmp=1ll*ca[i]*cb[i]%MOD;
		sum1=(sum1+1ll*tmp*f[i])%MOD;
		sum2=(sum2+1ll*tmp*phi[i])%MOD;
	}
	sum2=1ll*sum2*sum2%MOD;
	
	for(int i=1;i<=m;i++)//初始化 
	{
		gb[i]=1ll*phi[i]*cb[i]%MOD;
		ga[i]=1ll*phi[i]*ca[i]%MOD;
	}
	
	for(int i=1;i<=p[0];i++)
		for(int j=1;j<=m/p[i];j++)
			(ga[j*p[i]]+=ga[j])%=MOD,(gb[j*p[i]]+=gb[j])%=MOD;//高维前缀和！ 
	
	for(int i=1;i<=m;i++)
	{
		sum2=(sum2-1ll*gb[i]*gb[i]%MOD*a[i]%MOD+MOD)%MOD;
		sum2=(sum2-1ll*ga[i]*ga[i]%MOD*b[i]%MOD+MOD)%MOD;
	}
	
	sum2=(sum2+sum1)%MOD;
	
	sum1=1ll*sum1*mul[n-1]%MOD;
	sum2=1ll*sum2*mul[n-2]%MOD;
	sum2=(sum1+sum2)%MOD;
	ans=(sum1-1ll*invn*sum2%MOD+MOD)%MOD;
	ans=1ll*ans*invn%MOD;
	printf("%d",ans);
}