积性函数求和问题的一种筛法

前言

不知不觉鸽了这么久了,马上也要到今年的冬令营了… 还是早点把这个去年的想法写完吧。

这篇文章主要讨论曾经常在 OI 里出现的积性函数求和问题。

基本定义和性质

记号 ⌊ x ⌋ \lfloor x \rfloor x表示不超过 x x x的最大整数,记号 ⌈ x ⌉ \lceil x \rceil x表示不小于 x x x的最小整数。

记号 d ∣ n d | n dn表示 d d d整除 n n n gcd ⁡ ( x , y ) \gcd(x, y) gcd(x,y)表示 x x x y y y的最大公约数, [ x ] [x] [x] x x x为真时值为 1 1 1,否则值为 0 0 0

定义函数 f : Z + → F f : \mathbb{Z}^{+} \rightarrow \mathbb{F} f:Z+F为数论函数,其中 F \mathbb{F} F是一个域,称由全体数论函数构成的集合为 U ​ \mathcal{U} ​ U

若数论函数 g : Z + → F g : \mathbb{Z}^{+} \rightarrow \mathbb{F} g:Z+F满足 ∀ x , y ∈ Z + , gcd ⁡ ( x , y ) = 1 , ∃ g ( x y ) = g ( x ) g ( y ) \forall x, y \in \mathbb{Z}^{+}, \gcd(x, y) = 1, \exists g(xy) = g(x)g(y) x,yZ+,gcd(x,y)=1,g(xy)=g(x)g(y),则称 g g g为积性函数,称由全体积性函数构成的集合为 V ​ \mathcal{V} ​ V

若数论函数 h : Z + → F h : \mathbb{Z}^{+} \rightarrow \mathbb{F} h:Z+F满足 ∀ x , y ∈ Z + , ∃ h ( x y ) = h ( x ) h ( y ) \forall x, y \in \mathbb{Z}^{+}, \exists h(xy) = h(x)h(y) x,yZ+,h(xy)=h(x)h(y),则称 h ​ h ​ h为完全积性函数。

定义 P P P为素数集合,集合中的第 i i i个元素是第 i i i小的素数 p i p_i pi。对于一个正整数 x = ∏ p i e x i x = \prod p_i^{e_{x_i}} x=piexi,定义该正整数的坐标 e x = [ x ] P = ( e x , 1 , e x , 2 , e x , 3 , … ) e_x = [x]_P = (e_{x,1}, e_{x,2}, e_{x,3}, \ldots) ex=[x]P=(ex,1,ex,2,ex,3,)

p ∈ P p \in P pP,定义 V p = { 1 , p , p 2 , p 3 , … } V_p = \{1, p, p^2, p^3, \ldots\} Vp={1,p,p2,p3,},表示 p p p的非负次幂构成的集合。

定义 f p ( n ) = f ( n ) [ n ∈ V p ] ​ f_p(n) = f(n) [n \in V_p] ​ fp(n)=f(n)[nVp],即保留 f ​ f ​ f中下标在 V p ​ V_p ​ Vp中的项后的数论函数。

定义 U ​ \mathcal{U} ​ U上的 Dirichlet 卷积运算 ∗ ​ * ​ ,满足 ( f ∗ g ) ( n ) = ∑ d ∣ n f ( d ) g ( n d ) ​ (f * g)(n) = \sum_{d | n} f(d) g\left(\frac{n}{d} \right) ​ (fg)(n)=dnf(d)g(dn),用 ∏ ^ ​ \hat{\prod} ​ ^表示 Dirichlet 卷积的连续运算,Dirichlet 卷积有以下性质:

性质一: f ∈ V    ⟺    f = ∏ ^ p ∈ P f p ​ f \in \mathcal{V} \iff f = {\hat{\prod}}_{p \in P} f_p ​ fVf=^pPfp,其中 f p ​ f_p ​ fp只有下标在 V p ​ V_p ​ Vp中的位置有值。

证明: f ∈ V f \in \mathcal{V} fV,比较等式两端系数可得 f = ∏ ^ p ∈ P f p f = {\hat{\prod}}_{p \in P} f_p f=^pPfp;若 f = ∏ ^ p ∈ P f p f = {\hat{\prod}}_{p \in P} f_p f=^pPfp,容易验证 f ∈ V f \in \mathcal{V} fV

性质二: ∀ f , g ∈ V , f ∗ g ∈ V \forall f, g \in \mathcal{V}, f * g \in \mathcal{V} f,gV,fgV

证明: f ∗ g = ∏ ^ p ∈ P f p g p ​ f * g = {\hat{\prod}}_{p \in P} f_p g_p ​ fg=^pPfpgp,由性质一有 f ∗ g ∈ V ​ f * g \in \mathcal{V} ​ fgV

性质三: ∀ f ∈ V , ∃ g = f − 1   s . t . f ∗ g = ϵ \forall f \in \mathcal{V}, \exists g = f^{-1}~ \mathrm{s.t.} f*g=\epsilon fV,g=f1 s.t.fg=ϵ,其中 ϵ ( x ) = [ x = 1 ] \epsilon(x) = [x = 1] ϵ(x)=[x=1]

定义 n ∈ Z + ​ n \in \mathbb{Z}^{+} ​ nZ+的整除集合为 D n = { ⌊ n m ⌋ ∣ m ∈ Z + } ​ D_n = \left \{ \lfloor \frac{n}{m} \rfloor | m \in \mathbb{Z}^{+} \right \} ​ Dn={mnmZ+},整除集合有如下性质。

