首页 > 其他分享 >筛法

筛法

时间：2024-07-22 20:07:59浏览次数：13  
标签：lfloor right 筛法 dfrac sum aligned left

筛法

杜教筛，\(\textsf {Min-25}\) 筛

越来越边缘了啊，主要是板子，题很少

建议阅读 数论函数基础，数论分块，莫比乌斯函数，欧拉函数 章节，了解基本概念与先要知识


此处获取本节调试数据 / 代码包

全文 绝大多数 内容是对 [0] 中讲述的 粗略抄写 和 胡乱加工


1. 杜教筛

在 [0] 中被 跳过了，悲，所以这部分主要参考 [1] \(\color {black} \textsf {harryzhr}\) 老师 \(2023\) 年的一份课件

但是课件中 时间复杂度 部分是 伪证一例，故此处参考 [2] OI - Wiki 上的证明

概念

对于一个 数论函数 \(f\)，杜教筛 可以在 低于线性 的时间复杂度内计算 \(S (n) = \sum _ {i = 1} ^ n f (i)\)

也就是 更快的计算求和函数


做法

考虑找一个 积性函数 \(g\)，我们考虑 \(f, g\) 狄利克雷卷积 的 求和函数，即

\[\begin {aligned} \sum _ {i = 1} ^ n (f * g) (i) \end {aligned} \]

简单转化，有

