Min_26 筛学习笔记

· · 题解

前言

好像没啥要说的，就是没人写想写一写

Min_26 筛

当积性函数 f 在 p 处的点值为关于 p 的低次多项式，且 f(p^k) 的值容易求出时，我们可以使用 Min26 筛在 O\left(n^{\frac23}\right) 的时间复杂度内求出 \sum f 在所有 \left\lfloor\dfrac ni\right\rfloor 处的值。相较于 Min25 筛的 O\left(\dfrac{n^{0.75}}{\log n}\right) 快了不少。

先来约定一些记号：

该筛法共分为四步，总复杂度为 O\left(n^{\frac23}\right)。

Step 1

求出 F_{prime} 在 n 处的所有点值。

来看看 Min25 筛是怎么处理的。由于 f(p) 是关于 p 的一个低次多项式，我们可以把它拆开来，拆成 k 次方之和的形式。设 g(n,j)=\sum^n_i[i\in\mathbb P\lor\text{lpf}(i)\ge p_j]i^k。有转移公式：

g(n,j)=g(n,j-1)-p_j^k\left(g\left(\dfrac n{p_j},j-1\right)-g(p_{j-1},j-1)\right)

减去的部分是「最小质因子恰好为 p_j 的合数」。则 F_{prime}(x)=g(x,\pi(\sqrt n))。

如果直接计算 g 函数，时间复杂度为 O\left(\dfrac{n^{0.75}}{\log n}\right)。如何做到更优呢？

考虑 根 号 分 治。是的你没听错，根号分治。

首先对于 x\ge n^{\frac23} 的点值，我们可以直接用原版 Min25 筛的方式更新。剩下的部分，我们先将每个位置都初始化为 \sum^x_{i=2}i^k，再在之后筛去质数 p_j 时，将 F_p(x) 更新为 F_p(x)-\sum^x_{i=1}[\text{lpf(i)}=p_j\land i\not\in\mathbb P]i^k 即可。

看起来好像和原来完全没区别啊，还把好不容易搞出来的 dp 拆回去了。但实际上我们可以枚举 p_j，然后暴力搜索质因数分解，找出所有 \text{lpf(i)}=p_j 的 i。发现转移式实际上就是个区间修改单点查询，所以用 树 状 数 组 来优化。

但是复杂度依然爆炸。没事，对于 \text{lpf}(x)\le\sqrt[6]n 的点值，我们直接暴力计算。由于满足 x<n^{\frac23} 且 \text{lpf}(x)>\sqrt[6] n 的数有 O\left(\dfrac{n^{\frac 23}}{\log n}\right) 个，算上树状数组的复杂度，总时间复杂度为 O\left(n^{\frac 23}\right)。我们成功地优化了老版筛法的复杂度！

第一步可能把你震撼了三次，但后面就好很多了。

Step 2

求出 F_{\pi(\sqrt[3]n)+1} 在 n 处的所有点值。不妨设 t=\pi(\sqrt[3]n)+1。

因为能造成贡献的 f(x) 均满足 \text{lpf}(x)\ge p_t，所以 x 只可能是 1，p，或 p^2（p\ge p_t）。那么我们可以推出转移式：

F_{t}(x)=1+(F_{prime}(x)-F_{prime}(p_t-1))+\sum^{\sqrt x}_{p\ge p_t}f(p^2)+f(p)(F_{prime}\left(\frac xp\right)-F_{prime}(p))

容易发现，当 x<p_t^2 时，该式的转移都是 O(1) 的，这部分有 O(\sqrt n) 个数。当 x\ge p_t^2 时，后半部分求和的次数为 \pi(t,\sqrt x)。总复杂度为 O(\sqrt n+\sum^{\sqrt[3]n}_{i=1}\pi(\sqrt{\frac ni}))=O\left(\dfrac{n^{\frac 23}}{\log n}\right)。

Step 3

依次求出 F_{\pi(\sqrt[3]n)},F_{\pi(\sqrt[3]n)-1},F_{\pi(\sqrt[3]n)-2},\cdots,F_{\pi(\sqrt[6]n)+1} 在 n 处的所有点值。

