【题目】GCD of Divisors
【题意】给定f(n)=Σd|n gcd(d,n/d)的前缀和F(n)，n=10^15。
【算法】莫比乌斯反演
【题解】参考：任之洲数论函数.pdf
这个范围显然杜教筛也是做不了的，而且考虑直接化简f(n)也遇到了困难，所以考虑将前缀和的Σ一起化简。
$$F(n)=sum_{i=1}^{n}sum_{d|i}(d,frac{i}{d})$$
这一步很常见的是第一重改为枚举倍数，但这样化简后面就推不下去了。
这道题必须最后转成$sigma_0(n)$才能解出来。
所以直接枚举gcd值
$$F(n)=sum_{d=1}^{n}dsum_{i=1}^{n}sum_{g|i}[(g,frac{i}{g})=d]$$
这里gcd(g,i/g)=d，说明i中必须至少包含2个d，那么令i=i/d^2，g即可任取i的因子，最终的g和i/g各乘d即可，所以可以进行如下化简。（关键①）
$$F(n)=sum_{d=1}^{n}dsum_{i=1}^{frac{n}{d^2}}sum_{g|i}[(g,frac{i}{g})=1]$$
转化成φ希望不大，所以直接莫比乌斯反演。
$$F(n)=sum_{d=1}^{n}dsum_{i=1}^{frac{n}{d^2}}sum_{g|i}sum_{d'|g cap d'|frac{i}{g}}mu(d')$$
$$F(n)=sum_{d=1}^{n}sum_{d'=1}^{n}d*mu(d')sum_{i=1}^{frac{n}{d^2}}sum_{g|i}[d'|g cap d'|frac{i}{g}]$$
这里和上面一样，都是要求d'|g&&d'|i/g，因此从i中提取2个d'，即令i=i/d'^2。
$$F(n)=sum_{d=1}^{n}sum_{d'=1}^{n}d*mu(d')sum_{i=1}^{frac{n}{(dd')^2}}sum_{g|i}1$$
会发现后面是约数个数。（关键②）
$$F(n)=sum_{d=1}^{n}sum_{d'=1}^{n}d*mu(d')sum_{i=1}^{frac{n}{(dd')^2}}sigma_0(i)$$
前面部分发现d*μ(d')好像可以卷积到φ，考虑合并dd‘来构造卷积，令d=dd'。（关键③）
$$F(n)=sum_{d=1}^{sqrt{n}}sum_{g|d}g*mu(frac{n}{g})sum_{i=1}^{frac{n}{(dd')^2}}sigma_0(i)$$
这里d只枚举到√n，因为d>√n时后面的Σ=0，没有贡献。因此可以缩小实际枚举范围。（关键④）
然后根据狄利克雷卷积μ*id=φ可以化简
$$F(n)=sum_{d=1}^{sqrt{n}}varphi(d)sum_{i=1}^{frac{n}{(dd')^2}}sigma_0(i)$$
大功告成！
其中，约数个数的前缀和可以进行分块取值优化，如下
$$sum_{i=1}^{n}sigma_0(i)=sum_{i=1}^{n}sum_{d|i}1=sum_{i=1}^{n}sum_{j=1}^{frac{n}{i}}1$$
$$sum_{i=1}^{n}sigma_0(i)=sum_{i=1}^{n}left lfloor frac{n}{i} ight floor$$
然后线性筛φ的过程中求解即可。
复杂度分析：
$$sum_{i=1}^{sqrt{n}}O(sqrt{frac{n}{i^2}})=sum_{i=1}^{sqrt{n}}O(frac{sqrt{n}}{i})=O(sqrt{n} ln sqrt{n})$$
倒数第二步将√n提到最外面，Σ内就是调和数列，和近似为ln n。
#include<cstdio> #include<cmath> #define int long long const int maxn=32000000; int phi[maxn],n,prime[maxn],tot; int solve(int n){ int pos,sum=0; for(int i=1;i<=n;i=pos+1){ pos=n/(n/i); sum+=(pos-i+1)*(n/i); } return sum; } #undef int int main(){ #define int long long scanf("%lld",&n); int ans=1*solve(n),N=(int)sqrt(n)+1;phi[1]=1;//1 for(int i=2;i<=N;i++){ if(!phi[i])phi[prime[++tot]=i]=i-1; for(int j=1;j<=tot&&i*prime[j]<=N;j++){ if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;} phi[i*prime[j]]=phi[i]*(prime[j]-1); } ans+=phi[i]*solve(n/(i*i)); } printf("%lld",ans); return 0; }