Min_25 筛 题解

· · 题解

前言

由于其由 Min_25 发明并最早开始使用，故称「Min_25 筛」。

其可以在 O(\frac{n^{\frac{3}{4}}}{logn}) 的时间复杂度下解决一类 积性函数 的前缀和问题。

要求： f(p) 是一个关于 p 的项数较少的多项式或可以快速求值； f(p^{c}) 可以快速求值。

摘自 OI wiki

解析

例题

积性函数 f(p^k) = p^k(p^k - 1) 求：\sum_{i = 1} ^ n f(i)

准备 1

规定 p 代表质数, minp(i) 代表 i 的最小质因子（其中 minp(p) = p），则：

\sum_{i = 1} ^ n f(i) = \sum_{p \leq n} f(p) + \sum_{i \leq n \ and \ i \notin p} f(i)

这一步表示将答案拆成质数的答案和合数的答案。 那么对于合数部分：

根据积性函数的性质，则有：

提**尽** $i$ 的最小质因子 $p_e$，将 $i$ 分成两部分： $p_e ^ {c_e} \times \frac{i}{p_e ^ {c_e}}

则上式可变为：

f(i) = f(p_e ^ {c_e}) \times f(\frac{i}{p_e ^ {c_e})})

那么合数部分的答案可转化为：

\sum_{p ^ e \leq n} f(p ^ e) \sum_{1 \leq i \leq n \ and \ p ^ e \mid i \ and \ minp(\frac{i}{p ^ e}) > p} f(\frac{i}{p ^ e})

代表的意思是：

枚举质数 p 及它的指数 e，

在可以提的合数 i 中提出 p ^ e，将合数 i 转化为 p ^ e \times \frac{i}{p ^ e}，即上式乘法分配律展开的结果。

然后再判断如下条件：

这两个条件均可通过 minp(\frac{i}{p ^ e}) 是否大于 p 来判断：

对于 minp(\frac{i}{p ^ e}) < p 的 i，那么不满足条件 1。

对于 minp(\frac{i}{p ^ e}) = p 的 i，那么不满足条件 2。

这些均可通过对 i 分解质因数来理解。

那么对于 f(p ^ e) 由题意可以直接求，对于 f(\frac{i}{p ^ e}) 则是同样的求法，直到 \frac{i}{p ^ e} 是 p_k ^ {c_k} 可以直接求。

而上述的约束则保证了不会多求漏求。

准备 2

所以我们现在需要求出所有 f(p ^ k) 的值。

因为 n 的范围很大，所以我们没法线性筛求出 f 的前缀和。

考虑 DP。。。

设 h_k(i) = i ^ k，显然 h_k 是一个完全积性函数。

设： g_k(n,j) = \sum_{i = 1} ^ n [i \in Prime \ or \ minp(i) > p_j] h_k(i)

其中 p_j 表示第 j 个质数。

换句话说 g_k(n,j) 也就是 1 ~ n 中质数或者最小质因数大于 p_j 的 i 的 h_k(i) 的前缀和。

考虑转移。。。 对于 $g_k(n, j - 1)$ 根据定义，其中包含两部分： * $g_k(n, j)

所以可由 g_k(n, j - 1) 减去多余的部分得到 g_k(n, j)。

所以有状态转移方程如下：

g_k(n,j) = g_k(n, j - 1) - p_{j} ^ k (g_k(\frac{n}{p_{j}}, j - 1) - g_k(p_{j - 1}, j - 1))

从状态转移方程可以看出：

最小质因数等于 p_j 的 h_k(i) 的和即为：

原因如下： 根据定义不难得到 $g_k(p_{j - 1}, j - 1)$ 即为前 $j - 1$ 个质数的 $k$ 次方和。 而对于 $g_k(\frac{n}{p_j}, j - 1)$，包含了两部分： * $1$ ~ $p_{j - 1}$ 的质数和 * 最小质因数大于等于 $p_j$ 的合数和，及 $p_{j}$ ~ $\frac{n}{p_j}$ 的质数和 注意第二点，如果我们将这些最小质因数大于等于 $p_j$ 的合数或者后一部分的质数同时乘上 $p_j$ 那么我们得到的是不是都是最小质因数等于 $p_j$ 的合数，这点应该很好理解。 那也就是说得到了我们想要筛去的合数。 因为我们只知道上述两部分的和 $g_k(\frac{n}{p_j}, j - 1)$，所以在乘的过程中第一部分也会被乘上一个 $p_j$，那么我们应该把这一部分质数减去。 也就是 $g_k(p_j, j - 1)$。 又因为 $h_k$ 是完全积性函数，而我们谈论的和是 $h_k$ 的累加，根据乘法分配律，我们可以对 $g_k$ 运用完全积性函数的性质（$f(ab) = f(a)f(b)$）。 那么要得到最小质因数等于 $p_j$ 的 $h_k(i)$ 的和，我们可以直接给我们的第二部分乘上 $h_k(p_j)$ 的值，也就是 $p_j ^ k$，得到我们要求的函数值。 所以关于状态转移方程的解释就结束了。 关于初始化： 显然对于 $g_k(n, 0) = \sum_{i = 1} i ^ k

