SPOJ

目录

• 前言
• 题目
• 解法
  • $ ext{Part 1}$
  • $ ext{Part 2}$
• 总结
• 代码

前言

这样的函数 可以用线性筛，于是我就自己推了推，发现自己有好多问题没学透...

最后好不容易过了，但是常数巨大，实际上还不如 (mathcal O(nln n))...

题目

传送门

解法

( ext{Part 1})

[ ext{Ans}=n imes sum_{i=1}^n frac{i}{gcd(i,n)} ]

(i) 挺难搞的，但由于 (gcd(n,i)=gcd(n,n-i))，我们可以把 (i) 两两分组都凑成 (n)。

然后我试了试这个式子（这个式子到最后推不动，如果不想看的可以跳过）：

[=left(n imes sum_{i=1}^{(n-1)/2}frac{n}{gcd(i,n)} ight )+n+[nmod 2=0] imes n ]

具体就是 (n) 没有数进行配对，如果 (n) 是偶数那么 (frac{n}{2}) 也没有数配对，所以在后面再算。

然后想到枚举 (gcd)。

[=left(n^2 imes sum_{g=1}^{(n-1)/2}frac{1}{g}sum_{i=1}^{(n-1)/2g}[gcd(frac{n}{g},i)=1] ight )+n+[nmod 2=0] imes n ]

推到这里发现不能用欧拉函数搞，因为 ((n-1)/2g) 这个范围和 (frac{n}{g}) 不同。

发现了这个方法有问题的原因后，我们不妨换一种思路：将 (frac{i}{gcd(i,n)}) 看成 (frac{n}{gcd(i,n)})，这样需要将答案除以 (2)，但是我们发现 (i) 的值域可以达到 (n) 了！

[=left(n imes frac{1}{2}sum_{i=1}^{n}frac{n}{gcd(i,n)} ight )+n ]

有一个小细节，(frac{1}{2}sum_{i=1}^{n}frac{n}{gcd(i,n)}) 并没有算第 (n) 项，因为 (sum_{i=1}^{n}frac{n}{gcd(i,n)}) 是个奇数，而那个 (1) 是第 (n) 项贡献的，所以在除二的时候被除掉了。最后我又加回去了。

[=n+frac{n^2}{2}sum_{d|n}frac{1}{d} imes varphileft( frac{n}{d} ight) ]

[=n+frac{n}{2}sum_{d|n}d imes varphi(d) ]

在这个式子里，第 (n) 项就是 (d=1) 的情况，此时 (varphi(1)=0)，所以也没有算第 (n) 项。

( ext{Part 2})

然后我们发现 (g(i)=i imes varphi(i)) 肯定是个积性函数，只不过我们要求的是 (f(n)=sum_{d|n}d imes varphi(d))。此外我们还需要定义辅助函数 (pk(i))，假设 (i) 中最小质因子 (p) 在 (i) 中的幂次为 (k)，那么 (pk(i)=p^k)。

1. (i) 为质数。(f(i)=i(i-1),pk(i)=i)。

2. (i) 与 (p) 互质。

   考虑相对 (i)，(ip) 的因数增加了什么 —— 容易发现它们都是 (p) 的倍数。那么就有：

   [f(ip)=f(i)+sum_{kp|ip}g(kp) ]

   (k) 和 (p) 一定是互质的，可以用反证法，就不赘述。于是就有：

   [f(ip)=f(i)+(f(i)+1) imes g(p) ]
   [pk(ip)=p ]

   加一是因为 (f(i)) 没有计算约数为 (1) 的情况，但这在 (f(ip)) 里对应着 (p)。我就被这个坑了很久。

3. (i) 与 (p) 不互质。

   将 (i) 拆分成 (j imes p^k)。类似地，考虑增加的因数都是 (p^{k+1}) 的倍数。于是有：

   [f(ip)=f(i)+(f(j)+1) imes g(p^{k+1}) ]
   [pk(ip)=pk(i) imes p ]

   加一的理由同上。

总结

除了单纯的积性函数之外，应该还可以用线性筛搞搞积性函数的 “前缀和”。

我倒戈了，暴力加它不香吗？

代码

#include <cstdio>

#define print(x,y) write(x),putchar(y) 

template <class T> inline T read(const T sample) {
    T x=0; int f=1; char s;
    while((s=getchar())>'9'||s<'0') if(s=='-') f=-1;
    while(s>='0'&&s<='9') x=(x<<1)+(x<<3)+(s^48),s=getchar();
    return x*f;
}
template <class T> inline void write(const T x) {
    if(x<0) return (void) (putchar('-'),write(-x));
    if(x>9) write(x/10);
    putchar(x%10^48);
}

const int maxn=1e6+5;

int n,p[maxn],pc;
long long f[maxn],pk[maxn];
bool is[maxn];

void init() {
	f[1]=0;
	for(int i=2;i<=maxn-5;++i) {
		if(!is[i]) p[++pc]=i,f[i]=1ll*i*(i-1),pk[i]=i;
		for(int j=1;j<=pc and i*p[j]<=maxn-5;++j) {
			is[i*p[j]]=1;
			if(i%p[j]==0) {
				f[i*p[j]]=f[i]+(f[i/pk[i]]+1)*(pk[i]*p[j]-pk[i])*pk[i]*p[j];
				pk[i*p[j]]=pk[i]*p[j]; 
				break;
			}
			f[i*p[j]]=f[i]+(f[i]+1)*(p[j]-1)*p[j];
			pk[i*p[j]]=p[j];
		}
	}
}

signed main() {
	init();
	for(int T=read(9);T;--T) {
		int n=read(9);
		print(f[n]*n/2+n,'
');
	}
	return 0;
}