「学习笔记」Min25筛

前言

周指导今天模拟赛五分钟秒第一题，十分钟说第二题是 $\text{Min25}$ 筛板子题，要不是第三题出题人数据范围给错了，周指导十五分钟就 $\text{AK}$ 了，为了向 $\text{AK}$王学习，真诚的膜拜他，接受红太阳的指导，下午就学习了一下 $\text{Min25}$ 筛。

简介

如果 $f(n)$ 是一个积性函数，且 $f(n)$ 是一个关于 $n$ 的简单多项式，并可以快速算出 $f(p^k),\ p\text{ is prime, } k \geq 0$ 的值，那么大概可以用 $\text{Min25}$ 筛在 $\mathcal O(\frac{n^\frac{3}{4}}{log_n})$ 求出 $\sum_{i=1}^nf(i)$ 的值。

算法

首先令 $s$ 为所有 $i\in[1,n], \lfloor\frac{n}{i}\rfloor$ 的集合，并试图预处理出 $\forall x \in s \sum_{i=2}^x [i\text{ is prime}] f(i)$

可能当 $x$ 不是质数的时候，$f(x)$ 不太好求，我们先假装所有数都是质数，然后用类似埃式筛法的过程把不是质数的 $f(x)$ 值给筛掉。

具体的算法之前，先定义一些东西：

令 $P=\{prime_1\dots prime_m\}$ 表示前 $m$ 个质数的集合，并且 $prime_{m+1}>\sqrt n$ 。

令 $g(n,k) ={\sum_{i=2}^n[i \text{ is prime}\text{ or } minp(i)>prime_k]}f(i) $ 表示前 $n$ 个数中满足是质数或者最小质因子大于第 $k$ 个质数的数假装它为质数求出来的 $f(x)$ 之和。

$g(n, k)$ 的本质是筛掉了前 $n$ 个数筛去前 $k$ 个质数的倍数后剩下的那些数的 $f(x)$ 之和，显然，$g(S,m)$ 就是要预处理的东西。

可以通过这个东西求 $g$ 了

\[g(n,k) = g(n,k-1)\ \ \ \ \ \ \text{if } prime_k > \sqrt n \\
g(n,k) = g(n,k-1)-f(prime_k)\times [g(\lfloor\frac{n}{prime_k}\rfloor,k-1)-\sum_{i=1}^{k-1} f(prime_i)]\ \ \ \ \ \ \text{if } prime_k \leq \sqrt n
\]

第一个转移显然，$prime_k$ 筛不掉任何数了，第二个转移考虑把所有在前 $k-1$ 轮最小质因子不为 $prime_k$ 的数贡献减掉，因为 $f(x)$ 是积性的直接搞就好了，因为 $g(n,k)$ 仍然要保留前 $k-1$ 个质数的贡献，所以还要加回来。

现在已经预处理出了 $\forall x \in s \sum_{i=2}^x [i\text{ is prime}] f(i)$ 的值，也就是我们求出了质数的答案，接下来通过从小到大加入质因子来求出合数的答案。

令 $S(n,k)=\sum_{i=2}^n [minp(i)\geq k] f(i)$ ，这里的 $f(i)$ 就是真实的值了，不假装他是质数，那么答案就是 $f(1)+S(n,1)$ 。

质数的贡献之前已经算出来了，就是 $g(n,m)-\sum_{i=1}^{k-1}f(prime_i)$。

考虑枚举每个数的最小质因子以及它们幂次得到合数的转移

\[S(n,k)=g(n,m)-\sum_{i=1}^{k-1}f(prime_i)+\sum_{j=k}^m\sum_{t=1}^{prime_{j}^{t+1}\leq \ n}S(\lfloor\frac{n}{prime_j^t}\rfloor,j+1)\times f(prime_{j}^t)+f(prime_j^{t+1})
\]

其中 $S(\lfloor\frac{n}{prime_j^t}\rfloor,j+1)\times f(prime_{j}^t)$ 算的是后面还有别的质因子的情况，$f(prime_j^{t+1})$ 算的是后面没有别的质因子的情况，可以感性理解一下。

例题

「Loj6235」区间素数个数

求 $[1,n]$ 的素数个数，$n \leq10^{11}$

令 $f(x)$ 当 $x$ 是质数的时候为 $1$ ，$x$ 是其它数的时候为 $0$ ，虽然不满足这个东西是一个关于 $x$ 的简单多项式，但是只求 $f$ 是质数的时候的答案显然是对的，$g(n,m)$ 就是答案。

code

/*program by mangoyang*/

#include <bits/stdc++.h>

#define inf (0x7f7f7f7f)

#define Max(a, b) ((a) > (b) ? (a) : (b))

#define Min(a, b) ((a) < (b) ? (a) : (b))

