题解 LibreOJ #6539. 奇妙数论题

LibreOJ #6539. 奇妙数论题

同步更新在莫比乌斯反演做题记录。

---

题如其名，一道巧妙莫反题。

给定长为 \(n\) 的排列 \(a\) 求

\[\sum_{i=1}^n\sum_{j=1}^n\gcd(a_i,a_j)\times \gcd(i,j) \]

我们有部分分是 \(a_i=i\) 的，所以一直到 \(70\) 分的部分分就是求：

\[\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)^2 \]

就 Luogu P1390 公约数的和 稍微改一下就有了：

\[\sum_{k=1}^nk^2\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\lfloor\frac{n}{k\times d}\rfloor\lfloor\frac{n}{k\times d}\rfloor \]

正解的瓶颈其实就是排列变成无序后无从下手。

所以我们先将式子转换成：

\[\begin{aligned} \texttt{ans}&=\sum_{i=1}^n\sum_{j=1}^n\gcd(a_i,a_j)\times \gcd(i,j)\\ &=\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=d]\times \gcd(a_i,a_j)\\ &=\sum_{d=1}^nd\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)\sum_{i=1}^{\lfloor\frac{n}{d\times x}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d\times x}\rfloor}\gcd(a_{idx},a_{jdx}) \end{aligned} \]

这里就是根据上一题的套路，\(\gcd(i,j)\) 看成 \(\gcd(i,j)\times 1\) 于是可以在数论分块的部分直接带入对应下标 \(a_i\) 的 \(\gcd\)。

但是求 \(\gcd(a_i,a_j)\) 也就是 \(\gcd(a_{idx},a_{jdx})\) 的复杂度降不下来，考虑怎么做。

设：

\[f(x)=\sum_{i=1}^{\lfloor\frac{n}{x}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{x}\rfloor}\gcd(a_{ix},a_{jx}) \]

我们观察这个函数，其实就是要求在 \(a\) 这个排列中下标分别为 \(x\) 的倍数的元素彼此之间 \(\gcd\) 的和。

我们设 \(\mathbb S_x\) 是大小为 \(x\) 的倍数的下标的集合。

于是有：

\[f(x)=\sum_{i\in\mathbb S_x}\sum_{j\in\mathbb S_x}\gcd(a_i,a_j) \]

我们进一步抽象，我们将 \(i\in\mathbb S_x\) 映射到每一个 \(a_i\)，\(\forall i\in \mathbb S_x,\) 设 \(a_i\) 的集合为 \(\mathbf S_x\)。

现在有：

\[f(x)=\sum_{i\in\mathbf S_x}\sum_{j\in\mathbf S_x}\gcd(i,j)\Rightarrow \sum_d d\sum_{i\in\mathbf S_x}\sum_{j\in \mathbf S_x}[\gcd(i,j)=d] \]

\(i,j\) 的集合与范围发生了变化，显然我们不能再将 \([\gcd(i,j)=d]\) 化为 \([\gcd(i,j)=1]\) 去做。

我们有另一种方法求 \([x=d]\) 这一类的式子。

我们观察 \([x=d]\) 长得很像 \(\epsilon(x)=[x=1]\)。

我们有 \(\mu*1=\epsilon\) 所以有 \(\displaystyle \epsilon(x)=\sum_{d|x}\mu(d)\)。

那么我们的 \([x=d]\) 可以化为 \(\displaystyle[\frac{T}{d}=1\space\mathrm{ and } \space T|d]\)。

类似的我们展开得 \(\displaystyle [x=d]=\sum_{d|T,T|x}\mu(\frac{T}{d})\)

得到（这个地方看官方题解看了很久才懂的，感觉很 tricky）：

\[\begin{aligned} f(x)&=\sum_d d\sum_{i\in\mathbf S_x}\sum_{j\in \mathbf S_x}\sum_{d|T,T|\gcd(i,j)}\mu(\frac{T}{d})\\ &=\sum_d d\sum_{d|T}\mu(\frac{T}{d})(\sum_{i\in \mathbf S_x}[T|i])^2\\ &=\sum_T({\color{blue}{\sum_{d|T}d\times\mu(\frac{T}{d})}})(\sum_{i\in \mathbf S_x}[T|i])^2 \end{aligned} \]

我们知道有 \(\mu * \mathrm{id}=\varphi\)。

观察上面蓝色部分的 \(\displaystyle{\color{blue}{\sum_{d|T}d\times\mu(\frac{T}{d})}}\) 就是 \(\varphi(T)\)

最后得到：

\[f(x)=\sum_T \varphi(T)(\sum_{i\in \mathbf S_x}[T|i])^2 \]

我们回到 \(\texttt{ans}\) 的部分。

有：

\[\texttt{ans}=\sum_{d=1}^nd\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)f(dx) \]

具体实现方法：

• 首先我们筛出得到 \(\forall x\in[1,n]\) 每个 \(\mathbf S_x\)。
• 依次统计对于每个 \(\mathbf S_x\) 每个元素的的因子的贡献。
• 得到 \(f(x)\) 并筛出 \(\mu(x)\) 后直接求和即可。

复杂度分析的话，看到 loj 官方题解的评论区中有这样的结论：

我们设对于 \(x\)，\(x\) 的因子个数为 \(d(x)\)。

As Vaclav Kotesovec said Aug 30 2018, \(\displaystyle \sum_{k\leq n}d^2(k)\)is asymptotic to \(\Theta(\frac{1}{\pi^2}n\log^3 n+n\log^2 n)\).

代码（其实就是出题人的写法，好像多数采用了埃筛）：

#include<bits/stdc++.h>

#define ll long long
#define ull unsigned long long
#define INL inline
#define Re register

//Tosaka Rin Suki~

INL int read()
{
	int x=0,w=1;char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();if(ch=='-')w=-1,ch=getchar();
	while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+ch-48,ch=getchar();return x*w;
}

const int N=1e5+5,MOD=1e9+7;

std::vector <int> factors[N];
std::vector <int> pack;
int mu[N],g[N],f[N],d[N],n,a[N],ans;


int main()
{
	//freopen(".in","r",stdin);
	//freopen(".out","w",stdout);
	n=read();for(int i=1;i<=n;i++)a[i]=read();
	mu[1]=1;
	for(int i=1;i<=n;i++)
	{
		factors[i].push_back(i);
		for(int j=i*2;j<=n;j+=i)
			mu[j]-=mu[i],factors[j].push_back(i);
	}
	for(int i=1;i<=n;i++)
		for(int j=i;j<=n;j+=i)
			g[j]=(g[j]+mu[j/i]*i)%MOD;
	for(int i=1;i<=n;i++)
	{
		for(int j=i;j<=n;j+=i)
			for(int k=0;k<factors[a[j]].size();k++)
				++d[factors[a[j]][k]],pack.push_back(factors[a[j]][k]);
		for(int j=0;j<pack.size();j++)
			if(d[pack[j]])
				f[i]=(f[i]+1ll*d[pack[j]]*d[pack[j]]*g[pack[j]])%MOD,d[pack[j]]=0;
		pack.clear();
	}
	for(int d=1;d<=n;d++)
	{
		int s=0;
		for(int x=1;x<=n/d;x++)s=(s+1ll*mu[x]*f[d*x])%MOD;
		ans=(ans+1ll*s*d)%MOD;
	}
	printf("%d\n",ans);
	return 0;
}