min25筛学习笔记&模板详解

    xiaoxiao2022-07-07  151

    min25筛算法推导&模板详解

    问题引入

    引入函数 σ 0 ( n ) \sigma _0(n) σ0(n)为n的正因数的数量,求 S ( n , k ) = ∑ i = 1 n σ 0 ( i k ) S(n,k)=\sum_{i=1}^n\sigma_0(i^k) S(n,k)=i=1nσ0(ik), n , k ≤ 1 0 10 n,k\le 10^{10} n,k1010

    简化分析与引入

    对于 σ 0 ( n ) \sigma_0(n) σ0(n)函数,显然有: 当 p 为 素 数 , n 为 p k   σ 0 ( n ) = k + 1 当p为素数,n为p^k\ \sigma_0(n)=k+1\\ pnpk σ0(n)=k+1 且显然 σ 0 \sigma_0 σ0为积性函数

    n ≤ 1 0 7 n\le 10^7 n107时可以通过

    #include<bits/stdc++.h> const int size=1e7+5; bool prime[size]; int p[size]; int sigma[size]; int tot=0; int k; void init() { sigma[1]=1; for(int i=2;i<size;i++) prime[i]=true; for(int i=2;i<size;i++) { if(prime[i]) { p[tot++]=i; sigma[i]=k+1; } for(int j=0;j<tot&&p[j]*i<size;j++) { prime[i*p[j]]=false; if(i%p[j]==0) { sigma[i*p[j]]=sigma[i]*2-sigma[i/p[j]]; break; } else { sigma[i*p[j]]=sigma[i]*(k+1); } } } } int main() { int n; scanf("%d%d",&n,&k); init(); int ans=0; for(int i=1;i<=n;i++) { ans+=sigma[i]; } printf("%d\n",ans); }

    这样对于每个n,k我们都可以线性求解出答案。但是当n最大可以达到 1 0 10 10^{10} 1010时,这种做法显然是不行的为此就需要引入一种亚线性筛:“min25筛”

    虽然min25筛乃至杜教筛名字上都带有一个筛字但是实际上并不能称作一种筛,准确地说,其应该是一种可以亚线性地复杂度下求出积性函数前缀和的算法。其中min25筛求解积性函数前缀和时需要满足条件

    f ( x ) f(x) f(x)当x为质数时有一个多项式表示 f ( x c ) f(x^c) f(xc)在x为质数时可以快速计算

    这两条条件都将在下面的推导种起到重要的作用。

    前置知识

    埃拉托斯特尼筛法

    在线性求解一些数论问题时常常需要用到线性筛,但是有事也可以使用复杂度稍差一点,但是依然有效的埃拉托斯特尼筛法(下面简称埃氏筛法)。其主要的想法是每筛出来一个素数,就将其倍数也删去,这样不断向前递推,未被筛去的就是素数

    bool prime[N]; void init(){ for(int i=2;i<N;i++) prime[i]=true; for(int i=2;i<N;i++){ if(prime[i]){ for(int j=2*i;j<N;j+=i){ prime[j]=false; } } } }

    当然,min25筛其实并不用到埃氏筛,但是其借鉴了一点埃氏筛的思想

    引理1

    对于给定的正整数n和所有的正整数, ⌊ n m ⌋ \lfloor\frac{n}{m}\rfloor mn的取值只有 O ( n ) O(\sqrt{n}) O(n )

    证明. 当 m ≤ n ​ m\le \sqrt n​ mn 时满足条件的m只有 O ( n ) ​ O(\sqrt{n})​ O(n )种。当 m > n ​ m>\sqrt n​ m>n 时, 0 ≤ ⌊ n m ⌋ < n ​ 0\le \lfloor\frac{n}{m}\rfloor<\sqrt n​ 0mn<n 满足这一条件的 ⌊ n m ⌋ ​ \lfloor\frac{n}{m}\rfloor​ mn也只有 O ( n ) ​ O(\sqrt{n})​ O(n )种,故得证

    min25筛的分析与推导

    我们将问题泛化,令 f ( x ) = σ 0 ( x k ) f(x)=\sigma_0(x^k) f(x)=σ0(xk),由于 σ 0 ( x ) \sigma_0(x) σ0(x)函数为积性函数,易证*f(x)*也为积性函数。且因为 f ( x c ) = c ∗ k + 1 f(x^c)=c*k+1 f(xc)=ck+1,因此积性函数 f ( x ) f(x) f(x)满足上面所提到的min25筛求解积性函数前缀和的条件。现要求积性函数 f ( x ) ​ f(x)​ f(x)的前缀和。

    为了方便公式的书写,我们需要先引入一些定义 P : 素 数 集 合 l p f ( n ) : n 的 最 小 质 因 数 P i : 埃 氏 筛 法 筛 了 前 i 个 质 数 之 后 剩 下 的 自 然 数 集 合 p i : 第 i 个 质 数 P:素数集合\\ lpf(n):n的最小质因数\\ P_i:埃氏筛法筛了前i个质数之后剩下的自然数集合\\ p_i:第i个质数 P:lpf(n):nPi:ipi:i

    将原前缀和的求和方式进行转化 ∑ i = 1 n f ( i ) = 1 + ∑ 2 ≤ p ≤ n   p ∈ P ∑ 2 ≤ x ≤ n   l p f ( x ) = p f ( x ) \sum_{i=1}^nf(i)=1+\sum_{2\le p\le n\,p\in P}\sum_{2\le x\le n\,lpf(x)=p}f(x) i=1nf(i)=1+2pnpP2xnlpf(x)=pf(x) 由于*f(x)*为积性函数,故 1 + ∑ 2 ≤ p ≤ n   p ∈ P ∑ 2 ≤ x ≤ n   l p f ( x ) = p f ( x ) = 1 + ∑ 2 ≤ p e ≤ n , e ≥ 1 , p ∈ P f ( p e ) ( 1 + ∑ 2 ≤ x ≤ ⌊ n p e ⌋ , l p f ( x ) > p f ( x ) ) 1+\sum_{2\le p\le n\,p\in P}\sum_{2\le x\le n\,lpf(x)=p}f(x)=1+\sum_{2\le p^e\le n,e\ge1,p\in P}f(p^e)\left(1+\sum_{2\le x\le \lfloor\frac{n}{p^e}\rfloor,lpf(x)>p}f(x)\right) 1+2pnpP2xnlpf(x)=pf(x)=1+2pen,e1,pPf(pe)1+2xpen,lpf(x)>pf(x) 显然,对于一个合数来说,其最小的质数定然不会超过 n \sqrt{n} n ,因此上式就可以进行进一步的转化 ∑ i = 1 n f ( i ) = 1 + ∑ 2 ≤ p ≤ n   p ∈ P ∑ 2 ≤ x ≤ n   l p f ( x ) = p f ( x ) = 1 + ∑ 2 ≤ p e ≤ n , e ≥ 1 , p ∈ P f ( p e ) ( 1 + ∑ 2 ≤ x ≤ ⌊ n p e ⌋ , l p f ( x ) > p f ( x ) ) = 1 + ∑ 2 ≤ p ≤ n 2 ≤ p e ≤ n , e ≥ 1 , p ∈ P f ( p e ) ( 1 + ∑ 2 ≤ x ≤ ⌊ n p e ⌋ , l p f ( x ) > p f ( x ) ) + ∑ n < p ≤ n , p ∈ P f ( p ) \begin{aligned} \sum_{i=1}^nf(i)&=1+\sum_{2\le p\le n\,p\in P}\sum_{2\le x\le n\,lpf(x)=p}f(x)\\ &=1+\sum_{2\le p^e\le n,e\ge1,p\in P}f(p^e)\left(1+\sum_{2\le x\le \lfloor\frac{n}{p^e}\rfloor,lpf(x)>p}f(x)\right)\\&=1+\sum_{2\le p\le \sqrt{n}2\le p^e\le n,e\ge1,p\in P }f(p^e)\left(1+\sum_{2\le x\le \lfloor\frac{n}{p^e}\rfloor,lpf(x)>p}f(x)\right)+\sum_{\sqrt{n}<p\le n,p\in P}f(p) \end{aligned} i=1nf(i)=1+2pnpP2xnlpf(x)=pf(x)=1+2pen,e1,pPf(pe)1+2xpen,lpf(x)>pf(x)=1+2pn 2pen,e1,pPf(pe)1+2xpen,lpf(x)>pf(x)+n <pn,pPf(p) 但是直接求这个有一定的难度,为此引入两个新的函数 g ( n , m ) = ∑ 2 ≤ x ≤ n , l p f ( x ) > m f ( x ) h ( n ) = ∑ 2 ≤ p ≤ n , p ∈ P f ( p ) g(n,m)=\sum_{2\le x\le n,lpf(x)>m}f(x)\\ h(n)=\sum_{2\le p\le n,p\in P} f(p) g(n,m)=2xn,lpf(x)>mf(x)h(n)=2pn,pPf(p) 原式即是 g ( n , 0 ) g(n,0) g(n,0)

    为此我们对g(n,m)进行推导 g ( n , m ) = ∑ 2 ≤ x ≤ n , l p f ( x ) > m f ( x ) = ∑ m < p ≤ n , p e ≤ n , e ≥ 1 , p ∈ P f ( p e ) ( 1 + ∑ 2 ≤ x ≤ ⌊ n p e ⌋ , l p f ( x ) > p f ( x ) ) + ∑ n < p ≤ n , p ∈ P f ( p ) = ∑ m < p ≤ n , p e ≤ n , e ≥ 1 , p ∈ P f ( p e ) ( 1 + g ( ⌊ n p e , p ⌋ ) ) + h ( n ) − h ( n ) \begin{aligned} g(n,m)&=\sum_{2\le x\le n,lpf(x)>m}f(x)\\ &=\sum_{m<p\le \sqrt{n},p^e\le n,e\ge1,p\in P}f(p^e)\left(1+\sum_{2\le x\le \lfloor\frac{n}{p^e}\rfloor,lpf(x)> p}f(x)\right)+\sum_{\sqrt{n}<p\le n,p\in P}f(p)\\ &=\sum_{m<p\le \sqrt{n},p^e\le n,e\ge1,p\in P}f(p^e)\left(1+g(\lfloor\frac{n}{p^e},p\rfloor)\right)+h(n)-h(\sqrt{n}) \end{aligned} g(n,m)=2xn,lpf(x)>mf(x)=m<pn ,pen,e1,pPf(pe)1+2xpen,lpf(x)>pf(x)+n <pn,pPf(p)=m<pn ,pen,e1,pPf(pe)(1+g(pen,p))+h(n)h(n ) 但是这样不便于式子之间的转化,因此将原式中的素数部分的函数值的递推转化一下 g ( n , m ) = ∑ m < p ≤ n , p e ≤ n , e ≥ 1 , p ∈ P f ( p e ) ( [ e > 1 ] + g ( ⌊ n p e , p ⌋ ) ) + ∑ m < p ≤ n , p ∈ P f ( x ) + h ( n ) − h ( n ) = ∑ m < p ≤ n , p e ≤ n , e ≥ 1 , p ∈ P f ( p e ) ( [ e > 1 ] + g ( ⌊ n p e , p ⌋ ) ) + h ( n ) − h ( m ) \begin{aligned} g(n,m)&=\sum_{m<p\le \sqrt{n},p^e\le n,e\ge1,p\in P}f(p^e)\left([e>1]+g(\lfloor\frac{n}{p^e},p\rfloor)\right)+\sum_{m<p\le \sqrt{n},p\in P}f(x)+h(n)-h(\sqrt{n})\\ &=\sum_{m<p\le \sqrt{n},p^e\le n,e\ge1,p\in P}f(p^e)\left([e>1]+g(\lfloor\frac{n}{p^e},p\rfloor)\right)+h(n)-h(m) \end{aligned} g(n,m)=m<pn ,pen,e1,pPf(pe)([e>1]+g(pen,p))+m<pn ,pPf(x)+h(n)h(n )=m<pn ,pen,e1,pPf(pe)([e>1]+g(pen,p))+h(n)h(m) 对于我们想要得到的 g ( n , 0 ) g(n,0) g(n,0)我们需要枚举所有的小于等于 n \sqrt{n} n 的素数作为m并将其对应的 g ( n , m ) g(n,m) g(n,m)加上才可以最后得到答案。需要注意的是当 n ≤ m 2 n\le m^2 nm2时, g ( n , m ) = 0 g(n,m)=0 g(n,m)=0这是显然可以得知的,当 n ≤ m 2 n\le m^2 nm2时,1到n之间,不存在最小质因数还大于m的数字。

    由于前面的条件 f ( p e ) f(p^e) f(pe)是可以快速计算的,因此接下来只要预处理出h就可以进行递推了。

    接下来考虑如何求 h h h

    首先可以发现所有需要用到的h的标号均为0到 n ​ \sqrt{n}​ n 之间的质数,或为 ⌊ n p e ⌋ ​ \lfloor\frac{n}{p^e}\rfloor​ pen

    现在证明:对于所有需要用到的 h ( i ) h(i) h(i)其一定满足,存在m使得 ⌊ n m ⌋ = i \lfloor\frac{n}{m}\rfloor=i mn=i

    证 : 令 ⌊ n p ⌋ = k , p 为 任 意 的 小 于 等 于 n 的 数 则 有 n = p ∗ k + b ( b < p ) n ÷ k = p … … b ∵ p ≤ n ∴ k ≥ n ∴ k > b ∴ ⌊ n k ⌋ = p \begin{aligned} 证:&令\lfloor\frac{n}{p}\rfloor=k,p为任意的小于等于n的数\\ &则有n=p*k+b(b<p)\\ &n\div k=p……b\\ &\because p\le\sqrt{n}\\ &\therefore k\ge\sqrt{n}\\ &\therefore k>b\\ &\therefore\lfloor\frac{n}{k}\rfloor=p\\ \end{aligned} pn=kpnn=pk+b(b<p)n÷k=pbpn kn k>bkn=p 由此得出,所有小于等于 n \sqrt{n} n 的数字都是满足上述的可以通过 ⌊ n m ⌋ \lfloor\frac{n}{m}\rfloor mn得到的,而所有0到 n \sqrt{n} n 之间的数字自然也就满足这个条件,而 ⌊ n p e ⌋ \lfloor\frac{n}{p^e}\rfloor pen则是也显然满足这个条件,因此得证

    根据引理1,满足这种条件的数字共有 O ( n ) O(\sqrt n) O(n )个.由于h函数是对f函数求素数前缀和,且上面说到了,所有满足可以用min25筛求的积性函数都需要满足 f ( p ) f(p) f(p)可以通过多项式表示,也即f函数可以表示为 f ( p ) = ∑ a t p b t f(p)=\sum a_tp^{b_t} f(p)=atpbt,不妨令 f t ( p ) = a t p b t f_t(p)=a_tp^{b_t} ft(p)=atpbt, f t f_t ft函数的前缀和为 h t ( i ) h_t(i) ht(i)因此h函数也就可以表示为 h ( i ) = ∑ t h t ( i ) h(i)=\sum_{t}h_t(i) h(i)=tht(i)。进而,对于每个h(i)我们只需要求出所有的 h t ( i ) h_t(i) ht(i)即可。

    到了这里我们的问题就变成了:给定n和t,对所有满足 i = ⌊ n m ⌋ i=\lfloor\frac{n}{m}\rfloor i=mn,求: h t ( i ) = ∑ p ≤ i , p ∈ P p b t h_t(i)=\sum_{p\le i,p\in P}p^{b_t} ht(i)=pi,pPpbt 为了求这个函数我们需要引入一个新的函数 h t , i ′ ( j ) h'_{t,i}(j) ht,i(j),其含义为埃氏筛法筛出前i个数之后,前j个数字种剩余的数字的 b t b_t bt次方和 h t , i ′ ( j ) = ∑ 1 ≤ p ≤ j p = 1   o r   p ∈ P   o r   l p f ( p ) > p i p b t h'_{t,i}(j)=\sum_{1\le p \le j\\p=1\ or\ p\in P\ or\ lpf(p)>p_i}p^{b_t} ht,i(j)=1pjp=1 or pP or lpf(p)>pipbt 那么函数 h t ( j ) = h t , i ( j ) ( j < p i + 1 ∗ p i + 1 ) h_t(j)=h_{t,i}(j)(j<p_{i+1}*p_{i+1}) ht(j)=ht,i(j)(j<pi+1pi+1)

    对于 h t , i ′ ( j ) h'_{t,i}(j) ht,i(j)有转移式 h t , i ( j ) = h t , 0 ′ ( j ) − ∑ k = 0 k < i p k + 1 b t ( h t , k ′ ( ⌊ j p k + 1 ⌋ ) − h t , k ′ ( p k ) ) h_{t,i}(j)=h'_{t,0}(j)-\sum_{k=0}^{k<i}p_{k+1}^{b_t}(h'_{t,k}(\lfloor\frac{j}{p_{k+1}}\rfloor)-h'_{t,k}(p_k)) ht,i(j)=ht,0(j)k=0k<ipk+1bt(ht,k(pk+1j)ht,k(pk)) 由此,h函数的求值就结束了。

    由此就可以得到积性函数f的前缀和了。上述算法的复杂度为 O ( n 3 4 l o g ( n ) ) O(\frac{n^{\frac{3}{4}}}{log(n)}) O(log(n)n43),现在还有一种新版的min25筛可以优化到 O ( n 2 3 ) O(n^{\frac{2}{3}}) O(n32)但是实际应用上,由于常数巨大,实际上的运行时间并不优化多少,且代码冗长,故此处不予介绍。

    模板及解释

    //问题描述定义: //定义sigma(n)为n的正因数的数量,求f(n,k)=sum_{i=1}^n (sigma(i^k)) #include<stdio.h> #include<math.h> using namespace std; typedef unsigned long long ull; typedef long long ll; const int N=200005; int T,S,pr[N],pc; ll n,num[N],m,K; ull g[N]; bool fl[N]; // 给定一个数字X求出其为第几个可以得到有效的g的数字 inline int ID(ll x){return x<=S?x:m-n/x+1;} ull f(ll n,int i){ if(n<1||pr[i]>n)return 0; ull ret=g[ID(n)]-(i-1)*(K+1); while((ll)pr[i]*pr[i]<=n){ int p=pr[i]; ull e=K+1,t=n/p; while(t>=p)ret+=f(t,i+1)*e+e+K,t/=p,e+=K; // ret= sum{sigma(p^es)([n>1]+f(n/p^es,p)+g[n]-g[num[i-1]]} // 因为对于函数g(n,m)当n<=m^2,g(n,m)为0 // 且根据sigma函数的性质,对于质数p,有sigma(p^es)=es*k+1 // 所以 ret+=(es*k+1)*f(n/p^es,p)+((es+1)k+1)(1<=es&& n/p^es>p)即当前项的sigma(p^e)*f(n/p^e,p)加上下一项的sigma(p^e) // 这样的做法避免了e<1也就不必进行特判 i++; } return ret; } //g[i]即小于i的所有质数的sigma函数求和 //根绝前面的推导,只有当i可以通过[n/m]得到时,其才有用 ull solve(ll n){ int i,p,x;ull y; S=sqrt(n); while((ll)S*S<=n)S++; while((ll)S*S>n)S--; while(m)num[m--]=0; for(i=1;i<=S;i++)num[++m]=i; for(i=S;i>=1;i--)if(n/i>S)num[++m]=n/i; for(i=1;i<=m;i++)g[i]=num[i]-1; //此处g[i]为小于等于第i个满足可以通过[n/m]得到的数的素因子的数量 //故先减去“1”因为1一定不为素因子且无法被后续操作筛去 x=1;y=0; for(p=2;p<=S;p++)if(g[p]!=g[p-1]){ while(num[x]<(ll)p*p)x++; //令g'(i,j)为埃氏筛法筛出前j个素数后,1到j之间剩余数的数量 //则有g'(i,j)=g'(0,j)-sum_{k=0,k<i}(g'(k,[j/p[k+1]])-g'(k,p[k])) //其中p[0]=1,p[i](i>0)为第i个质数 //则g(j)即为g'(i,j)使得有p[i]<j&&p[i+1]>j //其中g'(k,p[k])即为k+1 //由于以上的递推式仅仅作用到num[x]>=p^2为止,故仅仅将更新到x //由于每次都要减去尚未去除当前素数的g以充当g'因此从大往小更新 for(i=m;i>=x;i--)g[i]-=g[ID(num[i]/p)]-y; y++; } for(i=1;i<=m;i++)g[i]*=K+1; return f(n,1)+1; } int main(){ int i,j; //素数筛 for(i=2;i<N;i++)if(!fl[i]){ pr[++pc]=i; for(j=i+i;j<N;j+=i)fl[j]=1; } for(scanf("%d",&T);T--;){ scanf("%lld%lld",&n,&K); printf("%llu\n",solve(n)); } return 0; }

    参考资料

    朱震霆的国家集训队论文: 《一些特殊的数论函数求和问题》, 国家集训队2018论文集新版min25筛(O(n^(2/3)))详解 https://zhuanlan.zhihu.com/p/60378354yyb的博客 https://www.cnblogs.com/cjyyb/p/9185093.html
    最新回复(0)