洛谷P5325 [模板] Min 25筛 题解
Link
数学部分
约定
$$ p\in\mathbb P,\quad p_i\text{is the}i^\text{th}\text{number}\in\mathbb P $$
第一部分($f(x)$)
min_25 筛的要求是 $f(p^k)$ 可以表示为 $a_0p_0^{b_0}+a_1p_1^{b_1}+\cdots$, 而此题中有:
$$ f(1)=1,\quad\mathop{f(p^k)=p^k(p^k-1)=p^{2k}-p^k}_\limits{k\ge1} $$
试着将 $\sum_{i=1}^nf(i)$ 拆开:
$$ \sum_{i=1}^nf(i){=f(1)+\sum_{2\le p\le n}\sum_{2\le j\le n,\mathrm{lpf}(x)=p}f(j)\\ =1+\sum_{2\le p^e\le n,e\ge1}f(p^e)\sum_{2\le j\le\lfloor\frac n{p^e}\rfloor,\mathrm{lpf}(j)>p}f(j) } $$
发现当 $\sqrt n<p^e\le n$ 时:
$$ \sum_{2\le i\le\lfloor\frac n{p^e}\rfloor,\mathrm{lpf}(i)>p}f(i)=\begin{cases}f(p),& e=1(\forall j\not\in\mathbb P,\mathrm{lcp}(j)>p\Rightarrow j>p^2>n)\\0,&e\not=1(p^e\ge p^2>n,\frac n{p^e}<1)\end{cases} $$
所以可以将柿子化为:
$$ \sum_{i=1}^nf(i)=1+\sum_{2\le p\le\sqrt n,2\le p^e\le n,e\ge1}f(p^e)\sum_{2\le j\le\lfloor\frac n{p^e}\rfloor,\mathrm{lpf}(j)>p}f(j)+\sum_{\sqrt n<p\le n}f(p) $$
第二部分 ($g(x)$)
记 $g_k(n,m)=\sum_{2\le i\le n,i\in\mathbb P~\text{or}~\mathrm{lpf}(i)>p_m}i^k$, 即 $g_k(n,m)$ 表示筛去前 $m$ 个质数后的 $\sum_{i\le n}i^k$; 又记 $g(n,m)=\sum_ka_kg_k(n,m)$. 则:
$$ \sum_{i=2}^nf(i)=g(n,0)=g_2(n,0)-g_1(n,0) $$
$$ g_k(n,m){=g_k(n,m-1)-\sum_{e\ge1}p_m^e(g_k(\lfloor\frac n{p_m}\rfloor,m-1)-g_k(p_{m-1},m-1))\\ =g_k(n,m-1)-\sum_{e\ge1}p_m^e(g_k(\lfloor\frac n{p_m}\rfloor,m-1)-\sum_{1\le i\le m}p_i^k) } $$
记 $g_k(n)=\sum_{2\le p\le n}p^k=g_k(n,\mathop p\limits_{\text{the largest prime}\le\sqrt n})$.
$\lfloor\frac{\lfloor\frac na\rfloor}b\rfloor=\lfloor\frac n{ab}\rfloor$, 所以只需要算出所有的 $g_k(\mathop{\lfloor\frac ni\rfloor}_\limits{1\le i\le\sqrt n})$, 显然 $\lfloor\frac ni\rfloor$ 共有 $\Theta(\sqrt n)$ 个.
记 $S(n,m)=\sum_{2\le i\le n,\mathrm{lpf}(i)>p_m}f(i)$, 则:
$$ S(n,m)=g(n)-g(m)+\sum_{p_i^e\le n,i>m}f(p_i^e)(S(\lfloor\frac n{p_i^e}\rfloor,i)+[e\ge 1]) $$
最终答案即为 $S(n,0)$.
Code
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
using namespace std;
const long long jly=1000000007,inv=333333336/*inv(3)*/;
namespace pr/*线性筛质数*/{
int s; bool f[1000005];
long long p[3][1000005]/*0表示质数,1表示质数前缀和,2表示质数平方前缀和*/;
inline void getp(int x) {
for (int i=2;i<=x;++i) {
if (!f[i]) {
p[0][++s]=i; p[1][s]=(p[1][s-1]+i)%jly;
p[2][s]=(p[2][s-1]+(long long)i*i)%jly;
}
for (int j=1;j<=s&&i*p[0][j]<=x;++j) {
f[i*p[0][j]]=true;
if (!(i%p[0][j])) break;
}
}
}
}
int d[2][1000005]/*0表示一次项,1表示二次项*/;
long long n,q/*sqrt(n)*/,w[1000005]/*n/i的可能取值*/,g[2][1000005];
long long s(long long a,int b) {
if (pr::p[0][b]>=a) return 0;
int c=a>q?d[1][n/a]:d[0][a];
/*求s(n,m)的前半部分*/
long long r=(g[1][c]-g[0][c]+pr::p[1][b]-pr::p[2][b]+jly*2)%jly,x;
for (int i=b+1;i<=pr::s&&pr::p[0][i]*pr::p[0][i]<=a;++i) {
x=pr::p[0][i]; for (int j=1;x<=a;++j,x*=pr::p[0][i])
/*求s(n,m)的后半部分*/
r=(r+x%jly*((x-1)%jly)%jly*(s(a/x,i)+(j==1?0:1)))%jly;
}
return r;
}
int main() {
int t;
scanf("%lld",&n);
pr::getp(q=sqrt(n));
for (long long i=1;i<=n;) /*求g_k(p)*/ {
long long j=(w[++t]=n/i)%jly;
g[0][t]=(j*(j+1)/2-1)%jly;
g[1][t]=(j*(j+1)/2%jly*(j*2+1)%jly*inv-1)%jly;
j=n/(n/i); w[t]>q?(d[1][j]=t):(d[0][n/i]=t);
i=j+1;
}
for (int i=1;i<=pr::s;++i) /*求剩下的g_k(w[])*/
for (int j=1;j<=t&&pr::p[0][i]*pr::p[0][i]<=w[j];++j) {
int k=w[j]/pr::p[0][i]>q?d[1][n/(w[j]/pr::p[0][i])]:d[0][w[j]/pr::p[0][i]];
g[0][j]=(g[0][j]-pr::p[0][i]*(g[0][k]-pr::p[1][i-1]+jly)%jly+jly)%jly;
g[1][j]=(g[1][j]-pr::p[0][i]*pr::p[0][i]%jly*(g[1][k]-pr::p[2][i-1]+jly)%jly+jly)%jly;
}
printf("%lld",(s(n,0)+1)%jly); return 0;
}