bzoj 4805 欧拉函数求和 杜教筛

Description

给出一个数字N，求sigma(phi(i)),1<=i<=N

Input

正整数N。N<=2*10^9

Output

输出答案。

Sample Input

10

Sample Output

32

HINT

传送门

太神了我竟然才会杜教筛……

大致写写推导过程吧（虽然都是按照网上的来的QAQ）

设phi(i)表示i的欧拉函数，那么有如下式子：

∑d|iphi(d)=i

这个东西怎么来的呢。。可以简单证明一下（可能不太严谨）：

显然phi(i)满足，1<=i<=4，

令n>4，且n=p∗m，p为一个质数，且phi(m)满足该性质

那么phi(n)=phi(m)∗p=phi(m)∗(phi(p)+1)

设m的约数为一个大小为Km的集合Sm

也就是phi(n)=∑i=1,i<=Kmphi(p)∗phi(Si)+phi(m)

=∑i=1,i<=Kmphi(p∗Si)+∑i=1,i<=Kmphi(Si)

右边的部分的原有的m的所有约数的phi的和，

而加入一个素数p，新产生的约数就是所有的p∗Si

答案加上了新约数的phi的和，故n的时候也满足，得证

现在令M[n]=∑phi[i]

对于上面得出的式子进行变形：

∑d|nphi(d)=n=∑i=1,i<=n∑d|iphi(d)=n∗(n+1)/2令i=kd，枚举k=∑k=1,k<=n∑d=1,d<=n/kphi(d)=∑k=1,k<=nM[n/k]=n∗(n+1)/2

于是做个差，就得到了：

M[n]=n∗(n+1)/2−∑k=2,k<=nM[n/k]

注意n/k下取整。

其实推到中由枚举d转化到枚举1~n/k是挺巧的……

可以看作对于k，1~n/k的每个约数次数肯定都+1了。

然后这个式子我们该怎么算呢？

这个式子很好地构成了一个”n\k”项，

而我们知道，n/k是只有O(n−√)种取值的，

也就是说，先枚举k=1~ n−√，求出这n−√种取值，

然后枚举n/k=i，i=1~n−√，

因为k更大的时候，答案种数不超过n−√个。

这就可以递归计算了，然后用一个记忆化记录，

先预处理出O(n2/3)内的M[]，

由于只有n−√种，所以记忆化的存储量不会很大。

记忆化我就用了个map来存了……（数组开不了）

时间复杂度怪怪的不会证。。说是O(n2/3)的。。

用了map的话……再乘个log吧。。。

神啊杜教筛……

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int
maxn=1600000;
ll M[maxn];
map<int,ll>f;
int pcnt,phi[maxn],prime[maxn/10];
bool notprime[maxn];
void Get_Phi(){
phi[1]=M[1]=1,pcnt=0;
notprime[1]=1;
for (int i=2;i<maxn;i++){
if (!notprime[i]) phi[i]=i-1,prime[++pcnt]=i;
for (int j=1;j<=pcnt;j++){
if (prime[j]*i>=maxn) break;
notprime[i*prime[j]]=1;
if (i%prime[j]==0){
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
}
for (int i=2;i<maxn;i++) M[i]=M[i-1]+phi[i];
}
ll solve(int n){
if (n<maxn) return M
;
if (f
) return f
;
int tmp=sqrt(n);ll t=0LL;
for (int i=2;i<=tmp;i++) t+=solve(n/i);
tmp=n/(tmp+1);
for (int i=1;i<=tmp;i++) t+=solve(i)*(n/i-n/(i+1));
f
=((ll)n*(n+1)>>1LL)-t;
return f
;
}
int main(){
Get_Phi();
int n;scanf("%d",&n);
printf("%lld\n",solve(n));
return 0;
}