Min_25 筛大学习

@liuzhangfeiabc 进行讲解，然后觉得似乎这个比较重要，那就记一下吧。

这个东西普适性很强，但是不好写。

依旧是求一个积性函数 \(f(n)\) 的前缀和，只不过这个 \(f(n)\) 不像 \(\varphi(n)\) 那样特殊。

这里我们要求 \(f\) 满足：

• 这个函数在质数 \(p\) 处的值 \(f(p)\) 是一个关于 \(p\) 的低阶多项式，例如：\(f(p)=p^2-2p+1\)。
• 这个函数在质数幂 \(p^c\) 处的值 \(f(p^c)\) 可以快速求出（不一定是低阶多项式，只要能 \(O(1)\) 单点计算即可）。

整个算法时间复杂度为 \(O(\frac{n^{\frac{3}{4}}}{\ln n})\)。

这个算法分为两大步：

• 第一步：求出所有质数处 \(f(p)\) 的和。

• 第二步：利用你算出的 \(f(p)\)，通过枚举“最小质因子”把合数的贡献加回来，得到完整的前缀和。

第一步

1.引入

我们先从更简单的问题开始：求 \(1\) 到 \(n\) 中质数的个数。\(n\le10^{11}\)。

显然不可能线性筛。但是我们可以利用埃氏筛的做法：用小于等于 \(\sqrt n\) 的质数把合数筛掉。

然后利用 DP 的思想，设 \(f(i,j)\) 表示不超过 \(i\) 的数中，本身就是质数或不含有前 \(j\) 个质数作为质因子的数有多少个。显然 \(f(i,0)=i\)，最终所求为 \(f(n,c)-1\)，\(c\) 是你筛出来的质数个数。

考虑转移。我们从小到大用每个质数 \(p_j\) 去筛。那么从 \(f(i,j-1)\) 转移到 \(f(i,j)\)，就意味着我们要把 \(i\) 以内最小质因子恰好是 \(p_j\) 的合数去掉。

设有一 \(i\) 以内的合数，它能被表示成 \(p_j\times m\) 的形式，其中 \(m\) 的最小质因子至少是 \(p_j\)。那么我们有 \(p_j\times m\le i\)。

于是可以得到 \(m\le\left\lfloor\frac{i}{p_j}\right\rfloor\)，且 \(m\) 要么是质数，要么是最小质因子大于等于 \(p_j\) 的合数。可以发现，这种 \(m\) 的取值一共有 \(f(\lfloor i/p_j\rfloor,j-1)\) 种。

所以我们是不是直接把这些 \(m\) 扔出去就行了？显然会出问题。

因为你发现，\(f(\lfloor i/p_j\rfloor,j-1)\) 中 \(p_1\) 到 \(p_{j-1}\) 这 \(j-1\) 个质数和 \(m=1\) 时对应的 \(p_j\) 被我们扔出去了。所以我们要把它们加回来，于是

\[f(i,j)=f(i,j-1)-f(\lfloor i/p_j\rfloor,j-1)+j \]

但是你发现又出问题了。因为当 \(\lfloor i/p_j\rfloor<p_j\)，\(\lfloor i/p_j\rfloor\) 中根本不可能会有 \(p_j\)，直接转移就是错的。然后又因为 \(i<p_j^2\) 的时候，\(f(i,j)\) 里面恰好代表所有不超过 \(i\) 的质数（以及 \(1\)），这个值不会随 \(j\) 增大而改变。所以转移的时候我们从后往前刷表，如果发现 \(i<p_j^2\)，说明更前面的值从此不再发生变化。那么我们此时直接停止刷表即可。

可以证明这样做复杂度是 \(O(\frac{n^{\frac{3}{4}}}{\ln n})\) 的。

2.深入

现在我们正式来看 \(f(p)\) 是一个低阶多项式的情况。

显然我们可以将 \(f(p)\) 拆分成若干幂次进行计算，然后再整合到一起。所以我们实际上要解决就是“求不超过 \(n\) 的质数的 \(k\) 次方和”这个问题。

