简单的数论函数及其求和方法
数论函数
记号与约定
称函数 \(f\) 为数论函数，当且仅当其定义域为 \(\mathbb Z\)，陪域为 \(\mathbb C\)。
称数论函数 \(f\) 为积性函数当且仅当对于任意 \(a,b \in \mathbb N_+,\gcd(a,b) = 1\)，有 \(f(ab) = f(a)f(b)\)。
称数论函数 \(f\) 为完全积性函数当且仅当对于任意 \(a,b \in \mathbb N_+\)，有 \(f(ab) = f(a)f(b)\)。
按照如下记号记录部分数论函数：
- \(I(n) = 1\)
- \(\epsilon(n) = [n = 1]\)
- \(id(n) = n\)
- \(id_k(n) = n^k\)
- \(d(n) = \sum \limits_{d \mid n} 1\)
- \(\sigma(n) = \sum \limits_{d \mid n} d\)
- \(\mu,\varphi\) 定义照常，不再赘述。
其中 \(I,\epsilon,id,id_k\) 完全积性。
另外，文中部分 \(a/b,\frac{a}{b}\) 型分数根据上下文可能有向下取整（\(\left \lfloor \frac ab \right \rfloor\)）的含义，请注意甄别。
线性筛
对于任意积性函数，我们可以线性筛求出。
我们只需要关心积性函数在素数幂 \(p^k\) 处的取值。为了检验一个数 \(i\) 是否是 \(p^k\) 型的，我们可以维护 \(low(i)\) 表示 \(i\) 中最小质因子的部分。例如 \(low(6) = 2,low(49) = 49,low(12)=4\)。当且仅当 \(i = low(i)\) 时 \(i\) 是 \(p^k\) 型的。
具体的 \(k\) 和 \(p\) 也是不难维护的。
下面的模板来自 Alex-Wei。
for(int i = 2; i < N; i++) {
if(!vis[i]) pr[++pcnt] = i, f[i] = ..., low[i] = i; // 单独算 f(p)
for(int j = 1; j <= pcnt && i * pr[j] < N; j++) {
vis[i * pr[j]] = 1;
if(i % pr[j] == 0) { // i 与 p 不互质
low[i * pr[j]] = low[i] * pr[j];
if(i == low[i]) f[i * pr[j]] = ...; // i = p ^ k，单独算 f(p ^ {k + 1})
else f[i * pr[j]] = f[i / low[i]] * f[low[i * pr[j]]];
break;
}
low[i * pr[j]] = pr[j];
f[i * pr[j]] = f[i] * f[pr[j]]; // i 与 p 互质，f(ip) = f(i)f(p)
}
}
两例变量分离
-
\(d(ij)\)
\[d(ij) = \sum \limits_{x \mid i} \sum \limits_{y \mid j} [\gcd(x,y) = 1] \]构造双射容易证明。
-
\(\varphi(ij)\)
\[\varphi(ij) =\frac{\varphi(i)\varphi(j)\gcd(i,j)}{\varphi(\gcd(i,j))} \]证明：
\[\varphi(i)\varphi(j)\gcd(i,j) = \gcd(i,j) \left(i \prod \limits_{p \mid i,p \in \mathbb P} \frac{p-1}{p}\right) \left(j \prod \limits_{q \mid j,q \in \mathbb P} \frac{q-1}{q}\right) \\ = ij\gcd(i,j) \left(\prod \limits_{p \mid ij,p \in \mathbb P} \frac{p-1}{p}\right) \left(\prod \limits_{q \mid \gcd(i,j),q \in \mathbb P} \frac{q-1}{q}\right) \\ = \varphi(ij)\varphi(\gcd(i,j)) \]
关于 \(\mu^2(i)\)
有性质：
当 \(\mu(i) \neq 0\) 时，必然不存在平方因子，因此恰有 \(d = 1\) 造成 \(1\) 的贡献。
当 \(\mu(i) = 0\) 时，设 \(p_1,p_2,\cdots,p_k\) 这些质因子在 \(i\) 中的次数至少为 \(2\)。\(d\) 要在右侧造成贡献，只能是这些质因子构成的集合的某个子集的乘积，根据二项式定理，贡献为 \(0\)。
狄利克雷卷积
狄利克雷卷积
狄利克雷卷积：是一个算子。对于数论函数 \(f,g\)，我们有
狄利克雷逆：若 \(f* g = \epsilon\) 则称 \(f,g\) 互为狄利克雷逆。
对于 \(f(1) \neq 0\) 的 \(f\)，\(f\) 的狄利克雷逆一定存在。显然这条件已经必要。如下推导得到 \(g\)：
基本性质
-
交换律
\[f * g= g* f \] -
结合律
\[(f * g) * h = f* (g* h) \] -
分配律
\[(f\pm g) * h = f * h \pm g * h \] -
消去律
\[f * h = g * h \Leftrightarrow f = g \ \ (h(1) \neq 0) \] -
两个积性函数的狄利克雷卷积是积性函数
-
积性函数的狄利克雷逆是积性函数
几个简单的狄利克雷卷积
-
\(\mu * I = \epsilon\)
这也是为什么我们可以把 \([\gcd(x,y) = 1]\) 化开变成 \(\sum \limits_{d \mid \gcd(x,y)} \mu(d)\)。
-
\(\varphi * I = \text{id}\)
这也是为什么我们可以把 \(\gcd(x,y)\) 化开变成 \(\sum \limits_{d \mid \gcd(x,y)} \varphi(d)\)。
-
\(I * I = d\)
根据 \(d\) 的定义。
-
\(\text{id} * I = \sigma\)
根据 \(\sigma\) 的定义。
狄利克雷卷积的求法 (特别鸣谢: pp_orange)
\(f,g\) 已知，每一项可以 \(O(1)\) 求出。设我们要计算 \(h = f*g\) 的前 \(n\) 项。
当 \(f,g\) 是普通函数时
采用定义式 \(O(n \log n)\) 计算。
Code：[V 我 50 我帮你写]。
当 \(f\) 为积性函数时
考察 \(f = I\) 这一特例。此时等价于狄利克雷前缀和，复杂度为 \(O(n \log \log n)\)。
当 \(f\) 为任意积性函数时，我们可以执行类似的操作。考虑枚举质数 \(p\)。初始我们的卷积中只有 \(h'(n) = f(1)g(n)\) 这一项。每枚举一个质数 \(p\)，我们就枚举不含因子 \(p\) 的正整数 \(n\)，把 \(f(p)h(n),f(p^2)h(n),f(p^3)h(n),\cdots\ (p \nmid n)\) 这样的项贡献到对应的新的 \(h'(np),h'(np^2),h'(np^3),\cdots\) 上。这样我们就相当于在 \(p\) 这个维度做了一遍前缀和，现在 \(h'(n)\) 里面包含了所有 \(f(x)g(n/x)\) 的项，其中 \(x\) 整除 \(n\)，且 \(x\) 的质因子集合中只有我们枚举过的质数。
因为高维前缀和是对的，那么这个操作也应当是对的。此时复杂度为 \(\sum \limits_{1 \leq p^k \leq n,k \ge 1,p \in \mathbb P} \left \lfloor \frac{n}{p^k} \right \rfloor = O(n \log \log n)\)。这是因为对于每个 \(p\)，后面的部分都可以看作余项（等比数列求和不超过常数倍的第一项），因此其复杂度和狄利克雷前缀和 / 埃氏筛法相同。
Code：[V 我 50 我帮你写]。
当 \(f,g\) 为积性函数时
时间复杂度为 \(O(n)\)。
我们只需要求出 \(h\) 在素数幂处的取值。对于每个素数直接暴力，小于等于 \(\sqrt n\) 的素数复杂度可以看成 \(O\left(\dfrac {\sqrt n}{\log \sqrt n} \times \log^2 n\right) = O(\sqrt n \log n)\)，而对于大于 \(\sqrt n\) 的素数只会有一项。因此，计算素数幂处取值的复杂度可以看作低阶项，不影响用 \(O()\) 计算的复杂度。
Code：[V 我 50 我帮你写]。
杜教筛
前置：\(\sum \limits_{i = 1}^{n} i^k \ \ (k\leq 4)\) 公式
- \(k = 1\)：\(n(n+1)/2\)。
- \(k = 2\)：\(n(n+1)(2n+1)/6\)。
- \(k = 3\)：\(n^2(n+1)^2/4\)。
- \(k = 4\)：\(n(n+1)(2n+1)(3n^2+3n-1)/30\)。已经很少用到了。
- 对于 \(k > 4\)，使用连续点值拉格朗日插值可以 \(O(k)\) 计算一次。
流程
考虑求某数论函数 \(f\) 的前缀和 \(S(n) = \sum \limits_{i = 1}^n f(i)\)。杜教筛的想法是通过一个比较好的数论函数 \(g\) 和比较好的 \(f*g\) 把 \(f\) 的前缀和化简。这和 \(f,g,f*g\) 是否是积性函数无关。
我们知道
如果 \(f*g\) 和 \(g\) 的前缀和比较简单，我们就可以快速求出 \(f\) 的前缀和了。
考虑如下伪代码：
inline int SumF(int n){ // 注意：这里计算的是 g(1)S(n). 但大部分时候 g 是积性函数或其点乘，因此 g(1)=1.
int ans=SumFG(n);
for(int l=2,r;l<=n;l=r+1){
r=n/(n/l);
ans-=(SumG(r)-SumG(L-1))*SumF(n/l);
}
return ans;
}
当使用记忆化的时候，根据整除分块经典结论，我们只会用到不超过 \(2 \sqrt n\) 个值。
那么 \(T(n) = O(\sqrt n) + \sum \limits_{i = 1}^{\sqrt n}\left(O(\sqrt i) +O\left(\sqrt {\dfrac{n}{i}}\right)\right)\)。求和改积分容易证得 \(T(n) = O(n^{3/4})\)。
否则时间复杂度不详（OI-Wiki 对不加以记忆化的时间复杂度证明作出了证伪：不能直接将第二层之后的 \(T\left(\sqrt{\dfrac{n}{ij}}\right)\) 看作余项）。
如果我们预处理部分 \(S(i) (i = 1,2,\cdots,m) \ \ (m \geq \sqrt n)\)，设预处理复杂度为 \(T_0(m)\)，则总复杂度为
那么一般 \(T_0(m) = O(m)\)，此时 \(m = \Theta(n^{2/3})\) 最优。
详见 OI-Wiki。
两例简单的杜教筛
-
\(f = \mu\)
根据 \(\mu * I = \epsilon\)，取 \(f = \mu,g = I, f*g = \epsilon\)。
-
\(f = \varphi\)
根据 \(\varphi * I = id\)，取 \(f = \varphi,g = I,f * g = id\)。
点乘
定义 \(f,g\) 点乘的结果 \(f\cdot g\) 是一个函数，满足 \((f\cdot g)(n) = f(n)g(n)\)。
对于数论函数 \(A,B,C\)，当 \(C\) 是完全积性函数的时候，我们有
那么当我们要求前缀和的函数 \(f\) 为两个简单函数点乘的时候，我们可以运用这个等式，从而寻求形式较好的 \(g\) 和 \(f*g\)，然后杜教筛解决。
下面是两例比较基本的点乘处理方法。
-
\(f = \mu \cdot id_k\)：取 \(g = id_k\)。
\[h = (\mu \cdot id_k )* id_k = (\mu \cdot id_k) * (I * id_k) = (\mu * I) * id_k = \epsilon \] -
\(f = \varphi \cdot id_k\)：取 \(g = id_k\)。
\[h = (\varphi \cdot id_k )* id_k = (\varphi \cdot id_k) * (I * id_k) = (\varphi * I) * id_k = id_{k+1} \]
下面是抄的一例点乘：
- \(f = \mu^2 * (\mu \cdot id)\)：取 \(g = id\)。\[h = \mu^2 * (\mu \cdot id) * id = \mu^2 * ((\mu * I) \cdot id) = \mu ^2 * \epsilon = \mu ^2 \]而 \(\mu^2\) 的前缀和可以通过公式 \(\mu^2(i) = \sum \limits_{d^2 \mid i} \mu(d)\) 计算：简单地交换求和号，前缀和变为 \(\sum \limits_{d = 1}^{\sqrt n} \mu(d) \left \lfloor \frac{n}{d^2} \right \rfloor\)。直接计算即可做到 \(O(\sqrt n)\)，求出 \(\mu(i)\) 的前缀和可以进一步做到 \(O(n^{1/3})\)。但是在低于 \(O(\sqrt n)\) 时间复杂度内求出 \(\mu(i)\) 的前缀和需要再写筛法，其实就不如直接 \(O(\sqrt n)\) 求了。
另一例略显奔放的点乘（出自某场联考的最后一部分）：
- \(f = id_k \cdot ((\mu \cdot id_k) * I)\)：取 \(g = id_{2k}\)。\[h = (id_k\cdot ((\mu \cdot id_k) * I) ) * (id_k \cdot id_k)= ({\color{red}(\mu \cdot id_k) * id_k}* I) \cdot id_k \\ = (\epsilon * I) \cdot id_k \\ = id_k \]红色部分是第一例点乘中出现的结果，它等于 \(\epsilon\)。
小结：欲求 \(f = A\cdot C\) 形的前缀和，我们希望 \(g\) 形如 \(B \cdot C\)，且 \(g\) 的前缀和易求。然后，运用上面的等式将 \(f * g\) 化作 \((A*B) \cdot C\)。适当地选取 \(B\) 可以将 \(A*B\) 化作简单的函数，从而 \(f*g\) 的前缀和易求。
简单运用
例 1
求
\[\sum \limits_{i = 1}^{n} \sum \limits_{j = 1}^{m} \varphi(\gcd(i,j)) \]对某个大质数取模
\(1 \leq n \leq m \leq 10^9\)
经过平凡的推导可以得出原式等于
想要整除分块，则需要求出 \(\mu * \varphi\) 的前缀和。
令 \(f = \mu * \varphi,g = I,f*g = \varphi\)。利用杜教筛求出 \(\varphi\) 和 \(\mu * \varphi\) 的前缀和即可。注意到复杂度仍然是 \(O(n^{2/3})\) 的，因为我们可以线性筛预处理 \(f\) 的前 \(O(n^{2/3})\) 项。
例 2 (Luogu P3768 简单的数学题)
求
\[\sum \limits_{i = 1}^{n} \sum \limits_{j = 1}^{n} ij\gcd(i,j) \]对某个大质数取模
\(1 \leq n \leq 10^9\)
经过平凡的推导可以得出原式等于
我们知道了怎么处理点乘，因此杜教筛也是容易的：\(f = \varphi \cdot id_2,g = id_2,(f*g) = id_3\)。
例 3
求
\[\sum \limits_{i = 1}^{n} \sum \limits_{j =1}^{n} \operatorname{lcm}(i,j)^k \]\(1 \leq n \leq 10^9,1 \leq k \leq 1000\)
这碟醋终于出现了！
平凡推导可得原式等于
按照上面点乘处的推导，令 \(f(n) = n^k \sum \limits_{d \mid n} \mu(d) \cdot d^k\)，请出 \(g(n) = n^{2k}\)，\(f * g = id_k\)。
一个粗劣的复杂度上界是 \(O((\sqrt n + n^{2/3})k)\)，但事实上整个杜教筛过程中也只会用到 \(O(\sqrt n)\) 个 \(S_k\) 的值，对它们记忆化可以分析出更好的复杂度 \(O(\sqrt n k+n^{2/3})\)。
[NOI2016] 循环之美
简要题面：
\[\sum \limits_{i = 1}^{n} \sum \limits_{j = 1}^m [\gcd(i,j) = 1][\gcd(j,k) = 1] \]\(1 \leq n,m \leq 10^9,1 \leq k \leq 10^3\)
这题处理技巧比较奇怪，我们不太能直接把两个艾弗森括号都拆成 \(\mu\)，而是先保留下来 \([\gcd(j,k) = 1]\) 的那部分。在这里，和一个给定常数互质这个条件是比较好的。
然后要求两个东西：
对于 \(f(n)\)，根据 \(\gcd(i,k) = \gcd(i+k,k)\)，艾弗森括号的取值是 \(k\) - 循环的，我们只需要预处理 \(f(1) \sim f(k)\) 即可。
对于 \(g(n,k)\) 则比较需要有想象力。
递归计算，直到递归到 \(d = 1\) 时用杜教筛求 \(\mu\) 前缀和。
注意到最终要算 \(O(\sqrt n)\) 个 \(g(n,...)\)，\(...\) 这一维是 \(k\) 的因子，因此总复杂度为 \(O(d(k)\sqrt n+ n^{2/3} + k)\)。
BZOJ3512: DZY Loves Math IV
求
\[\left(\sum \limits_{i = 1}^{n} \sum \limits_{j = 1}^{m} \varphi(ij) \right)\bmod 998244353 \]\(1 \leq n \leq 10^5,1 \leq m \leq 10^9,n \leq m\)
根据《毒瘤之神的考验》一题快进到
记 \(S(a,b) = \sum \limits_{i = 1}^{a} \varphi(ib),f(n) = \sum \limits_{d \mid n} \frac{d\mu\left(\frac Td\right)}{\varphi(d)}\)，则上式可写作
在《毒瘤之神的考验》一题中，要用到的全体 \(S(x,y)\) 可 \(O(n \log n)\) 预处理，从而最终复杂度也是 \(O(n \log n)\) 的。
本题需要神秘的观察。我们考虑类似扰动法的策略，
注意到 \(k\) 这一维是小于等于给定的 \(n\) 的（和此处的变量 \(n\) 作区分），于是对于 \(S(m/i,i)\) 这样的状态来说，计算它们的复杂度不太大。我们可以采用记忆化 + 预处理小的部分来减小递归的复杂度。
Min_25 筛
综述
考虑求一类积性函数的前缀和，其中 \(f\) 是积性函数，满足：
- 对于任意素数 \(p\)，\(f(p)\) 处取值是关于 \(p\) 的低次多项式。
- \(f(p^c)\) 处取值容易求得。
则 Min_25 筛可以在 \(O(\frac{n^{3/4}}{\log n})\) 或 \(O(n^{1 - \epsilon})\) 的复杂度内求出 \(f\) 的前缀和。
流程
Min_25 筛的求解分为两步：
- 对于 \(W(n) = \left\{\left\lfloor \dfrac{n}{x}\right\rfloor \mid 1 \leq x \leq n\right\}\)，考虑其中每一个元素 \(m\)，求出 \(\sum \limits_{1 \leq i \leq m,i \in \mathbb P} f(i)\)，即前缀素数点处取值之和；
- 根据第一步求得的 \(O(\sqrt n)\) 个答案进而求解出整个前缀和。
第一步
因为 \(f(p)\) 的取值是关于 \(p\) 的低次多项式，因此，在这一步中，我们只需要考虑 \(f(p) = p^k\) 的情形如何求解（对于项数更多的情形，分别求解后求和即可）。
我们通常找一个 \(f'\) 来替代 \(f\)。这样的 \(f'\) 具有完全积性，且在素数 \(p\) 处和 \(f\) 相同。这里考虑函数 \(f'(x) = x^k\)。
Min_25 筛的第一步，我们要对于每个 \(m\)，求出数组 \(g\)。其中 \(g\) 的定义如下：
其中 \(L(i)\) 表示 \(i\) 的最小质因子，\(P_j\) 表示第 \(j\) 个质数，\(P_0 = 0\)。上式即前缀素数点处或最小质因子大于 \(P_j\) 处 \(f'\) 取值之和（注意是 \(f'\) 而非 \(f\)。我们用 \(f'\) 暂时替代了 \(f\)）。
对于 \(1\)，其不存在最小质因子，故 \([L(1) > P_j]\) 总是 \(0\)。此处可改作 \(2 \leq i \leq m\)，但为了推导方便保留可能的 \(i = 1\) 取值，虽然 \(f'(1)\) 永远不能被统计进贡献。
设 \(q(m)\) 表示最大的 \(j\) 使得 \(P_j^2 \leq m\)，那么 \(g(m,q(m))\) 即为 \(\sum \limits_{1 \leq i \leq m,i \in \mathbb P} f(i)\)。
考虑递推求出 \(g(m,j)\)。初始有 \(g(m,0) = \sum \limits_{i = 2}^{m} i^k\)。对于 \(j \geq 1\)，我们会在 \(g(m,j-1)\) 的基础上删除若干个点。质数是不会被删除的，因此我们会恰好删除最小质因子是 \(P_j\) 的非质数点。
若 \(P_j^2 > m\)，\(g(m,j)\) 不会有变化。
若 \(P_j^2 \leq m\)，形如 \(P_jx(P_jx \leq m)\)，且 \(P_jx\) 的最小质因子恰好是 \(P_j\) 的点会造成贡献。这可以表示为下面的式子：
来解释一下这个式子：首先，我们可以根据完全积性把 \(f'(P_jx)\) 拆开。现在考虑 \(x\) 的取值，很显然是来源于 \(g\left(\left \lfloor \dfrac{m}{P_j} \right \rfloor,j-1\right)\)。但是既小于等于 \(\left \lfloor \dfrac{m}{P_j} \right \rfloor\) 又小于等于 \(P_{j-1}\) 的质数是不能作为 \(x\) 的，因为我们有 \(P_{j-1} < P_j \leq \sqrt m\)，所以 \(\left \lfloor \dfrac{m}{P_j} \right \rfloor \geq P_{j-1}\)，只需要考虑 \(g(P_{j-1},j-1)\) 就可以了。
因此，我们得到了下式：
复杂度分析 · 一
现在考虑复杂度。根据经典结论 \(\left \lfloor \dfrac{n}{ab} \right \rfloor=\left \lfloor \dfrac{\frac{n}{a}}{b} \right \rfloor\)，我们对于 \(m \in W(n)\) 有 \(W(m) \subseteq W(n)\)。同时可以证明 \(\forall 1 \leq i \leq \sqrt n\)，\(i \in W(n)\)，因此对于任意 \(m \in W(m)\) 我们求 \(g(m,j)\) 的时候，只会遍历 \(g\) 数组中第一维在 \(W(n)\) 中的元素。同时，求出单个状态的复杂度是 \(O(1)\) 的，因此状态数即为复杂度。
根据范围素数近似估计，复杂度可写作
注意到 \(i \leq \sqrt n\)，因此 \(\ln \sqrt \frac ni\) 只和 \(\ln n\) 差常数倍。直接将这部分提出来，上面的部分积分，得到
第二步
综述
现在我们已经知道所有 \(\sum \limits_{1 \leq i \leq m,i \in \mathbb P} f(i),\ \text{where}\ m \in W(n)\) 了。将其简记为 \(g(m)\)。考虑使用这些 \(g(m)\) 求出原来的解。
设 \(S(m,j)\) 表示 \(\sum \limits_{1 \leq i \leq m}[L(i) > P_j] f(i)\)，即最小质因子大于 \(P_j\) 的 \(f(i)\) 之和（我们现在抛弃 \(f'\) 了，同时也不再限定 \(i\) 是质数；只不过 \(1\) 仍然被抛弃在外）。这样再调用 \(S(n,0)\) 就得到答案了。
方法一
这种方法复杂度为 \(O(n^{1 - \epsilon})\)，在 \(n \leq 10^{13}\) 的时候也可以看作 \(O\left(\frac{n^{3/4}}{\log n}\right)\)，够用。实际测试略快于第二种。这里着重介绍。
仍然考虑如何求解。分 \(i\) 是质数或合数两种情况讨论：
-
\(i\) 是质数
根据 \(g\) 的定义，这部分的贡献是 \(g(n) - g(P_j)\)。
-
\(i\) 是合数
枚举最小质因子 \(P_k (k > j)\)，枚举指数 \(e\)，利用积性函数的性质进行变量分离，得到
\[\sum \limits_{k > j} \sum \limits_{e \geq 1,P_{k}^e \leq m} f(P_{k}^e)\left(S\left(\left\lfloor\frac{m}{P_k^e}\right\rfloor,k\right)+ [e>1]\right) \]这里 \([e>1]\) 有些耐人寻味。考虑前面的 \(S\left(\left\lfloor\frac{m}{P_k^e}\right\rfloor,k\right)\) 没有考虑 \(P_{k}^e\) 自己的贡献。如果 \(e = 1\)，那么 \(P_k^e\) 就是 \(P_k\)，应该归到 \(i\) 是质数的分类里，没有贡献。否则，\(P_{k}^e\) 会造成 \(f(P_{k}^e)\) 的贡献。
则 \(S(n,0) + 1\) 就是答案（我们之前没有考虑 \(f(1)\) 的贡献）。
代码 (For Luogu P4213)
略好于之前的版本。
# include <bits/stdc++.h>
const int N=1000010,INF=0x3f3f3f3f,mod=1e9+7,inv6=166666668,inv2=500000004;
typedef long long ll;
int pr[N],pc;
bool vis[N];
int g1[N],g2[N];
ll w[N],idx1[N],idx2[N];
int wtot;
ll MAXN;
ll n;
inline int read(void){
int res,f=1;
char c;
while((c=getchar())<'0'||c>'9')
if(c=='-')f=-1;
res=c-48;
while((c=getchar())>='0'&&c<='9')
res=res*10+c-48;
return res*f;
}
inline void init(void){
for(int i=2;i<=MAXN;++i){
if(!vis[i]) pr[++pc]=i;
for(int j=1;i*pr[j]<=MAXN&&j<=pc;++j){
vis[i*pr[j]]=true;
if(i%pr[j]==0) break;
}
}
return;
}
inline int sum1(ll x){
x%=mod;
return x*(x+1)%mod*inv2%mod;
}
inline int sum2(ll x){
x%=mod;
return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
}
inline int idx(ll x){
return (x<=MAXN)?idx1[x]:idx2[n/x];
}
inline int adc(int a,int b){
return (a+b>=mod)?(a+b-mod):(a+b);
}
inline int dec(int a,int b){
return (a<b)?(a-b+mod):(a-b);
}
inline void add(int &a,int b){
a=adc(a,b);
}
inline void del(int &a,int b){
a=dec(a,b);
}
inline int mul(int a,int b){
return 1ll*a*b%mod;
}
ll S(ll x,int i){
if(pr[i]>=x) return 0;
int xid=idx(x),pid=idx(pr[i]);
ll res=dec(dec(g2[xid],g1[xid]),dec(g2[pid],g1[pid]));
for(int j=i+1;j<=pc&&1ll*pr[j]*pr[j]<=x;++j){
for(ll e=1,p=pr[j],_p;p<=x;++e,p*=pr[j]){
_p=p%mod;
res=adc(res,mul(mul(dec(_p,1),_p),S(x/p,j)+(e>1)));
}
}
return res;
}
int main(void){
scanf("%lld",&n);
MAXN=1ll*sqrt(n);
init();
for(ll l=1,r;l<=n;l=r+1){
ll val=n/l;
r=n/val,w[++wtot]=val;
g1[wtot]=dec(sum1(val),1),g2[wtot]=dec(sum2(val),1);
if(val<=MAXN) idx1[val]=wtot;
else idx2[n/val]=wtot;
}
for(int i=1;i<=pc;++i){
for(int j=1;1ll*pr[i]*pr[i]<=w[j]&&j<=wtot;++j){
int mid=idx(w[j]/pr[i]),pid=idx(pr[i-1]);
del(g1[j],mul(pr[i],dec(g1[mid],g1[pid])));
del(g2[j],mul(mul(pr[i],pr[i]),dec(g2[mid],g2[pid])));
}
}
printf("%d",adc(S(n,0),1));
return 0;
}
方法二
再度审视 \(S\) 的定义：\(S(m,j)\) 表示 \(\sum \limits_{1 \leq i \leq m}[L(i) > P_j] f(i)\)，即最小质因子大于 \(P_j\) 的 \(f(i)\) 之和。
我们直接套用第一步中的递推方式。
这里的取 \(\min\) 是不好的。我们注意到当 \(P_{j+1}^{e+1} > m\) 的时候，\(S() - g()\) 没有贡献，因为这意味着 \(P_{j+1} > \dfrac{m}{P_{j+1}^e}\)。于是可以改为
使用和第一步相同的方法进行递推即可。
至此可以使用与第一步类似的方法分析出 \(O\left(\frac{n^{3/4}}{\log n}\right)\)。
方法二的优势是，我们可以求出每个 \(m \in W(n)\) 的前缀 \(m\) 的答案。这无疑在整除分块需要计算 \(O(\sqrt n)\) 个前缀的函数值之和的时候给我们提供了便利。
代码 (For LibreOJ 6784)
# include <bits/stdc++.h>
const int N=2000010,INF=0x3f3f3f3f,mod=1e9+7,inv6=166666668,inv2=500000004;
typedef long long ll;
int pr[N],pc;
bool vis[N];
int g1[N],g0[N];
ll w[N],idx1[N],idx2[N];
int s[N];
int wtot;
ll MAXN;
ll n;
inline int read(void){
int res,f=1;
char c;
while((c=getchar())<'0'||c>'9')
if(c=='-')f=-1;
res=c-48;
while((c=getchar())>='0'&&c<='9')
res=res*10+c-48;
return res*f;
}
inline void init(void){
for(int i=2;i<=MAXN;++i){
if(!vis[i]) pr[++pc]=i;
for(int j=1;i*pr[j]<=MAXN&&j<=pc;++j){
vis[i*pr[j]]=true;
if(i%pr[j]==0) break;
}
}
return;
}
inline int sum1(ll x){
x%=mod;
return x*(x+1)%mod*inv2%mod;
}
inline int sum0(ll x){
return x%mod;
}
inline int idx(ll x){
return (x<=MAXN)?idx1[x]:idx2[n/x];
}
inline int adc(int a,int b){
return (a+b>=mod)?(a+b-mod):(a+b);
}
inline int dec(int a,int b){
return (a<b)?(a-b+mod):(a-b);
}
inline void add(int &a,int b){
a=adc(a,b);
}
inline void del(int &a,int b){
a=dec(a,b);
}
inline int mul(int a,int b){
return 1ll*a*b%mod;
}
inline int f(int x,int k){
return x^k;
}
std::vector <int> vec;
int main(void){
scanf("%lld",&n);
MAXN=1ll*sqrt(n);
init();
for(ll l=1,r;l<=n;l=r+1){
ll val=n/l;
r=n/val,w[++wtot]=val;
g1[wtot]=dec(sum1(val),1),g0[wtot]=dec(sum0(val),1);
if(val<=MAXN) idx1[val]=wtot;
else idx2[n/val]=wtot;
}
for(int i=1;i<=pc;++i){
for(int j=1;1ll*pr[i]*pr[i]<=w[j]&&j<=wtot;++j){
int mid=idx(w[j]/pr[i]),pid=idx(pr[i-1]);
del(g1[j],mul(pr[i],dec(g1[mid],g1[pid])));
del(g0[j],dec(g0[mid],g0[pid]));
}
}
for(int i=1;i<=wtot;++i) s[i]=g1[i]=adc(dec(g1[i],g0[i]),(w[i]>=2)*2);
for(int i=pc;i;--i){
for(int j=1;1ll*pr[i]*pr[i]<=w[j]&&j<=wtot;++j){
for(ll e=1,p=pr[i];p<=w[j]/pr[i];++e,p*=pr[i]){
add(s[j],mul(f(pr[i],e),dec(s[idx(w[j]/p)],g1[idx(pr[i])])));
add(s[j],f(pr[i],e+1));
}
}
}
vec.resize(wtot);
for(int i=0;i<wtot;++i) vec[i]=adc(s[i+1],1);
std::sort(vec.begin(),vec.end());
vec.erase(std::unique(vec.begin(),vec.end()),vec.end());
int ans=0;
for(auto v:vec) ans^=v;
printf("%d",ans);
return 0;
}
小试牛刀
例 1
(杜教筛 例 3)
求
\[\sum \limits_{i = 1}^{n} \sum \limits_{j =1}^{n} \operatorname{lcm}(i,j)^k \]\(1 \leq n \leq 10^9,1 \leq k \leq 1000\)
快进到这里。
那么我们只要瞪眼法看出后面是积性函数，直接上 Min_25 筛就完事了，不需要观察到可以卷 \(id_{2k}\)。
Typora 崩掉了，不写了。
浙公网安备 33010602011771号