性质一: D n = { 0 , 1 , 2 , … , ⌊ n ⌋ , ⌊ n ⌊ n ⌋ ⌋ , ⌊ n ⌊ n ⌋ − 1 ⌋ , … , n } D_n = \left \{ 0, 1, 2, \ldots, \lfloor\sqrt{n}\rfloor, \left \lfloor \frac{n}{\lfloor \sqrt{n} \rfloor} \right \rfloor, \left \lfloor \frac{n}{\lfloor \sqrt n\rfloor - 1} \right \rfloor, \ldots, n \right \} Dn={0,1,2,,n ,n n,n 1n,,n}

证明: m ≤ n ​ m \le \sqrt n ​ mn 时, ⌊ n / m ⌋ ​ \lfloor n/m \rfloor ​ n/m的取值已经被包含在集合中,否则 ⌊ n m ⌋ ≤ ⌊ n ⌋ ​ \left \lfloor \frac{n}{m} \right \rfloor \le \lfloor \sqrt n \rfloor ​ mnn 。而此时 m > n ​ m > \sqrt n ​ m>n ,若 ⌊ n m ⌋ ≠ ⌊ n m + 1 ⌋ ​ \lfloor \frac{n}{m} \rfloor \ne \lfloor \frac{n}{m + 1} \rfloor ​ mn=m+1n,则有 n m + 1 ≤ ⌊ n m ⌋ ≤ n m ​ \frac{n}{m + 1} \le \left \lfloor \frac{n}{m} \right \rfloor \le \frac{n}{m} ​ m+1nmnmn,但 n m − n m + 1 < 1 ​ \frac{n}{m} - \frac{n}{m + 1} < 1 ​ mnm+1n<1,因此 ⌊ n m ⌋ = ⌊ n m + 1 ⌋ + 1 ​ \lfloor \frac{n}{m} \rfloor = \lfloor \frac{n}{m + 1} \rfloor + 1 ​ mn=m+1n+1,于是结论成立。

性质二: D ( ⌊ n / m ⌋ ) ⊆ D ( n ) , n , m ∈ Z + ​ D\left( \left \lfloor n/m \right \rfloor \right) \subseteq D(n), n,m \in \mathbb{Z}^{+} ​ D(n/m)D(n),n,mZ+

证明: n = k a b + r b + c , r ∈ [ 0 , a ) ∩ Z , c ∈ [ 0 , b ) ∩ Z ​ n = kab + rb + c, r \in [0, a) \cap \mathbb{Z}, c \in [0, b) \cap \mathbb{Z} ​ n=kab+rb+c,r[0,a)Z,c[0,b)Z,则 ⌊ n / a b ⌋ = k = ⌊ ( k a + r ) / a ⌋ = ⌊ ⌊ n / b ⌋ / a ⌋ ​ \left \lfloor n/ab\right \rfloor = k = \left \lfloor (ka+r)/a \right \rfloor = \left \lfloor \lfloor n/b \rfloor /a \right \rfloor ​ n/ab=k=(ka+r)/a=n/b/a,取 b = m ​ b = m ​ b=m可得结论成立。

性质三: U n ( x ) = min ⁡ y ∈ D n , y ≥ x y ​ U_n(x) = \min_{y \in D_n, y \ge x} y ​ Un(x)=minyDn,yxy,其中 1 ≤ x ≤ n ​ 1 \le x \le n ​ 1xn,则 ⌊ n / U n ( x ) ⌋ = ⌊ n / x ⌋ ​ \lfloor n / U_n(x) \rfloor = \lfloor n/x \rfloor ​ n/Un(x)=n/x

证明: 不妨设 y = ⌊ n / x ⌋ y = \lfloor n / x \rfloor y=n/x,那么 x y ≤ n x y \le n xyn,于是 x ≤ n / y x \le n / y xn/y,因此这样的 x x x的上界为 ⌊ n / y ⌋ ∈ D n \lfloor n / y \rfloor \in D_n n/yDn,因此 ⌊ n / U n ( x ) ⌋ ≤ ⌊ n / x ⌋ \lfloor n / U_n(x) \rfloor \le \lfloor n / x \rfloor n/Un(x)n/x,而另一个方向是显然的,于是得证。

推论: ∀ x , y ∈ [ 1 , n ] ∩ Z , x y ≤ n    ⟺    U n ( x ) U n ( y ) ≤ n ​ \forall x, y \in [1, n] \cap \mathbb{Z}, xy \le n \iff U_n(x) U_n(y) \le n ​ x,y[1,n]Z,xynUn(x)Un(y)n

f ∈ U , n ∈ Z + ​ f \in \mathcal{U}, n \in \mathbb{Z}^{+} ​ fU,nZ+定义 f ​ f ​ f定义在 n ​ n ​ n的整除集合上的前缀和函数为 S n , f ( x ) = ∑ y ≤ U n ( x ) f ( y ) , 1 ≤ x ≤ n ​ S_{n,f}(x) = \sum_{y \le U_n(x)} f(y), 1 \le x \le n ​ Sn,f(x)=yUn(x)f(y),1xn,那么在存储 S n , f ​ S_{n, f} ​ Sn,f时我们只要存储 D n ​ D_n ​ Dn中所有下标的值即可。

定理: 给定 S n , f ​ S_{n, f} ​ Sn,f S n , g ​ S_{n, g} ​ Sn,g,可以在 O ( n 2 / 3 log ⁡ 1 / 3 n ) ​ O(n^{2/3} \log^{1/3} n) ​ O(n2/3log1/3n)的时间内求得 S n , f ∗ g ​ S_{n, f * g} ​ Sn,fg