typedef long long ll;

using namespace std;

template <class T>

inline void read(T &x){

	int ch = 0, f = 0; x = 0;

	for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = 1;

	for(; isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;

	if(f) x = -x;

}

#define ID(x) ((x) <= T ? id[0][(x)] : id[1][n/(x)])

const int N = 1000005;

ll prime[N], id[2][N], a[N], g[N], n, m, tot, T;

int b[N];

int main(){

	read(n), T = sqrt(n);

	for(ll l = 1; l <= n; l = n / (n / l) + 1){

		a[++m] = n / l, g[m] = a[m] - 1;

		id[a[m]<=T?0:1][a[m]<=T?a[m]:n/a[m]] = m;

	}

	for(int i = 2; i <= T; i++){

		if(!b[i]) prime[++tot] = i;

		for(int j = 1; j <= tot && i * prime[j] <= T; j++){

			b[i*prime[j]] = 1;

			if(i % prime[j] == 0) break;

		}

	}

	for(int i = 1; i <= tot; i++)

		for(int j = 1; j <= m && prime[i] * prime[i] <= a[j]; j++)

			g[j] -= g[ID(a[j]/prime[i])] - (i - 1);

	cout << g[ID(n)] << endl;

	return 0;

}

「51NOD1222」最小公倍数计数

令 $f(n)=\sum_{i=1}^n\sum_{j=i}^n [\text{lcm}(i,j)=n]$，求 $\sum_{i=a}^b f(i),a\leq b\leq 10^{11}$ 。

首先先令 $f'(n)=\sum_{i=1}^n\sum_{j=1}^n[\text{lcm}(i,j)=n]$ ，显然 $f(n)=\frac{f'(n)+n}{2}$ 。

\[\sum_{i=1}^n\sum_{j=1}^n[\text{lcm}(i,j)=n]\\
=\sum_{i|n}^n\sum_{j|n}^n[\frac{ij}{\gcd(i,j)}=n] \\
=\sum_{i|n}^n\sum_{j|n}^n[ij=\gcd(ni,nj)] \\
=\sum_{i|n}^n\sum_{j|n}^n[\gcd(\frac{n}{i},\frac{n}{j})=1] \\
=\sum_{i|n}^n\sum_{j|n}^n[\gcd(i,j)=1] \\
=\sum_{d|n}2^\omega(d)
\]

最后一步考虑组合意义，每一种质因子，要么全部归 $i$ 要么全部归 $j$ ，然后你会发现，$f'(n)$ 是一个积性函数，且满足 $f(prime_i^k)=2k+1$ ，直接 $\text{Min25}$ 筛即可。

code

#include<bits/stdc++.h>

typedef long long ll;

using namespace std;

const int N = 1000005;

namespace Min25{

	#define id(x) ((x) <= T ? id1[x] : id2[n/(x)])

	ll prime[N];

	int b[N], tot, id1[N], id2[N], m;

	ll a[N], g[N], T, n;

	inline void init(){

		m = tot = 0, T = sqrt(n + 0.5);

		for(int i = 2; i <= T; i++){

			if(!b[i]) prime[++tot] = i;

			for(int j = 1; j <= tot && i * prime[j] <= T; j++){

				b[i*prime[j]] = 1;

				if(i % prime[j] == 0) break;

			}

		}

		for(ll l = 1; l <= n; l = n / (n / l) + 1){

			a[++m] = n / l;

			if(a[m] <= T) id1[a[m]] = m; else id2[n/a[m]] = m;

			g[m] = 3ll * (a[m] - 1);

		}

		for(int j = 1; j <= tot; j++)

			for(int i = 1; i <= m && prime[j] * prime[j] <= a[i]; i++){

					g[i] -= g[id(a[i]/prime[j])] - 3ll * (j - 1);

				}

	}

	inline ll solve(ll a, ll b){

		if(a < prime[b]) return 0;

		ll ans = g[id(a)] - 3 * (b - 1);

		for(int j = b; j <= tot && prime[j] * prime[j] <= a; j++)

			for(ll p = prime[j], t = 1; p * prime[j] <= a; t++, p = p * prime[j])

				ans += solve(a / p, j + 1) * (2 * t + 1) + 2 * t + 3;

		return ans;

	}

	inline ll gao(ll x){

		if(x == 0) return 0;

		if(x == 1) return 1;

		n = x, init();

		return (solve(n, 1) + 1 + n) / 2ll;

	}

}

int main(){

	ll a, b;

	cin >> a >> b;

	cout << Min25::gao(b) - Min25::gao(a - 1) << endl;

	return 0;

}

巴特西

「学习笔记」Min25筛

「学习笔记」Min25筛

前言

简介

算法

例题

最新文章

热门文章