洛谷P5325 [模板] Min 25筛 题解

Jun 08th, 2020
  • 在其它设备中阅读本文章

Link

洛谷 P5325 [模板] Min 25 筛

数学部分

约定

$$ 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;
}

owo

mo-ha