首页 > 技术文章 > BZOJ4916 神犇和蒟蒻 【欧拉函数 + 杜教筛】

Mychael 原文

题目

很久很久以前,有一只神犇叫yzy;
很久很久之后,有一只蒟蒻叫lty;

输入格式

请你读入一个整数N;1<=N<=1E9,A、B模1E9+7;

输出格式

请你输出一个整数A=sum_{i=1}^N{mu (i^2)};
请你输出一个整数B=sum_{i=1}^N{varphi (i^2)};

输入样例

1

输出样例

1

1

题解

首先很明显= =

[ans1 = 1 ]

然后重点是(ans2)
我们会发现这样一个性质:

[varphi(i^2) = i*varphi(i) ]

证明:
(i^2)(i)拥有相同的质因子,所以与(i)互质的同样与(i^2)互质,与(i)不互质的同样与(i^2)不互质
由于(gcd(i,j) = gcd(i + i,j)),且(i^2 = i * i),所以(varphi(i^2) = i * varphi(i))
【ps:突然发现我傻逼了,直接上(varphi)的计算公式就得证了QAQ】

然后就可以筛了:
我们令

[f(i) = i * varphi(i) ]

拿出杜教筛的套式:

[g(1) * S(n) = sumlimits_{i=1}^{n}(f * g)(i) - sumlimits_{i=2}^{n} g(i) * S(lfloor frac{n}{i} floor) ]

我们企图找出这样一个好算的(g(i))
考虑卷积:

[(f * g)(n) = sumlimits_{d|n} g(frac{n}{d}) * d * varphi(d) ]

我们想把(varphi(d))外的变为常数,这样可以提出来,成为(varphi * 1 = Id)
很容易发现,如果(g = Id),那么就刚好可以做到:

[(Id * f)(n) = sumlimits_{d|n} frac{n}{d} * d * varphi(d) = n * sumlimits_{d|n} varphi(d) = n^2 ]

再由:

[sumlimits_{i=1}^{n} i^2 = frac{n * (n + 1) * (2 * n + 1)}{6} ]

那么杜教筛的式子就化成了:

[S(n) = frac{n * (n + 1) * (2 * n + 1)}{6} - sumlimits_{i=2}^{n} i * S(lfloor frac{n}{i} floor) ]

就可以做了

#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
#include<cstring>
#include<algorithm>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<' '; puts("");
using namespace std;
const int maxn = 1000005,maxm = 100005,INF = 1000000000,P = 1000000007;
typedef map<LL,LL> Map;
Map _phi;
Map::iterator it;
LL p[maxn],pi,N,n,phi[maxn],v6;
int isn[maxn];
LL qpow(LL a,LL b){
	LL ans = 1;
	for (; b; b >>= 1,a = a * a % P)
		if (b & 1) ans = ans * a % P;
	return ans;
}
void init(){
	v6 = qpow(6,P - 2);
	N = (LL)pow(n,2.0 / 3.0);
	phi[1] = 1;
	for (LL i = 2; i < N; i++){
		if (!isn[i]) p[++pi] = i,phi[i] = i - 1;
		for (int j = 1; j <= pi && i * p[j] < N; j++){
			isn[i * p[j]] = true;
			if (i % p[j] == 0){
				phi[i * p[j]] = phi[i] * p[j] % P;
				break;
			}
			phi[i * p[j]] = phi[i] * (p[j] - 1) % P;
		}
	}
	for (LL i = 1; i < N; i++) phi[i] = (i * phi[i] % P + phi[i - 1]) % P;
}
LL S(LL n){
	if (n < N) return phi[n];
	if ((it = _phi.find(n)) != _phi.end())
		return it->second;
	LL ans = n * (n + 1) % P * (2 * n + 1) % P * v6 % P;
	for (LL i = 2,nxt; i <= n; i = nxt + 1){
		nxt = n / (n / i);
		ans = (ans - (nxt + i) * (nxt - i + 1) / 2 % P * S(n / i) % P) % P;
	}
	ans = (ans % P + P) % P;
	return _phi[n] = ans;
}
int main(){
	cin >> n;
	init();
	cout << 1 << endl << S(n) << endl;
	return 0;
}

推荐阅读