LOJ #6202. 叶氏筛法(min_25 筛)

题意

求 $[L, R]$ 之间的素数之和 . $L≤10^{10},2×10^{10} \le R \le 10^{11}$

题解

一个有点裸的 min_25筛 ? 现在我只会筛素数的前缀和 , 合数的过几天再学吧 .

首先推荐一波 yyb大佬博客 ~~这个人很强 , 别那么fake就好啦~~

令 $F(x) = x$ 显然此处 $F(x)$ 是完全积性函数 .

我们需要求的就是 $\displaystyle \sum_{i=1}^{n} [i \in Prime] F(i)$ .

这个就是 min_25筛的预处理部分 啦.

由 yyb博客可得 .

对于下面这个表达式 ,

g (n, j) = \sum_{i = 1}^{n} i [i \in P o r min (p) > P_{j}, p | i, p \in P]

$\displaystyle g(n,j)=\sum_{i=1}^ni[i\in P\ or \min(p)>P_j,p|i,p\in P\ ]$

有如下一个递推式 :

g (n, j) = {\begin{cases} g (n, j - 1) & P_{j}^{2} > n \\ g (n, j - 1) - P_{j} [g (\frac{n}{P_{j}}, j - 1) - g (P_{j} - 1, j - 1)] & P_{j}^{2} \leq n \end{cases}

$g(n,j)=\begin{cases} g(n,j-1)&P_j^2\gt n\\ g(n,j-1)-P_{j}[g(\frac{n}{P_j},j-1)-g(P_j-1,j-1)]&P_j^2\le n\end{cases}$

这个直接实现是 $\displaystyle O(\frac{n^{\frac{3}{4}}}{\log n})$ 的复杂度 qwq

然后考虑代码实现 .

首先预处理前 $\sqrt n$ 的素数 , (假设处理到了 $lim$ )

以及 $i=1 \sim lim$ 的 $F(x)$ 前缀和 . (此处可适当处理多一点 , 时间效率会提高)

然后假设我们当前考虑的是 $g(n, m)$ .

直观上共有三步 .

$P_{m+1}^2 > n$ 且 $n \le lim$ , 那么此时可以直接用之前预处理的答案 , 因为此时存在有贡献的数只可能为素数 .
将 $m$ 一直缩小到 $n < P_m^2$ . 这个利用了递推式 $P_j^2 > n$ 时 $g(n,j) = g(n, j - 1)$ 的那个定理 .
然后不断将 $m$ 下降, 直至下降到 $0$ , 此间需要要递归计算 $P_{j}[g(\frac{n}{P_j},j-1)+g(P_j-1,j-1)]$ 的值 , 其中后者 $g(P_j-1,j-1)$ 可以用前面预处理 , 递归下去也只会有一层 . 而前者会不断递归计算 . 其中 $m$ 到 $0$ 的时候 , 就是边界情况 : $\displaystyle g(n, 0)=\sum_{i=2}^{n} i$ . 这个是很显然的 , 因为此时所有数都计入了贡献 , 但 $1$ 是无法给予这个贡献的 . (一开始此处调试许久... )

代码

递归版本

直接按前面的式子模拟即可。

#include <bits/stdc++.h>
#define For(i, l, r) for(register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for(register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Set(a, v) memset(a, v, sizeof(a))
#define debug(x) cout << #x << ':' << x << endl
using namespace std;

void File() {
#ifdef zjp_shadow
	freopen ("6202.in", "r", stdin);
	freopen ("6202.out", "w", stdout);
#endif
}

const int N = 1e7 + 1e3, Lim = 1e7;

typedef __int128 ll;

inline ll read() {
    ll x = 0, fh = 1; char ch = getchar();
    for (; !isdigit(ch); ch = getchar()) if (ch == '-') fh = -1;
    for (; isdigit(ch); ch = getchar()) x = (x << 1) + (x << 3) + (ch ^ 48);
    return x * fh;
}

inline void Out(ll x, bool fir = true) { 
	if (!x) { if (fir) puts("0"); return ; }Out(x / 10, false); putchar (x % 10 + 48); if (fir) putchar ('\n'); 
}

int prime[N], cnt = 0; bitset<N> is_prime;

ll sump[N];

void Init(int maxn) {
	is_prime.set(); is_prime[0] = is_prime[1] = false;
	For (i, 2, maxn) {
		sump[i] = sump[i - 1];
		if (is_prime[i])
			prime[++ cnt] = i, sump[i] += i;
		For (j, 1, cnt) {
			register int res = i * prime[j];
			if (res > maxn) break;
			is_prime[res] = false;
			if (!(i % prime[j])) break;
		}
	}
}

inline ll Sum(ll x) { return x * (x + 1) / 2 - 1; }

int tot = 0;
ll Sump(ll n, int m) {
	if (n <= Lim && n < (ll)prime[m + 1] * prime[m + 1]) return sump[n];

	for (; n < (ll)prime[m] * prime[m]; -- m);

	ll res = Sum(n);
	for (; m; -- m)
		res -= (Sump(n / prime[m], m - 1) - /*Sump(prime[m] - 1, m - 1)*/ sump[prime[m] - 1]) * prime[m];

	return res;
}

inline ll Calc(ll x) { return Sump(x, cnt - 1); }

int main () {

	File();

	Init(Lim);

	ll l = read(), r = read(); Out(Calc(r) - Calc(l - 1));

	return 0;
}

非递归版本

如何做非递归呢？

考虑每次的 $j$ 都是不断减小的，就是我们会依次会除掉一个个质因子 $p_j$ 。

我们从小到大枚举每个质因子 $p_j \le \lfloor \sqrt n \rfloor$ ，然后再从小到大枚举每个 $n$ 的块 $\displaystyle \lfloor \frac{n}{x} \rfloor$ （只有 $\sqrt n$ 个）。

然后利用上面的式子直接递推即可。

具体实现的时候由于只有 $\sqrt n$ 个，我们对于 $\le \sqrt n$ 的和 $> \sqrt n$ 的进行分块即可。

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl

using namespace std;

//typedef long long ll;
typedef __int128 ll;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

inline ll read() {
    ll x(0), sgn(1); char ch(getchar());
    for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
    for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
    return x * sgn;
}

inline void Out(ll x) {
	if (!x) { puts("0"); }
	char stk[25]; int top = 0;
	for (; x; x /= 10) stk[++ top] = (x % 10) + 48;
	while (top) putchar (stk[top --]); putchar ('\n');
}

void File() {
#ifdef zjp_shadow
	freopen ("6202.in", "r", stdin);
	freopen ("6202.out", "w", stdout);
#endif
}

const int N = sqrt(1e11) * 2 + 10;

bitset<N> is_prime; 
int prime[N / 10], pcnt; ll sump[N / 10];

void Linear_Sieve(int maxn) {
	is_prime.set();
	is_prime[0] = is_prime[1] = false;

	For (i, 2, maxn) {
		if (is_prime[i]) {
			prime[++ pcnt] = i;
			sump[pcnt] = sump[pcnt - 1] + i;
		}
		for (int j = 1; j <= pcnt && 1ll * prime[j] * i <= maxn; ++ j) {
			is_prime[prime[j] * i] = false; if (!(i % prime[j])) break;
		}
	}
}

int d, id1[N], id2[N]; 

ll val[N], res[N];

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

inline ll Sum(ll n) {
	return n * (n + 1) / 2 - 1;
}

ll Min25_Sieve(ll n) {
	int cnt = 0; d = sqrt((long long)(n));
	for (ll i = 1; i <= n; i = n / (n / i) + 1)
		val[id(n / i) = ++ cnt] = n / i, res[cnt] = Sum(n / i);

	For (i, 1, pcnt)
		for (int j = 1; j <= cnt && 1ll * prime[i] * prime[i] <= val[j]; ++ j)
			res[j] -= (res[id(val[j] / prime[i])] - sump[i - 1]) * prime[i];
	return res[id(n)];
}

int main () {

	File();

	ll l = read(), r = read();

	Linear_Sieve(sqrt(r + .5));

	Out(Min25_Sieve(r) - Min25_Sieve(l - 1));

	return 0;

}

__EOF__

本文作者：zjp_shadow
本文链接：https://www.cnblogs.com/zjp-shadow/p/9190865.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！