Loading... ## 摘要 有人也叫做扩展埃氏筛?其实确实挺像的,对于一类积性函数,函数在质数上的取值容易求得,整个函数在任意一点的值、求和、乘积都有办法通过该类方法求得。对于函数有如下要求: - 积性函数 - 设 $p$ 为质数,那么 $f(p)$ 是一个低阶多项式 - $f(p^k)$ 的值可以快速求出 ## 定义 以 Luogu 模板题为例 *def.* 积性函数 $f(x)$ 且满足 $f(p^k) = p^{k}(p^{k} - 1)$,求 $\sum\limits_{i=1}^{n} f(i)$ ## 步骤一:按照质数分类 整个函数的求和答案可以分为函数在质数上的值的和与合数上的值的和。也就是:(为了方便下面均使用 $P$ 表示质数集,$p$ 表示枚举的质数) $$ \sum\limits_{i=1}^{n} = \sum\limits_{p \in [1,n]} f(p) + \sum\limits_{i=1 \& i \notin P}^{n} f(i) $$ 这一步还是非常简单 ## 步骤二:将合数分解质因数以计算函数值 枚举后面合数的最小质因子和次数,注意所有质因子都因 $\leq \sqrt{n}$ 。注意下式中使用 $minp$ 代表最小质因子。 $$ \sum\limits_{p \in [1,n]} f(p) + \sum\limits_{p^e \in [1,n] \& p \in [1, \sqrt{n}]} (\sum_\limits{i = 1 \& minp > p}^{\frac{p^e}{n}} f(i)) $$ 看一下这个看起来非常牛逼的式子。其实也很好理解。前面的项没有变,仍然是函数在质数上的值的和。后面的项首先枚举 $[1,n]$ 中的每一个质因子 $p$ 和它的指数 $e$,并且 $p \leq \sqrt{n}$,后面再单独乘上函数在 $minp > \sqrt{n}$ 处的取值。 ## 步骤三:通过第一个 dp 计算后面项的值 我们知道函数在质数上的取值是一个低阶的多项式,我们通过一个 dp 来计算。整个 Min_25 筛的过程中运用了两次 dp,这是第一次。这个 dp 确实挺智慧的。 设 $g(n,i)$ 满足:($\lor$ 符号表示逻辑或运算) $$ g(n,i) = \sum\limits_{i=1}^{n} [i\in p \lor minp > p_j] i^k $$ 因为是一个低阶多项式,所以 $k$ 很小,不会超过二次幂(应该没有恶心题有三次幂吧...) 这个 dp 的作用就是求 $[1,n]$ 中所有满足条件:质数或 $minp > p_j$ 的数的 $k$ 次幂和。 考虑转移。 - 若 $n < p_{j}^{2}$ ,那么 $g(n,i) = g(n,i - 1)$ - 若 $p_j^{2} \leq n$, 则被筛掉的质数一定有 $p_j$ ,那么要 $-g(\frac{n}{p_j}, j - 1)$ - 注意上面减去的项中,因为 $p_j^2 \leq n \Longleftrightarrow p_j \leq n / p_j$, 所以有 $minp < p_j$ 的 $n’$ 被减去,这部分要加回来:$g(p_{j-1}, j-1)$ 总的转移式: $$ g(n,j) = g(n,j-1) - p_j^k(g(\frac{n}{p_j}, j-1) - g(p_{j-1},j-1)) $$ 注意若 $n \in P$ ,就无需减去后面那一坨了,直接 $g(n,j) = g(n,j-1)$ 性质来了,$g(p_{j-1},j-1)$ 这一项就是前 $j-1$ 个质数的共 $k$ 次幂和。我们设 $sp_n = \sum\limits_{i=1}^{n} p^k$ *claim:* $\lfloor \frac{\lfloor \frac{n}{a} \rfloor} {b} \rfloor = \lfloor \frac{n}{ab} \rfloor$ 利用整除分块可以计算这一部分的答案。 ## 步骤四:求解最终答案 现在仅仅需要将上面的计算整合,算出最终答案。 再次使用 dp,这个 dp 我们叫第二个 dp。设 $s(n,x)$ 表示求 $[1,n]$ 中所有 $minp > p_x$ 的 $f(x)$ 之和。 $$ s(n,x) = g(n,p_x) - sp_x + \sum\limits_{p^e \leq n \& k>x} f(p_k^e) (s(\frac{n}{p^e}, k + [e \neq 1])) $$ 上式分为两部分,前一部分是计算 $ > p_x$ 的数,后一部分是 $minp > p_x$ 的合数,枚举质因子计算即可。 我们可以递归求解 $s$,最后答案在 $s(n,0)$ ## 实现 我估计自己写还是写不出,还是要看代码。先讲一些要注意的地方。以模板题为例。 首先你要预处理处质数,还有他们的一次幂前缀和二次幂前缀和。由于 $n$ 很大,所以计算 $n$ 的时候对于每一个 $\lfloor \frac{n}{x} \rfloor$ 的下表需要用一个 `map` 来存,这样多一个 $\log$ ,但是可以开两个数组分别存 $x$ 和 $n / x$ 这两个数的下标,这样不会超过 $\sqrt{n}$ ,可以接受。 难点可能主要在求 $g$ 上难实现,可能求和的时候还需要用到一些差分技巧。求 $s$ 确实是递归处理,里面有一个整除分块。 文末放一下代码吧。 ## 时间复杂度 并不会证,但是可以知道 $n$ 在 $10^{10}$ 级别的时候大概是 1s, 常数很小。可能接受度高一点的时间复杂度是 $\Theta (\frac{n^{\frac{3}{4}}}{\log n})$ 有一个曲线吧,出处就是 Luogu 模板题的第一篇题解。 ![](http://47.98.232.80/usr/uploads/2021/08/950998997.png) ## Code 我知道讲了这么多你都不一定会,但是有了这个说不定就会了。 实现的并不是很认真(你看我都 `define int long long` 了,所以一共跑了 1.70s) ```cpp #include <bits/stdc++.h> #define int long long using namespace std; const int SIZE = 1e6 + 5; const int mod = 1e9 + 7; int n, cnt, tot, sqrn; int p[SIZE], isPrme[SIZE], sum1[SIZE], sum2[SIZE], g1[SIZE], g2[SIZE], w[SIZE]; int vis1[SIZE], vis2[SIZE], vis[SIZE]; namespace GTR { const int bufl = 1 << 15; char buf[bufl], *s = buf, *t = buf; inline int fetch() { if (s == t) { t = (s = buf) + fread(buf, 1, bufl, stdin); if (s == t) return EOF; } return *s++; } inline int read() { int a = 0, b = 1, c = fetch(); while (c < 48 || c > 57) b ^= c == '-', c = fetch(); while (c >= 48 && c <= 57) a = (a << 1) + (a << 3) + c - 48, c = fetch(); return b ? a : -a; } } using GTR::read; void init(int N) { isPrme[1] = 1; for (int i = 2; i <= N; ++ i) { if (!isPrme[i]) p[++ cnt] = i, sum1[cnt] = (sum1[cnt - 1] + i) % mod, sum2[cnt] = (sum2[cnt - 1] + i * i % mod) % mod; for (int j = 1; j <= cnt && i * p[j] <= N; ++ j) { isPrme[i * p[j]] = 1; if (i % p[j] == 0) break; } } } int qPow(int a, int b) { int ans = 1ll; for ( ; b; b >>= 1, a = a * a % mod) { if (b & 1) ans = ans * a % mod; } return ans % mod; } const int inv2 = qPow(2ll, mod - 2ll), inv3 = qPow(3ll, mod - 2ll); int calc(int x, int y) { if (p[y] >= x) return 0; int k = x <= sqrn ? vis1[x] : vis2[n / x], ans = (g2[k] - g1[k] + mod - (sum2[y] - sum1[y] + mod) % mod) % mod; for (int i = y + 1; i <= cnt && p[i] * p[i] <= x; ++ i) { int pw = p[i]; for (int c = 1; pw <= x; ++ c, pw = pw * p[i]) { int t = pw % mod; ans = (ans + t * (t - 1) % mod * (calc(x / pw, i) + (c != 1)) % mod) % mod; } } return ans % mod; } signed main() { n = read(), sqrn = sqrt(n); init(sqrn); for (int i = 1; i <= n; ) { int j = n / (n / i); w[++ tot] = n / i; g1[tot] = w[tot] % mod, g2[tot] = g1[tot] * (g1[tot] + 1) / 2ll % mod * (2 * g1[tot] + 1) % mod * inv3 % mod; -- g2[tot], g1[tot] = g1[tot] * (g1[tot] + 1) / 2ll % mod - 1ll; if (n / i <= sqrn) vis1[n / i] = tot; else vis2[n / (n / i)] = tot; i = j + 1; } for (int i = 1; i <= cnt; ++ i) { for (int j = 1; j <= tot && p[i] * p[i] <= w[j]; ++ j) { int k = w[j] / p[i] <= sqrn ? vis1[w[j] / p[i]] : vis2[n / (w[j] / p[i])]; g1[j] -= p[i] * (g1[k] - sum1[i - 1] + mod) % mod; g2[j] -= p[i] * p[i] % mod * (g2[k] - sum2[i - 1] + mod) % mod; g1[j] %= mod, g2[j] %= mod; if (g1[j] < 0) g1[j] += mod; if (g2[j] < 0) g2[j] += mod; } } printf("%lld\n", (calc(n, 0) + 1ll) % mod); return 0; } ``` 最后修改:2021 年 08 月 29 日 © 允许规范转载 打赏 赞赏作者 支付宝微信 赞 如果觉得我的文章对你有用,请随意赞赏