依葫芦画瓢，设 \(f(i,j)\) 表示不超过 \(i\) 的数中，本身就是质数或不含有前 \(j\) 个质数作为质因子的数的 \(k\) 次方和是多少。显然 \(f(i,0)=\sum_{j=1}^ij^k\)。因为 \(k\) 比较小，可以光速幂 \(O(1)\) 求出每一项。

再设 \(sp_j\) 表示前 \(j\) 个质数的 \(k\) 次方和。转移的时候，我们要把对应的 \(m\) 的贡献扔出去（一共是 \(p_j^kf(\lfloor i/p_j\rfloor,j-1)\)），然后把多扔出去的贡献加回来（一共是 \(p_j^ksp_{j-1}+p_j^k\)），那么有转移

\[f(i,j)=f(i,j-1)-p_j^kf(\lfloor i/p_j\rfloor,j-1)+p_j^ksp_{j-1}+p_j^k \]

或许你会产生疑问：为什么我们一定要把 \(f(p)\) 拆成 \(p\) 的幂次分别计算最后整合到一起，而不是直接计算呢？其实你可以发现，对于你要扔出去的那些 \(p_j\times m\)，它们的 \(k\) 次方是完全积性函数，但是你这个低阶多项式本身可能并不是一个完全积性函数。

第二步

现在我们要利用刚才求出的 \(f\) 来求积性函数 \(g(n)\) 的前缀和，满足对于质数 \(p\)， \(g(p)\) 是关于 \(p\) 的低阶多项式，\(p^c\) 处的值 \(g(p^c)\) 容易求出。

仍然依葫芦画瓢，设 \(s(i,j)\) 表示不超过 \(i\) 的数中不包含前 \(j\) 个质数作为质因子的数的 \(g\) 的和。则答案为 \(s(n,0)\)。设 \(sp_j\) 为前 \(j\) 个质数的低阶多项式之和。

容易得到 \(s(i,c)=f(i,c)-sp_c\)，\(c\) 还是表示你筛出来的质数个数。

求解 \(s\) 的过程与 \(f\) 刚好相反：要从大到小枚举质数 \(p_j\)，从 \(s(i,j)\) 转移到 \(s(i,j-1)\)。整个过程中要把一些数“加进去”。

然后你发现 \(s(i,j-1)\) 相比 \(s(i,j)\)，多出来了最小质因子是 \(p_j\) 的数。所以我们枚举 \(p_j\) 的次数 \(t\)，则有

\[s(i,j-1)=s(i,j)+\sum_{t\ge1}g({p_j}^t)s\left(\left\lfloor\frac{i}{{p_j}^t}\right\rfloor,j\right) \]

刷表法实现复杂度仍为 \(O(\frac{n^{\frac{3}{4}}}{\ln n})\)。因为 \(i<p_j^2\) 的时候 \(s(i,j)\) 中只含有质数，那么就可以通过 \(f\) 和 \(sp\) 求出。

但是显然你都写成递归形式了那你写递归何乐而不为呢。

方便起见，我们把 \(s(i,j)\) 的定义中去掉 \(1\) 这一项，则所求即为 \(s(n,0)+1\)。那么 \(s(i,j)\) 中包含了质数和最小质因子大于等于 \(p_{j+1}\) 的合数。

对于质数部分，贡献是 \(f(i,c)-sp_{j}-1\)。对于合数部分，枚举最小质因子 \(p_k\)（\(k>j\)）的次数 \(t\)，然后把 \(g({p_k}^t)s\left(\left\lfloor\frac{i}{{p_k}^t}\right\rfloor,k\right)+g(p_k^{t+1})\) 累加到 \(s(i,j)\) 上即可。多出来的一项是因为我们把 \(1\) 的贡献去掉了，要加回来。

写递归的复杂度是 \(O(n^{1-\varepsilon})\)，虽然没有刷表法更优，但是实际运行起来比刷表法更快。

