【算法】min_25 筛（扩展埃氏筛）笔记

我们考虑这样一个问题。

$ans=\sum_{i=1}^nf(i)$
其中

$1\leq n\leq 10^{10}$

其中 $f(i)$ 是一个非常奇怪的函数，并不像 $\mu(i),\varphi(i),i\varphi(i)$ 那样具有那么好的性质。但是满足以下条件：

1. 若 $p$ 为质数，则 $f(p)$ 是一个关于 $p$ 的多项式，比如 $\mu(p)=-1,\varphi(p)=p-1$ .

2. 若 $p$ 为质数， $e$ 为正整数，则 $f(p^e)$ 可以很快求出。（通常是 $O(1)$ ）

3. $f(n)$ 为积性函数。

然后就可以使用 min_25 筛了。（顾名思义是 min_25 发明的）

首先，我们需要知道 min_25 筛是埃氏筛的一个拓展，它的思想很大一部分借助于埃氏筛。

回想一下埃氏筛，我们是每次将最小质因子为 $P_i$ 的合数筛去，剩下的就是质数。

我们知道这些最小质因子至多为 $\sqrt{n}$ ，所以合数可以通过枚举最小质因子来计算，质数我们则使用另外的方法。

首先我们看质数的怎么做。

$\sum_{i=1}^n[i\in Prime]f(i)$

根据条件 1，我们知道 $f(i)$ 是一个多项式，这样的话我们可以按照次数将 $f(i)$ 拆成 $i^k$ 之和，因为 $i^k$ 是一个完全积性函数（很快就有用的）。

$\sum_{i=1}^n[i\in Prime]i^k$

为了计算这个，我们需要引入一个辅助数组 $g(n,j)$ 。（鬼知道是怎么想到的）

$g(n,j)=\sum_{i=1}^n[i\in Prime \ or \ minp(i)>P_j]i^k$

其中 $minp(i)$ 表示 $i$ 的最小质因子，所以

$\sum_{i=1}^n[i\in Prime]i^k=g(n,|P|)$

既然我们要使用质数，所以我们可以先用欧拉筛把所有 $\leq \sqrt{n}$ 的质数筛出来，同时还要预处理 $\sum_{i=1}^j[i\in Prime]i^k$ .

我们考虑 dp 计算。既然是埃氏筛，我们就要在 $g(n,j-1)$ 中最小质因子为 $P_j$ 的合数筛去。

我们假设 $P_j^2\leq n$ ，否则肯定不行。

首先，由于它是完全积性函数，所以 $P_j$ 可以直接提出来，剩下的减去 $i\leq \lfloor\frac{n}{P_j}\rfloor$ 中的数就可以了。

这些数中要求质因子 $\geq P_j$ ，所以是 $g(\lfloor\frac{n}{P_j}\rfloor,j-1)$ ，但是这里面质数被重复计算了，所以要减去里面的质数。

$g(n,j)=g(n,j-1)-P_j^k(g(\lfloor\frac{n}{P_j}\rfloor,j-1)-\sum_{i=1}^{j-1}[i\in Prime]i^k)$

但是这样的话是 $O(n*|P|)$ 的，时间和空间都承受不了。但是我们发现我们可以使用一个优化。

我们发现 $g(n,j)$ 中 $n$ 只有 $O(\sqrt{n})$ 种取值，因为每次递归的时候是 $n$ 变为 $\lfloor\frac{n}{P_j}\rfloor$ ，而我们发现

$\lfloor\frac{\lfloor\frac{a}{b}\rfloor}{c}\rfloor=\lfloor\frac{a}{bc}\rfloor$

所以 $n$ 只会变为 $\lfloor\frac{n}{x}\rfloor$ ，于是我们就直接 “手动” 离散化，这个可以看代码。

然后 $g(n,j)$ 的第二维也可以滚动数组滚掉。所以时间 $O(\sqrt{n}*|P|)$ ，空间 $O(\sqrt{n})$ .

预处理部分终于结束了，接下来我们考虑计算答案，首先我们还是需要一个辅助数组。

$S(n,j)=\sum_{i=1}^n[minp(i)>P_j]f(i)$

