题解 LibreOJ #6539. 奇妙数论题
LibreOJ #6539. 奇妙数论题
同步更新在莫比乌斯反演做题记录。
题如其名，一道巧妙莫反题。
给定长为 \(n\) 的排列 \(a\) 求
我们有部分分是 \(a_i=i\) 的，所以一直到 \(70\) 分的部分分就是求：
就 Luogu P1390 公约数的和 稍微改一下就有了：
正解的瓶颈其实就是排列变成无序后无从下手。
所以我们先将式子转换成：
这里就是根据上一题的套路，\(\gcd(i,j)\) 看成 \(\gcd(i,j)\times 1\) 于是可以在数论分块的部分直接带入对应下标 \(a_i\) 的 \(\gcd\)。
但是求 \(\gcd(a_i,a_j)\) 也就是 \(\gcd(a_{idx},a_{jdx})\) 的复杂度降不下来，考虑怎么做。
设：
我们观察这个函数，其实就是要求在 \(a\) 这个排列中下标分别为 \(x\) 的倍数的元素彼此之间 \(\gcd\) 的和。
我们设 \(\mathbb S_x\) 是大小为 \(x\) 的倍数的下标的集合。
于是有：
我们进一步抽象，我们将 \(i\in\mathbb S_x\) 映射到每一个 \(a_i\)，\(\forall i\in \mathbb S_x,\) 设 \(a_i\) 的集合为 \(\mathbf S_x\)。
现在有：
\(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）：
我们知道有 \(\mu * \mathrm{id}=\varphi\)。
观察上面蓝色部分的 \(\displaystyle{\color{blue}{\sum_{d|T}d\times\mu(\frac{T}{d})}}\) 就是 \(\varphi(T)\)
最后得到：
我们回到 \(\texttt{ans}\) 的部分。
有：
具体实现方法：
- 首先我们筛出得到 \(\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;
}
