Min_25 筛学习笔记

质数前缀统计

求出 \(n\) 以内所有质数的 \(c\) 次方之和。

  考虑埃氏筛,每次用一个质数枚举其的所有倍数,筛去所有不合法的数。

一些定义以及解释

  • \(p_{\min}(x)\):\(x\) 的最小质因子,\(p_{\min}(1) = +\infty\),如果 \(x\) 非质数,那么 \(p_{\min}(x) \le \sqrt{x}\)。
  • \(m:\left\lfloor \sqrt{n}\right\rfloor\) 以内质数的个数。
  • \(p_k\):第 \(k\) 个质数。
  • \(P_k\):前 \(k\) 个质数构成的集合。
  • \(S(n, k)\):\(n\) 以内的数用前 \(k\) 个质数删除了不合法的数后,剩下的数的集合,也就是 \(S(n, k) = \{x \in \mathbb{N}^* | x \le n \wedge p_{\min}(x) \ge p_k\}\)。同时,我们注意到,对于 \(S(n, m)\) 就是所有还没有遍历的质数。
  • \(D(n, k)\):\(n\) 以内的数中用第 \(k\) 个质数删去的数的集合,即 \(D(n, k) = S(n, k - 1) - S(n, k)\)。
  • \(f(n, k)\):\(\displaystyle \sum_{x \in S(n, k)} x^c\),于是题目要求的答案就是 \(f(n, m) - 1 + \displaystyle \sum_{x \in P_m} x^c\)。

推导过程

\(S(n, k) = S(n, k - 1) - p_k S(\left\lfloor \frac{n}{p_k}\right\rfloor , k - 1)\)

  也就是考虑被 \(p_k\) 删除的数,这个可以根据 \(p_k S(\left\lfloor \frac{n}{p_k}\right\rfloor , k - 1)\) 与 \(D(n, k)\) 是个双射来证明。

  于是,我们可以得到 \(f\) 的计算式子:

\(f(n, k) = f(n, k - 1) - p_k^c f(\left\lfloor\frac{n}{p_k}\right\rfloor, k - 1)\)

  如果我们对于每个质数都直接用上面的式子计算一遍,那么总共有 \(\frac{\sqrt{n}}{\log n}\) 个质数,而根据整除分块,我们实际上要计算的值有 \(\sqrt{n}\) 个,于是复杂度就是 \(\mathcal O(\frac{n}{\log n})\) 了!

  我们可以考虑优化,我们会注意到一个这样的事情,对于 \(k > m\),\(S(n, k) = S(n, k - 1) - \{p_k\}\),因为 \(S\) 还在改变,于是我们还要计算那些 \(k > m\) 的值,如果我们让 \(S\) 不改变,那么是否会有一些优化呢?

  我们可以定义一些新的东西:

  • \(T(n, k) = S(n, k) +P_k\)
  • \(g(n, k) = f(n, k) + \displaystyle \sum_{x \in P_m} x^c\),最后答案就是 \(g(n, m) - 1\)

  于是我们可以利用之前 \(S\) 的式子得到:

\(T(n, k) = T(n, k - 1) - p_k \left(T(\left\lfloor \frac{n}{p_k}\right\rfloor , k - 1) - P_{k - 1}\right) + \{p_k\}\)

  又我们注意到,\(T(p_k - 1, k - 1)\) 就是 \(P_{k - 1} + {1}\),因此我们也可以得到 \(g\) 的式子:

\(g(n, k) = g(n, k - 1) - \left(g(\left\lfloor \frac{n}{p_k}\right\rfloor , k - 1) - g(p_k - 1, k - 1)\right)\)

  至于 \(g(n, 0) = \sum_{i \le n} i^c\),可以在 \(\mathcal O(c)\) 的时间里用拉格朗日插值求出。

  然后复杂度就是 \(\mathcal O(\dfrac{n^{3 / 4}}{\log n})\),不太会证。

  模板 tea:Luogu5493 质数前缀统计