\[\begin {aligned} \sum _ {i = 1} ^ n (f * g) (i) &= \sum _ {i = 1} ^ n { \sum _ {d \mid i} { g (d) f \left ( \dfrac i d \right ) } } \\ &= \sum _ {i = 1} ^ n { g (i) S \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

由于 \(g\) 为 积性函数，故有 \(g (1) = 1\)，于是 \(S (n) = g(1)S(n)\)，有下式

\[\begin {aligned} g(1)S(n) &= \sum _ {i = 1} ^ n { g (i) S \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } - \sum _ {i = 2} ^ n { g (i) S \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \\ &= \sum _ {i = 1} ^ n { (f * g) (i) } - \sum _ {i = 2} ^ n { g (i) S \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

设 \(f * g = h\)，那么 如果 \(h, g\) 的 前缀和 是容易求的，那么使用 数论分块 就可递归求解 \(\sum f\)

这个东西简单来讲其实就是 直到 \(g, f * g\) 的 求和函数，快速求解 \(f\) 的 求和函数


时间复杂度

若 \(h, g\) 的 前缀和 均可 \(O (1)\) 取得，由 数论分块 性质，我们知道 \(\left \lfloor \dfrac n i \right \rfloor\) 有 \(O (\sqrt n)\) 种取值

我们不妨设取值集合为 \(R (n) = \left \{ \left \lfloor \dfrac n k \right \rfloor : k \in [1, n] \right \}\)，存在如下引理

对于任意的 \(m \in R(n)\)，有 \(R (m) \subseteq R (n)\)

证明是容易的，考虑若 \(m = \left \lfloor \dfrac n x \right \rfloor\)，则任意 \(\left \lfloor \dfrac m y \right \rfloor = \left \lfloor \dfrac n {xy} \right \rfloor \in R (n)\)

这告诉我们，所有递归过程中涉及到的 \(S (x)\) 中的 \(x\)，取值均被包含在 \(R (n)\) 中

而我们若 使用记忆化，则对于每个 \(S (x)\)，我们至多计算一次

换言之，我们至多对于 \(R (n)\) 中每个值 \(x\) 计算 \(S (x)\) 就能得到最终答案，于是

\[\begin {aligned} T (n) &= \sum _ {k \in R (n)} T (k) \\ &= O (\sqrt n) + \sum _ {k = 1} ^ { \left \lfloor \sqrt n \right \rfloor } O (\sqrt k) + \sum _ {k = 2} ^ { \left \lfloor \sqrt n \right \rfloor } O \left ( \sqrt { \dfrac n k } \right ) \\ &= O \left ( \int _ 0 ^ {\sqrt n} { \left ( \sqrt x + \sqrt {\dfrac n x} \right ) dx } \right ) \\ &= O (n ^ {0.75}) \end {aligned} \]


而如果我们能 在合理的时间 \((T _ 0 (m))\) 内，预处理出一部分 \(S (k) ~ (k \in [1, m], m \ge \left \lfloor \sqrt n \right \rfloor)\)

则时间复杂度将变成

\[\begin {aligned} T (n) &= T _ 0 (m) + \sum _ {k \in R (n), k > m} T (k) \\ &= T _ 0 (m) + \sum _ {k = 1} ^ { \left \lfloor \frac n m \right \rfloor } O \left ( \sqrt { \dfrac n k } \right ) \\ &= O \left ( T_0 (m) + \int _ {0} ^ {\frac n m} { \sqrt {\dfrac n x} dx } \right ) \\ &= O \left ( T_0 (m) + \dfrac n {\sqrt m} \right ) \end {aligned} \]

若使用线性筛，即 \(T _ 0 (m) = O (m)\)，则由 均值不等式，当 \(m = n ^ {\frac 2 3}\) 时， \(T (n) _ {\min} = O (n ^ \frac 2 3)\)


Luogu P4213 【模板】杜教筛

求 \(S_1 (n) = \sum _ {i = 1} ^ n \varphi (i)\)，\(S_2 (n) = \sum _ {i = 1} ^ n \mu (i)\)

对于 \(S_1\)，考虑有 \(\varphi * 1 = \operatorname {id}\)，故有

\[\begin {aligned} S_1 (n) &= \sum _ {i = 1} ^ n { \operatorname {id} (i) } - \sum _ {i = 2} ^ n { S_1 \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \\ &= \dfrac {n (n + 1)} 2 - \sum _ {i = 2} ^ n { S_1 \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

对于 \(S_2\)，考虑有 \(\mu * 1 = \epsilon\)，故有

\[\begin {aligned} S_2 (n) &= \sum _ {i = 1} ^ n { \epsilon (i) } - \sum _ {i = 2} ^ n { S_2 \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \\ &= 1 - \sum _ {i = 2} ^ n { S_2 \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

然后递归求解即可，记得需要 记忆化

悲惨的消息是 std::unordered_map 常数太大，所以这个题 \(O (n ^ {0.75})\) 的方法是 过不去的

于是我们只能 线性筛预处理部分前缀和，然后做到 \(O (n ^ {\frac 2 3})\) 的复杂度，才能通过

#include <bits/stdc++.h>

const int MAXN = 1500005;
const int maxn = 1500000;

std::unordered_map <int32_t, int64_t> sum1, sum2;
int32_t vis[MAXN], pri[MAXN], phi[MAXN], mu[MAXN], cnt;
int64_t s1[MAXN], s2[MAXN];

inline void Sieve (const int n = maxn) {
	phi[1] = mu[1] = 1;
	for (int i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ cnt] = i, mu[i] = -1, phi[i] = i - 1;
		for (int k = 1; k <= cnt && pri[k] * i <= n; ++ k) {
			vis[i * pri[k]] = 1;
			if (i % pri[k] == 0) {
				phi[i * pri[k]] = pri[k] * phi[i];
				break ;
			}
			mu[i * pri[k]] = - mu[i];
			phi[i * pri[k]] = (pri[k] - 1) * phi[i]; 
		}
	}
	for (int i = 1; i <= n; ++ i) s1[i] = s1[i - 1] + phi[i];
	for (int i = 1; i <= n; ++ i) s2[i] = s2[i - 1] + mu [i];
}

inline int64_t Sieve1 (const int64_t n) {
	if (n <= maxn) return s1[n];
	if (sum1.count (n)) return sum1[n];
	int64_t ans = n * (n + 1) / 2;
	for (int64_t l = 2, r; l <= n; l = r + 1) 
		r = n / (n / l), ans = ans - (r - l + 1) * Sieve1 (n / l);
	return sum1[n] = ans;
}

inline int64_t Sieve2 (const int64_t n) {
	if (n <= maxn) return s2[n];
	if (sum2.count (n)) return sum2[n];
	int64_t ans = 1;
	for (int64_t l = 2, r; l <= n; l = r + 1) 
		r = n / (n / l), ans = ans - (r - l + 1) * Sieve2 (n / l);
	return sum2[n] = ans;
}

int32_t N;
int64_t x;

int main () {
	
	Sieve (), sum1[0] = sum2[0] = 0;
	
	std::cin >> N;
	
	for (int i = 1; i <= N; ++ i) {
		std::cin >> x;
		std::cout << Sieve1 (x) << ' ' << Sieve2 (x) << '\n';
	}
	
	return 0;
}

Luogu P3768 简单的数学题

并不困难，我们考虑先对式子进行转化，惯例枚举 \(\gcd\)

\[\begin {aligned} \sum _ {d = 1} ^ n { \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ n { [\gcd (i, j) = d] ij } } } \end {aligned} \]

只枚举 \(d\) 的倍数

\[\begin {aligned} \sum _ {d = 1} ^ n { d \sum _ {i = 1} ^ { \left \lfloor \frac n d \right \rfloor } { \sum _ {j = 1} ^ { \left \lfloor \frac n d \right \rfloor } { [\gcd (i, j) = 1] ijd ^ 2 } } } \end {aligned} \]

\(d ^ 2\) 放到前面去，然后使用 莫比乌斯函数 - Trick 3 - ex，有

\[\begin {aligned} \sum _ {d = 1} ^ n { d ^ 3 \sum _ {i = 1} ^ { \left \lfloor \frac n d \right \rfloor } i ^ 2 \varphi (i) } \end {aligned} \]

这时候已经结束了，考虑如果我们能快速求得 后半部分，那么整除分块就可以做

于是考虑 杜教筛，筛出 \(S (n) = \sum _ {i = 1} ^ n i ^ 2 \varphi (i)\) 的值，注意到 欧拉反演 \(\sum _ {d \mid n} \varphi (d) = n\)

我们有 \(f (i) = i ^ 2 \varphi (i)\)，可以想到设 \(g (i) = i ^ 2\)，于是

\[\begin {aligned} h (n) &= (f * g) (n) \\ &= \sum _ {d \mid n} { f (d) g \left ( \dfrac n d \right ) } \\ &= \sum _ {d \mid n} { d ^ 2 \varphi (d) \left ( \dfrac n d \right ) ^ 2 } \\ &= n ^ 3 \end {aligned} \]

牛得很，于是就可以做了

\[\begin {aligned} S (n) &= \sum _ {i = 1} ^ n h (i) - \sum _ {i = 2} ^ n { g (i) S \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \\ &= \sum _ {i = 1} ^ n i ^ 3 - \sum _ {i = 2} ^ n { i ^ 2 S \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

需求 \(\sum i ^ 2\) 与 \(\sum i ^ 3\) 的公式

但是杜教筛一次不是只能求 求和函数单点值 吗？

如果这样那么时间复杂度就是 \(O (n ^ {0.25} n ^ {\frac 2 3})\)，爆掉了

实质上，筛一次我们求出了 \(R (n)\) 这个集合内所有点的值，而这正是我们这个题需要的

故时间复杂度 \(O (n ^ {\frac 2 3} + \sqrt n)\)，可以通过本题

#include <bits/stdc++.h>

const int MAXN = 5000005;
const int PREN = 5000000;

int64_t MOD = 998244353;

std::unordered_map <int64_t, int64_t> sum;

int64_t s[MAXN], I2, I6;
int64_t vis[MAXN], pri[MAXN], phi[MAXN], cnt;

inline int64_t qpow (int64_t a, int64_t b = MOD - 2) {
	int64_t res = 1;
	a = a % MOD;
	while (b) {
		if (b & 1) res = res * a % MOD;
		a = a * a % MOD, b >>= 1;
	}
	return res;
}

inline void Sieve (const int64_t n = PREN) {
	phi[1] = 1, I2 = qpow (2), I6 = qpow (6);
	for (int64_t i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ cnt] = i, phi[i] = i - 1;
		for (int64_t k = 1; k <= cnt && pri[k] * i <= n; ++ k) {
			vis[i * pri[k]] = 1;
			if (i % pri[k] == 0) {
				phi[i * pri[k]] = pri[k] * phi[i];
				break ;
			}
			phi[i * pri[k]] = (pri[k] - 1) * phi[i]; 
		}
	}
	for (int64_t i = 1; i <= n; ++ i) s[i] = 1ll * i * i % MOD * phi[i] % MOD;
	for (int64_t i = 1; i <= n; ++ i) s[i] = (s[i - 1] + s[i]) % MOD;
} 

inline int64_t Sum1 (int64_t x) {
	x = x % MOD;
	return x * (x + 1) % MOD * I2 % MOD;
}

inline int64_t Sum2 (int64_t x) {
	x = x % MOD;
	return x % MOD * (x + 1) % MOD * (x + x + 1) % MOD * I6 % MOD;
}

inline int64_t Sum3 (const int64_t x) {
	return Sum1 (x) * Sum1 (x) % MOD;
}

inline int64_t Sieve1 (const int64_t n) {
	if (n <= PREN) return s[n];
	if (sum.count (n)) return sum[n];
	int64_t ans = Sum3 (n);
	for (int64_t l = 2, r; l <= n; l = r + 1) 
		r = n / (n / l), ans = (ans - (Sum2 (r) - Sum2 (l - 1) + MOD) % MOD * Sieve1 (n / l) % MOD + MOD) % MOD; 
	return sum[n] = ans;
}

inline int64_t Solve (const int64_t n) {
	int64_t ans = 0;
	for (int64_t l = 1, r; l <= n; l = r + 1) 
		r = n / (n / l), ans = (ans + (((Sum3 (r) - Sum3 (l - 1)) % MOD + MOD) % MOD * Sieve1 (n / l) % MOD + MOD) % MOD) % MOD;
	return ans;
}

int64_t n;

int main () {
	
	std::cin >> MOD >> n;
	
	Sieve (), Sieve1 (n);
	
	std::cout << Solve (n) << '\n';
	
	return 0;
}

Luogu P1587 [NOI2016] 循环之美

题目大意就是计数 在 \(k\) 进制下 纯循环小数 的 取值 个数

我们来分析题目中给出的条件，首先 每个取值只计算一次，我们可以转化成 计数最简分数

即对于分数 \(\dfrac x y\)，要求 \(x \perp y\)

然后考虑刻画 纯循环小数 这个性质，注意到其满足 小数部分左 / 右移 \(l\) 位后不变

即 \(x \equiv x k ^ l \pmod y\)（其中 \(k ^ l\) 表示 左移 \(l\) 位，\(\bmod y\) 表示取 小数部分）

于是有 \(k ^ l \equiv 1 \pmod y\)，我们可以证明这个条件等价于 \(k \perp y\)

设 \(\gcd (k, y) = d\)，由于 \(k ^ l \equiv 1 \pmod y\)，故有 \(k ^ l = my + 1\)（\(m \in \Z\)）

由于 \(d \mid k \and d \mid y\)，显然需要有 \(d \mid my + 1\)，即 \(d \mid 1\)，故 \(d = 1\)，\(k \perp y\)

所以这个题就等价于让我们计数满足 \(x \perp y\) 且 \(k \perp y\) 的 二元组个数，即下式

\[\begin {aligned} \sum _ {x = 1} ^ n { \sum _ {y = 1} ^ m { [\gcd (x, y) = 1] [\gcd (y, k) = 1] } } \end {aligned} \]

容易利用 莫比乌斯反演，转换枚举顺序 后得到

\[\begin {aligned} \sum _ {d = 1} ^ {\min (n, m)} { \sum _ {x = 1} ^ { \left \lfloor \frac n d \right \rfloor } { \sum _ {y = 1} ^ { \left \lfloor \frac m d \right \rfloor } { \mu (d) [\gcd (yd, k) = 1] } } } \end {aligned} \]

将无用的 \(x\) 提出，分离 \(y, d\)，于是

\[\sum _ {d = 1} ^ {\min (n, m)} { \mu (d) \left \lfloor \dfrac n d \right \rfloor [\gcd (d, k) = 1] } \sum _ {y = 1} ^ { \left \lfloor \frac m d \right \rfloor } { [\gcd (y, k) = 1] } \]

前面部分的处理没有头绪，但是后面一个 \(\sum\) 部分是简单的，设 \(g (n) = \sum _ {y = 1} ^ n [\gcd (i, k) = 1]\)

显然这东西 \([\gcd (i, k) = 1]\) 这东西有循环节长 \(k\)，容易知道其通项公式

\[\begin {aligned} g (n) &= \left \lfloor \dfrac n k \right \rfloor \varphi (k) + g (n \bmod k) \\ &= \left \lfloor \dfrac n k \right \rfloor g (k) + g (n \bmod k) \end {aligned} \]

于是我们只需要 \(O (k)\) 预处理 \(1 \sim k\) 的答案即可 \(O (1)\) 求解

然后此时我们可以将要求的式子转化成下面的形式

\[\begin {aligned} \sum _ {d = 1} ^ {\min (n, m)} { \left \lfloor \dfrac n d \right \rfloor g \left ( \left \lfloor \dfrac m d \right \rfloor \right ) \mu (d) [\gcd (d, k) = 1] } \end {aligned} \]

显然前两项可以 整除分块 处理掉，于是我们只在意 \(f (n) = \sum _ {i = 1} ^ n \mu (i) [\gcd (i, k) = 1]\) 这部分

用类似 杜教筛 的式子，我们可以把它表示成 两和式相减 的形式

\[\begin {aligned} f (n) = \sum _ {i = 1} ^ n { [\gcd (i, k) = 1] f \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } - \sum _ {i = 2} ^ n { [\gcd (i, k) = 1] f \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

前面部分展开就变成了（下式不包含后面部分）

\[\begin {aligned} \sum _ {i = 1} ^ n { [\gcd (i, k) = 1] \sum _ {j = 1} ^ { \left \lfloor \frac n i \right \rfloor } { \mu (j) [\gcd (j, k) = 1] } } \end {aligned} \]

把 \(\sum _ {j}\) 放到前面，合并 \(\gcd\) 条件，有

\[\begin {aligned} \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ { \left \lfloor \frac n i \right \rfloor } { \mu (j) [\gcd (ij, k) = 1] } } \end {aligned} \]

经典套路枚举 \(d = ij\)，\(j\) 作为因数

\[\begin {aligned} \sum _ {d = 1} ^ n { \sum _ {j \mid d} { \mu (j) [\gcd (d, k) = 1] } } \end {aligned} \]

我们在其中发现了 \(\sum _ {j \mid d} \mu (j)\)，好得很

\[\begin {aligned} \sum _ {d = 1} ^ n { [d = 1] [\gcd (d, k) = 1] } \end {aligned} \]

于是这部分就是 \(1\)，那么 \(f (n)\) 就可以表示为

\[\begin {aligned} 1 - \sum _ {i = 2} ^ n { [\gcd (i, k) = 1] f \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

这就是一个 典型的杜教筛 形式了，我们筛出 \(f (n)\)，套上前面 整除分块 即可

时间复杂度似乎带点玄学，因为 整除分块 致使我们需要 \(O (\sqrt n)\) 次 杜教筛

显然 \(O (n ^ {\frac 1 2 + \frac 2 3})\) 的复杂度 怎么想都爆了

但显然这个题仍然具有上一个题 筛一次得到多个有效值 的特性

只是在这个题中，具体 跑满筛了多少次 似乎较难刻画，所以说是带点玄学的

#include <bits/stdc++.h>

const int MAXN = 1000005;
const int PERN = 1000000;

int vis[MAXN], pri[MAXN], mu[MAXN], g[MAXN], f[MAXN], cnt;
int N, M, K;

std::unordered_map <int32_t, int64_t> sum;

inline void Sieve (const int n = PERN) {
	for (int i = 1; i <= n; ++ i) g[i] = g[i - 1] + (std::__gcd (i, K) == 1);
	mu[1] = f[1] = 1;
	for (int i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ cnt] = i, mu[i] = - 1;
		for (int k = 1; k <= cnt && pri[k] * i <= n; ++ k) {
			vis[i * pri[k]] = 1;
			if (i % pri[k] == 0) break ;
			mu[pri[k] * i] = - mu[i];
		}
		f[i] = f[i - 1] + mu[i] * (g[i] - g[i - 1]);
	}
}

inline int64_t G (const int x) {
	return x / K * g[K] + g[x % K];
}

inline int64_t F (const int n) {
	if (n <= PERN) return f[n];
	if (sum.count (n)) return sum[n];
	int64_t ans = 1;
	for (int l = 2, r; l <= n; l = r + 1) 
		r = n / (n / l), ans = ans - F (n / l) * (G (r) - G (l - 1));
	return sum[n] = ans;
}

inline int64_t Solve (const int n, const int m) {
	int64_t ans = 0;
	for (int l = 1, r; l <= std::min (n, m); l = r + 1) 
		r = std::min (n / (n / l), m / (m / l)), ans = ans + (F (r) - F (l - 1)) * (n / l) * G (m / l);
	return ans;
}

int main () {
	
	std::cin >> N >> M >> K, Sieve ();
	
	std::cout << Solve (N, M) << '\n';
	
	return 0;
}

2. \(\textsf {Min-25}\) 筛

概念

快速筛 积性函数 \(f\) 前缀和 的 算法，要求 \(f (p ^ k) ~ (p \in \mathbb P)\) 为 关于 \(p ^ k\) 的 常数次多项式

也就是 \(f (p ^ k)\) 容易求得


做法

算法流程简单来说就是 先 构造完全积性函数 \(f'\) 使得 \(\forall p ^ k, p \in \mathbb P, k \in \Z, f (p ^ k) = f' (p ^ k)\)

通过函数 \(f'\)，我们可以推出 \(f\) 在所有 质数幂 处的取值

然后构造函数 \(S\)，利用 质数幂处的取值 以及 \(g\) 函数的值求解 原函数前缀和


时间复杂度

有说 \(O \left (\dfrac {n ^ {0.75}} {\log n} \right )\) 的，但事实上在 \(OI\) 数据范围里可能比 杜教筛 的 \(O (n ^ {\frac 2 3})\) 更快一些

常数较小，可以通过 \(10 ^ {10}\) 的数据，具体证明笔者并不会，可以参考 [2] OI Wiki 上的分析


Luogu P5325 【模板】Min_25 筛

此处参考 [3]

设 完全积性函数 \(f'\) 满足上述条件，我们钦定一个新函数 \(g (n, p)\)，有

\[\begin {aligned} g (n, j) = \sum _ {i = 1} ^ n f'(i) (i \in \mathbb P \or \min _ {p \mid i} p > P _ j) \end {aligned} \]

\(P_j\) 表示第 \(j\) 个质数

那这玩意儿怎么求呢？考虑通过 \(g (n, j - 1)\) 来 转移得出

考虑 \(g (n, j)\) 与 \(g (n, j - 1)\) 的差距

容易发现，\(g (n, j - 1)\) 比 \(g (n, j)\) 多的就是 最小质因数 为 \(P_j\) 的 \(i\) 对应的 $f' (i) $ 的值

我们可以 强行钦定 这些数一定有质因子 \(P_j\)，于是把这些数 除以 \(P_j\)

而这些数剩下部分的 最小质因子 需要 大于等于 \(P_j\)，于是这些数贡献就是 \(g \left (\left \lfloor \dfrac n {P_j} \right \rfloor, j - 1 \right)\)

然后我们还需要 去掉小于 \(P_j\) 的质数 的贡献，于是需要减去 \(\sum _ {i = 1} ^ {j - 1} f' (P _ i)\)

\[\begin {aligned} g (n, j) = g (n, j - 1) - f' (P _ j) \left ( g \left ( \dfrac n {P _ j} , j - 1 \right ) - \sum _ {i = 1} ^ {j - 1} f' (P_i) \right ) \end {aligned} \]

考虑求值，设 \(S (n, j) = \sum _ {i = 1} ^ n f (i) (\min _ {p \mid i} p > P _ j)\)，我们分成 质数 和 合数 来计算贡献

在 质数 时是简单的，考虑 \(f'\) 与 \(f\) 在 质数幂 处取值相同，故质数时答案就是

\[\begin {aligned} g \left ( n , |P| \right ) - \sum _ {i = 1} ^ j f' (P _ i) \end {aligned} \]

合数时我们可以 枚举每个数的最小质因数以及次数，根据积性函数性质转移，即

\[\begin {aligned} \sum _ {k > j} { \sum _ {c = 1} ^ {P _ k ^ c \le n} { f (P _ k ^ c) \left ( S \left ( \dfrac n {P _ k ^ c} , k \right ) + [c > 1] \right ) } } \end {aligned} \]

这里 \(c > 1\) 的意义是显然的，我们 不要重复统计质数的贡献

但是我们需要统计 质数幂（但不是质数，也就是说不是一次） 的贡献

于是有 \(S\) 的转移式

\[\begin {aligned} S (n, j) = g (n, |P|) - \sum _ {i = 1} ^ j f' (P_i) + \sum _ {k > j} { \sum _ {c = 1} ^ {P _ k ^ c \le n} { f (P _ k ^ c) \left ( S \left ( \dfrac n {P _ k ^ c} , k \right ) + [c > 1] \right ) } } \end {aligned} \]

最终答案就是 \(S (n, 0)\)，但是加上 \(f (1)\)，因为 \(S\) 函数未统计 \(1\) 处的贡献

#include <bits/stdc++.h>

const int MAXN = 200005;
const int MOD  = 1000000007;
const int I2   = 500000004;
const int I6   = 166666668;

int64_t vis[MAXN], pri[MAXN], cnt = 0;
int64_t pos[MAXN], id1[MAXN], id2[MAXN], tot = 0, sqrtn;
int64_t N;
int64_t h1[MAXN], h2[MAXN], g1[MAXN], g2[MAXN];

inline void Sieve (const int n = sqrtn) {
	for (int i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ cnt] = i;
		for (int k = 1; k <= cnt && pri[k] * i <= n; ++ k) {
			vis[pri[k] * i] = 1;
			if (i % pri[k] == 0) break ;
		}
	}
	for (int i = 1; i <= cnt; ++ i) h1[i] = (h1[i - 1] + pri[i]) % MOD;
	for (int i = 1; i <= cnt; ++ i) h2[i] = (h2[i - 1] + pri[i] * pri[i] % MOD) % MOD;
}

inline int32_t Id (const int64_t x) {
	return x <= sqrtn ? id1[x] : id2[N / x];
}

inline int64_t S1 (int64_t x) {
	x = x % MOD;
	return x * (x + 1) % MOD * I2 % MOD;
}

inline int64_t S2 (int64_t x) {
	x = x % MOD;
	return x * (x + 1) % MOD * (x + x + 1) % MOD * I6 % MOD;
}

inline int64_t Fp (int64_t n) {
	n = n % MOD;
	return n * (n - 1) % MOD;
}

inline int64_t G (const int64_t n) {
	return (g2[Id (n)] - g1[Id (n)] + MOD) % MOD;
}

inline int64_t H (const int64_t n) {
	return (h2[n] - h1[n] + MOD) % MOD;
}

inline int64_t S (const int64_t n, const int k) {
	if (pri[k] > n) return 0;
	int64_t ans = (G (n) - H (k) + MOD) % MOD;
	for (int i = k + 1; i <= cnt && pri[i] * pri[i] <= n; ++ i)
		for (int64_t p = pri[i]; p <= n; p *= pri[i])
			ans = (ans + Fp (p) * (S (n / p, i) + (p > pri[i])) % MOD) % MOD;
	return ans;
}

int main () {
	
	std::cin >> N, sqrtn = sqrt (N), Sieve ();
	
	for (int64_t l = 1, r; l <= N; l = r + 1) {
		r = N / (N / l), pos[++ tot] = N / l;
		g1[tot] = S1 (pos[tot]) - 1, g2[tot] = S2 (pos[tot]) - 1;
		pos[tot] <= sqrtn ? id1[pos[tot]] = tot : id2[N / pos[tot]] = tot;
	}
	
	for (int i = 1; i <= cnt; ++ i) {
		for (int j = 1; j <= tot && pri[i] * pri[i] <= pos[j]; ++ j) {
			g1[j] = (g1[j] - pri[i] * (g1[Id (pos[j] / pri[i])] - h1[i - 1]) % MOD + MOD) % MOD;
			g2[j] = (g2[j] - pri[i] * pri[i] % MOD * (g2[Id (pos[j] / pri[i])] - h2[i - 1]) % MOD + MOD) % MOD;
		}
	}
	
	pri[cnt + 1] = 1e18;
	
	std::cout << (S (N, 0) + 1) % MOD << '\n';
	
	return 0;
}

3. 引用资料

[0] Number Theory —— H_W_Y

[1] 筛法 —— harryzhr

[2] 杜教筛 —— OI Wiki

[3] Min_25 筛学习笔记 —— zhoukangyang

标签：lfloor,right,筛法,dfrac,sum,aligned,left
From： https://www.cnblogs.com/FAKUMARER/p/18316794

相关文章

  • 筛法学习笔记
    0.更新upd2023.5.21更新了关于powerfulnumber数量的证明upd2023.5.25更新了关于杜教筛的时间复杂度证明正文1.筛质数筛法其实就是判断质数的一个算法，但是是解决\([1,n]\)这一段区间的算法筛质数是最简单的一个用法1.1暴力最简单的方式就是对于每一个数去判断......
  • Min25 筛法
    之前学习的的确是太浅薄了，于是重新学习一下。可以做什么？对于满足条件的积性函数\(f(n)\)，求其前\(n\)项和\(\sum_{i=1}^nf(i)\)。需要准备些什么？设\(N=\{x|x=\lfloor\frac{n}{i}\rfloor\}\)，即为所有整除的不同点值。设\(B=\sqrtn,p_{k}\)代表\(\leB\)的所有质......
  • 取素数优化——埃拉托斯特尼筛法（Sieve of Eratosthenes）
    埃拉托斯特尼筛法（SieveofEratosthenes）是一种用来生成一定范围内所有素数的算法。其基本思想是从小到大遍历每个数，如果当前数是素数，则将其所有的倍数标记为非素数。这个过程中，所有未被标记为非素数的数即为素数。下面是使用埃拉托斯特尼筛法来计算区间[x,y]内的素数个数的修......
  • 算法课程笔记——素数朴素判定&埃氏筛法
    算法课程笔记——素数朴素判定&埃氏筛法sqrt返回浮点数，而且这样可防溢出优化i*i会更快......
  • 筛法学习笔记
    埃氏筛枚举质数\(p_i\)，每次去除所有是\(p_i\)倍数的数，总效率大概是\(O(n\log\logn)\)。int_prm=0,prm[M];boolisprm[M];voidGet_phi(intn){ for(inti=2;i<=n;i++){ if(isprm[i])continue; prm[++_prm]=i; for(intj=i*i;j<=n;j+=i)isprm[j]=1; }}值......
  • 《算法笔记》系列----质数的判断（埃氏筛法）
    目录一、朴素算法二、埃氏筛法1、与朴素算法对比2、算法介绍   3、例题即代码实现一、朴素算法 从素数的定义中可以知道，一个整数n要被判断为素数，需要判断n是否能被2.3.n-1中的一个整除。只2，3..n-1都不能整除n，n才能判定为素数，而只要有一个能整除n的数出现，n就......
  • C++实现欧拉筛法
    Euler筛法介绍以筛出100以内（含100）的所有素数为例来说明一下欧拉筛法的原理。和Eratosthenes筛法一样，Euler筛法也从2开始筛，但Eratosthenes筛法会把2的倍数一批全部筛掉，而Euler筛法用2筛时仅仅把2*2（即4）筛掉，而把其它偶数留到后面再筛掉，从而避免了一个偶数被多次筛除带来的性能开销，......
  • 浅谈筛法
    浅谈筛法Euler筛Eratosthenes筛可以证明（不是“不难证明”），Eratosthenes筛的复杂度为\(\Theta(n\log\logn)\)。Eratosthenes筛的复杂度证明Dirichlet前缀和以P5495【模板】Dirichlet前缀和为例。给定\(a_1,a_2,\cdots,a_n\)，求\(b_1,b_2,\cdots,b_n\)，满足\[b_i......
  • 筛法思想的题目
    这道题目比较经典，或者说这种思想比较经典。这种筛法的思想。我们正着想对于每一个\(n、n-1、n-2、...、2、1\)都分解一遍质因数显然是来不及的时间复杂度达到\(O(n\sqrt{n})\)我们考虑对于每一个1e6以内的质因数的个数跑了一下程序是\(78498\)个素数定理告诉我们不超过x......
  • 质数基础筛法
    目录埃氏筛线性筛埃氏筛埃氏筛是一种筛素数的方法，埃氏筛的思想很重要，主要是时间复杂度朴素的埃氏筛的时间复杂度是\(O(nlogn)\)这个复杂度是调和级数vector<int>p;intvis[N];voidsolve(){ rep(i,2,n){ if(!vis[i]) p.pb(i); for(intj=i+i;j<=n;j+=i) vis[j]=1; ......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99