洲阁筛学习笔记

引言

为啥机房人均会 \(\rm min\_25\) 筛！

就我只会杜教筛和 PN 筛吗……

这样甚至不能解决某些稍微困难的积性函数前缀和。

\(\rm min\_25\) 还是太复杂，学不动，所以学个洲阁筛好了。常数似乎会劣一点，但是复杂度差不多。

参考了 2016 年候选队论文《积性函数求和的几种方法》，即任之洲洲阁筛原始文献。

第二部分采用了论文中“另一种实现”的方法。这样子显然更好推了吧（

值得注意的是，洲阁筛是块筛！

算法流程

贡献拆分

设 \(P\) 为质数集。

设 \(f\) 为一积性函数，使得 \(f(p)\) 为关于 \(p\) 的低次多项式，\(f(p^c)\) 可以快速计算，其中 \(p\in P,c\in{\mathbb N}^*\)。

现在要求

\[S(n)=\sum_{i=1}^nf(i) \]

则由于 \(n\) 以内的数最多只有一个大于 \(\sqrt n\) 的质因子，我们分类讨论之，得到

\[S(n)=\sum_{i=1}^n[\exist p>\sqrt n\land p\in P,p|i]f(i)+\sum_{i=1}^n[\forall p>\sqrt n\land p\in P,p\nmid i]f(i) \\=\sum_{i=1}^{\lfloor\sqrt n\rfloor}f(i)\sum_{\sqrt n<p\le\lfloor\frac ni\rfloor\land p\in P}f(p)+\sum_{1\le i\le n\land\forall p>\sqrt n\land p\in P,p\nmid i}f(i) \]

设 \(T=\lfloor\sqrt n\rfloor\)。

答案即

\[S(n)=\sum_{i=1}^Tf(i)\sum_{T<p\le\lfloor\frac ni\rfloor\land p\in P}f(p)+\sum_{1\le i\le n\land\forall p>T\land p\in P,p\nmid i}f(i) \]

于是我们应当快速计算如下两部分贡献。

\[\sum_{T<p\le y\land p\in P}f(p) \]

\[\sum_{1\le i\le n\land\forall p>T\land p\in P,p\nmid i}f(i) \]

接下来的推导将说明，此两部分均可以使用质数前缀统计的技巧继续优化。

第一部分

\[\sum_{T<p\le y\land p\in P}f(p) \]

设不超过 \(T\) 的 \(m\) 个质数分别依次为 \(p_1,p_2,p_3,\dots,p_m\)。

设

\[g_k(j,y)=\sum_{v=1}^yv^k[\forall1\le i\le j,p_i\nmid v] \]

显然 \(\sum_{T<p\le y\land p\in P}f(p)\) 可以被 \(g_k(m,y)-1\) 线性表出。（\(f(p)\) 为关于 \(p\) 的低次多项式）

因此考虑怎么快速计算 \(g_k\)。

对于 \(g_k(0,y)\)，

\[g_k(0,y)=\sum_{v=1}^yv^k \]

这个显然可以 Lagrange 插值，在 \(O(k)\) 时间内计算单项。当 \(k\) 是常数时，其可视作 \(O(1)\)。

考虑对 \(g_k(j,y)\) 做递推。

考虑 \(j\) 增量 \(1\) 时，哪部分贡献会被除去。

显然是是 \(p_j\) 倍数的部分。这部分的贡献可以与 \(g_k(j-1,\lfloor\frac y{p_j}\rfloor)\) 的贡献构成双射，即

\[g_k(j,y)=g_k(j-1,y)-p_j^kg_k(j-1,\lfloor\frac y{p_j}\rfloor) \]

对第一部分直接递推做的复杂度是 \(O(\frac n{\log n})\) 的，需要优化。

考虑如下事实：

\(p_{j+1}>y\) 时，\(g_k(j,y)=1\)。

因此当 \(p_j\le y<p_j^2\) 时，

\[g_k(j,y)=g_k(j-1,y)-p_j^kg_k(j-1,\lfloor\frac y{p_j}\rfloor) \\=g_k(j-1,y)-p_j^k \]

也即，只用计算到最后一个使得 \(p_j^2\ge y\) 的部分为止，剩下的 \(p_j^k\) 的贡献可以后缀和预处理。

这样第一部分复杂度就被优化到了 \(O(\frac{n^{\frac34}}{\log n})\)。

第二部分

\[\sum_{1\le i\le n\land\forall p>T\land p\in P,p\nmid i}f(i) \]

设

\[f(j,y)=\sum_{v=1}^nf(v)[\forall p>T\land p\in P,p\nmid v][\forall1\le i\le j,p_i\nmid v] \]

显然

\[\sum_{1\le i\le n\land\forall p>T\land p\in P,p\nmid i}f(i)=f(0,n) \]

容易发现，

\[f(m,y)=1 \]

仍然考虑递推，使 \(j\) 减量 \(1\)，则

\[f(j-1,y)=\sum_{c\ge0}f(p_j^c)f(j,\lfloor\frac y{p_j^c}\rfloor) \]

朴素实现仍为 \(O(\frac n{\log n})\) 的，考虑优化。

考虑如下事实：

\(p_{j+1}>y\) 时，\(f(j,y)=1\)。