#define in read<int>()
#define lin read<ll>()
#define rep(i, x, y) for(int i = (x); i <= (y); i++)
#define per(i, x, y) for(int i = (x); i >= (y); i--)

using namespace std;

using ll = long long;

template < typename T > void chkmax(T &x, const T &y) { x = x > y ? x : y; }
template < typename T > void chkmin(T &x, const T &y) { x = x < y ? x : y; }

const int N = 1e6 + 10;
const int inf = 0x7fffffff;

int mod, C;
ll lim, n;

ll g0[N], g1[N];
// g0[i] 是 g(i), g1[i] 是 g(n / i)
ll qp(ll x, int t) { x %= mod; ll res = 1; for(; t; t = t >> 1, x = x * x % mod) if(t & 1) res = res * x % mod; return res; }

namespace calcg {
	const int N = 20;
	int y[N], fac[N], ifac[N], pre[N], suf[N], m;
	void init() {
		m = C + 2;
		fac[0] = 1; rep(i, 1, m) fac[i] = 1ll * fac[i - 1] * i % mod; 
                ifac[m] = qp(fac[m], mod - 2);
		per(i, m - 1, 0) ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
		rep(i, 1, m) y[i] = qp(i, C), y[i] = (y[i] + y[i - 1]) % mod;
	}
	int calc(ll v) {
		v %= mod;
		pre[0] = 1; rep(i, 1, m) pre[i] = 1ll * pre[i - 1] * (v - i + mod) % mod;
		suf[m + 1] = 1; per(i, m, 1) suf[i] = 1ll * suf[i + 1] * (v - i + mod) % mod;
		ll res = 0;
		rep(i, 1, m) {
			res += 1ll * pre[i - 1] * suf[i + 1] % mod * ifac[i - 1] % mod
                         * ifac[m - i] % mod * y[i] % mod * (m - i & 1 ? mod - 1 : 1) % mod;
			res %= mod;
		} return res;
	}
}

using calcg :: calc;

int main() {
	n = lin, C = in, mod = in, lim = sqrtl(n);
	calcg :: init();
	rep(i, 1, lim) g0[i] = (calc(i) - 1) % mod, g1[i] = (calc(n / i) - 1) % mod;
	rep(i, 2, lim) {
		g0[i] = (g0[i] + mod) % mod; if(g0[i] == g0[i - 1]) continue;
		ll pc = qp(i, C), pg = g0[i - 1];
		int l1 = lim / i, l2 = min(n / i / i, lim); chkmin(l1, l2);
		rep(j, 1, l1) g1[j] = (g1[j] - pc * (g1[j * i] - pg)) % mod;
		ll tt = n / i;
		if(tt < inf) { // 卡常
			int t = tt;
			rep(j, l1 + 1, l2) g1[j] = (g1[j] - pc * (g0[t / j] - pg)) % mod;
		} else {
			ll t = tt;
			rep(j, l1 + 1, l2) g1[j] = (g1[j] - pc * (g0[t / j] - pg)) % mod;
		}
		ll st = 1ll * i * i;
		per(j, lim, st) g0[j] = (g0[j] - pc * (g0[j / i] - pg)) % mod;
	}
	ll ans = 0;
	rep(i, 1, lim) {
		g1[i] = (g1[i] + mod) % mod;
		ans += 1ll * i * i % mod * g1[i] % mod; ans %= mod;
	} printf("%lld\n", ans);
	return 0;
}

  当然,这个算法可以进行一波魔改,然后就可以做 LOJ6028「from CommonAnts」质数计数 II 了,代码

Min_25 筛