求解

而我们得到的 g_k 有什么用呢？

设 g(n, j) = \sum_{i = 1} ^ n [i \in Prime \ or \ minp(i) > p_j] f(i)。

因为 f(p) = p(p - 1) = p ^ 2 - p。

那么将满足条件的 f(i) 累加，就得到两部分：

如果令上文的 k 分别等于 1, 2。

那么就有： g(n, j) = g_2(n, j) - g_1(n, j)。

如果设 p_x 为小于等于 \sqrt{n} 的最大的质数，那么 1 ~ n 的质数的 f 的和 g(n, x) 记为 g(n)。

设 S(n, j) = \sum_{minp(i) > p_j} f(i)。

显然我们也需要 DP 。。。

同样的，考虑 S(n, j) 的组成：

设 sum_k(x) 表示小于等于 p_x 的质数的 k 次方的和。

那么质数部分就等于：

g(n) - (sum_{2}(j) - sum_{1}(j))

对于合数部分：

我们模仿 g_k 的求法

因为 f 不是完全积性函数，而是积性函数，我们不能像求 g_k 那样只提一个 p_j。

所以我们要把 p_j 提尽。

得到：

\sum_{k > j \ and \ p_k ^ e \leq n} f(p_k ^ e)(S(\frac{n}{p_k ^ e} + [e \neq 1]))

利用 S(\frac{n}{p_k ^ e} 的思想和求 g_k 的时候完全一样，这里就不在赘述了，实在不懂的同学可以手推一下。

后面加上一个 [e \neq 1] 是因为如果 e > 1，那么p ^ e 在我们提 p_k ^ e 时没有算进去，而 e = 1 时，在质数部分就算了，就不用了算进去。

下标的离散化

因为 n = 1e10，直接开肯定不行。 由性质 \lfloor \frac{\lfloor \frac{x}{a} \rfloor}{b} \rfloor = \lfloor \frac{x}{ab} \rfloor，在递归的时候，无论我们把 n 除以几，得到的都是 n 除以某个数，所以我们可以直接只保存 n / x 的值。 来自 wucstdio 大佬

代码

#include<cstdio>
#include<cmath>
#include<cstring>

using namespace std;

typedef long long LL;

const LL mod = 1e9 + 7, inv2 = 500000004, inv3 = 333333336;
const int N = 1e5 + 5;

LL n, Prime[N + 5], num, sp1[N + 5], sp2[N + 5], tot =  0, w[3 * N], g1[3 * N], g2[3 * N], sqr;
int id1[N + 5], id2[N + 5];
bool notPrime[N + 5];

void init() {
    memset(notPrime, false, sizeof(notPrime));
    num = 0;
    for(int i = 2; i <= N; i++) {
        if(!notPrime[i]) {
            Prime[++num] = i;
            sp1[num] = (sp1[num - 1] + i) % mod;
            sp2[num] = (sp2[num - 1] + 1ll * i * i % mod) % mod;
        }
        for(int j = 1; j <= num && Prime[j] * i <= N; j++) {
            notPrime[i * Prime[j]] = true;
            if(i % Prime[j] == 0) break;
        }
    }
}

LL s(LL x, int y) {
    if(Prime[y] >= x) return 0;
    LL k = x <= sqr ? id1[x] : id2[n / x];
    LL ans = (g2[k] - g1[k] + mod - (sp2[y] - sp1[y]) + mod) % mod;
    for(int i = y + 1; i <= num && Prime[i] * Prime[i] <= x; i++) {
        LL P = Prime[i];
        for(int e = 1; P <= x; e++, P = P * Prime[i]) {
            LL xx = P % mod;
            ans = (ans + xx * (xx - 1) % mod * (s(x / P, i) + (e != 1))) % mod;
        }
    }
    return ans % mod;
}

int main() {
    init();
    scanf("%lld", &n);
    sqr = sqrt(n);
    for(LL l = 1, r; l <= n; l = r + 1) {
        r = (n / (n / l));
        w[++tot] = n / l;
        g1[tot] = w[tot] % mod;
        g2[tot] = g1[tot] * (g1[tot] + 1) / 2 % mod * (2 * g1[tot] + 1) % mod * inv3 % mod - 1;
        if(g2[tot] < 0) g2[tot] += mod;
        g1[tot] = g1[tot] * (g1[tot] + 1) / 2 % mod - 1;
        if(g1[tot] < 0) g1[tot] += mod;
        if(w[tot] <= sqr) id1[w[tot]] = tot;
            else id2[n / w[tot]] = tot;
    }
    for(int i = 1; i <= num; i++)
        for(int j = 1; j <= tot && Prime[i] * Prime[i] <= w[j]; j++) {
            LL k = w[j] / Prime[i];
            k = k <= sqr ? id1[k] : id2[n / k];
            g1[j] -= Prime[i] * (g1[k] - sp1[i - 1] + mod) % mod;
            g2[j] -= Prime[i] * Prime[i] % mod * (g2[k] - sp2[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", (s(n, 0) + 1) % mod); 
    return 0;
}