【GDSOI 2016】第一题 互补约数

Description

简化题意：求
$$\sum _{i=1}^n\sum _{j|i}gcd(j,\frac{i}{j})$$

Input

输入包含一行，一个正整数 n。

Output

输出只有一行， F(n)。

Sample Input

Sample Input1：
10

Sample Input2：
1000

Sample Output

Sample Output1：
32

Sample Output2：
12776

Data Constraint

$n<=10^{11}$

Solution

两个sigma外加gcd，直接考虑反演
$$\sum _{i=1}^n\sum _{j|i}gcd(j,\frac{i}{j})$$
这里是枚举i和i的约数，我们转换为枚举i和i的倍数
$$\sum _{i=1}^n\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }gcd(i,j)$$
套路反演

设$$f[d]=\sum _{i=1}^n\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }gcd(i,j)=d$$

设$$g[d]=\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }f[id]$$

反演得$$f[d]=\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }g[id]\ast \mu (i)$$
$$ans=\sum _{d=1}^{n}d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }g[id]\ast \mu (i)$$

---

设 $T=id$，我们发现gcd(i,j)的值不可能大于$\sqrt{n}$，所以
$$ans=\sum _{T=1}^{\sqrt{n}}g[T]\ast \sum _{i|T}\frac{T}{i}\ast \mu (i)$$
将${\sum }_{T=1}^{\sqrt{n}}{\sum }_{i|T}\frac{T}{i}\ast \mu (i)$预处理，设为$a(i)$
为了方便预处理，设$i\ast j=T$变形为
$$a(i)=\sum _{i=1}^{\sqrt{n}}\sum _{i=1}^{\lfloor \frac{\sqrt{n}}{d}\rfloor }j\ast \mu (i)$$

我们知道，如果gcd(x,y)为T的倍数，那么x,y也是T的倍数
那么$x=x^{\prime }\ast T,y=y^{\prime }\ast T$
因为$x\ast y<=n$
所以$x^{\prime }$
枚举$x^{\prime }\ast T\ast y^{\prime }\ast T<=n$
所以$$g[T]=\sum _{x=1}^{\lfloor \frac{n}{T^2}\rfloor }\lfloor \frac{n}{T^{2}x}\rfloor$$

---

于是，你会惊奇的发现，答案变成了
$$ans=\sum _{T=1}^{\sqrt{n}}g[T]\ast a(t)$$

很简单对吧
但还是过不了
求g[t]的时候用分块优化一下，而t本身不能分块（为什么？）

时间复杂度$O(n^{\frac{3}{4}})$

Code

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define ll long long
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define N 401000
using namespace std;
ll n,mu[N],a[N],bz[N],s[N];
int main()
{
	freopen("gcd.in","r",stdin);freopen("gcd.out","w",stdout);
	scanf("%lld",&n);
	ll m=sqrt(n);
	mu[1]=1;
	fo(i,2,m)
	{
		if(!bz[i]) s[++s[0]]=i,mu[i]=-1;
		fo(j,1,s[0])
		{
			int k=i*s[j];
			if(k>m) break;
			bz[k]=1;
			if(i%s[j]==0) {mu[k]=0;break;}
			mu[k]=-mu[i];
		}
	}
	fo(i,1,m) fo(j,1,m/i) a[i*j]+=j*mu[i];
	ll ans=0;
	for(ll t=1;t<=m;t++)
	{
		ll jy=0,q=n/(t*t);
		for(ll x=1,y=1;x<=q;y=x+1,x=min((q/(q/y)),q)) 
		{
			jy+=(x-y+1)*(q/x);
			if(x>=q) break;
		}
		ans+=(jy*a[t]);
	}
	printf("%lld",ans);
}