证明: 即对 m ∈ D n m \in D_n mDn S n , f ∗ g ( m ) = ∑ x ∈ D m ( S n , f ( ⌊ m / x ⌋ ) − S n , f ( ⌊ m / U m ( x + 1 ) ⌋ ) ) S n , g ( x ) S_{n, f * g}(m) = \sum_{x \in D_m} \left(S_{n, f}(\lfloor m / x \rfloor) - S_{n, f}(\lfloor m / U_m(x + 1) \rfloor)\right) S_{n, g}(x) Sn,fg(m)=xDm(Sn,f(m/x)Sn,f(m/Um(x+1)))Sn,g(x),由此可以在 O ( m ) O(\sqrt m) O(m ) 的时间内对特定的 m m m求得该值。另一方面,满足 x , y ∈ D n , x y ≤ S x, y \in D_n, xy \le S x,yDn,xyS ( x , y ) (x, y) (x,y)对数是 O ( S log ⁡ S ) O(S \log S) O(SlogS)的,因此由性质二和性质三的推论,对所有 m ≤ S m \le S mS可以在 O ( S log ⁡ S ) O(S \log S) O(SlogS)内求得答案。而对超过 S S S m m m消耗的总时间为 O ( ∫ 1 n / S ( n / x ) 1 / 2 ) = O ( n / S 1 / 2 ) O(\int_1^{n/S} (n/x)^{1/2}) = O(n / S^{1/2}) O(1n/S(n/x)1/2)=O(n/S1/2),取 S = n 2 / 3 log ⁡ − 2 / 3 n S = n^{2/3} \log^{-2/3} n S=n2/3log2/3n时得到 O ( n 2 / 3 log ⁡ 1 / 3 n ) O(n^{2/3} \log^{1/3} n) O(n2/3log1/3n)的时间复杂度,故命题正确。

另外注意到,若 f f f g g g是积性的,我们可以在 O ( S ) O(S) O(S)的时间内求得前 S S S项的值,于是复杂度可以降到 O ( n 2 / 3 ) O(n^{2/3} ) O(n2/3)

特别地,如果 S n , f ( x ) = c S_{n, f}(x) = c Sn,f(x)=c对所有 x ≤ n x \le \sqrt n xn 成立,合并可以用暴力在 O ( n log ⁡ n ) O(\sqrt n \log n) O(n logn)的时间内完成。

定义 Powerful Number 为所有的满足 ∀ k ∈ Z + , e x , k ≠ 1 \forall k \in \mathbb{Z}^{+}, e_{x, k} \ne 1 kZ+,ex,k=1 x x x,令 W W W为全体 Powerful Number 组成的集合。

性质: n ​ n ​ n以内的 Powerful Number 的数目是 Θ ( n ) ​ \Theta(\sqrt n) ​ Θ(n )的。

证明: 一方面,所有的平方数均 powerful,因此数目是 Ω ( n ) ​ \Omega(\sqrt{n}) ​ Ω(n )的。另一方面,所有的 powerful number 总能表示为 a 2 b 3 ​ a^2 b^3 ​ a2b3的形式,这样的数的数目为 O ( ∫ 1 n 1 / 3 ( n / x 3 ) 1 / 2 d x ) + O ( n 1 / 3 ) = O ( n 1 / 2 ) ​ O(\int_1^{n^{1/3}} (n/x^3)^{1/2} \mathrm{d} x) + O(n^{1/3}) = O(n^{1/2}) ​ O(1n1/3(n/x3)1/2dx)+O(n1/3)=O(n1/2),因此原命题得证。

定理: f f f中只有 Powerful Number 处有值,且我们有 S n , g S_{n, g} Sn,g,我们可以在 O ~ ( n 0.6 ) \tilde{O}(n^{0.6}) O~(n0.6)的时间内求得 S n , f ∗ g S_{n, f * g} Sn,fg

证明: 由 Powerful Number 的性质,按照 m 2 / 3 m^{2/3} m2/3分段可知, S m , f S_{m, f} Sm,f中值不同的项的数目是 O ( m 1 / 3 ) O(m^{1/3}) O(m1/3)的。同样采用类似的方法,平衡较小部分的复杂度 O ~ ( S ) \tilde{O}(S) O~(S)和较大部分的复杂度 O ( n / S 2 / 3 ) O(n / S^{2/3}) O(n/S2/3)可在 S = O ~ ( n 0.6 ) ​ S = \tilde{O}(n^{0.6}) ​ S=O~(n0.6)时取得平衡。于是命题成立。

一种筛法

OI 中出现的积性函数求和问题往往分为两种:给定 f p f_p fp的生成规则,求 S n , f ( n ) S_{n, f}(n) Sn,f(n)或是求出 S n , f S_{n, f} Sn,f

由于对于 n < p ≤ n \sqrt n < p \le n n <pn p p p O ( n log ⁡ − 1 n ) O(n \log^{-1} n) O(nlog1n)的,为了给出生成规则我们一般对于这一部分的 p p p一定有 f ( p ) = y ( p ) f(p) = y(p) f(p)=y(p),其中 y ( x ) = ∑ i ∈ T a i x i , T ⊂ Z + y(x) = \sum_{i \in T} a_i x^i, T \subset \mathbb{Z}^{+} y(x)=iTaixi,TZ+,并且 ∣ T ∣ |T| T是一个常数,且其中元素的大小不太大。

