CVBB ACM Team

用户工具

站点工具

您在这里: front_page » 2020-2021 » teams » i_dont_know_png » potassium » sieve
2020-2021:teams:i_dont_know_png:potassium:sieve

差别

这里会显示出您选择的修订版和当前版本之间的差别。

到此差别页面的链接

两侧同时换到之前的修订记录 前一修订版
后一修订版 		前一修订版
2020-2021:teams:i_dont_know_png:potassium:sieve [2020/05/26 11:32]
potassium [实例] 		2020-2021:teams:i_dont_know_png:potassium:sieve [2020/06/07 12:46] (当前版本)
nikkukun fix some typos
行 618: 行 618:
  
 ===== min_25 筛 ===== ===== min_25 筛 =====
 +
 +
 +
 +min_25 筛可以用来解决具有下列性质的积性函数 $f$ 在单点处的前缀和：
 +
 +  - 存在完全积性函数 $f'​$，使得质数 $p$ 处有 $f'(p) = f(p)$；
 +  - $f'​(p)$ 的前缀和很好计算；
 +  - $f'​(p^k)$ 的值很好算（是多项式或可以很快求出）。
 +
 +为了方便，我们约定：
 +
 +  * 若没有特别说明，下文中的所有 $p$ 都取质数；
 +  * $\mathrm{pri}_j$ 表示第 $j$ 小的质数；
 +  * $\mathrm{minp}(n)$ 表示 $n$ 的最小质因子。
 +
 +==== 计算质数贡献 ====
 +
 +我们先考虑解决求所有分块点处质数的贡献，即 $G(m) = \sum _{i=1}^m [i \text{ is prime}] f'​(i)$，其中 $m$ 是 $n$ 的数论分块点。
 +
 +设状态 $g(j, m)$ 表示 $\sum _{i=1}^m [i \text{ is prime or } \mathrm{minp}(i) > \mathrm{pri}_j] f'​(i)$，则我们要求的就是 $G(m) = g(\sum _{j \land \mathrm{pri}_j \leq m} 1, m)$ 了。这个状态可以理解为，Eratosthenes 筛法中用 $\mathrm{pri}_j$ 筛掉它的倍数后，剩余没被筛掉部分的和。
 +
 +显然，边界条件是 $g(0, m) = \sum _{i=2}^m f'​(i)$。请注意一定要去掉 $f'​(1)$，因为它不会被筛掉，我们只要最后加上即可。
 +
 +考虑转移。若 $\mathrm{pri}_j > \sqrt m$，此时筛不掉任何数，故 $g(j, m) = g(j-1, m)$，故考虑 $\mathrm{pri}_j \leq \sqrt m$ 的情况。此时 $\mathrm{pri}_j$ 能筛掉的**合数**都能表示为 $\mathrm{pri}_j$ 乘一个不超过 $\left\lfloor \dfrac m{\mathrm{pri}_j} \right\rfloor$ 的、最小质因子**大于等于** $\mathrm{pri}_j$ 的数，即
 +
 +$$
 +g(j, m) = g(j - 1, m) - f'​(\mathrm{pri}_j) \times [g(j - 1, \left\lfloor \dfrac m{\mathrm{pri}_j} \right\rfloor) - \sum _{i=1}^{j-1} f'​(\mathrm{pri}_i)]
 +$$
 +
 +因为 $f'$ 是完全积性函数，所以即使分解的东西不互质也可以直接乘在一起。另外，式子最后的部分用于减去质数的贡献（参考 $g$ 的定义）。
 +
 +具体实现时，可以用滚动数组存放 $g$，滚掉它的第一维 $j$。可以证明这一部分的时间复杂度是 $O\left(\dfrac {n^{3/​4}}{\log n} \right)$ 的。
 +
 +
 +==== 计算合数贡献 ====
 +
 +设 $S(n, j) = \sum _{i=1}^m [\mathrm{minp}(i) \geq \mathrm{pri}_j] f(i)$，注意中间对 $\mathrm{minp}(i)$ 的要求和 $G$ 有所不同。根据定义，我们最后要求的就是 $S(n, 1) + f(1)$。
 +
 +对于 $S(n, j)$，质数的贡献为 $G(n) - \sum _{i=1}^{j-1} f'​(\mathrm{pri}_i)$，而合数的贡献需要枚举其最小质因子及其次数，分为两个互质的部分。转移如下：
 +
 +$$
 +S(n, j) = G(n) - \sum _{i=1}^{j-1} f'​(\mathrm{pri}_i) + \sum _{i = j}^{\mathrm{pri}_i \leq n} \sum _{e=1} ^{\mathrm{pri}_i ^{e+1} \leq n} [f(\mathrm{pri}_i ^e) S(\left\lfloor \dfrac n{\mathrm{pri}_i^e} \right\rfloor,​ i + 1) + f(\mathrm{pri}_i^{e+1})]
 +$$
 +
 +这部分不需要记忆化，直接递归的时间复杂度是 $O(n^{1 - \varepsilon})$ 的，而且跑 $n = 10^{10}$ 的数据也很飞快。
  
  
