Posted 算法9 minutes read (About 1384 words)

Min_25筛

前言

大家都会啊，还是一次考试题的部分分，听 ftq 讲了

咕了几天就来学了

好像比杜教筛少很多思维难度吧

update:

好难啊，重新理解了好多次，还改了文章结构

简介

如果一个积性函数 \(f(i)\) 在 \(i\) 是质数时是一个关于 \(i\) 的低次多项式，并且在质数的幂处的能快速求，那么大概可以用Min_25筛来求

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

对于 \(f(1)\) 特判

时间复杂度 \(\mathcal O(\frac{n^\frac{3}{4}}{\log n})\)，空间复杂度\(\mathcal O(\sqrt{n})\)

并且同时我们可以求出每个 \(\lfloor \frac{n}{x} \rfloor\) 的 \(\sum\limits_{i=2}^{\lfloor \frac{n}{x} \rfloor}[\text{$i$ 是质数}]f(i)\)

前置技能：埃拉托斯特尼筛法

埃拉托斯特尼筛法的第 \(i\) 轮找出第 \(i\) 个质数，并把所有 \(i\) 的倍数筛去

对所有 \(\le \sqrt{n}\) 的质数执行完之后就可以结束

事实上被 \(i\) 第一次筛去的 \(i\) 的倍数 \(x\) 必然满足 \(x \ge i^2\)

算法流程

第一部分 处理质数

我们先来处理每个 \(\lfloor \frac{n}{x} \rfloor\) 的

\[\sum_{i=2}^{\lfloor \frac{n}{x} \rfloor}[\text{$i$ 是质数}]f(i)\]

对于质数 \(p\)，\(f(p)\) 可以被表示为若干 \(p^k\) 的和，所以我们只考虑 \(f(p)=p^k\) 的情况

令 \[g(a,b)=\sum_{i=2}^a [\text{$i$ 是质数 或 $pmin_i>prime_b$}] * i^k\]

其中 \(pmin_i\) 表示 \(i\) 的最小质因子，\(prime_b\) 表示第 \(b\) 个质数

直观地说，就是对 \(2..a\) 做埃拉托斯特尼筛法 \(b\) 轮后剩下的所有数的 \(k\) 次幂和，包括 \(prime_1,..,prime_b\) 这些处理过的质数

需要求的就是

\[\sum_{i=2}^n [\text{$i$ 是质数}]f(i)=g(n,\infty)\]

我们从小到大枚举质数，从 \(g(* ,b-1)\) 推出 \(g(* ,b)\)

\[g(a,0)=\sum_{i=2}^a i^k\]

\[ g(a,b)= \begin{cases} g(a,b-1), & a<prime_b^2\\ g(a,b-1)-prime_b^k\left(g(\lfloor \frac{a}{prime_b} \rfloor,b-1)-g(prime_{b-1},b-1)\right), & a\ge prime_b^2 \end{cases} \]

  • 若 \(a<prime_b^2\)，所有 \(a\) 以内的合数都被筛掉了，不会再产生影响

  • 若 \(a\ge prime_b^2\)，此时我们考虑一轮埃拉托斯特尼筛法会筛掉的数

    每一个 \(\le \lfloor \frac{a}{prime_b} \rfloor\) 的且不含有 \(<prime_b\) 的质因子的数的 \(prime_b\) 倍，都会被筛掉

    但是 \(g(\lfloor \frac{a}{prime_b} \rfloor,b-1)\) 包含了 \(<prime_b\) 的质数，把这些影响消去即可

    事实上做完 \(b-1\) 轮的时候 \(\le prime_{b-1}\) 的只剩下质数了，于是可以直接取 \(g(prime_{b-1},b-1)\)

可以发现 \(a\) 的取值都形如 \(\lfloor \frac{n}{x} \rfloor\)，其中 \(prime_{b-1}<\sqrt{n}\)，所以没有影响

我们记录 \(\mathcal O(\sqrt{n})\) 个值，递推即可，复杂度 \(\mathcal O(\frac{n^\frac{3}{4}}{\log n})\)

第二部分 计算答案

令

\[S(a,b)=\sum_{i=2}^a [\text{$i$ 是质数 或 $pmin_i\ge prime_b$}]f(i)\]

显然所求

\[\sum_{i=1}^n f(i)=S(n,1)+f(1)\]

和上面类似地，但是需要枚举指数