只求单一的 S n , f ( n ) S_{n, f}(n) Sn,f(n)的值可以用 Min-25 筛,虽然是 Ω ( n c )   ∀ c < 1 \Omega(n^{c}) \ \forall c < 1 Ω(nc) c<1的但在现在的范围 ( ∼ 1 0 11 ) (\sim 10^{11}) (1011)里还是非常 work 的,求 S n , f S_{n, f} Sn,f的传统方法有:

  1. 杜教筛:已知 S n , f S_{n, f} Sn,f S n , g S_{n, g} Sn,g,可以在 O ( n 2 / 3 ) O(n^{2/3}) O(n2/3)内求 S n , f ∗ g − 1 S_{n, f * g^{-1}} Sn,fg1
  2. 洲阁筛(扩展埃筛):若 ∣ S ∣ |S| S是常数,可以在 O ( n 3 / 4 log ⁡ − 1 n ) O(n^{3/4} \log^{-1} n) O(n3/4log1n)内求 S n , f S_{n, f} Sn,f

以上资料可以在《2016/2018年候选队论文集》以及 Min-25的博客 中找到。在博客中还有另一个比本文方法复杂度稍劣的方法,并且和本文方法一样没有很大的效率提升

下面我们将介绍一种 O ( n 2 / 3 log ⁡ − 1 n ) O(n^{2/3} \log^{-1} n) O(n2/3log1n)的求 S n , f S_{n, f} Sn,f的方法。设 f = ∏ ^ p ∈ P f p f = \hat{\prod}_{p \in P} f_p f=^pPfp,并定义 f D = ∏ ^ p ∈ D f p f_{\mathcal{D}} = \hat{\prod}_{p \in \mathcal{D}} f_p fD=^pDfp,那么:

定理一: 给定 f p f_p fp,可以在 O ( n log ⁡ p n ) O(\sqrt n \log_p n) O(n logpn)的时间内由 S n , g S_{n, g} Sn,g推得 S n , g ∗ f p S_{n, g * f_p } Sn,gfp

证明: 暴力枚举 f p f_p fp中的每一项并更新即可。

因此,若我们对于 p ≤ n 1 / α ​ p \le n^{1/\alpha} ​ pn1/α使用该方法更新,其复杂度为 O ( α n 1 / α + 1 / 2 ) ​ O(\alpha n^{1/\alpha + 1/2}) ​ O(αn1/α+1/2)

接下来我们将 f = ∏ ^ f p ​ f = \hat{\prod} f_p ​ f=^fp按照 p ​ p ​ p分成三个部分: D 1 = P ∩ [ 1 , n 1 / α ] ​ \mathcal{D}_1 = P \cap [1, n^{1/\alpha}]​ D1=P[1,n1/α], D 2 = P ∩ ( n 1 / α , n 1 / 2 ] ​ \mathcal{D}_2 = P \cap (n^{1/\alpha}, n^{1/2}] ​ D2=P(n1/α,n1/2], D 3 = P ∩ ( n 1 / 2 , n ] ​ \mathcal{D}_3 = P \cap (n^{1/2}, n] ​ D3=P(n1/2,n]。若我们求出了 S n , f D 2 ​ S_{n, f_{\mathcal{D}_2}} ​ Sn,fD2 S n , f D 3 ​ S_{n, f_{\mathcal{D}_3}} ​ Sn,fD3,由之前的性质,我们可以在 O ( α n 1 / α + 1 / 2 ) ​ O(\alpha n^{1/\alpha + 1/2}) ​ O(αn1/α+1/2)的时间内将其合并为 S n , f ​ S_{n, f} ​ Sn,f。首先考虑如何求 S n , f D 2 ​ S_{n, f_{\mathcal{D}_2}} ​ Sn,fD2

定理二: ∀ 1 ≤ k ≤ log ⁡ 2 n \forall 1 \le k \le \log_2 n 1klog2n,恰含 k k k个素因子的素因子大小至少为 n 1 / α n^{1/ \alpha} n1/α的数的个数是 O ( ( α + 1 ) k n log ⁡ − 1 n ) O((\alpha + 1)^k n \log^{-1} n) O((α+1)knlog1n)的。

证明: k = 1 k = 1 k=1时该式成立。若当 k = m k = m k=m时该式成立,考虑 k = m + 1 k = m + 1 k=m+1时这样的数的个数显然为
O ( ( α + 1 ) m ∫ n 1 / α n d x ln ⁡ x n x ln ⁡ x + ( α + 1 ) m ∫ n 1 / α n d x ln ⁡ x ) = O ( ( α + 1 ) m + 1 n log ⁡ − 1 n ) O\left((\alpha + 1)^m\int_{n^{1/\alpha}}^n \frac{\mathrm{d} x}{\ln x} \frac{n}{x \ln x} + (\alpha + 1)^m\int_{n^{1/\alpha}}^n \frac{\mathrm{d} x}{\ln x} \right) = O((\alpha + 1)^{m + 1} n \log^{-1} n) O((α+1)mn1/αnlnxdxxlnxn+(α+1)mn1/αnlnxdx)=O((α+1)m+1nlog1n)
传统的计算 S n , f D 2 ​ S_{n, f_{\mathcal{D}_2}} ​ Sn,fD2的方法总要考虑素因子的顺序,这导致了状态的两个维度的产生。即便每个状态的转移已经最优,复杂度也难以优化。因此我们尝试减少这一部分的状态数:如果我们不考虑素因子的顺序会发生什么呢?也就是说,令 g k ( x ) = f D 2 ( x ) [ ∑ e x , i = k ] ​ g_k(x) = f_{\mathcal{D}_2}(x) [\sum e_{x,i} = k] ​ gk(x)=fD2(x)[ex,i=k],我们希望通过计算 g k − 1 ​ g_{k - 1} ​ gk1 ∑ p ∈ D 2 f p ​ \sum_{p \in \mathcal{D}_2} f_p ​ pD2fp的卷积来推得 S n , g k ​ S_{n, g_k} ​ Sn,gk

