神犇和蒟蒻

题目

传送门 to VJ

题目描述
给定 $n(1\le n\le 10^9)$ ，求下列两个值取模 $10^9+7$ 后的结果。

$$\sum _{i=1}^n\mu (i^2),\sum _{i=1}^n\phi (i^2)$$

思路

对于第一个，直接输出 $1$ 就好了……定义告诉我们，指数为 $2$ 就会凉。只有 $\mu (1^2)=1$ 可以挽回一些局面。

对于第二个，不难发现 $\phi (i^2)=i\cdot \phi (i)$ 。证明方法很多，我给两个比较好理解的，否则这篇博客就太短啦。

等价转化

$$\begin{array}{rl}\phi (x^2)& =\sum _{i=1}^{x^2}[\gcd (i,x^2)=1]\\ & =\sum _{i=1}^{x^2}[\gcd (i,x)=1]\\ & =\sum _{i=1}^{x^2}[\gcd (i\mathrm{m}\mathrm{o}\mathrm{d}x,x)=1]\\ & =\sum _{d=1}^x[\gcd (d,x)=1]\sum _{i\equiv d(\mathrm{m}\mathrm{o}\mathrm{d}x)}1\\ & =\sum _{d=1}^x[\gcd (d,x)=1]\cdot x\\ & =x\cdot \phi (x)\end{array}$$

狄利克雷

$$\begin{array}{rl}\phi (x^2)& =\sum _{d|x^2}\mu (d)\cdot \frac{x^2}{d}\\ & =x\sum _{d|x}\mu (d)\cdot \frac{x}{d}\\ & =x\cdot \phi (x)\end{array}$$

用上 $\phi =\mu \ast I_d$ 就非常好推呀！并且这个方向也比较自然。

当然，$d|x^2$ 变为 $d|x$ 也是比较重要的操作。

最后一步，直接杜教筛即可。用 $I_d$ 去筛。

代码

#include <cstdio>
#include <iostream>
#include <vector>
#include <cstring>
using namespace std;
typedef long long int_;
inline int readint(){
	int a = 0; char c = getchar(), f = 1;
	for(; c<'0'||c>'9'; c=getchar())
		if(c == '-') f = -f;
	for(; '0'<=c&&c<='9'; c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}

const int MaxN = 5e6;
const int Mod = 1e9+7;
const int inv2 = (Mod+1)>>1;
const int inv3 = (Mod+1)/3;
vector< int > primes;
bool isPrime[MaxN];
int phi[MaxN], f[MaxN];
void sievePrime(){
	memset(isPrime+2,1,MaxN-2);
	f[1] = phi[1] = 1;
	for(int i=2,len=0; i<MaxN; ++i){
		if(isPrime[i]){
			primes.push_back(i);
			++ len, phi[i] = i-1;
		}
		f[i] = 1ll*i*phi[i]%Mod;
		for(int j=0,x; j<len; ++j){
			if((x = primes[j]*i) >= MaxN)
				break; // good try
			isPrime[x] = false;
			if(i%primes[j] == 0){
				phi[x] = phi[i]*primes[j];
				break;
			}
			phi[x] = phi[i]*(primes[j]-1);
		}
	}
	for(int i=2; i<MaxN; ++i)
		f[i] = (f[i]+f[i-1])%Mod;
}

const int SqrtN = 200005;
int haxi[2][SqrtN];
# define index_(x) (((x) >= SqrtN) ? haxi[1][n/(x)] : haxi[0][x])
int w[SqrtN<<1]; // f(xez[x])
int_ xez[SqrtN<<1]; // from haxi to real
int djs(int_ n){ // Dc.Du sieve
	int id = 0; // allocate index
	for(int_ i=1; i<=n; i=n/(n/i)+1){
		index_(n/i) = ++ id;
		xez[id] = n/i; // real value
	}
	int_ x; // for convenience
	for(int i=id; i>=1; --i){
		if((x = xez[i]) < MaxN){
			w[i] = f[x]; continue;
		}
		w[i] = ((x<<1)+1)*x%Mod
			*(x+1)%Mod*inv2%Mod
			*inv3%Mod; // sum x^2
		int lst = 1, now; // init 1
		for(int_ l=2,r; l<=x; l=r+1){
			r = x/(x/l);
			now = (r*(r+1)>>1)%Mod;
			w[i] -= 1ll*(now-lst)
				*w[index_(x/l)]%Mod;
			w[i] = (w[i]%Mod+Mod)%Mod;
			lst = now; // move on
		}
	}
//	for(int i=id; i>=1; --i)
//		printf("f[%lld] = %d\n",xez[i],w[i]);
	return w[1]; // maximum = n
}

int main(){
	int_ n; scanf("%lld",&n);
	sievePrime();
	printf("1\n%d\n",djs(n));
	return 0;
}