给出一个积性函数 \(F\):

  • \(F(p)\) 为关于 \(p\) 的简单多项式。
  • \(F(p^k)\) 能够较快速求出。
  • 求出 \(\sum_{i = 1}^n F(i)\) 。

  按照 \(k\) 从大到小的顺序,考虑 \(p_{\min} = k\) 的数,那么记 \(S(n, k) = \{x \in \mathbb{N}^* | x \le n \wedge p_{\min}(x) \ge p_k\}\),记 \(f(n, k) = \sum_{x \in S(n, k)} F(x)\),最后答案就是 \(f(n, 1)\)。

  于是有 \(\displaystyle f(n, k) = \sum_{i = k}^{p_i \le n} \sum_{v} F(p_i ^ v)f\left(\left\lfloor\frac{n}{p_i ^ v}\right\rfloor, k+ 1\right)\),当 \(p_i \ge \sqrt{n}\) 的时候,后面这个式子的值就是 \(F(p_i)\),于是可以用之前预处理的质数前缀和解决。

  复杂度为 \(\mathcal O(\dfrac{n^{3 / 4}}{\log n})\),不太会证。

  比上面的更好的写法 \(\downarrow\)

#include <bits/stdc++.h>

#define in read<int>()
#define lin read<ll>()
#define rep(i, x, y) for(int i = (x); i <= (y); i++)
#define per(i, x, y) for(int i = (x); i >= (y); i--)

using namespace std;

using ll = long long;

const int N = 1e6 + 10;
const int mod = 1e9 + 7;
const int inv2 = (mod + 1) >> 1;
const int inv3 = (mod + 1) / 3;

ll n, lim, sp1[N], sp2[N], g1[N], g2[N], w[N];
int ind1[N], ind2[N], tot;
int pri[N / 10], pnum;
bool v[N];

void shai(int l) {
	rep(i, 2, l) {
		if(!v[i]) {
			pri[++pnum] = i;
			sp1[pnum] = (sp1[pnum - 1] + i) % mod;
			sp2[pnum] = (sp2[pnum - 1] + 1ll * i * i) % mod;
		}
		rep(j, 1, pnum) {
			if(pri[j] * i > l) break;
			v[pri[j] * i] = true; if(i % pri[j] == 0) break;
		}
	}
}

int S(ll x, int y) {
	if(x < pri[y]) return 0;
	ll res = 0;
	res += g2[x <= lim ? ind1[x] : ind2[n / x]] - sp2[y - 1];
	res -= g1[x <= lim ? ind1[x] : ind2[n / x]] - sp1[y - 1];
	res = (res % mod + mod) % mod;
	rep(k, y, pnum) {
		if(1ll * pri[k] * pri[k] > x) break;
		ll cur = pri[k];
		for(int c = 1; cur <= x; cur *= pri[k], c++) {
			res += cur % mod * ((cur - 1) % mod) % mod * (S(x / cur, k + 1) + (c != 1)) % mod;
			res %= mod;
		} // 枚举计算的次数,当 c!=1 时,应该包括 1,c = 1 的时候,我们已经计算了。
	} return res % mod;
}

int main() {
	n = lin, lim = sqrtl(n);
	shai(lim * 2); // 稍微多筛一点
	for(ll i = 1; i <= n; i++) {
		++tot, w[tot] = (n / i); ll v = w[tot] % mod;
		g1[tot] = (v * (v + 1) / 2 - 1) % mod;
		g2[tot] = (v * (v + 1) % mod * (v * 2 + 1) % mod * inv2 % mod * inv3 % mod - 1) % mod;
		if(w[tot] <= lim) ind1[w[tot]] = tot; else ind2[n / w[tot]] = tot;
		i = n / (n / i);
	}
	rep(i, 1, pnum) { 
		int p = pri[i]; ll x1 = sp1[i - 1], x2 = sp2[i - 1];
		rep(x, 1, tot) {
			if(1ll * p * p > w[x]) break;
			int y = w[x] / p <= lim ? ind1[w[x] / p] : ind2[n / (w[x] / p)];
			g1[x] = (g1[x] - 1ll * p * (g1[y] - x1) % mod + mod) % mod;
			g2[x] = (g2[x] - 1ll * p * p % mod * (g2[y] - x2) % mod + mod) % mod;
		}
	}
	printf("%d\n", (S(n, 1) + 1) % mod);
	return 0;
}
上一篇:【论文笔记】Adversarial Training for Weakly Supervised Event Detection


下一篇:关系