行 623: 行 668:
 ==== 实例 ==== ==== 实例 ====
  
 +
 +=== Luogu 5325 ===
 +
 +**题意**：积性函数 $f(n)$ 的质数幂次取值 $f(p^k) = p^k(p^k - 1)$，求前缀和。
 +
 +**题解**：$f(p^k) = (p^2)^k - p^k$，拆成两部分计算质数处贡献，最后合并即可。
 +
 +
 +
 +<hidden 参考实现>​
 +<​code:​cpp>​
 +#include <​bits/​stdc++.h>​
 +using namespace std;
 +
 +#define fi first
 +#define se second
 +#define lch (o << 1)
 +#define rch (o << 1 | 1)
 +
 +typedef double db;
 +typedef long long ll;
 +typedef unsigned int ui;
 +typedef pair<​int,​ int> pint;
 +
 +const int N = 2e5 + 5;
 +const int MOD = 1e9 + 7;
 +const int INV2 = 500000004;
 +const int INV6 = 166666668;
 +const int INF = 0x3f3f3f3f;
 +const ll INF_LL = 0x3f3f3f3f3f3f3f3f;​
 +
 +// 注意两倍空间
 +int idx[N];
 +ll val[N];
 +ll _n, sqrtN;
 +int nIdx;
 +
 +int GetId(ll x) {
 + if (x <= sqrtN)
 + return idx[x];
 + else
 + return idx[sqrtN + _n / x];
 +}
 +
 +void InitId(ll n) {
 + // x =  n / i 是该分块点中能取到最大 / 最右侧的值
 + for (ll i = 1; i <= n; i++) {
 + ll x = n / i;
 + val[++nIdx] = x;
 + if (x <= sqrtN)
 + idx[x] = nIdx;
 + else
 + idx[sqrtN + n / x] = nIdx;
 + i = n / (n / i);
 + }
 +}
 +
 +int pri[N];
 +ll s1[N], s2[N];
 +bool notPri[N];
 +
 +void InitPri() {
 + for (int i = 2; i < N; i++) {
 + if (!notPri[i])
 + pri[++pri[0]] = i;
 + for (int j = 1; j <= pri[0]; j++) {
 + int p = pri[j];
 + if (i * p >= N) break;
 + notPri[i * p] = 1;
 + if (i % p == 0) break;
 + }
 + }
 + for (int i = 1; i <= pri[0]; i++) {
 + s1[i] = (s1[i - 1] + 1LL * pri[i] * pri[i]) % MOD;
 + s2[i] = (s2[i - 1] + pri[i]) % MOD;
 + }
 +}
 +
 +ll f1[N], f2[N];
 +
 +// Sieve1(n) = sum [i = 1 -> n] [i is isprime] * f(i)
 +// 传入的 n 是原始值
 +void Sieve1(ll n) {
 + // 预处理初值
 + // 需要注意这里一般不计算 f(1)
 + for (int i = 1; i <= nIdx; i++) {
 + ll x = val[i] % MOD;
 + f1[i] = x * (x + 1) % MOD * (2 * x + 1) % MOD * INV6 % MOD - 1;
 + f2[i] = (1 + x) * x % MOD * INV2 % MOD - 1;
 + }
 + for (int j = 1; j <= pri[0]; j++) {
 + ll p = pri[j];
 + if (p * p > n) break;
 + for (int i = 1; i <= nIdx; i++) {
 + ll x = val[i];
 + if (p * p > x) break;
 + // f(m) = f(idx[m / pri[j]]) - sum[j - 1];
 + f1[i] -= (p * p) % MOD * (f1[GetId(x / p)] - s1[j - 1]);
 + f2[i] -= p * (f2[GetId(x / p)] - s2[j - 1]);
 + f1[i] %= MOD;
 + f2[i] %= MOD;
 + }
 + }
 +}
 +
 +// f(p ^ e)
 +ll F(ll p, int e) {
 + ll t = 1;
 + while (e--) t *= p;
 + t %= MOD;
 + return t * (t - 1) % MOD;
 +}
 +
 +// Sieve2(n, j) = sum [i = 1 -> n] [minp(i) >= pri[j]] * f(i)
 +// 传入的 n 是原始值
 +ll Sieve2(ll n, int j) {
 + if (n <= 1 || pri[j] > n) return 0;
 + int id = GetId(n);
 + // ret = Sieve1(n) - sum [i = 1 -> j - 1] f(i) + ...
 + ll ret = (f1[id] - s1[j - 1]) - (f2[id] - s2[j - 1]);
 +
 + // ret = ... + 合数部分
 + for (int i = j; i <= pri[0]; i++) {
 + ll p = pri[i], t = p; // t = p ^ e
 + if (p * p > n) break;
 + for (int e = 1; p * t <= n; t *= p, e++) {
 + ll tmp = F(p, e) * Sieve2(n / t, i + 1) + F(p, e + 1);
 + ret = (ret + tmp) % MOD;
 + }
 + }
 +
 + return ret;
 +}
 +
 +int main() {
 + ios::​sync_with_stdio(0);​
 +
 + ll n;
 + cin >> n;
 + sqrtN = floor(sqrt(n));​
 + _n = n;
 +
 + InitId(n);
 + InitPri();
 + Sieve1(n);
 + ll ans = (Sieve2(n, 1) + 1) % MOD;
 + cout << (ans + MOD) % MOD;
 +
 + return 0;
 +}
 +</​code>​
 +</​hidden>​
 +\\
  
  
2020-2021/teams/i_dont_know_png/potassium/sieve.1590463940.txt.gz · 最后更改: 2020/05/26 11:32 由 potassium

页面工具

除额外注明的地方外，本维基上的内容按下列许可协议发布： CC Attribution-Share Alike 4.0 International
CC Attribution-Share Alike 4.0 International Donate Powered by PHP Valid HTML5 Valid CSS Driven by DokuWiki
京ICP备20018844号-3