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

$$\sum_{d|i} phi(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的约数为一个大小为K_m的集合S_m$
$也就是phi(n)=\sum_{i=1,i<=K_m} phi(p)*phi(S_i)+phi(m)$
$=\sum_{i=1,i<=K_m} phi(p*S_i)+\sum_{i=1,i<=K_m}phi(S_i)$
$右边的部分的原有的m的所有约数的phi的和，$
$而加入一个素数p，新产生的约数就是所有的p*S_i$
$答案加上了新约数的phi的和，故n的时候也满足，得证$

现在令$M[n]=\sum phi[i]$
对于上面得出的式子进行变形：
$$\sum_{d|n} phi(d)=n\\ =\\ \sum_{i=1,i<=n}\sum_{d|i}phi(d)=n*(n+1)/2\\ 令i=kd，枚举k\\ =\\ \sum_{k=1,k<=n}\sum_{d=1,d<=n/k}phi(d)\\ =\\ \sum_{k=1,k<=n}M[n/k]=n*(n+1)/2$$
于是做个差，就得到了：
$$M[n]=n*(n+1)/2-\sum_{k=2,k<=n}M[n/k]$$
注意n/k下取整。
其实推到中由枚举d转化到枚举1~n/k是挺巧的……
可以看作对于k，1~n/k的每个约数次数肯定都+1了。
然后这个式子我们该怎么算呢？
这个式子很好地构成了一个”n\k”项，
而我们知道，n/k是只有O($\sqrt n$)种取值的，
也就是说，先枚举k=1~ $\sqrt n$，求出这$\sqrt n$种取值，
然后枚举n/k=i，i=1~$\sqrt n$，
因为k更大的时候，答案种数不超过$\sqrt n$个。
这就可以递归计算了，然后用一个记忆化记录，
先预处理出O($n^{2/3}$)内的M[]，
由于只有$\sqrt n$种，所以记忆化的存储量不会很大。
记忆化我就用了个map来存了……（数组开不了）
时间复杂度怪怪的不会证。。说是O($n^{2/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[n];
    if (f[n]) return f[n];
    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[n]=((ll)n*(n+1)>>1LL)-t;
    return f[n];
}
int main(){
    Get_Phi();
    int n;scanf("%d",&n);
    printf("%lld\n",solve(n));
    return 0;
}