像上面说的一样，分质数和合数两类计算。

$S(x,y)=g(x,|P|)-\sum_{i=1}^{y}f(P_i)+\sum_{P_k\leq \sqrt{x},k>j}\sum_{e>0,P_k^e\leq x}f(p_k^e)(S(\lfloor\frac{x}{P_k^e}\rfloor,k)+[e>1])$

前面两项指的是质数的部分，后面的和式是枚举合数的最小质因子 $P_k$ 和它的次数 $e$ ，

这个跟 $g$ 不同， $S$ 是要按照第二维倒着计算的，但是我们也可以使用递归的方法来计算。

$S(n,0)$ 就是最终答案。

要注意的是， $1$ 既不是质数也不是合数，所以最后要加上。

至于上面说的那个手动离散化，我们要开两个数组 $id1$ 和 $id2$ ，分别记录 $\leq \sqrt{n}$ 和 $>\sqrt{n}$ 的部分的数值的编号。这样就不用 map 了，可以省掉一个 log。

至于时间复杂度？我也不知道，总之跟杜教筛差不多，甚至有时候比杜教筛还要快。

有人说是什么 $O(\frac{n^{\frac{3}{4}}}{\log n})$ ，或者什么 $O(n^{1-\epsilon})$ 。总之差不多就可以了。

#include<bits/stdc++.h>
#define Rint register LL
using namespace std;
typedef long long LL;
const int N = 1000003, mod = 1e9 + 7, inv2 = 500000004, inv3 = 333333336;
LL n, Sqr, pri[N], tot, pre1[N], pre2[N], ind1[N], ind2[N], g1[N], g2[N], w[N], cnt;
bool notp[N];
inline void init(LL m){
    notp[0] = notp[1] = true;
    for(Rint i = 2;i <= m;i ++){
        if(!notp[i]){
            pri[++ tot] = i;
            pre1[tot] = (pre1[tot - 1] + i) % mod;
            pre2[tot] = (pre2[tot - 1] + i * i) % mod;
        }
        for(Rint j = 1;j <= tot && i * pri[j] <= m;j ++){
            notp[i * pri[j]] = true;
            if(!(i % pri[j])) break;
        }
    }
}
inline LL S(LL x, int y){
    if(pri[y] >= x) return 0;
    LL k = (x <= Sqr) ? ind1[x] : ind2[n / x];
    LL ans = (g2[k] - g1[k] + pre1[y] - pre2[y] + mod + mod) % mod;
    for(Rint i = y + 1;i <= tot && pri[i] * pri[i] <= x;i ++){
        LL pe = pri[i];
        for(Rint e = 1;pe <= x;e ++, pe *= pri[i]){
            LL xx = pe % mod;
            ans = (ans + xx * (xx - 1) % mod * (S(x / pe, i) + (e > 1))) % mod;
        }
    }
    return ans % mod;
}
int main(){
    scanf("%lld", &n);
    Sqr = sqrt(n);
    init(Sqr);
    for(Rint i = 1, last;i <= n;i = last + 1){
        last = n / (n / i);
        w[++ cnt] = n / i;
        LL xx = w[cnt] % mod;
        g1[cnt] = (xx * (xx + 1) / 2 + mod - 1) % mod;
        g2[cnt] = (xx * (xx + 1) / 2 % mod * (2 * xx + 1) % mod * inv3 % mod + mod - 1) % mod;
        if(n / i <= Sqr) ind1[w[cnt]] = cnt;
        else ind2[last] = cnt;
    }
    for(Rint i = 1;i <= tot;i ++)
        for(Rint j = 1;j <= cnt && pri[i] * pri[i] <= w[j];j ++){
            LL k = (w[j] / pri[i] <= Sqr) ? ind1[w[j] / pri[i]] : ind2[n / (w[j] / pri[i])];
            g1[j] -= pri[i] * (g1[k] - pre1[i - 1] + mod) % mod;
            g2[j] -= pri[i] * pri[i] % mod * (g2[k] - pre2[i - 1] + mod) % mod;
            if(g1[j] < 0) g1[j] += mod;
            if(g2[j] < 0) g2[j] += mod;
        }
    printf("%lld", (S(n, 0) + 1) % mod);
}