因此当 \(p_j\le y<p_j^2\) 时，

\[f(j-1,y)=f(j,y)+f(p_j) \]

预处理 \(f(p_j)\) 后缀和即可优化到 \(O(\frac{n^{\frac34}}{\log n})\)。

Code

摆。

块筛代码类似，这里鸽了。

以下是 min-25 板子代码，能不能过要看评测机波动……

封装了起来，之后直接用即可。

const ullt Mod=1e9+7;
typedef ConstMod::mod_ullt<Mod>modint;
typedef std::vector<modint>modvec;
namespace ZhougeSieve{
    ullt SQRT(ullt v){
        ullt ans=sqrt(v);
        while(ans*ans<=v)ans++;
        while(ans*ans>v)ans--;
        return ans;
    }
    const uint Lim=200000;
    const uint KLim=30;
    ullt V[Lim+5],Prime[Lim+5];
    uint End[Lim+5],End2[Lim+5];
    modint F[Lim+5],SufSum[Lim+5],M[Lim+5],Val[Lim+5];
    modint P[KLim+5],Q[KLim+5],B[KLim+5];
    modint solve(ullt n,modvec Ppoly,std::function<modint(ullt,uint,ullt)>Func){
        modint ans;
        uint tp=0,tp2=0;for(ullt l=1,r;l<=n;l=r+1)r=n/(n/l),V[tp2++]=n/l;
        std::unordered_map<ullt,uint>Find;for(uint i=0;i<tp2;i++)Find.insert({V[i],i});
        uint T=SQRT(n);F[1]=1;
        {
            static bol User[Lim+5];for(uint i=2;i<=T;i++)User[i]=0;
            static uint Pre[Lim+5];
            for(ullt i=2;i<=T;i++){
                if(!User[i]){
                    Prime[++tp]=i;uint k=1;Pre[i]=Pre[i-1]+1;
                    for(ullt j=i;j<=T;j*=i,k++)F[j]=Func(i,k,j),User[j]=1;
                }
                else Pre[i]=Pre[i-1];
                for(uint q=1;;q++){
                    ullt p=Prime[q];
                    if(i%p&&i*p<=T)for(ullt t=p;i*t<=T;t*=p)F[i*t]=F[i]*F[t],User[i*t]=1;
                    else break;
                }
            }
            for(uint i=0;i<tp2;i++)End[i]=Pre[SQRT(V[i])];
            for(uint i=0;i<tp2;i++)End2[i]=V[i]>T?tp:Pre[V[i]];
        }
        SufSum[tp]=0;
        {
            uint K=Ppoly.size()-1;
            for(uint i=1;i<=tp;i++)Val[i]=1;
            P[0]=1;for(uint i=1;i<=K+1;i++)P[i]=P[i-1]*i;
            Q[K+1]=P[K+1].inv();for(uint i=K+1;i;i--)Q[i-1]=Q[i]*i;
            B[0]=1;for(uint i=1;i<=K;i++){B[i]=0;for(uint j=0;j<i;j++)B[i]-=P[i]*Q[j]*Q[i-j+1]*B[j];}
            for(uint k=0;k<=K;k++){
                if(Ppoly[k]()){
                    for(uint i=tp;i;i--)SufSum[i-1]=SufSum[i]+Val[i];
                    for(uint i=0;i<tp2;i++){
                        modint v(1),b(V[i]+(bol)k);M[i]=0;
                        for(uint j=k;~j;j--)M[i]+=P[k]*Q[j]*Q[k+1-j]*B[j]*(v*=b);
                    }
                    for(uint j=1;j<=tp;j++)for(uint i=0;j<=End[i];i++){
                        uint a=Find[V[i]/Prime[j]];
                        M[i]-=Val[j]*(End[a]<j?M[a]-SufSum[End[a]]+SufSum[j-1]:M[a]);
                    }
                    for(uint i=1,a=0;i<=T;i++){
                        while(V[a]!=n/i)a++;
                        ans+=F[i]*Ppoly[k]*(M[a]-SufSum[End[a]]-1);
                    }
                }
                for(uint i=0;i<=tp;i++)Val[i]*=Prime[i];
            }
        }
        {
            for(uint i=tp;i;i--)SufSum[i-1]=SufSum[i]+F[Prime[i]];
            for(uint i=0;i<tp2;i++)M[i]=1+SufSum[End[i]]-SufSum[End2[i]];
            for(uint j=tp;j;j--)for(uint i=0;j<=End[i];i++){
                uint k=1;
                for(ullt v=Prime[j];v<=V[i];v*=Prime[j],k++){
                    uint a=Find[V[i]/v];
                    M[i]+=Func(Prime[j],k,v)*(j>End[a]?(j>End2[a]?1:1+SufSum[j]-SufSum[End2[a]]):M[a]);
                }
            }
            ans+=M[0];
        }
        return ans;
    }
}
int main()
{
#ifdef MYEE
    freopen("QAQ.in","r",stdin);
    // freopen("QAQ.out","w",stdout);
#endif
    ullt n;scanf("%llu",&n);
    ZhougeSieve::solve(n,{0,-modint(1),1},[](ullt p,uint c,ullt v){return v*modint(v-1);}).println();
    return 0;
}