可以发现前后区别就是 \text{lpf}(x)=p_k 的这些数。得出有递推：

F_k(x)=\sum_{p_k^i\le x}f(p^i_k)F_{k+1}\left(\dfrac m{p_k}\right)

直接暴力转移复杂度爆炸，仍然考虑根号分治：对于 x\ge n^{\frac23} 处的点值暴力更新，由于有 F_k(x)=F_{k+1}(x)+\sum^x_{i=1}[\text{lpf}(i)=p_k]f(i)，剩下的部分按第一步的方法更新即可，不再赘述。

Step 4

依次求出 F_{\pi(\sqrt[6]n)},F_{\pi(\sqrt[6]n)-1},F_{\pi(\sqrt[6]n)-2},\cdots,F_{1} 在 n 处的所有点值。

这里直接用上一步的式子暴力做就可以了。时间复杂度 O\left(\dfrac{n^{\frac 23}}{\log n}\right)。

参考代码（P5325 【模板】Min_25 筛）

有 f(p^k)=(p^k)^2-p^k，直接计算。

没有特意卡常的情况下跑出了总和 1.57s 的好成绩。

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int MAXN = 1e6 + 10;
const int mod = 1e9 + 7;
const int inv = (mod + 1) / 6;

inline 
ll qpow(ll b, ll p) {
    ll res = 1;
    for (; p; p >>= 1, b = b * b % mod) if (p & 1) res = res * b % mod;
    return res;
}

int p[MAXN], tot, t3 = 1, t6 = 1; bool vis[MAXN];

ll s[MAXN][3], n2, n3, n6, n23; // sqrt(n), cbrt(n), sqrt[6](n), n^2/3

inline 
void init(int n) {
    vis[1] = 1, n += 10;
    for (int i = 1; i <= n; i++) {
        if (!vis[i]) p[++tot] = i;
        for (int j = 1; j <= tot; j++) {
            if (i * p[j] > n) break;
            vis[i * p[j]] = 1;
            if (i % p[j] == 0) break;
        }
        if (i == n3) t3 += tot;
        if (i == n6) t6 += tot;
    }
    for (int i = 1; i <= tot; i++) {
        s[i][1] = (s[i][1] + p[i]) % mod;
        s[i][2] = (s[i][2] + (ll)p[i] * p[i] % mod) % mod;
    }
}

inline 
ll calc(ll n, int k) {
    n %= mod;
    if (!k) return n;
    if (k == 1) return (n * (n + 1) >> 1) % mod;
    if (k == 2) return n * (n + 1) % mod * (2 * n + 1) % mod * inv % mod;
}

inline 
ll f(ll n) {
    n %= mod;
    return (n * n - n) % mod;
}

ll F[MAXN << 1], g[MAXN << 1][3], val[MAXN << 1]; // F(n / i)

ll n; int m;

inline 
ll get(ll k) {
    return k <= n2 ? k : m - n / k + 1;
}

ll c[MAXN];

inline 
void add(int k, ll x) {
    for (int i = k; i <= m; i += i & -i) c[i] = (c[i] + x) % mod;
}

inline 
ll ask(int k) {
    ll res = 0;
    for (int i = k; i; i &= i - 1) res = (res + c[i]) % mod;
    return res;
}

void dfs(int k, int tp, ll v, ll fv) { // 搜索找 lpf(x)=p_k
    if (!~tp || v != p[k]) add(get(n / (n / v)), ~tp ? mod - fv : fv);
    for (int i = k + 1; p[i] * v < n23 && i <= tot; i++) {
        for (ll x = p[i]; v * x < n23; x *= p[i]) {
            dfs(i, tp, v * x, fv * (~tp ? qpow(x, tp) : f(x)) % mod);
        }
    }
}

inline 
void upd_lpf(int k, int tp) { // 搜索入口 
    for (ll i = p[k]; i < n23; i *= p[k]) {
        dfs(k, tp, i, ~tp ? qpow(i, tp) : f(i));
    }
}

