admin管理员组文章数量:1794759
Min
为了表述方便,下用 P \mathbb P P 表示素数集,并用 p 1 < p 2 < ⋯ < p c p_1<p_2<\cdots<p_c p1<p2<⋯<pc 表示 1 1 1 到 n n n 内的所有质数, t t t 表示第一个大于 n \sqrt n n 的质数下标.
特别的,我们令 p 0 = 1 p_0=1 p0=1.
算法介绍
一、简介
Min_25
筛是一种能快速求出某些满足特殊条件的积性函数前缀和的亚线性筛法,在通常的递归实现中复杂度为 Θ ( n 1 − ϵ ) \Theta(n^{1-\epsilon}) Θ(n1−ϵ),而经过略复杂的后缀和优化后复杂度甚至能达到 O ( n 3 4 log 2 n ) \Omicron(\frac{n^{\frac34}}{\log_2 n}) O(log2nn43) 级别.
由于后者实现较困难且前者在数据范围较小时效率更高,这里只介绍前者.
二、适用范围
要求所求函数 f ( x ) f(x) f(x) 为积性函数,同时对任意 p ∈ P , e ∈ N + p\in \mathbb{P},e\in\mathbb{N}^+ p∈P,e∈N+ 要能够快速计算 f ( p e ) f(p^e) f(pe) 的值且 f ( p e ) f(p^e) f(pe) 最好能表示成关于 p p p 的低阶多项式.
三、算法流程
为了能快速转移出所有数的和,我们设
S ( m , k ) = ∑ i = 2 m f ( i ) [ min p ∈ P , p ∣ i { p } > p k ] S(m,k)=\sum_{i = 2}^mf(i)\Big[\min_{p\in \mathbb{P},p|i}\{p\} > p_k\Big] S(m,k)=i=2∑mf(i)[p∈P,p∣imin{p}>pk]
则所求答案为 S ( n , 0 ) S(n,0) S(n,0). 注意到 i i i 是从 2 2 2 开始循环的,这是因为 1 1 1 没有最小质因子,计算时容易引起麻烦,所以我们最后再加上 f ( 1 ) f(1) f(1) 即可.
我们考虑对其递归计算,将其分为质数与合数分别计算贡献。对于合数,我们枚举剩下的数的最小质因子与其次数进行递归转移,且因为我们只算合数,只枚举枚举小于等于 n \sqrt n n 的质数即可;而质数部分似乎无法直接转移,故我们考虑再设出一个类似的函数 h ( m ) h(m) h(m) 表示 f f f 在小于等于 m m m 的所有质数处的和,那么 S S S 的转移式就非常显然了:
S ( m , k ) = h ( m ) − h ( p k ) + ∑ k < j < t ∑ p j e ⩽ m f ( p j e ) ( S ( ⌊ m p j e ⌋ , j ) + [ e ≠ 1 ] ) S(m,k) = h(m) - h(p_k) + \sum_{k<j<t}\sum_{p_j^e\leqslant m}f(p_j^e)\bigg(S\bigg(\bigg\lfloor\frac{m}{p_j^e}\bigg\rfloor,j\bigg) + [e\not=1]\bigg) S(m,k)=h(m)−h(pk)+k<j<t∑pje⩽m∑f(pje)(S(⌊pjem⌋,j)+[e=1])
观察转移式,容易发现 k k k 恒小于 t t t, 而 t t t 很小,故我们可以直接线筛对所有 k < t k< t k<t 预处理出 s k = h ( p k ) s_k=h(p_k) sk=h(pk),式子化为:
S ( m , k ) = h ( m ) − s k + ∑ k < j < t ∑ p j e ⩽ m f ( p j e ) ( S ( ⌊ m p j e ⌋ , j ) + [ e ≠ 1 ] ) S(m,k) = h(m) - s_k + \sum_{k<j<t}\sum_{p_j^e\leqslant m}f(p_j^e)\bigg(S\bigg(\bigg\lfloor\frac{m}{p_j^e}\bigg\rfloor,j\bigg) + [e\not=1]\bigg) S(m,k)=h(m)−sk+k<j<t∑pje⩽m∑f(pje)(S(⌊pjem⌋,j)+[e=1])
可以证明如果能够 O ( 1 ) \Omicron(1) O(1) 计算 h ( m ) h(m) h(m) 的值,那么剩下部分直接递归计算复杂度是 Θ ( n 1 − ϵ ) \Theta(n^{1-\epsilon}) Θ(n1−ϵ)的。
于是问题变成如何快速求 h ( m ) h(m) h(m),注意到当 S ( m , k ) S(m,k) S(m,k) 被调用时,只会递归到 S ( ⌊ m p j e ⌋ , j ) S\bigg(\bigg\lfloor\frac{m}{p_j^e}\bigg\rfloor,j\bigg) S(⌊pjem⌋,j),又根据 ⌊ ⌊ a b ⌋ c ⌋ = ⌊ a b c ⌋ \left\lfloor\frac{\left\lfloor\frac{a}{b}\right\rfloor}{c}\right\rfloor=\left\lfloor\frac{a}{bc}\right\rfloor ⌊c⌊ba⌋⌋=⌊bca⌋ 可知调用 S S S 所导致的所有递归中的 m m m 值均只与最开始的 m m m 即 n n n 有关,且均为形如 ⌊ n d ⌋ \left\lfloor\frac{n}{d}\right\rfloor ⌊dn⌋ 的数(其中 d ∈ N + , d ⩽ n d\in \mathbb N^+,d\leqslant n d∈N+,d⩽n)。设这些数组成的集合为 A A A,下面我们用 a a a 表示遍历 A A A 中所有数。
易证 ∣ A ∣ ⩽ 2 n |A| \leqslant 2\sqrt n ∣A∣⩽2n ,故考虑能否直接预处理出所有的 h ( a ) h\left(a\right) h(a) 。
尝试模仿上文中的方法,我们设一个辅助函数 g g g 来转移出需要的 h h h ,设
g ( m , k ) = ∑ i = 1 m f ( i ) [ ( i ∈ P ) ∨ ( min p ∈ P , p ∣ i { p } > p k ) ] g(m,k)=\sum_{i = 1}^mf(i)\Big[(i\in \mathbb{P}) \; \vee \;\Big(\min_{p\in \mathbb{P},p|i}\{p\} > p_k\Big)\Big] g(m,k)=i=1∑mf(i)[(i∈P)∨(p∈P,p∣imin{p}>pk)]
即用 g ( m , k ) g(m,k) g(m,k) 表示 f f f 在 1 1 1 到 m m m 中最小质因子大于 p k p_k pk 的数与所有质数处的值之和,则 h ( n ) = g ( n , t ) h(n)=g(n,t) h(n)=g(n,t) 。 考虑对 g g g 进行转移,因为 f f f 为积性函数,所以有:
g ( m , k ) = g ( m , k − 1 ) − f ( p k ) ( g ( ⌊ m p k ⌋ , k − 1 ) − g ( p k − 1 , k − 1 ) ) g(m,k)=g(m,k - 1) - f(p_k)\bigg(g\bigg(\left\lfloor\frac{m}{p_k}\right\rfloor,k - 1\bigg) - g(p_{k - 1},k - 1)\bigg) g(m,k)=g(m,k−1)−f(pk)(g(⌊pkm⌋,k−1)−g(pk−1,k−1))
显然 k ⩽ t k\leqslant t k⩽t 仍然成立,可以将 g ( p k − 1 , k − 1 ) g(p_{k - 1},k - 1) g(pk−1,k−1) 用 s k − 1 s_{k - 1} sk−1 换掉,式子变为:
g ( m , k ) = g ( m , k − 1 ) − f ( p k ) ( g ( ⌊ m p k ⌋ , k − 1 ) − s k − 1 ) g(m,k)=g(m,k - 1) - f(p_k)\bigg(g\bigg(\left\lfloor\frac{m}{p_k}\right\rfloor,k - 1\bigg) - s_{k - 1}\bigg) g(m,k)=g(m,k−1)−f(pk)(g(⌊pkm⌋,k−1)−sk−1)
观察这个式子,发现所有对 g g g 的调用仍然满足 m ∈ A m\in A m∈A,故只用 m ∈ A m\in A m∈A 处的 g g g 值即可转移出需要的 h ( a ) = g ( a , t ) h(a)=g(a,t) h(a)=g(a,t) . 于是我们只要将 O ( n ) \Omicron\left(\sqrt n\right) O(n ) 个 A A A 内的数离散化并通过转移预处理出它们的值即可,可以证明这样做时间复杂度是 O ( n 3 4 log 2 n ) \Omicron(\frac{n^{\frac34}}{\log_2 n}) O(log2nn43) 的.
这里讲一个关于离散化的小技巧:我们开两个长度为 n \sqrt n n 的数组,对于 t ⩽ n t \leqslant \sqrt n t⩽n ,我们直接以 t t t 为下标放进数组一;而对于 t > n t>n t>n,我们以 ⌊ n t ⌋ \left\lfloor\frac{n}{t}\right\rfloor ⌊tn⌋ 为下标放进数组二,查询同理,这样就避免了开长度为 n n n 的数组.
于是每次递归 S S S 时要用到的 h ( m ) h(m) h(m) 就处理完了,剩下部分直接递归即可,时间复杂度为 Θ ( n 1 − ϵ ) \Theta(n^{1-\epsilon}) Θ(n1−ϵ) 。
别忘了最后加 f ( 1 ) f(1) f(1) 。
代码实现
下面是洛谷Min25筛模板的代码
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const int maxm = 1e5 + 50,mod = 1e9 + 7,inv3 = 333333336;
int m,cnt,tot,p[maxm],vis[maxm],ind1[maxm],ind2[maxm],s1[maxm],s2[maxm],g1[2 * maxm],g2[2 * maxm];
long long n,w[2 * maxm];
inline int add(int x,int y){if(x + y < mod) return x + y;else return x + y - mod;
}
inline int dec(int x,int y){if(x - y >= 0) return x - y;else return x - y + mod;
}
inline int pos(long long x){return x <= m ? ind1[x] : ind2[n / x];
}
int S(int d,int k){if(p[k] >= w[d]) return 0;int s = dec(dec(g2[d],g1[d]),dec(s2[k],s1[k]));for(int i = k + 1; i <= tot && 1ll * p[i] * p[i] <= w[d]; i ++){long long r = p[i];for(int j = 1; r <= w[d]; j ++,r *= p[i]){int t = r % mod;s = add(s,1ll * t * dec(t,1) % mod * (S(pos(w[d] / r),i) + (j != 1)) % mod);}}return s;
}
int main(){scanf("%lld",&n);m = sqrt(n);for(int i = 2; i <= m; i ++){if(!vis[i]) p[++tot] = i,s1[tot] = add(s1[tot - 1],i),s2[tot] = add(s2[tot - 1],1ll * i * i % mod);for(int j = 1; j <= tot && i * p[j] <= m; j ++){vis[i * p[j]] = 1;if(i % p[j] == 0) continue;}}for(long long i = 1; i <= n; i = n / (n / i) + 1){long long t = n / i;w[++cnt] = t;if(t <= m) ind1[t] = cnt;else ind2[i] = cnt;t %= mod;g1[cnt] = dec(t * (t + 1) / 2 % mod,1);g2[cnt] = dec(t * (t + 1) / 2 % mod * add(2 * t,1) % mod * inv3 % mod,1);//为了方便,直接将原函数拆成两部分算,这也是通常的做法,直接将函数每项拆开且不考虑系数,最后再相加//这也是为什么要求f(x)要能表示成关于p的低阶多项式} for(int j = 1; j <= tot; j ++){for(int i = 1; i <= cnt && 1ll * p[j] * p[j] <= w[i]; i ++){int t = pos(w[i] / p[j]);g1[i] = dec(g1[i],1ll * p[j] * dec(g1[t],s1[j - 1]) % mod);g2[i] = dec(g2[i],1ll * p[j] * p[j] % mod * dec(g2[t],s2[j - 1]) % mod);}}printf("%d\n",add(S(pos(n),0),1));//加上f(1)return 0;
}
本文标签: min
版权声明:本文标题:Min 内容由林淑君副主任自发贡献,该文观点仅代表作者本人, 转载请联系作者并注明出处:http://www.xiehuijuan.com/baike/1693256544a252546.html, 本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容,一经查实,本站将立刻删除。
发表评论