【SPOJ-GCDEX】GCD Extreme（欧拉函数）

题目：

分析：

求：

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

这道题给同届新生讲过，由于种种原因只讲了 (O(n)) 预处理欧拉函数 (O(n)) 查询的暴力做法，顺带提了一句 “这题能根号查询” 被教练嘴了 QAQ 。以及小恐龙给我说有 (O(nlog n)) 预处理 (O(1)) 查询的另一种写法。

重点是前几天某学长讲课讲这道题，才知道有 (O(n)) 预处理 (O(1)) 查询的神仙做法，并且据说因为 SPOJ 上时限 0.237s ，复杂度比这个高的都过不去（只交了最后一种，所以我也不知道是真的假的……）。下面分别介绍这三种做法（在无特殊说明的情况下，本文中所有除法均为向下取整）：

方法一

为了方便计算，先给和式中添上 (sum_{i=1}^{n}gcd(i,i)=frac{n(n+1)}{2})，再换一下枚举顺序，得到：

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

然后先枚举 (gcd(i,j)) 的值，并把 (i) 变成 (icdot d)，(j) 变成 (jcdot d) ：

[sum_{d=1}^{n}dsum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{i}{d}}[gcd(i,j)=1]-frac{n(n+1)}{2} ]

（方括号表示其中式子成立时值为 (1) ，否则为 (0) ）

根据欧拉函数的定义（ (varphi(x)) 表示不超过 (x) 的与 (x) 互质的正整数的个数），得到：

[sum_{d=1}^{n}dsum_{i=1}^{frac{n}{d}}varphi(frac{i}{d})-frac{n(n+1)}{2}=sum_{d=1}^{n}dcdot S_varphi(frac{n}{d})-frac{n(n+1)}{2} ]

（其中 (S_varphi) 表示 (varphi) 的前缀和，即 (S_varphi(x)=sum_{i=1}^{x}varphi(x)) ）

用线性筛处理出欧拉函数前缀和后数论分块即可。预处理 (O(n)) ，查询 (O(sqrt{n})) 。

方法二

设：

[F(i)=sum_{j=1}^{i}gcd(i,j) ]

则答案为：

[sum_{i=1}^{n}F(i)-frac{n(n+1)}{2} ]

考虑如何计算 (F(i)) 。同样先枚举 (gcd(i,j)) 的值：

[F(i)=sum_{d|i}dsum_{j=1}^{frac{i}{d}}[gcd(i,j)=1]=sum_{d|i}dvarphi(frac{i}{d}) ]

于是先线性筛 (varphi) ，然后枚举 (d) ，每个 (d) 向所有 (F(icdot d) (0<ileq frac{n}{d})) 贡献 (dcdot varphi(i)) 。最后处理 (F) 函数的前缀和。预处理时间复杂度为 (n) 乘上 (sum_{i=1}^{n}frac{1}{i})，约为 (O(nlog n)) 。可以 (O(1)) 查询。

方法三：

接着方法二：

[F(n)=sum_{d|n}dvarphi(frac{n}{d})=sum_{icdot d=n}dcdot varphi(i) ]

则答案为：

[sum_{i=1}^{n}F(i)-frac{n(n+1)}{2} ]

考虑如果已经求出了 (F(n)) ，则有（其中 (p) 是质数且与 (n) 互质）：

[F(np)=pF(n) ]

对于和式的每一项，这个 (p) 要么乘到 (d) 上，要么乘到 (varphi(i)) 上，而欧拉函数是积性函数所以 (varphi(ip)=pcdot varphi(i)) 。

这个性质很好，如果能快速计算 (n=p^k) ((p) 是质数）的情况就能愉快地线性筛了。（为什么？请看【知识总结】线性筛_杜教筛_Min25筛的线性筛部分）

考虑当 (n=p^k) ，(F(n)=sum_{i=0}^{k}p^{k-i}varphi(p^i)) 。根据欧拉函数的公式：

[F(n)=n+sum_{i=1}^{k}p^{k-i}p^icdot frac{p-1}{p}=kp^{k-1}(p-1)+p^k ]

（注意当 (i=0) 时 (p^0=1) 中没有因子 (p) ，要特殊处理）

所以当 (k>1) ， (F(p^{k-1})=(k-1)p^{k-2}(p-1)+p^{k-1}) ，给它乘上 (p) 再加上 (p^{k-1}(p-1)) 就是 (F(p^k)) 。 (k=1) （即 (n) 是质数）的情况直接手算一下，等于 (2n-1) 。

于是可以线性筛出 (F) 函数然后处理前缀和， (O(n)) 预处理 (O(1)) 查询。

代码：

（方法三）

#include <cstdio>
#include <cctype>
#include <cstring>
#include <algorithm>
using namespace std;

namespace zyt
{
	template<typename T>
	inline bool read(T &x)
	{
		char c;
		bool f = false;
		x = 0;
		do
			c = getchar();
		while (c != EOF && c != '-' && !isdigit(c));
		if (c == EOF)
			return false;
		if (c == '-')
			f = true, c = getchar();
		do
			x = x * 10 + c - '0', c = getchar();
		while (isdigit(c));
		if (f)
			x = -x;
		return true;
	}
	template<typename T>
	inline void write(T x)
	{
		static char buf[20];
		char *pos = buf;
		if (x < 0)
			putchar('-'), x = -x;
		do
			*pos++ = x % 10 + '0';
		while (x /= 10);
		while (pos > buf)
			putchar(*--pos);
	}
	typedef long long ll;
	const int N = 1e6 + 10;
	ll f[N];
	int prime[N], last[N], cnt;
	bool mark[N];
	void init()
	{
		f[1] = 1;
		for (int i = 2; i < N; i++)
		{
			if (!mark[i])
				prime[cnt++] = last[i] = i, f[i] = i * 2 - 1;
			for (int j = 0; j < cnt && (ll)i * prime[j] < N; j++)
			{
				int k = i * prime[j];
				mark[k] = true;
				if (i % prime[j] == 0)
				{
					last[k] = last[i] * prime[j];
					if (last[k] == k)
						f[k] = f[i] * prime[j] + (ll)i * (prime[j] - 1);
					else
						f[k] = f[i / last[i]] * f[last[k]];
				}
				else
				{
					last[k] = prime[j];
					f[k] = f[i] * f[prime[j]];
				}
			}
		}
		for (int i = 1; i < N; i++)
			f[i] += f[i - 1];
	}
	int work()
	{
		int n;
		init();
		while (read(n) && n)
			write(f[n] - (ll)(n + 1) * n / 2), putchar('
');
		return 0;
	}
}
int main()
{
	return zyt::work();
}