Min_25 筛

发表于 分类于 ，

概述:

Min_25 筛 (aka. 扩展埃氏筛) 可以做到在渐进亚线性的时间复杂度下求出大部分数论函数的前缀和.

与杜教筛相比其最大的优点是适用范围广， 符合如下要求的函数 \(f\) 均可使用:

  • 在质数 \(p\) 处的取值是个与 \(p\) 有关的项数较小的多项式.
  • \(f(p^e)\) 易求.

实现:

抛开 \(f(1)\), 可以将 \(S_f(n)=\sum_{i=2}^{n}f(i)\) 分为质数部分与合数部分:

\[ S_f(n) = \sum_{i \in \mathbb{P}}f(i) + \sum_{i\notin \mathbb{P}}f(i) \]

其中 \(\mathbb{P}\) 表示全体素数集合.

Part 1: 质数前缀

首先思考如何求得质数部分的和.

由于前面的条件 1, 可以先将多项式拆开， 求得每个幂次的质数部分和再加起来即可.

定义一个函数 \(g(n,m)\):

\[ g(n,m) = \sum_{i=1}^n[i\in\mathbb{P} \lor \mathrm{LPF}(i)>P_m]i^k \]

其中 \(\mathrm{LPF}(i)\) 表示 \(i\) 的最小质因子 (Least Prime Factor), \(P_m\) 为从小到大的第 m 个质数，\(k\) 即为分离出来的幂次.

这个函数是整个筛法最精华的部分，同时也是另一个名字扩展埃氏筛的由来: 这个函数就像是进行了 m 次埃筛后剩下的部分.

考虑进行 dp 转移，有:

其中 \(S_{P^k}(m-1)\) 为前 \(m-1\) 个质数的 \(k\) 次方和，可以使用线性筛处理.

当 \(P_m \le \sqrt{n}\) 时， 此时需要减去最小质因子为 \(P_m\) 的数的影响。由于 \(i^k\) 是完全积性函数 , 所以可以提出一个 \(P_m^k\), 同时 \(g(\lfloor n/P_m \rfloor , m-1)\) 中还有小于 \(P_m\) 的质数部分，需要减去这部分，即为 \(S_{P^k}(m-1)\).

同时根据一个结论:

\[ \left\lfloor \frac{\left\lfloor \frac{n}{a} \right\rfloor}{b} \right\rfloor = \left\lfloor \frac{n}{ab} \right\rfloor \]

所以我们只需要求出所有的 \(g(\lfloor n/i \rfloor)\) 即可.

存储可以使用根号分治，使用 \(idl_x\) 表示小于等于 \(\sqrt{n}\) 的 \(g(x)\) 下标，使用 \(idh_{\lfloor n/x \rfloor}\) 表示大于 \(\sqrt{n}\) 的 \(g(x)\) 下标.

一点点证明：对于大于 \(\sqrt{n}\) 的 \(x=\lfloor n/i\rfloor\), \(\lfloor n/x \rfloor\) 是唯一的。由数论分块， 所证即为 \(i=\lfloor n/ \lfloor n/i \rfloor \rfloor\) \[ \begin{aligned} \left\lfloor \frac{n}{\lfloor \frac ni \rfloor} \right\rfloor & = \left\lfloor \frac{ni}{n-n \mod i} \right\rfloor \\ & \le \left\lfloor \frac{ni}{n-i+1} \right\rfloor \\ \because \lfloor n/i\rfloor > \sqrt{n} & \Rightarrow i < \sqrt{n} \\ \therefore \left\lfloor \frac{ni}{n-i+1} \right\rfloor & < \left\lfloor \frac{ni+n-i^2+1}{n-i+1} \right\rfloor = i+1 \\ \Rightarrow \left\lfloor \frac{n}{\lfloor \frac ni \rfloor} \right\rfloor & < i+1 \\ \because \left\lfloor \frac{n}{\lfloor \frac ni \rfloor} \right\rfloor & \ge i \\ \therefore \left\lfloor \frac{n}{\lfloor \frac ni \rfloor} \right\rfloor & = i \end{aligned} \]

Part 2: 合数部分

\([1,n]\) 中合数的最小质因子一定小于等于 \(\sqrt{n}\), 对合数部分枚举最小质因子的幂次， 有:

\[ \sum_{i\notin \mathbb{P}}f(i) \\ = \sum_{p\in\mathbb{P}, p\le\sqrt{n}, p^e\le n} f(p^e) \sum_{i\le n/p^e, \mathrm{LPF}(i)>p} f(i) + [e>1] \]

与前面的第一部分类似，定义函数 \(S(n,m)\):

\[ S(n,m) = \sum_{i=2}^n[\mathrm{LPF}(i)>P_m]f(i) \]

则有:

