先定义几个符号：
[]：若方括号内为一个值，则向下取整，否则为布尔判断
集合P：素数集合。
题目分析：
题目是一个积性函数。做法之一是洲阁筛，也可以采用Min_25筛。
对于一个可以进行Min_25筛法的积性函数，它需要满足与洲阁筛相同的条件，即：
对于$f(p), p in P$，它可以多项式表出。对于$f(p^k),p in P$可以被快速计算出。
这道题中$f(p) = p-1$再对$2$进行修正即可。
对于1的情况我们单独考虑，现在我们对答案进行一些变换。
$$sum_{i=2}^{n}f(i) = sum_{2<=p^e<=n, p in P} f(p^e)(1+sum_{2<=i<=[frac{n}{p^e}],i质因子大于p}f(i)) = sum_{2<=p<=sqrt n,p^e<=n, p in P} f(p^e)(1+sum_{2<=i<=[frac{n}{p^e}],i质因子大于p}f(i))+sum _{p in P,sqrt n<p<=n}f(p)$$
我们注意到中间刺眼的中文，不妨以此为根据设置dp状态。
令$$g_{n,m}=sum_{i=2,i的质因子大于m}^{n}f(i)$$
它的递推式不难写出来。根据上面那串式子:
$$g_{n,m}=sum_{i=m,i in P}^{sqrt n}[p^e leq n]([e eq 1]+sum_{i=p+1,i的质因子大于p}^{[frac{n}{p^e}]}f(i))+sum_{i=m+1,i in P}^{n}f(i) = sum_{i=m,i in P}^{sqrt n}[p^e leq n]([e eq 1]+g_{[frac{n}{p^e}],p})+sum_{i=m+1,i in P}^{n}f(i)$$ 之后我们可以发现后面那一部分等于$m+1$到$n$的所有素数.
朱老大的论文似乎在这个式子上有一些问题，这个式子是我加了修正的，可以对比原式找到区别。
现在我们重点考虑后面的部分，质数前缀和的求解。 我们可以用DP模拟埃拉托斯特尼筛法来完成。现在对f(p)这个多项式的每个单项式分开考虑。
令$$h_{n,m}=sum_{i=2,i为素数或i质因子大于m}f(i)$$
由埃氏筛可知$$h_{n,m}=h_{n,m-1}-p_m^s*(h_{[frac{n}{p}],m-1}-h_{p_m,m-1})$$起始状态将所有数看做质数，做k次幂和。
由于n以下的合数至多的因子大小是$sqrt n$,所以对于一个n只用筛前$sqrt n$个素数。
g数组同理。
这样我们就完成整个过程了，现在我们来思考时间复杂度上的问题。这一部分我的证明存在争议，与论文上的有出入。希望有人可以帮助我确认出入之处。
先给出两个引理吧。
素数定理:$pi(x) approx frac{x}{lnx} approx li(x)$
数论分块:对于一个数n，它对1到n整除的结果只有$2sqrt n$个，分别是1到$sqrt n$与$n/sqrt n$到$n/1$
首先是dp模拟埃氏筛。对于数论分块中的每一个数，每个小于它的根号的素数都被考虑，所以时间为，
$$O(sum_{i=2}^{sqrt n}frac{sqrt i}{lni})+O(sum_{i=1}^{sqrt n}frac{sqrt frac{n}{i}}{ln{frac{n}{i}}})$$
对于前面的大O,将$sqrt i$放缩到$sqrt n$,$lni$放缩到$lnn$,因为根号的增长远大于对数，所以放缩成立。前面被证明是$O(frac{n^{0.75}}{lnn})$的。
后面的比较复杂，将它看做积分，即$$O(int_1^{sqrt n}frac{sqrt frac{n}{x}}{lnfrac{n}{x}}dx)$$
由于对数函数增长极慢，不妨将分母用同样的放缩。当然也可以查积分表得到与我相同的结论。这样再对单项式函数做积分就简单了，得到$$O(frac{sqrt n}{lnn}*[2sqrt x]_1^{sqrt n})=O(frac{n^{0.75}}{lnn})$$。
埃氏筛的复杂度得证，我们再来看看Min_25筛的表现。论文中通过证明得出一个与运行效率不相符的时间复杂度，我试着用其它方法得到它的复杂度。
观察Min_25的特点，它与DP埃氏筛的区别在于多枚举了质数的次幂。这样我们不妨对于质数的次幂分开考虑。
首先是质数本身，如果只枚举它，时间与埃氏筛DP相同。实际上它代表的是二次方，否则不被考虑。
然后是质数的平方，它要小于等于n。接着是质数的三次方，它也要小于等于n，然后是四次方，五次方。我们对第x次方小于等于n做一些转化使得他们统一。
$p^x leq n leftrightarrow p leq n^{1/x} leftrightarrow p^2 leq n^{2/x}$
这时候我们的复杂度相当于对$n^{2/2},n^{2/3},n^{2/4}$做埃氏筛DP。每个的时间在上面分析了，所以下面直接转化成积分的形式。
$$O(int_{2}^{log_{2}n}frac{n^{frac{3}{2x}}}{ln{n^{frac{2}{x}}}}dx)$$
$$=O(frac{1}{lnn}int_{2}^{log_{2}n}frac{n^{frac{3}{2x}}}{frac{2}{x}}dx)$$
$$=O(frac{1}{2lnn}int_{2}^{log_{2}n}{x*n^{frac{3}{2x}}}dx)$$
$$=O(frac{1}{2lnn}int_{2}^{log_{2}n}{x*a^{frac{1}{x}}}dx),a = n^{frac{3}{2}}$$
接下俩数学好的可以换元，令$y=frac{1}{x}$,再令$z = a^y$，求$$int frac{1}{ln^3{z}dz}$$即可。
现在跳过这些步骤，得到结论是
$$O([frac{x(x+ln a)a^frac{1}{x}-ln^2 aoperatorname{li}(a^{frac{1}{x}})}{4lnn}]_{2}^{log_{2}n})$$
这个答案是负的，证明我可能某个步骤写错了。但是代入2可以得到时间是$O(n^{frac{3}{4}})$，递归的过程不会重复计算。
继续写下去我也写不动了。到此打止。
代码：
1 #include<bits/stdc++.h> 2 using namespace std; 3 4 const int SQR = 100000; 5 const int mod = 1000000007; 6 7 long long n; 8 9 int prime[SQR],flag[SQR+5],num; 10 long long h[SQR*2+5],sqr; 11 long long fh[SQR*2+5]; 12 13 void getprime(int N){ 14 for(int i=2;i<=N;i++){ 15 if(!flag[i]) prime[++num] = i; 16 for(int j=1;j<=num&&i*prime[j] <= N;j++){ 17 flag[i*prime[j]] = 1; 18 if(i%prime[j] == 0) break; 19 } 20 } 21 } 22 23 void geth(int ft){ 24 for(int i=1;i<=num;i++){int zt = (ft > 0?prime[i]:1); 25 for(int j=1;j<=sqr&&(1ll*prime[i]*prime[i]<=(n/j));j++){ 26 long long nxt = 1ll*j*prime[i]; int plu = 2*sqr-j+1; 27 if(nxt > sqr) h[plu] -= 1ll*zt*(h[n/nxt]-h[prime[i]-1])%mod; 28 else h[plu] -= 1ll*zt*(h[2*sqr-nxt+1]-h[prime[i]-1])%mod; 29 h[plu]<0?h[plu]+=mod:1; 30 } 31 for(int j=sqr;j>=2&&(1ll*prime[i]*prime[i]<=j);j--){ 32 h[j] -= 1ll*zt*(h[j/prime[i]]-h[prime[i]-1])%mod; 33 h[j]<0?h[j]+=mod:1; 34 } 35 } 36 } 37 38 int f(int now,int zt){ 39 if(now <= sqr){ 40 if(now <= 1 || prime[zt] >= now) return 0; 41 int ans = fh[now]-fh[prime[zt]]; if(ans < 0) ans += mod; 42 if(1ll*prime[zt]*prime[zt] >= now) return ans; 43 for(int i=zt+1;1ll*prime[i]*prime[i] <= now;i++){ 44 long long dt = prime[i],k=1; 45 while(dt <= now){ 46 ans += ((prime[i]^k)*(f(now/dt,i)+1))%mod; ans %= mod; 47 dt = dt*prime[i];k++; 48 } 49 ans-=(prime[i]^1); if(ans < 0) ans += mod; 50 } 51 return ans; 52 }else{ 53 long long rd = n/(2*sqr-now+1); 54 int ans = fh[now]-fh[prime[zt]]; if(ans < 0) ans += mod; 55 if(1ll*prime[zt]*prime[zt] >= rd) return ans; 56 for(int i=zt+1;1ll*prime[i]*prime[i] <= rd&&i<=num;i++){ 57 long long dt = prime[i],k=1; 58 while(dt <= rd){ 59 if(rd/dt <= sqr) ans += ((prime[i]^k)*(f(rd/dt,i)+1))%mod; 60 else ans += ((prime[i]^k)*(f(2*sqr-(2*sqr-now+1)*dt+1,i)+1))%mod; 61 ans %= mod; dt = dt*prime[i];k++; 62 } 63 ans-=(prime[i]^1); if(ans < 0) ans += mod; 64 } 65 return ans; 66 } 67 } 68 69 void work(){ 70 for(int i=2;i<=sqr;i++){h[i] = (h[i-1] + i)%mod;} 71 for(int i=1;i<=sqr;i++){ 72 long long fak=(n/i)%mod;h[sqr*2-i+1]=((2+fak)*(fak-1)/2)%mod; 73 } 74 geth(1); for(int i=2;i<=2*sqr;i++) fh[i] += h[i]; 75 for(int i=2;i<=sqr;i++) h[i] = i-1; for(int i=1;i<=sqr;i++) h[2*sqr-i+1]=(n/i-1)%mod; 76 geth(0); for(int i=2;i<=2*sqr;i++) fh[i] = fh[i]-h[i]+2,fh[i] =(fh[i]+mod)%mod; 77 printf("%d",f(2*sqr,0)+1); 78 } 79 80 int main(){ 81 scanf("%lld",&n); sqr = sqrt(n); 82 if(n == 1){puts("1");return 0;} 83 getprime(sqr);work(); 84 return 0; 85 }