【模板】Min_25 筛

嗯对，板子。直接给代码吧。

Code

i64 n;
i64 pri[N], tot;
bool prime[N];
i64 f1[N], g1[N];//二次项的贡献
i64 f2[N], g2[N];//一次项的贡献
i64 sp[N];
i64 t;

inline void init () {
    pri[0] = 1;
    for (int i = 2; i <= N - 15; ++ i) {
        if (!prime[i]) pri[++ tot] = i;
        for (int j = 1; j <= tot; ++ j) {
            if (i * pri[j] > N - 15) break;
            prime[i * pri[j]] = true;
            if (i % pri[j] == 0) break;
        }
    }
    for (int i = 1; i <= tot; ++ i) {
        sp[i] = (sp[i - 1] + pri[i] * (pri[i] - 1) % mod) % mod;
    }
}//预处理

inline i64 sum1(i64 x) {
    x %= mod;
    return x * (x + 1) % mod * inv2 % mod;
}//等差数列求和

inline i64 sum2 (i64 x) {
    x %= mod;
    return x * (x + 1) % mod * (2 * x + 1) % mod * inv6 % mod;
}//平方和

inline void solve () {
    t = sqrt (n) + 10;
    for (i64 l = 1, r; l <= n; l = r + 1) {
        i64 j = n / l;
        r = n / j;
        if (j <= t) {
            f1[j] = (sum2(j) - 1 + mod) % mod;
            g1[j] = (sum1(j) - 1 + mod) % mod;
        } else {
            f2[r] = (sum2(j) - 1 + mod) % mod;
            g2[r] = (sum1(j) - 1 + mod) % mod;
        }
    }
    for (int i = 1; i <= tot; ++ i) {
        i64 x = pri[i];
        for (i64 l = 1, r; l <= n; l = r + 1) {
            i64 j = n / l;
            r = n / j;
            if (j < x * x) break;//保证复杂度的灵魂所在
            i64 v1, v2;
            if (j / x <= t) v1 = f1[j / x], v2 = g1[j / x];
            else v1 = f2[n / (j / x)], v2 = g2[n / (j / x)];
            i64 t1, t2;
            if (x - 1 <= t) t1 = f1[x - 1], t2 = g1[x - 1];
            else t1 = f2[n / (x - 1)], t2 = g2[n / (x - 1)];
            if (j <= t) {
                f1[j] = (f1[j] - x * x % mod * (v1 - t1 + mod) % mod + mod) % mod;
                g1[j] = (g1[j] - x % mod * (v2 - t2 + mod) % mod + mod) % mod;
            } else {
                f2[r] = (f2[r] - x * x % mod * (v1 - t1 + mod) % mod + mod) % mod;
                g2[r] = (g2[r] - x % mod * (v2 - t2 + mod) % mod + mod) % mod;
            }
        }
    }
}//第一步

inline i64 get (i64 x) {
    if (x <= 1) return 0;
    if (x <= t) return (f1[x] - g1[x] + mod) % mod;
    return (f2[n / x] - g2[n / x] + mod) % mod;
}

i64 S (i64 x, int y) {
    if (x <= pri[y]) return 0;
    i64 ans = (get(x) - sp[y] + mod) % mod;
    for (int i = y + 1; i <= tot && 1ll * pri[i] * pri[i] <= x; ++ i) {
        i64 p = pri[i];
        i64 p1 = p, p2 = p * p;
        for (int c = 1; p2 <= x; ++ c, p1 *= p, p2 *= p) {
            i64 v1 = (p1 % mod) * ((p1 - 1) % mod) % mod;
            i64 v2 = (p2 % mod) * ((p2 - 1) % mod) % mod;
            ans = (ans + v1 * S (x / p1, i) % mod + v2) % mod;
        }
    }
    return ans;
}//第二步

int main () {
    init ();
    cin >> n;
    solve ();
    cout << (S(n, 0) + 1) % mod << endl;
    return 0;
}