\[ \begin{aligned} & S(n,m) \\ = & G(n,m') - \sum_{i=1}^{m-1}f(P_i) \\ & + \sum_{P_m<P_i\le\sqrt{n}, P_i^e\le n} f(P_i^e) (S(\left\lfloor \frac{n}{P_i^e} \right\rfloor, i) + [e>1]) \end{aligned} \]

其中 \(G\) 即为第一部分求得的各幂次项之和，\(m'\) 是使得 \(P_{m'}^2\le n\) 的最大值.

最后所求的前缀和即为 \(S(n,0)+f(1)\).

时间复杂度:

事实上 Min_25 筛的时间复杂度比较玄学， 第一部分的时间复杂度使用类似于杜教筛的分析方法可以得到是 \(O(n^{3/4}/\log n)\). 而第二部分的复杂度是 \(O(n^{1-\epsilon})\), 且 \(\epsilon\) 随 \(n\) 增长无限趋近于 0 (所以是渐进亚线性). 事实上 Min_25 筛在 \(10^{11}\) 内的数据表现良好.

示例代码:

luogu5325 【模板】Min_25 筛

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
#include <cstdio>
#include <cmath>
using LL = long long;
const LL MP(1000000007);
const LL I2((MP+1)/2), I3((MP+1)/3);    //2和3的逆元
const int Msn(100500);

inline constexpr LL madd(const LL& a, const LL& b) {    //模加法
    return a+b>=MP ? a-MP+b : a+b;
}
inline constexpr LL mmns(const LL& a, const LL& b) {    //模减法
    return a-b<0 ? a+MP-b : a-b;
}

LL n,sqn;
int pc(0);
bool isv[Msn];
LL prm[Msn/10]; //质数P
LL sp1[Msn/10], sp2[Msn/10];
int dc(0);  //分块数量
int idl[Msn], idh[Msn];
LL v[Msn<<1];
LL g1[Msn<<1], g2[Msn<<1];  //g使用滚动数组
void init() {   //线性筛质数
    sqn = sqrt(n);
    for(int i(2);i<=sqn;++i) {
        if(!isv[i]) {
            prm[++pc] = i;
            sp1[pc] = madd(sp1[pc-1],i);    //s_{p^1}
            sp2[pc] = madd(sp2[pc-1],1ll*i*i%MP);   //s_{p^2}
        }
        for(int j(1);j<=pc && i*prm[j]<=sqn;++j) {
            isv[i*prm[j]] = true;
            if(i%prm[j]==0) {
                break;
            }
        }
    }
    LL l(1),r;  //使用数论分块寻找所有的n/i
    for(; l <= n; l= r + 1) {
        r = n/(n/l);
        v[++dc] = n / l;
        if(v[dc] <= sqn) {
            idl[v[dc]] = dc;
        } else {
            idh[r] = dc;
        }
        /*g1, g2的初值为除1外的前缀和与前缀平方和*/
        g2[dc] = v[dc] % MP;
        g1[dc] = g2[dc] * (g2[dc] + 1) % MP * I2 % MP;
        g2[dc] = g1[dc] * (g2[dc] * 2 + 1) % MP * I3 % MP;
        g1[dc] = mmns(g1[dc], 1);
        g2[dc] = mmns(g2[dc], 1);
    }
}

void calc_g() { //第一部分
    for(int i(1);i<=pc;++i) {
        for(int j(1);v[j]>=prm[i]*prm[i];++j) {
            LL tx = v[j]/prm[i];
            LL tp = tx<=sqn ? idl[tx] : idh[n/tx];  //tx对应的g下标
            g1[j] = mmns(g1[j], prm[i]*mmns(g1[tp], sp1[i-1])%MP);
            g2[j] = mmns(g2[j], prm[i]*prm[i]%MP*mmns(g2[tp], sp2[i-1])%MP);
        }
    }
}

LL S(LL x, int k) { //第二部分
    if(prm[k]>=x) {
        return 0;
    }
    int pg = x<=sqn ? idl[x] : idh[n/x];    //x对应的g下标
    LL ret = mmns(mmns(g2[pg],g1[pg]), mmns(sp2[k],sp1[k]));
    for(int y=k+1;y<=pc && prm[y]*prm[y]<=x;++y) {
        ret = madd(ret, mmns(prm[y]*prm[y]%MP,prm[y]) * S(x/prm[y],y) % MP);
        for(LL pe=prm[y]*prm[y];pe<=x;pe*=prm[y]) {
            LL mpe = pe%MP;
            ret = madd(ret, mmns(mpe*mpe%MP, mpe) * (S(x/pe,y) + 1) % MP);
        }
    }
    return ret;
}

int main() {
    scanf("%lld",&n);
    init();
    calc_g();
    printf("%lld",S(n,0)+1);
    return 0;
}