UOJ188#【UR #13】Sanrd

题目链接：UOJ

这道题，也算是 min_25 的一个基础应用吧。。。

我们要求

$\sum_{i=1}^nf(i)$ ，其中

$f(i)$ 表示

$i$ 的次大质因子。

按照套路，我们设

$S(n,j)=\sum_{i=1}^n[minp(i)>P_j]f(i)$

所以

$S(n,j)=\sum_{k>j,P_k\leq \sqrt{n}}\sum_{e>0,P_k^{e+1}\leq n}(S(\lfloor\frac{n}{P_k^e}\rfloor,k)+P_k\sum_{i=P_k}^{\lfloor\frac{n}{P_k}\rfloor}[i\in Prime])$

素数个数用 min_25 筛可以很快求出来。

#include<bits/stdc++.h>
#define Rint register LL
using namespace std;
typedef long long LL;
const int N = 1000003;
LL Sqr, pri[N], tot, g[N], w[N], id1[N], id2[N], cnt;
bool notp[N];
inline void init(int m){
    notp[0] = notp[1] = true;
    for(Rint i = 2;i <= m;i ++){
        if(!notp[i]) pri[++ tot] = i;
        for(Rint j = 1;j <= tot && i * pri[j] <= m;j ++){
            notp[i * pri[j]] = true;
            if(!(i % pri[j])) break;
        }
    }
}
inline LL solve(LL n, LL x, int y){
    if(x <= 1) return 0;
    LL ans = 0;
    for(Rint k = y + 1;k <= tot && pri[k] * pri[k] <= x;k ++){
        for(Rint pe = pri[k];pe * pri[k] <= x;pe *= pri[k]){
            LL kk = (x / pe <= Sqr) ? id1[x / pe] : id2[n / (x / pe)];
            ans += solve(n, x / pe, k) + pri[k] * (g[kk] - k + 1);
        }
    }
    return ans;
}
inline LL solve(LL n){
    Sqr = sqrt(n);
    tot = cnt = 0;
    init(Sqr);
    for(Rint i = 1, last;i <= n;i = last + 1){
        w[++ cnt] = n / i;
        last = n / w[cnt];
        g[cnt] = w[cnt] - 1;
        if(w[cnt] <= Sqr) id1[w[cnt]] = cnt;
        else id2[last] = cnt;
    }
    for(Rint i = 1;i <= tot;i ++)
        for(Rint j = 1;j <= cnt && pri[i] * pri[i] <= w[j];j ++){
            LL k = (w[j] / pri[i] <= Sqr) ? id1[w[j] / pri[i]] : id2[n / (w[j] / pri[i])];
            g[j] -= g[k] - i + 1;
        }
    return solve(n, n, 0);
}
int main(){
    LL l, r;
    scanf("%lld%lld", &l, &r);
    printf("%lld\n", solve(r) - solve(l - 1));
}

6 条评论

Icontofig · 2020年3月24日 10:17 下午

那个 $g(n,j)$ 应该是埃氏筛法的扩展的思路吧，感觉应该这就是说的 min25 筛是埃氏筛法的扩展的意思

Remmina · 2020年3月25日 12:16 上午

额，我当时学这玩意的时候似乎是听说这两者是同一个东西，不太记得清了。

AThousandMoon · 2019年5月30日 1:10 下午

https://www.cnblogs.com/AThousandMoons/p/10827172.html

AThousandMoon · 2019年5月27日 3:58 下午

对不起， $S(n,j)$ 的计算公式中后面是 $k$ 而不是 $k+1$ 。
本人能力有限，请各位指出错误。（或者各位看看代码）

Remmina · 2019年5月30日 12:57 下午

OvO…
才看懂您的意思，很抱歉耽误这么久，现在我就修改

Remmina · 2019年5月30日 1:10 下午

已经提升了 “投稿者” 权限，可以修改和删除自己已发布的文章啦，您可以自行修改了（我没找到您所描述的问题之处）