\[ S(a,b)= \begin{cases} S(a,b+1), & a<prime_b^2 \\ \\ S(a,b+1) + \\ \sum\limits_{prime_b^{i+1} \le a} \left(f(prime_b^i)* \left(S(\lfloor \frac{a}{prime_b^i} \rfloor, b+1) - g(prime_b, \infty)\right) + f(prime_b^{i+1}) \right), & a\ge prime_b^2 \end{cases} \]

相当于每次把最小质因数为 \(prime_b\) 的合数加进来

复杂度也是 \(\mathcal O(\frac{n^{\frac{3}{4}}}{\log n})\)

第二部分的另一种实现

不记忆化，暴力递归求一个值

重新定义，令

\[S(a,b)=\sum_{i=2}^a [pmin_i\ge prime_b]f(i)\]

区别就是这里不包含额外的质数

  • 计算 \(S(a,b)\) 中质数的贡献

    即\(\sum_{i=2}^a [\text{$i\ge prime_b$ 且 $i$ 是质数}]f(i)\)

    这可以通过处理的 \(g\) 快速得到，就是\(g(a,\infty)-g(prime_{b-1},\infty)\)

    实际上这里的 \(g(prime_{b-1},\infty)\) 由于 \(prime_{b-1}\) 是很小的，下面有些代码中在筛质数的时候预处理了这部分值，本质是一样的

  • 计算合数的贡献

    我们枚举最小的质因子 \(prime_i\) 和次数 \(t\)

    有贡献

    \[\sum_{i=b}^{\infty} \sum_{t\ge 1, prime_i^{t+1}\le a}\left( S(\lfloor \frac{a}{prime_i^t} \rfloor,i+1) * f(prime_i^t) + f(prime_i^{t+1}) \right)\]

    其中后一部分是计算该质数幂的贡献，前一部分是其它的合数

    \(S( * ,i+1)\) 保证了没有 \(\le prime_i\) 的质因子，\(prime_i^{t+1}\le a\) 保证了不会越界

所以两部分加起来就好了

\[ S(a,b)= \begin{cases} 0, & a<prime_b\\ \\ g(a,\infty)-g(prime_{b-1},\infty)+ \\ \sum\limits_{i=b}^{\infty} \sum\limits_{t\ge 1, prime_i^{t+1}\le a}\left( S(\lfloor \frac{a}{prime_i^t} \rfloor,i+1) * f(prime_i^t) + f(prime_i^{t+1}) \right), &a\ge prime_b \end{cases} \]

这部分不加记忆化也跑得飞快，好像也被证明是 \(\mathcal O(\frac{n^{\frac{3}{4}}}{\log n})\) 的

但是这种方法在有些求多个值的时候不适用

实现

不是很会优化常数

  • 我们可以预处理所有的 \(\lfloor \frac{n}{x} \rfloor\) 的值，并从小到大地记为 \(a_1,..,a_m\)

    可以发现

    • 对于 \(1\le i \le \sqrt{n}\) 有 \(a_i=i\)

    • 对于其他的 \(i\) 有 \(\lfloor \frac{n}{a_i} \rfloor=m-i+1\)

    这样我们可以实现 \(\mathcal O(1)\) 从数字到编号的转换

    下面有一些代码是相反的顺序，并且使用了两个数组来实现转换

  • 在处理 \(g\) 的时候直接用转换之后的值计算，滚动第二维，空间复杂度 \(O(\sqrt{n})\)

  • 一种写法是线性筛 \(\le \sqrt{n}\) 的质数

    但是如果我们需要筛 \(f(i)=1\) 在质数处的和时，可以在第一部分中顺便处理，具体可以参考例题一的代码

例题

例题一 区间素数个数

LOJ #6235. 区间素数个数

题意

求 \(1,..,n\) 中的质数数量

\(n\le 10^{11}\)

实现

并不需要上述的第二部分求 \(S\)，令 \(f(i)=1\)，求出 \(g(n,\infty)\) 就是答案

代码中使用了两个数组 \(id1[x]\) 和 \(id2[x]\) 分别存 \(\le \sqrt{n}\) 和 \(> \sqrt{n}\) 的数字对应的编号

注意这里和上面描述的不同，\(a_1,..,a_n\) 是从大到小记录对应的数字

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#include<cstdio>
#include<math.h>

#define ll long long