定理三: α \alpha α为常数,则计算一次这样的卷积的时间复杂度为 O ( n 2 / 3 log ⁡ − 1 n ) O(n^{2/3} \log^{-1} n) O(n2/3log1n)

证明: 由于只在素数及其幂次处有值,可以注意到 S n , ∑ p ∈ D 2 f p S_{n, \sum_{p \in \mathcal{D}_2} f_p} Sn,pD2fp中不同的值的期望个数是 O ( n 1 / 2 log ⁡ − 1 n ) O(n^{1/2} \log^{-1} n) O(n1/2log1n)的。因次在计算 S n , g k − 1 ∗ ∑ p ∈ D 2 f p S_{n, g_{k - 1} * \sum_{p \in \mathcal{D}_2} f_p} Sn,gk1pD2fp时,对一个 m m m单次枚举的复杂度可以做到 O ( m 1 / 2 log ⁡ − 1 m ) O(m^{1/2} \log^{-1} m) O(m1/2log1m):我们仍然按照 m \sqrt m m 分段,较大部分时我们的素数个数是除了一个 log 的,而较小时值的个数也有限。而由定理二,若 α \alpha α为常数,则计算较小部分的复杂度也可以做到 O ( S log ⁡ − 1 n ) O(S \log^{-1} n) O(Slog1n),于是取 S = O ( n 2 / 3 ) S = O(n^{2/3}) S=O(n2/3)可以达到定理中的时间复杂度。

由定理三,当 α \alpha α为常数时,我们总的计算时间复杂度为 O ( n 2 / 3 log ⁡ − 1 n ) O(n^{2/3} \log^{-1} n) O(n2/3log1n)

但我们究竟算了什么呢?容易发现,对一个数 x ​ x ​ x来说,其对答案的贡献为 ( ∑ e x , i ) ! ∏ f D 2 ( p i ) e x , i ( e x , i ) ! − 1 ​ (\sum e_{x, i})! \prod f_{\mathcal{D}_2}(p_i)^{e_{x,i}} (e_{x,i})!^{-1} ​ (ex,i)!fD2(pi)ex,i(ex,i)!1,因此若 f D 2 ( p i e ) = f D 2 ( p i ) e ( e ! ) − 1 ​ f_{\mathcal{D}_2} (p_i^{e}) = f_{\mathcal{D}_2}(p_i)^e (e!)^{-1} ​ fD2(pie)=fD2(pi)e(e!)1,我们算出的值恰好就是 S n , g k ​ S_{n, g_k} ​ Sn,gk k ! ​ k! ​ k!倍,只要在最后除去即可。

注意到若我们构造另一个积性函数 f 0 f_0 f0满足 f 0 ( p i e ) = f D 2 ( p i ) e ( e ! ) − 1 f_0(p_i^e) = f_{\mathcal{D}_2}(p_i)^e (e!)^{-1} f0(pie)=fD2(pi)e(e!)1,那么 f D 2 ∗ f 0 − 1 f_{\mathcal{D}_2} * f_0^{-1} fD2f01不为 0 0 0的位置一定是 Powerful Number 集合 W W W的子集。于是我们可以先求出 f 0 f_0 f0的前缀和,再和一个 Powerful 的函数做卷积,这样我们就在 O ( n 2 / 3 log ⁡ − 1 n ) O(n^{2/3} \log^{-1} n) O(n2/3log1n)的时间内求出了 S n , f D 2 S_{n, f_{\mathcal{D}_2}} Sn,fD2

那么如何求 S n , f D 3 S_{n, f_{\mathcal{D}_3}} Sn,fD3呢,我们将问题转化为 ∣ T ∣ |T| T个求所有最小素因子超过 n \sqrt n n 的数的 k k k次方之和的问题,令 u k ( x ) = x k [ x ∈ P ] [ x > n ] u_k(x) = x^k [x \in P] [x > \sqrt n] uk(x)=xk[xP][x>n ],即求 S n , u k S_{n, u_k} Sn,uk。注意到超过 n \sqrt n n 的因子至多只有一个,令 v k ( x ) = x k [ ∀ p i > n , ∃ e x , i = 0 ] v_k(x) = x^k [\forall p_i > \sqrt n, \exists e_{x, i} = 0] vk(x)=xk[pi>n ,ex,i=0],则有: S n , u k ∗ ( v k + ϵ ) = S n , i d k S_{n, u_k * (v_k + \epsilon)} = S_{n, id^k} Sn,uk(vk+ϵ)=Sn,idk,而 S n , u k S_{n, u_k} Sn,uk也只有 n \sqrt n n 以上的项有值,因此有 S n , v k + ϵ S_{n, v_k + \epsilon} Sn,vk+ϵ S n , i d k S_{n, id^k} Sn,idk后可以用一个简单的容斥做到 O ( n log ⁡ n ) O(\sqrt n \log n) O(n logn) S n , u k S_{n, u_k} Sn,uk。我们可以利用之前的方法求出 S n , v k S_{n, v_k} Sn,vk,这样我们这一部分的时间复杂度也是 O ( n 2 / 3 log ⁡ − 1 n ) O(n^{2/3} \log^{-1} n) O(n2/3log1n)

