闲话 22.12.17
闲话
好那杜教筛的内容到这里就告一段落了。
如果有想学习而且不嫌我写的垃圾的同学可以经下方链接阅读。
前传 上篇 下篇
贝尔级数
写到这突然发现之前浅学了一点贝尔级数。
贝尔级数总之就是描述一个积性函数在质数幂次处的取值,可以用来构造杜教筛需要的函数。
仍然记 Dirichlet 卷积为 \(*\),积性函数的数乘为 \(\times\),质数集为 \(\mathbb P\)。
记某小写字母表示的积性函数的贝尔级数为该字母大写对应的相关符号。对于特殊数论函数的贝尔级数用该函数本身替代 \(F\)。若该函数有上下标则在该函数全体外加括号。
我们定义一个积性函数 \(f\) 关于质数 \(p\) 的贝尔级数 \(F_p(x)\) 按如下方式给出:
然后我们就能在某个质数幂次的取值处研究积性函数了。
\(\text{Theorem 1}\)
两个数论函数相卷积,其贝尔级数相乘。
证明:
令 \(h = f * g\),其中 \(f,g\) 为积性函数。
\(\text{Theorem 2}\)
对于完全积性函数 \(f\),有 \(F_p(x) = \frac 1 {1-f(p)x}\)。
证明:
\(\text{Theorem 3}\)
对于积性函数 \(f,g\),若 \(\forall p \in \mathbb P,\ F_p(x) = G_p(x)\),则 \(f = g\)。
证明:
令 \(n = \prod p_i ^{c_i}\)。有
\(\text{Theorem 5}\)
对于积性函数 \(f\) 和数论函数 \(g\),若 \(\forall p,n \ge 1, \ f(p^{n+1}) = f(p) f(p^n) - g(p) f(p^{n-1})\),则
\[F_p(x) = \frac{1}{1 - f(p)x + g(p)x^2} \]
证明:
\(f(p^{n+1})x^n = f(p) f(p^n)^n - g(p) f(p^{n-1})^n\)
两边对 \(n > 0\) 求和,有
整理可以得到
一些贝尔级数:
- \(e\)
\(e_p(x) = \sum_{i=0} [p^i=1] x^i = 1\) - \(\text I\)
\(\text I_p(x) = \sum_{i=0} x^i = \frac 1 {1-x}\) - \(\text{id}_ k\)
\((\text{id}_ k)_p(x) = \sum_{i=0} p^{ki}x^i = \frac 1 {1 - p^k x}\) - \(\mu = \text I^{-1}\)
\(\mu_p(x) = 1-x\) - \(\mu^2\)
\((\mu^2)_p(x) = 1+x\) - \(\varphi = \text{id} * \mu\)
\(\varphi_p(x) = \frac 1 {1 - p x} \times (1-x) = \frac {1-x} {1 - p x}\) - \(\text d = \text I * \text I\)
\(\text d_p(x) = \frac 1 {(1-x)^2}\) - \(\sigma_k = \text I * \text{id}_k\)
\((\sigma_k)_p(x) = \frac 1 {1-x}\times \frac 1 {1 - p^k x} = \frac 1 {(1-x)(1 - p^kx)}\) - \(\lambda\)
\(\lambda_p(x) = \frac 1 {1+x}\)
就是这一类的东西吧。
用在杜教筛里面?首先记住上面这些贝尔级数,这些级数都是好求前缀和的。然后把题目给你的式子写成贝尔级数,随后用上面的级数乘题目的式子,化简成上面的形式之积/分母的形式简单的函数。
看几个例子吧。
例 \(1\):求 \(F(n) = \sum_{d|n} \mu^2 (d) \frac nd \mu\left(\frac nd\right)\) 的前缀和。
容易发现这个东西是 \(\mu^2\) 和 \(\text{id}\times \mu\) 的卷积。前面一个好做,但是后面一个咋办呢?不难发现 \((\text{id}\times \mu)_p(x) = 1 - px\)。
因此 \(F_p(x) = (1 + x)(1 - px)\)。注意到 \(\text{id}_p(x) = \frac {1}{1 - px}\),因此构造 \(g = \text{id}\) 即可。
例 \(2\):\(f(1) = 1, \ f(p^c) = p^c + (-1)^c\),且 \(f\) 为积性函数。求 \(f\) 的前缀和。
先别急着用 \(\text{Min_25}\) 筛(
不难发现,\(F_p(x) = \text{id}_p(x) + \lambda_p(x) - e_p(x)\),这是由于 \(c^k\) 会为倍数贡献 \(k\) 个 \(-1\),在 \(1\) 处特判即可。
这样可以写出 \(F_p(x) = \frac 1 {1 - p x} + \frac 1 {1+x} - 1 = \frac{1 - px^2}{(1 - px)(1 + x)}\)。
发现分母是例 \(1\) 中的贝尔级数,因此卷上例 \(1\) 中的函数即可。我们构造 \(g(n) = \sum_{d|n} \mu^2 (d) \frac nd \mu\left(\frac nd\right), h(n) = [n = 1] + -p \times [n = p^2]\) 即可。
说了这么多,有什么实际作用吗?不知道但是挺帅的
例题
\(T\) 组数据,给定 \(n,m\),求
\[\sum_{i=1}^n \sum_{j=1}^m \varphi(ij) \pmod{998244353} \]\(1 \le T, n, m, \le 10^5\)。
这题和杜教筛没啥关系。
考虑用来拆开 \(\varphi\) 的第二个式子:
带入能得到
我们设 \(f(n) = \sum_{d|n} \frac{d}{\varphi(d)} \mu\left(\frac nd\right)\),这个东西可以枚举倍数模拟 Dirichlet 卷积,在 \(O(n\log n)\) 的复杂度内得到 \(1\sim n\) 的值。
这个东西没有什么好办法,只能暴力求。但是暴力铁定 T 掉。
咋办呢?观察咱们能用的方法。
第一种是对每个 \(n / T\) 和 \(m / T\) 处理前缀和,时间挺好的,但是空间铁炸。
第二种是枚举 \(T\),三部分分别预处理。这样空间随随便便压下来,但是没法数论分块,时间炸掉。
欸你看到这俩做法,想到了什么?
对,阈值分治后即可解决本题。
\(\text{Bonus : }\) 你是否能写出 \(F_p(x)\) 的封闭形式?
\(\text{Bonus2 : }\) 本题又名毒瘤之神的考验。
给定 \(n\)。求
\[\sum_{i=1}^n \sum_{j=1}^i \frac{\text{lcm}(i, j)}{\gcd(i, j)} \pmod{998244353} \]\(n \le 10^9\)。
如题面。我们不妨先考虑两个 \(\sum\) 上界都是 \(n\) 的情况,得到 \((ans + n) / 2\) 就是答案。设 \(S(n) = n(n+1) / 2\),直接化简:
然后考虑 \(f(n) = \sum_{d | n} \mu(d)d^2\) 如何求解。容易发现 \(f = (\mu \times \text{id}_2) * \text I\)。
应用贝尔级数,\(F_p(x) = \frac{1 - p^2x}{1 - x}\)。这样我们只需要乘上 \(\frac 1{1 - p^2x}\) 就能得到 \(\frac{1}{1 - x}\) 了。也就是说 \(f * \text{id}_2 = \text I\)。
欸好简单就构造出来了呢!
我就说吧我就说吧!
于是应用杜教筛与数论分块即可解决问题。\(f\) 的计算可以采用 Dirichlet 卷积。
给定 \(n\)。求
\[\sum_{i=1}^n \gcd\left(\lfloor i^{1/3} \rfloor, i\right) \pmod{998244353} \]\(n \le 10^{30}\)。
为啥这题不叫 problemprovidercreep 啊?
1e30,这你肯定不能枚举 \([gcd(i, j) = x]\) 的套路了。如果有数论巨佬用枚举 \(\gcd\) 的方法爆踩了我那受我一拜顺便教教我。
考虑一个 \(\lfloor i^{1/3} \rfloor\) 对应的 \(i\) 其实是很长一段。这些段只有最后一段是残的,因此我们可以将答案表示为
分开求解。
先出一个结论吧。
这个简单的,不展开了。
然后第一部分可以直接套上去。设 \(m = \lfloor n^{1/3} \rfloor - 1\),
采用两种杜教筛(\(\sum \varphi, \sum i\times\varphi\))可完成这部分。总时间复杂度为 \(O(n^{2/9})\)。
然后第二部分可以直接套上去。设 \(k = \lfloor n^{1/3} \rfloor\)。
这个可以直接朴素求解。具体地,我们枚举每个 \(d | n\)。对于每个 \(d\),它的质因子必定是 \(n\) 的质因子。因此预处理 \(n\) 的质因子后枚举计算 \(\varphi\) 即可。
跑得挺快。听说预处理前 \(10^7\) 的 \(\varphi\) 后总时间复杂度就是 \(O(n^{2/9})\) 的了。没证。
\(T\) 组询问,给定 \(n,a,b\),求
\[f(n,a,b) = \sum_{i=1}^n\sum_{j=1}^i \gcd(i^a - j^a, i^b - j^b)\times [\gcd(i,j) = 1] \pmod{1000000007} \]\(1 \le n, a, b \le 10^9\)。保证 \(a,b\) 互质,且 \(n > 10^6\) 的询问不会超过 \(10\) 次。
首先需要证明一个结论:
对于 \(i < j, \gcd(i, j) = 1\),有 \(\gcd(i^a - j^a, i^b - j^b) = i^{\gcd(a, b)} - j^{\gcd(a, b)}\)。
证明:
不妨假设 \(a \ge b\),令 \(a = kb + r\)。
提取 \((i^b - j^b)\) 项,我们有 \(i^a - j^a = (i^b - j^b)(i^{a-b} + i^{a - 2b}j^b + \cdots + i^{r}j^{a - b - r}) + i^rj^{a-r} - j^a\)。
原式 \(= \gcd(i^rj^{a-r} - j^a, i^b - j^b) = \gcd(j^{a - r}(i^r - j^r), i^b - j^b)\)
先不管 \(i^r - j^r\)。我们着眼于 \(\gcd(j^{a - r}, i^b - j^b) = \gcd(j^kb, i^b - j^b)\)。
仍然提取 \((i^b - j^b)\) 项,我们有 \(j^{kb} = -(i^b - j^b)(j^{(k-1)b} + i^bj^{(k-2)b} + \cdots + a^{(k-1)b}) + a^{kb}\)。
因此 \(\gcd(j^kb, i^b - j^b) = \gcd(i^kb, i^b - j^b) = d\)。由于 \(\gcd(i, j) = 1\) 且 \(d|i, d|j\),\(d\) 必等于 \(1\)。因此这部分不需要在意。
原式 \(= \gcd(j^{a - r}(i^r - j^r), i^b - j^b) = \gcd(i^r - j^r, i^b - j^b)\)。
根据更相减损术,结论成立。
因此可以使用杜教筛解决。
HDU 看不了代码
#include <bits/stdc++.h>
using namespace std;
template <typename T1, typename T2> T1 min(T1 a, T2 b) { return a < b ? a : b; }
template <typename T1, typename T2> T1 max(T1 a, T2 b) { return a > b ? a : b; }
#define rep(i,s,t) for (register int i = (s), i##_ = (t) + 1; i < i##_; ++ i)
#define pre(i,s,t) for (register int i = (s), i##_ = (t) - 1; i > i##_; -- i)
const int N = 1e6 + 10;
int T, n, a, b, L = 1e6;
#define S(n) (((n) * (n) * (n) - (n)) * 166666668)
modint _f[N], mu[N];
int prime[N], cnt;
bool vis[N];
void sieve() {
mu[1] = 1;
rep(i, 2, L) {
if (!vis[i]) prime[++cnt] = i, mu[i] = mod - 1;
rep(j, 1, cnt) {
int t = i * prime[j];
if (t > L) break;
vis[t] = 1;
if (i % prime[j] == 0) break;
mu[t] -= mu[i];
}
}
rep(i,1,L) _f[i] = mu[i] * i + _f[i - 1];
}
modint mp[N];
int stk[N], top;
modint getv(modint x) {
if (x <= L) return _f[*x];
if (mp[n / (*x)] != -1) return mp[n / (*x)];
modint ret = 1;
for (modint l = 2, r; l <= x; l = r + 1) {
r = x / (x / l);
ret -= (r + l) * (r - l + 1) * 500000004 * getv(x / l);
} stk[++ top] = n / (*x);
return mp[n / (*x)] = ret;
}
modint f(int n, int a, int b) {
rep(i,1,top) mp[stk[i]] = -1;
top = 0;
modint ans = 0;
for (modint l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
ans += S(n / l) * (getv(r) - getv(l - 1));
} return ans;
}
signed main() {
cin >> T;
sieve(); rep(i,0,100000) mp[i] = -1;
while (T --) {
cin >> n >> a >> b;
cout << f(n, a, b) << '\n';
}
}
以下是博客签名,与正文无关。
请按如下方式引用此页:
本文作者 joke3579,原文链接:https://www.cnblogs.com/joke3579/p/chitchat221217.html。
遵循 CC BY-NC-SA 4.0 协议。
请读者尽量不要在评论区发布与博客内文完全无关的评论,视情况可能删除。