const int N = 316300;
ll n, g[N<<1], a[N<<1];
int id, cnt, sn, prime[N];
inline int Id(ll x){ return x<=sn?x:id-n/x+1;}
int main() {
	scanf("%lld", &n), sn=sqrt(n);
	for(ll i=1; i<=n; i=a[id]+1) a[++id]=n/(n/i), g[id]=a[id]-1;
	for(int i=2; i<=sn; ++i) if(g[i]!=g[i-1]){
        // 这里 i 必然是质数，因为 g[] 是前缀质数个数
        // 当 <i 的质数的倍数都被筛去，让 g[] 发生改变的位置只能是下一个质数
		prime[++cnt]=i;
		ll sq=(ll)i*i;
		for(int j=id; a[j]>=sq; --j) g[j]-=g[Id(a[j]/i)]-(cnt-1);
	}
	return printf("%lld\n", g[id]), 0;
}

例题二 Sum

Luogu P4213 【模板】杜教筛（Sum）

题意

给定正整数 \(n<2^{31}\)，求

\[ \sum_{i=1}^n \varphi(i) \\ \sum_{i=1}^n \mu(i) \]

实现

令 \(g(i)\) 表示 \(i\) 以内质数的和，\(h(i)\) 表示 \(i\) 以内质数的数量

那么求 \(\varphi(i)\) 的前缀和就用 \(g(i)-h(i)\)，\(\mu(i)\) 的前缀和用 \(-h(i)\) 做就好了

直接不记忆化递归实现

求 \(\mu\) 的时候有些系数为 \(0\) 的可以不搜下去，应该会快不少

目前跑的最快

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include<cstdio>
#include<math.h>

using namespace std;
#define ll long long

const int N = 46345;
ll g[N<<1];
int T, id, sum[N], h[N<<1];
unsigned cnt, sn, n, id1[N], id2[N], prime[N], a[N<<1];
bool p[N];
inline unsigned Id(unsigned x){ return x<=sn?x:id-n/x+1;}
ll SolvePhi(unsigned a, int b){
	if(a<prime[b]) return 0;
	ll ans=g[Id(a)]-(sum[b-1]-(b-1));
	for(unsigned i=b; i<=cnt && prime[i]*prime[i]<=a; ++i){
        // 这里是强行展开了一层，可能会快一点，因为条件必然满足，事实上可以和下面的写在一起
		ans+=(SolvePhi(a/prime[i], i+1)+prime[i])*(prime[i]-1);
		for(unsigned x=prime[i]*prime[i], f=x-prime[i]; (ll)x*prime[i]<=a; x=x*prime[i], f*=prime[i])
			ans+=(SolvePhi(a/x, i+1)+prime[i])*f;
	}
	return ans;
}
int SolveMu(unsigned a, int b){
	if(a<prime[b]) return 0;
	int ans=h[Id(a)]+(b-1);
	for(unsigned i=b; i<=cnt && prime[i]*prime[i]<=a; ++i)
		ans-=SolveMu(a/prime[i], i+1);
	return ans;
}
int main() {
	scanf("%d", &T);
	while(T--){
		scanf("%u", &n), sn=sqrt(n);
		if(!n){ puts("0 0"); continue;}
		cnt=id=0;
		for(unsigned i=1; i<=n; i=a[id]+1)
			a[++id]=n/(n/i), g[id]=(ll)a[id]*(a[id]+1)/2-1, h[id]=a[id]-1;
		for(unsigned i=2; i<=sn; ++i) if(h[i]!=h[i-1]){
			prime[++cnt]=i, sum[cnt]=sum[cnt-1]+i;
			unsigned sq=i*i;
			for(int j=id; a[j]>=sq; --j){
				unsigned t=Id(a[j]/i);
				g[j]-=i*(g[t]-sum[cnt-1]);
				h[j]-=h[t]-(cnt-1);
			}
		}
		for(int i=1; i<=id; ++i) g[i]-=h[i], h[i]*=-1;
        // 上面的计算都是不考虑 1 的函数值的，要加上去
		printf("%lld %d\n", SolvePhi(n, 1)+1, SolveMu(n, 1)+1);
	}
	return 0;
}

例题三 简单的函数

LOJ #6053. 简单的函数

题意

  1. \(f(1)=1\)

  2. 对于质数 \(p\)，\(f(p^c)=p \oplus c\)，其中 \(\oplus\) 表示按位异或运算

  3. 对于互质的 \(a,b\)，\(f(ab)=f(a)f(b)\)

求 \[\sum_{i=1}^n f(i) \mod (10^9+7)\]

\(n\le 10^{10}\)

实现

一步步代进去就好了

对于质数 \(p\) 的函数值，

\[ f(p)= \begin{cases} p+1, & p=2 \\ p-1, & p\ne 2 \end{cases} \]

除了 \(2\) 的情况，这是和 \(\varphi(p)\) 一样的

质数的幂直接在枚举的时候计算

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include<cstdio>
#include<math.h>