int main() {
    // 初始化 
    scanf("%lld", &n);
    n2 = sqrtl(n), n3 = cbrtl(n), n6 = pow(n, 1.l / 6);
    init(n2), n23 = n / n3;
    for (ll i = 1; i < n; i = n / (n / (i + 1))) val[++m] = i; val[++m] = n;

    // step 1
    for (int i = 1; i <= m; i++) g[i][1] = calc(val[i], 1) - 1, g[i][2] = calc(val[i], 2) - 1;
    // 1. 对 lpf(x) <= sqrt[6](n) 暴力。 
    for (int i = 1; i < t6; i++) {
        for (int j = m, x; j && p[i] <= val[j] / p[i]; j--) {
            x = get(val[j] / p[i]);
            g[j][1] = (g[j][1] - p[i] * (g[x][1] - g[p[i - 1]][1]) % mod + mod) % mod;
            g[j][2] = (g[j][2] - p[i] * p[i] % mod * (g[x][2] - g[p[i - 1]][2]) % mod + mod) % mod;
        }
    }

    for (int k = 1; k <= 2; k++) { // 枚举次数 
        for (int i = 1; val[i] < n23; i++) c[i] = (g[i][k] - g[i & i - 1][k] + mod) % mod;
        // 2. 对 x <= n^2/3 使用树状数组优化原来的递推。
        for (int i = t6; i < t3; i++) {
            for (int j = m; j && val[j] >= n23 && p[i] <= val[j] / p[i]; j--) {
                ll x = val[j] / p[i];
                if (x < n23) g[j][k] = (g[j][k] - qpow(p[i], k) * (ask(get(x)) - ask(p[i - 1])) % mod + mod) % mod;
                else g[j][k] = (g[j][k] - qpow(p[i], k) * (g[get(x)][k] - ask(p[i - 1])) % mod + mod) % mod;
            }
            upd_lpf(i, k);
        }
        for (int i = 1; i <= m && val[i] < n23; i++) g[i][k] = ask(i);
        // 3. 其余暴力。
        for (int i = t3; i <= tot; i++) {
            for (int j = m; j && val[j] >= n23 && p[i] <= val[j] / p[i]; j--) {
                g[j][k] = (g[j][k] - qpow(p[i], k) * (g[get(val[j] / p[i])][k] - ask(p[i - 1])) % mod + mod) % mod;
            }
        }
    }
    for (int i = 1; i <= m; i++) F[i] = (g[i][2] - g[i][1] + mod) % mod;

    // step 2
    for (int i = m; i; i--) {
        if (val[i] < p[t3]) { F[i] = 1; continue; }
        F[i] = (F[i] + 1 - F[p[t3] - 1] + mod) % mod;
        for (int j = t3; j <= tot && p[j] <= val[i] / p[j]; j++) {
            ll x = f(p[j]) * (F[get(val[i] / p[j])] - F[p[j]] + mod) % mod;
            F[i] = (F[i] + f((ll)p[j] * p[j] % mod) + x) % mod;
        }
    }

    // step 3
    for (int i = 1; i <= m && val[i] < n23; i++) c[i] = (F[i] - F[i & i - 1] + mod) % mod;
    for (int i = t3 - 1; i >= t6; i--) {
        for (int j = m; j && val[j] >= n23; j--) {
            for (ll x = p[i]; x <= val[j]; x *= p[i]) {
                if (val[j] / x >= n23) F[j] = (F[j] + F[get(val[j] / x)] * f(x) % mod) % mod;
                else F[j] = (F[j] + ask(get(val[j] / x)) * f(x) % mod) % mod;
            }
        }
        upd_lpf(i, -1);
    }
    for (int i = 1; i <= m && val[i] < n23; i++) F[i] = ask(i);

    // step 4
    for (int i = t6 - 1; i; i--) {
        for (int j = m; j; j--) {
            for (ll x = p[i]; x <= val[j]; x *= p[i]) {
                F[j] = (F[j] + F[get(val[j] / x)] * f(x) % mod) % mod;
            }
        }
    }

    printf("%lld\n", F[m]);
}
参考资料