JZOJ 4496. 【GDSOI 2016】第一题 互补约数

( ext{Problem})

求

[sum_{i=1}^n sum_{d|n} gcd(d, frac{i}{d}) ]

有 (n le 10^{11})

( ext{Analysis})

一眼就知道要欧拉反演（虽然考场写了莫反）
那么就要套路套路地推式子了
先给出欧拉反演的一般形式

[n = sum_{d|n} varphi(d) ]

然后对题目中的式子推导一波

[egin{aligned} sum_{i=1}^n sum_{d|i} gcd(d, frac{i}{d}) &= sum_{xy le n} gcd(x,y) \ &= sum_{xy le n} sum_{d|gcd(x,y)} varphi(d) \ &= sum_{d=1}^{sqrt n} varphi(d) sum_{d|x} sum_{d|y} 1 \ &= sum_{d=1}^{sqrt n} varphi(d) sum_{x=1}^{lfloor frac{n}{d^2} floor} lfloor frac{n}{xd^2} floor\ end{aligned} ]

发现这个式子后面可以数论分块
然后直接计算复杂度是 (O(sqrt n log sqrt n))
可以通过 (10^{11})

( ext{Code})

#include<cstdio>
#define LL long long
using namespace std;

LL n;

int vis[1000005], prime[600005], phi[1000005], totp;
inline void Euler()
{
	vis[1] = 1, phi[1] = 1;
	for(register int i = 2; i <= 1e6; i++)
	{
		if (!vis[i]) prime[++totp] = i, phi[i] = i - 1;
		for(register int j = 1; j <= totp && prime[j] * i <= 1e6; j++)
		{
			vis[prime[j] * i] = 1;
			if (i % prime[j]) phi[prime[j] * i] = (prime[j] - 1) * phi[i];
			else{phi[prime[j] * i] = prime[j] * phi[i]; break;}
		}
	}
}

inline LL solve()
{
	LL sum = 0;
	for(register LL d = 1; d * d <= n; d++)
	{
		LL s = 0, m = n / d / d, j;
		for(register LL i = 1; i * d * d <= n; i = j + 1) 
		{
			j = m / (m / i);
			s += (j - i + 1) * (m / i);
		}
		sum += phi[d] * s;
	}
	return sum;
}

int main()
{
	freopen("gcd.in", "r", stdin), freopen("gcd.out", "w", stdout);
	Euler(), scanf("%lld", &n), printf("%lld
", solve());
}

因为原式是一个 (gcd) 的形式，一个数算入贡献
所以我们走欧拉反演的路
但考场看到 (gcd) 直接走莫反了
于是推了这么一个东西

[sum_{d=1}^{sqrt n} sum_{k=1}^{sqrt lfloor frac{n}{d^2} floor} mu(k) sum_{i=1}^{lfloor frac{n}{d^2k^2} floor} lfloor frac{n}{id^2k^2} floor ]

单单数论分快暴力跑显然过不了
发现可以预处理后面一坨式子，因为 (d^2k^2) 的取值只有 (O(sqrt n)) 种