Min_25筛目害学

带着理解抄了一遍模板题解来加深印象

求解积性函数前缀和一般可以使用杜教筛。但当这种积性函数不常见且难于卷出已于求前缀和的函数时,杜教筛是乏力的。这时可以考虑使用Min_25筛。

以下内容中约定 \(p\in\mathbb{P}\)\(p_i\) 为第 \(i\) 大的质数。

使用Min_25筛的条件:该函数(不妨令其为 \(F\) ) 是积性函数,且 \(F(p)\) 为关于 \(p\) 的低阶多项式, \(F(p^k)\) 可快速求值。

低次多项式的限制是因为Min_25筛主要用来解决求 \(f(p)=p^k\) 的前缀和,把多项式拆成若干单项式,分别求出前缀和后相加就能得到答案。

Min_25筛求解主要分两个步骤。

Step1 分类

我们按质数与合数将 \(\sum_{i=1}^nf(i)\) 分类。

\[\sum_{i=1}^nf(i)=\sum_{p\leq n}f(p)+\sum_{i=1}^nf(i)[i\notin\mathbb{P}] \]

对后一部分枚举最小质因子来满足限制。设 \(lpf_i\)\(i\) 的最小质因子,

\[\sum_{i=1}^nf(i)=\sum_{p\leq n}f(p)+\sum_{p\leq n}\sum_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor}f(i)[lpf_i\geq p] \]

于是对两部分分别求解。

Step2 求解

第一部分

线性筛出所有 \(f(p)\) 肯定是不现实的。我们定义 \(g_{n,i}\)\(n\) 以内所有质数与最小质因子大于 \(p_i\) 的数的 \(k\) 次方和,我们要求的东西就是 \(g(x,\sqrt x)\),考虑DP求这个东西。

\(g(n,0)\) 可以 \(O(1)\) 求,考虑如何从 \(g(n,j-1)\) 转移到 \(g(n,j)\) 。这步操作去掉的部分为最小质因子为 \(p_j\) 且不大于 \(n\) 的合数 \(k\) 次方值。

\[g(n,j)=g(n,j-1)-p_j^k(g(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1)-g(p_{j-1},j-1)) \]

后一个 \(g\) 去掉了所有小于 \(p_j\) 的质数。

实际上我们需要的只是 \(g(x,p_y)\) ,其中 \(p_y\) 为最后一个小于 \(x\) 的质数。称它为 \(pw_x\)

求解答案

继续设 \(S_{n,x}\)\([1,n]\) 内最小质因子大于 \(p_x\) 的数函数值之和。继续分类考虑转移。第一部分是大于 \(p_x\) 的质数,可写作 \(pw_n-pw_x\) 。另一部分是剩余合数,可以递归求解。

\[S_{n,x}=pw_n-pw_x+\sum_{p_k^e\leq n\wedge k>x}f(p_k^e)(S(\left\lfloor\frac{n}{p_k^e}\right\rfloor,k)+[e>1]) \]

然后发现我们实际上只需要 \(\frac{n}{x}\) 处的 \(pw\) 值,因为合数的最小质因子不超过根号,所以 \(pw_x\) 都可以直接预处理出来。于是求解 \(g\) 时只需求 \(O(\sqrt n)\) 个值。

细节

离散化时可以开两个数组,分别记录 \(\frac{n}{x}\)\(\frac{n}{\frac{n}{x}}\) 离散化后的位置,存的时候哪个小就存哪个,可以保证数组大小都是 \(O(\sqrt n)\) 的。

1的贡献单算,最后加一。


\(code:\)

P5325[模板]Min_25筛
#include<bits/stdc++.h>
using namespace std;

namespace IO{
    #define int LL
    typedef long long LL;
    int read(){
        int x=0,f=0; char ch=getchar();
        while(ch<'0'||ch>'9'){ f|=ch=='-'; ch=getchar(); }
        while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
        return f?-x:x;
    } char outp[50];
    void write(int x,char sp){
        int len=0;
        if(x<0) x=-x, putchar('-');
        do{ outp[len++]=x%10+'0'; x/=10; }while(x);
        for(int i=len-1;~i;i--) putchar(outp[i]); putchar(sp);
    }
} using namespace IO;

const int NN=300010,mod=1e9+7,inv3=333333336;
int n,sqr;
void pls(int& x,int y){ (x+=y)>=mod?(x-=mod):(x<0?(x+=mod):0); }

namespace Min_25{
    int cnt,pri[NN],p1[NN],p2[NN];
    int tot,g1[NN],g2[NN],has[NN],pos1[NN],pos2[NN];
    bool vis[NN];
    void sieve(int n){
        for(int i=2;i<=n;i++){
            if(!vis[i]){ pri[++cnt]=i; pls(p1[cnt]=p1[cnt-1],i); pls(p2[cnt]=p2[cnt-1],i*i%mod); }
            for(int j=1;j<=cnt&&i*pri[j]<=n;j++){
                vis[i*pri[j]]=1;
                if(i%pri[j]==0) break;
            }
        }
    }
    void divblock(){
        for(int l=1,r;l<=n;l=r+1){
            r=n/(n/l);
            has[++tot]=n/l;
            g1[tot]=has[tot]%mod;
            pls(g2[tot]=g1[tot]*(g1[tot]+1)/2%mod*(g1[tot]*2+1)%mod*inv3%mod,-1);
            pls(g1[tot]=g1[tot]*(g1[tot]+1)/2%mod,-1);
            n/l<=sqr?pos1[n/l]=tot:pos2[r]=tot;
        }
    }
    void getg(){
        for(int i=1;i<=cnt;i++)
            for(int j=1;j<=tot&&pri[i]*pri[i]<=has[j];j++){
                int k=has[j]/pri[i]<=sqr?pos1[has[j]/pri[i]]:pos2[n/(has[j]/pri[i])];
                pls(g1[j],-pri[i]*(g1[k]-p1[i-1])%mod);
                pls(g2[j],-pri[i]*pri[i]%mod*(g2[k]-p2[i-1])%mod);
            }
    }
    int calc(int x,int y,int res=0){
        if(pri[y]>=x) return 0;
        int k=x<=sqr?pos1[x]:pos2[n/x];
        pls(res,(g2[k]-g1[k]+p1[y]-p2[y])%mod);
        for(int i=y+1;i<=cnt&&pri[i]*pri[i]<=x;i++)
            for(int j=1,p=pri[i];p<=x;j++,p=p*pri[i]){
                int pp=p%mod;
                pls(res,pp*(pp-1)%mod*(calc(x/p,i)+(j!=1))%mod);
            }
        return res;
    }
} using namespace Min_25;

signed main(){
    n=read(); sqr=sqrt(n);
    sieve(sqr); divblock(); getg();
    write((calc(n,0)+1)%mod,'\n');
    return 0;
}
posted @ 2022-01-16 19:21  keen_z  阅读(43)  评论(0编辑  收藏  举报