不妨取 α = 6 \alpha = 6 α=6,总的时间复杂度即为 O ( n 2 / 3 log ⁡ − 1 n ) O(n^{2/3} \log^{-1} n) O(n2/3log1n)

下面是一份示例代码,不保证没有错误、也希望有人能写一个常数更小一点的版本。

const int mod = 1e9 + 7;
inline void inc(int &a, int b){ (a += b) < mod ?: a -= mod; }
inline void inc(int &a, i64 b){ a = (a + b) % mod; }
inline void dec(int &a, int b){ (a -= b) >= 0 ?: a += mod; }
inline void dec(int &a, i64 b){ inc(a, mod - b % mod); }
int fpow(int a, int t){
	static int r;
	for (r = 1; t; t >>= 1, a = (i64)a * a % mod) if (t & 1) r = (i64)r * a % mod;
	return r;
}
const int unit = 7;
i64 x;
int B, S, RS, size, count, fac[unit], ifac[unit], inv[unit];
vector<int> primes, pi;
struct array_t {
	int f[unit];
	array_t(){ memset(f, 0, unit * 4); }
	int &operator[](int i){ return f[i]; }
	int to_val(){
		int sum = 0;
		for (int i = 0; i < unit; ++i) inc(sum, (i64)f[i] * ifac[i]);
		return sum;
	}
} ;
inline void inc(array_t &a, array_t &b){
	for(int i = 0; i < unit; ++i) if (b[i]) inc(a[i], b[i]);
}
inline void inc(array_t &a, array_t &b, int v){
	for(int i = 0; i + 1 < unit; ++i) if (b[i]) inc(a[i + 1], (i64)b[i] * v);
}
template<class value_t>
struct data_t {
	vector<value_t> lsum, ssum;
	data_t(){ lsum.resize(size), ssum.resize(size); }
	void update(int p, vector<int> g){
		int i, j, e, l = S / p;
		for (i = 1; i <= S; ++i) dec(ssum[i], lsum[S]);
		for (i = S; i >= 2; --i) dec(lsum[i], lsum[i - 1]);
		for (i = 1; i <= S; ++i) for (j = i, e = 0; j <= l; inc(ssum[i], (i64)ssum[j *= p] * g[++e]));
		for (i = 2; i <= S; ++i) dec(ssum[i - 1], ssum[i]);
		for (i = S; i >= 1; --i) if (lsum[i]) {
			for (j = i, e = 0; j <= l; inc(lsum[j *= p], (i64)lsum[i] * g[++e]));
			for (j = x / p / j; j; j /= p) inc(ssum[j], (i64)lsum[i] * g[++e]);
		}
		for (i = 2; i <= S; ++i) inc(lsum[i], lsum[i - 1]);
		inc(ssum[S], lsum[S]);
		for (i = S; i >= 2; --i) inc(ssum[i - 1], ssum[i]);
	}
	void gsum(){
		int i;
		for (i = 2; i <= S; ++i) inc(lsum[i], lsum[i - 1]);
		inc(ssum[S], lsum[S]); if (S == RS) lsum[S] = ssum[S];
		for (i = S; i >= 2; --i) inc(ssum[i - 1], ssum[i]);
	}
	void gdif(){
		int i;
		for (i = 2; i <= S; ++i) dec(ssum[i - 1], ssum[i]);
		dec(ssum[S], lsum[S]);
		for (i = S; i >= 2; --i) dec(lsum[i], lsum[i - 1]);
	}
	void add(i64 a, value_t b){ inc(a <= S ? lsum[a] : ssum[x / a], b); }
	void print(){
		debug("function address %llu: \n", u64(this));
		for (int i = 1; i <= S; ++i) debug("S(%d) = %d\n", i, lsum[i]);
		for (int i = S; i >= 1; --i) debug("S(%llu) = %d\n", x / i, ssum[i]);
	}
} ;
void init(){
	int i, j;
	S = sqrtl(x), RS = x / S, size = S + 1, B = powl(x, 1.0 / unit) + 1e-7, pi.resize(size);
	vector<char> sieve(size);
	for (fac[0] = i = 1; i < unit; ++i) fac[i] = (i64)fac[i - 1] * i;
	for (inv[1] = 1, i = 2; i < unit; ++i) inv[i] = mod - (i64)(mod / i) * inv[mod % i] % mod;
	for (ifac[0] = i = 1; i < unit; ++i) ifac[i] = (i64)ifac[i - 1] * inv[i] % mod;
	for (i = 4; i <= S; i += 2) sieve[i] = true;
	for (i = 9; i <= S; i += 6) sieve[i] = true;
	for (i = 5; i * i <= S; i += 2) if (!sieve[i]) {
		for (j = i * i; j <= S; j += 6 * i) sieve[j] = true;
		for (j = i * (i + 6 - 2 * (i % 3)); j <= S; j += 6 * i) sieve[j] = true;
	}
	primes.pb(0);
	for (i = 2; i <= S; ++i) {
		pi[i] = pi[i - 1] + !sieve[i];
		if (!sieve[i]) primes.pb(i), ++count;
	}
	primes.pb(size);
}
data_t<int> solve(vector<vector<int>> f){
	vector<vector<int>> g(count + 1), h(count + 1);
	int i, j, k, l, v, m;
	for (i = 1; i <= count; ++i) {
		m = f[i].size(), g[i].resize(m), h[i].resize(m);
		if (i > pi[B]) {
			for (g[i][0] = j = 1; j < m; ++j) g[i][j] = (i64)g[i][j - 1] * f[i][1] % mod;
			for (j = 0; j < m; ++j) g[i][j] = (i64)g[i][j] * ifac[j] % mod;
		} else g[i] = f[i];
		for (j = 0; j < m; ++j) for (h[i][j] = f[i][j], k = 0; k < j; ++k) dec(h[i][j], (i64)h[i][k] * g[i][j - k]);
	}
	data_t<int> F, G;
	vector<int> fval(count + 1);
	for (i = 1; i <= count; ++i) fval[i] = g[i][1];
	for (i = 1; i <= count; ++i) inc(fval[i], fval[i - 1]);
	{
		data_t<array_t> ts;
		i64 block = std::max((i64)S, std::min(x, i64(powl(x, 2.0 / 3) * 9)));
		block = x / (x / block);
		vector<i64> nxt(size);
		function<void(int, i64, int, int)> dfs = [&](int pos, i64 cur, int cnt, int now){
			int i, u, l, r;
			for (i = pos; cur * primes[i] * primes[i] <= block; ++i) {
				i64 t = cur, a;
				int e = 0, b;
				while (t * primes[i] * primes[i] <= block) {
					t *= primes[i], ++e;
					dfs(i + 1, t, cnt + e, (i64)now * g[i][e] % mod);
					a = t * primes[i], b = (i64)now * g[i][e + 1] % mod * fac[cnt + e + 1] % mod;
					if (a <= S) nxt[a] = a, inc(ts.lsum[a][cnt + e + 1], b);
					else inc(ts.ssum[x / a][cnt + e + 1], b);
				}
			}
			now = (i64)now * fac[cnt + 1] % mod;
			for (i = pos; cur * primes[i] <= S; ++i) {
				u = cur * primes[i], nxt[u] = u;
				inc(ts.lsum[u][cnt + 1], (i64)now * (fval[i] + mod - fval[i - 1]));
			}
			i64 m = x / cur;
			l = i, r = count;
			if (cur * S > block) r = pi[block / cur];
			for (u = l; u <= r; ++u) inc(ts.ssum[m / primes[u]][cnt + 1], (i64)now * (fval[u] + mod - fval[u - 1]));
		} ;
		dfs(pi[B] + 1, 1, 0, 1);
		++ts.lsum[nxt[1] = 1][0];
		ts.gsum();
		for (i = 1; i <= S; ++i) if (!nxt[i]) nxt[i] = nxt[i - 1];
		l = x / block - 1;
		for (i = l; i >= 1; --i) {
			i64 m = x / i;
			array_t sum;
			sum[0] = 1;
			for (j = pi[B] + 1; (i64)i * primes[j] <= S && j <= count; ++j) inc(sum, ts.ssum[i * primes[j]], fval[j] + mod - fval[j - 1]);
			for (; (i64)primes[j] * primes[j] <= m && j <= count; ++j) inc(sum, ts.lsum[m / primes[j]], fval[j] + mod - fval[j - 1]);
			for (; j <= count; j = v + 1) {
				k = nxt[m / primes[j]];
				v = (i64)k * S > m ? pi[m / k] : count;
				inc(sum, ts.lsum[k], fval[v] + mod - fval[j - 1]);
			}
			ts.ssum[i] = sum;
		}
		for (i = 1; i <= S; ++i) G.lsum[i] = ts.lsum[i].to_val(), G.ssum[i] = ts.ssum[i].to_val();
	}
	for (i = pi[B]; i >= 1; --i) G.update(primes[i], g[i]);
	{
		vector<pair<i64, int>> vh;
		function<void(int, i64, int)> search = [&](int pos, i64 rem, int now){
			vh.emplace_back(x / rem, now);
			while ((i64)primes[pos] * primes[pos] <= rem) {
				i64 t = rem / primes[pos] / primes[pos];
				int e = 2;
				for (; t; t /= primes[pos], ++e) if (h[pos][e]) search(pos + 1, t, (i64)now * h[pos][e] % mod);
				++pos;
			}
		} ;
		search(1, x, 1);
		sort(vh.begin(), vh.end());
		for (i = 0; i + 1 < vh.size(); ++i) if (vh[i].first == vh[i + 1].first) inc(vh[i + 1].second, vh[i].second), vh[i].second = 0;
		vector<pair<i64, int>> vv;
		for (i = 0; i < vh.size(); ++i) if (vh[i].second) vv.pb(vh[i]);
		vh.swap(vv), vv.clear();
		vh.emplace_back(x + 1, 0);
		vector<int> hsum(RS + 1), hcnt(size), nxt(size);
		for (auto it : vh) inc(it.first > S ? hsum[x / it.first] : hcnt[it.first], it.second);
		for (i = RS - 1; i >= 1; --i) inc(hsum[i], hsum[i + 1]);
		for (i = 1; i <= S; ++i) inc(hcnt[i], hcnt[i - 1]); 
		nxt[S] = S;
		for (i = S - 1; i >= 1; --i) nxt[i] = hcnt[i] == hcnt[i + 1] ? nxt[i + 1] : i;
		i64 low = pow(x, 0.6);
		low = x / (x / low);
		for (j = 0; vh[j].first <= low; ++j) {
			for (i = 1; i * vh[j].first <= S; ++i) inc(F.lsum[i * vh[j].first], (i64)vh[j].second * (G.lsum[i] + mod - G.lsum[i - 1]));
			for (i64 it = x / vh[j].first; i * vh[j].first <= low && i <= S; ++i) inc(F.ssum[it / i], (i64)vh[j].second * (G.lsum[i] + mod - G.lsum[i - 1]));
		}
		F.gsum(), low = std::min(x / low, S + 1LL);
		for (i = S; i >= low; --i) for (j = 0; i * vh[j].first <= S; ++j) inc(F.ssum[i], (i64)vh[j].second * (G.ssum[i * vh[j].first] + mod - G.lsum[S]));
		for (i = low - 1; i >= 1; --i) {
			int w = 0;
			i64 m = x / i;
			for (j = 0; i * vh[j].first <= S; ++j) inc(w, (i64)vh[j].second * G.ssum[i * vh[j].first]);
			for (; vh[j].first <= S; ++j) inc(w, (i64)vh[j].second * G.lsum[m / vh[j].first]);
			for (j = 1; i * (j + 1) <= RS && j <= S; ++j) inc(w, (i64)G.lsum[j] * (hsum[i * j] + mod - hsum[i * (j + 1)]));
			if (j <= S) inc(w, (i64)G.lsum[j] * hsum[i * j]);
			F.ssum[i] = w;
		}
	}
	return F;
}
struct lagrange{
	vector<int> init, coef;
	int n;
	lagrange(vector<int> init): init(init){
		int i;
		coef.resize(n = init.size());
		vector<int> ifac(n), inv(n);
		for (ifac[0] = i = 1; i < n; ++i) {
			inv[i] = i > 1 ? mod - (i64)inv[mod % i] * (mod / i) % mod : 1;
			ifac[i] = (i64)ifac[i - 1] * inv[i] % mod;
		}
		for (i = 0; i < n; ++i) {
			coef[i] = (i64)init[i] * ifac[n - i - 1] % mod * ifac[i] % mod;
			if ((n & 1) == (i & 1)) coef[i] = mod - coef[i];
		}
	}
	int get(i64 w){
		int x = w % mod; 
		if (x < n) return init[x];
		int res = 0, a = 1, i;
		for (i = 0; i < n; ++i) {
			res = ((i64)res * (x - i) + (i64)a * coef[i]) % mod;
			a = (i64)a * (x - i) % mod;
		}
		return (res + mod) % mod;
	}
};
vector<vector<int>> generate(function<int(int, int)> func){
	int i, e;
	i64 t;
	vector<vector<int>> res(count + 1);
	for (i = 1; i <= count; ++i) for (res[i].pb(t = e = 1); t * primes[i] <= x; t *= primes[i], ++e) res[i].pb(func(primes[i], e));
	return res;
}
vector<vector<int>> conv(vector<vector<int>> a, vector<vector<int>> b){
	int i, j, k, l;
	decltype(a) c(count + 1);
	for (i = 1; i <= count; ++i) for (c[i].resize(l = a[i].size()), j = 0; j < l; ++j) for (k = 0; j + k < l; ++k) inc(c[i][j + k], (i64)a[i][j] * b[i][k]);
	return c;
}
data_t<int> primesum(int s){
	int i, j, e, L = std::max(s + 2, S);
	i64 t;
	vector<int> coef(L + 1);
	for (i = 1; i <= L; ++i) coef[i] = 1;
	for (i = 2; i <= L; ++i) if (coef[i] == 1) for (coef[i] = fpow(i, s), j = 2; i * j <= L; ++j) coef[i * j] = (i64)coef[i] * coef[j] % mod;
	vector<int> psum(coef.begin(), coef.begin() + s + 3);
	for (i = 1; i < psum.size(); ++i) inc(psum[i], psum[i - 1]);
	lagrange sol(psum);
	data_t<int> fw = solve(generate([=](int p, int e){ return fpow(p, e * s); }));
	for (i = 1; i <= S; ++i) fw.ssum[i] = (sol.get(x / i) + mod - fw.ssum[i]) % mod, fw.lsum[i] = 0;
	for (i = S; i >= 1; --i) for (j = 2; i * j <= S; ++j) dec(fw.ssum[i], (i64)coef[j] * fw.ssum[i * j]);
	return fw;
};
template<class value_t>
data_t<value_t> work(data_t<value_t> a, const function<int(int)> &op){
	for (int i = 1; i <= S; ++i) a.ssum[i] = op(a.ssum[i]), a.lsum[i] = op(a.lsum[i]);
	return a;
}
template<class value_t>
data_t<value_t> work(data_t<value_t> a, data_t<value_t> b, const function<int(int, int)> &op){
	data_t<value_t> c;
	for (int i = 1; i <= S; ++i) {
		c.ssum[i] = op(a.ssum[i], b.ssum[i]);
		c.lsum[i] = op(a.lsum[i], b.lsum[i]);
	}
	return c;
}
data_t<int> polysum(vector<int> poly){
	int n = poly.size(), i, j;
	vector<int> num(n + 2);
	for (i = 1; i <= n + 1; ++i) {
		for (j = n - 1; j >= 0; --j) num[i] = ((i64)num[i] * i + poly[j]) % mod;
		num[i] = (num[i] + num[i - 1]) % mod;
	}
	lagrange sol(num);
	data_t<int> res;
	for (i = 1; i <= S; ++i) res.lsum[i] = sol.get(i), res.ssum[i] = sol.get(x / i);
	return res;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值