引入函数 σ 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,k≤1010
对于 σ 0 ( n ) \sigma_0(n) σ0(n)函数,显然有: 当 p 为 素 数 , n 为 p k σ 0 ( n ) = k + 1 当p为素数,n为p^k\ \sigma_0(n)=k+1\\ 当p为素数,n为pk σ0(n)=k+1 且显然 σ 0 \sigma_0 σ0为积性函数
当 n ≤ 1 0 7 n\le 10^7 n≤107时可以通过
#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筛其实并不用到埃氏筛,但是其借鉴了一点埃氏筛的思想
对于给定的正整数n和所有的正整数, ⌊ n m ⌋ \lfloor\frac{n}{m}\rfloor ⌊mn⌋的取值只有 O ( n ) O(\sqrt{n}) O(n )种
证明. 当 m ≤ n m\le \sqrt n m≤n 时满足条件的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 0≤⌊mn⌋<n 满足这一条件的 ⌊ n m ⌋ \lfloor\frac{n}{m}\rfloor ⌊mn⌋也只有 O ( n ) O(\sqrt{n}) O(n )种,故得证
我们将问题泛化,令 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)=c∗k+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):n的最小质因数Pi:埃氏筛法筛了前i个质数之后剩下的自然数集合pi:第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=1∑nf(i)=1+2≤p≤np∈P∑2≤x≤nlpf(x)=p∑f(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+2≤p≤np∈P∑2≤x≤nlpf(x)=p∑f(x)=1+2≤pe≤n,e≥1,p∈P∑f(pe)⎝⎛1+2≤x≤⌊pen⌋,lpf(x)>p∑f(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=1∑nf(i)=1+2≤p≤np∈P∑2≤x≤nlpf(x)=p∑f(x)=1+2≤pe≤n,e≥1,p∈P∑f(pe)⎝⎛1+2≤x≤⌊pen⌋,lpf(x)>p∑f(x)⎠⎞=1+2≤p≤n 2≤pe≤n,e≥1,p∈P∑f(pe)⎝⎛1+2≤x≤⌊pen⌋,lpf(x)>p∑f(x)⎠⎞+n <p≤n,p∈P∑f(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)=2≤x≤n,lpf(x)>m∑f(x)h(n)=2≤p≤n,p∈P∑f(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)=2≤x≤n,lpf(x)>m∑f(x)=m<p≤n ,pe≤n,e≥1,p∈P∑f(pe)⎝⎛1+2≤x≤⌊pen⌋,lpf(x)>p∑f(x)⎠⎞+n <p≤n,p∈P∑f(p)=m<p≤n ,pe≤n,e≥1,p∈P∑f(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<p≤n ,pe≤n,e≥1,p∈P∑f(pe)([e>1]+g(⌊pen,p⌋))+m<p≤n ,p∈P∑f(x)+h(n)−h(n )=m<p≤n ,pe≤n,e≥1,p∈P∑f(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 n≤m2时, g ( n , m ) = 0 g(n,m)=0 g(n,m)=0这是显然可以得知的,当 n ≤ m 2 n\le m^2 n≤m2时,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⌋=k,p为任意的小于等于n的数则有n=p∗k+b(b<p)n÷k=p……b∵p≤n ∴k≥n ∴k>b∴⌊kn⌋=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)=p≤i,p∈P∑pbt 为了求这个函数我们需要引入一个新的函数 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)=1≤p≤jp=1 or p∈P or lpf(p)>pi∑pbt 那么函数 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+1∗pi+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=0∑k<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)但是实际应用上,由于常数巨大,实际上的运行时间并不优化多少,且代码冗长,故此处不予介绍。