using namespace std;
#define ll long long

const int N = 100005, P = 1000000007;
ll n, g[N<<1], h[N<<1], a[N<<1];
int sn, id, cnt, prime[N];
inline int Id(ll x){ return x<=sn?x:id-n/x+1;}
inline int S(ll a, int b){
	if(a<prime[b]) return 0;
	ll ans=g[Id(a)]-g[prime[b-1]]+P;
	for(int i=b, lim=sqrt(a); i<=cnt && prime[i]<=lim; ++i){
		ans+=(ll)S(a/prime[i], i+1)*(prime[i]^1)+(prime[i]^2);
		for(ll x=(ll)prime[i]*prime[i], j=2; x*prime[i]<=a; x*=prime[i], ++j)
			ans+=(ll)S(a/x, i+1)*(prime[i]^j)+(prime[i]^(j+1));
	}
	return ans%P;
}
int main() {
	scanf("%lld", &n), sn=sqrt(n);
	for(ll i=1; i<=n; i=a[id]+1){
		a[++id]=n/(n/i);
		ll tmp=a[id]%P;
		g[id]=tmp*(tmp+1)/2%P-1, h[id]=a[id]-1;
	}
	for(int i=2; i<=sn; ++i) if(h[i]!=h[i-1]){
		prime[++cnt]=i;
		ll sq=(ll)i*i;
		for(int j=id; a[j]>=sq; --j){
			int t=Id(a[j]/i);
			(g[j]-=i*(g[t]-g[i-1]))%=P; // 这里 g[i-1] 和 g[prime[cnt-1]] 没有区别
			h[j]-=h[t]-(cnt-1);
		}
	}
	for(int i=1; i<=id; ++i) (g[i]+=(i>1)*2+P-h[i])%=P;
	return printf("%d\n", (S(n, 1)+1)%P), 0;
}

例题四 Counting Divisors

SPOJ DIVCNT2 - Counting Divisors (square)

SPOJ DIVCNT3 - Counting Divisors (cube)

SPOJ DIVCNTK - Counting Divisors (general)

题意

令 \(\sigma_0(n)\) 表示 \(n\) 的约数个数

求 \[\sum_{i=1}^n \sigma_0(i^k)\]

对 \(2^{64}\) 取模

前两题即 \(k=2,3\) 的情况，数据范围略有不同

\(n\le 10^{11}\)

实现

差不多了

这里只贴DIVCNTK的代码

另外两题代码参考

「SPOJ DIVCNT2」Counting Divisors (square)

「SPOJ DIVCNT3」Counting Divisors (cube)

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include<cstdio>
#include<math.h>

using namespace std;
#define ll long long
#define ull unsigned ll

const int N = 100005;
int T, sn, id, cnt, prime[N];
ll n, k, a[N<<1];
ull g[N<<1];
inline int Id(ll x){ return x<=sn?x:id-(n/x)+1;}
ull S(ll a, int b){
	if(a<prime[b]) return 0;
	ull ans=g[Id(a)]-(b-1)*(k+1);
	for(int i=b, lim=sqrt(a); i<=cnt && prime[i]<=lim; ++i){
		ans+=S(a/prime[i], i+1)*(k+1)+(k*2+1);
		for(ll x=(ll)prime[i]*prime[i], j=(k*2+1); x*prime[i]<=a; x*=prime[i], j+=k)
			ans+=S(a/x, i+1)*j+j+k;
	}
	return ans;
}
int main() {
	scanf("%d", &T);
	while(T--){
		scanf("%llu%llu", &n, &k), sn=sqrt(n);
		cnt=id=0;
		for(ll i=1; i<=n; i=a[id]+1) a[++id]=n/(n/i), g[id]=(ull)(a[id]-1)*(k+1);
		for(int i=2; i<=sn; ++i) if(g[i]!=g[i-1]){
			prime[++cnt]=i;
			ll lim=(ll)i*i;
			for(int j=id; a[j]>=lim; --j) g[j]-=g[Id(a[j]/i)]-(ull)(cnt-1)*(k+1);
		}
		printf("%llu\n", S(n, 1)+1);
	}
	return 0;
}

例题五 SPOJ APS2 Amazing Prime Sequence (hard)

参考这里

例题六 51nod 1847 奇怪的数学题

参考这里

例题七 51nod 1965 奇怪的式子

参考这里

Min_25筛

https://cekavis.github.io/min-25/

Author

Cekavis

Posted on

2018-12-03

Updated on

2022-06-16

Licensed under

#

Comments

Follow
×