Min_25筛思想及例题集

难得一批!!!

Min_25是针对积性函数f(i)求和问题，可处理范围可达到$n^{11}-n^{12}$。 所以复杂度是低于线性的。大概是$\frac{n^{\frac{3}{4}}}{\log \sqrt n}$。可以说是完爆洲阁筛甚至所有筛法。 但同时，虽然有这些优点，可理解起来相当不容易，所以很难。

粗略思路： 先将质数分为小于等于$\sqrt{n}$和大于$\sqrt{n}$两种。 再预处理所有质数函数的前缀和，利用小于等于 $\sqrt{n}$的质数进行转移，推出所有函数的和。

思想

筛质数

问题：假设我们需要求出n以内所有质数和。

质数筛求出所有小于等于$\sqrt{n}$的质数，$p_i$表示第i个质数，一共有$p_0$个质数，ps为p的前缀和。

设函数$S(x,j)=\sum_{i=2}^xi*[i为质数或i的最小质因数大于p_j]$

$设Sum(n)=\frac{n(n+1)}{2}-1.$

$显然对于任意x，S(x,0)=Sum(x)-1，我们要求的ans=S(x,p_0)$

考虑转移，如果$p_j^2>x,S(x,j)=S(x,j-1).$ 否则: $$S(x,j)=S(x,j-1)-p_j*(S(\frac{x}{p_j},j-1)-ps_{j-1})$$

一直转移、递归，就可以得到$S(n,p_0).$

可以看出，这种筛法可以筛很多的积性函数，但应该不是所有积性函数都能筛吧？

积性函数求和

通过上面方法，可以求出很多积性函数f(x)，当x为质数时的和，现在要求： $$\sum_{i=1}^nf(i)$$

考虑把上面的式子倒过来推，定义跟前面的一样：(可以用上面的S预处理$G(x,p_0))$ $设fs_x=\sum_{i=1}^xf(p_i)$

$$G(x,j)=G(x,j+1)+\sum_{k=1}(G(\frac{x}{p_j^k},j+1)-fs_j)*f(p_j^k)+\sum_{k=2}f(p_j^k)$$

则ans=G(n,0).

Code

由于递归消耗的时间会过长，所以采用循环代替递归.

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
typedef long long ll;

const int N = 1e5 + 10;

namespace Min25 {
    int prime[N], id1[N], id2[N], flag[N], ncnt, m;

    ll g[N], sum[N], a[N], T, n;

    inline int ID(ll x) {
        return x <= T ? id1[x] : id2[n / x];
    }

    inline ll calc(ll x) {
        return x * (x + 1) / 2 - 1;
    }

    inline ll f(ll x) {
        return x;
    }

    inline void init() {
        ncnt = 0, m = 0;
        T = sqrt(n + 0.5);
        for (int i = 2; i <= T; i++) {
            if (!flag[i]) prime[++ncnt] = i, sum[ncnt] = sum[ncnt - 1] + i;
            for (int j = 1; j <= ncnt && i * prime[j] <= T; j++) {
                flag[i * prime[j]] = 1;
                if (i % prime[j] == 0) break;
            }
        }
        for (ll l = 1; l <= n; l = n / (n / l) + 1) {
            a[++m] = n / l;
            if (a[m] <= T) id1[a[m]] = m; else id2[n / a[m]] = m;
            g[m] = calc(a[m]);
        }
        for (int i = 1; i <= ncnt; i++)
            for (int j = 1; j <= m && (ll)prime[i] * prime[i] <= a[j]; j++)
                g[j] = g[j] - (ll)prime[i] * (g[ID(a[j] / prime[i])] - sum[i - 1]);
    }

    inline ll Solve(ll x) {
        if (x <= 1) return x;
        return n = x, init(), g[ID(n)];
    }
}

总结

第一部分

考虑一个递推公式g(j,n)，表示n以内筛过前j个质数的倍数的合数和。 $g(j,n)=\sum_{i\epsilon f(i)最小质因子>p_j}f(i)=\sum_{i\epsilon f(i)最小质因子>p_j}i^k$

状态转移：

$$g(j,n)=g(j-1,n)-f(p_j)[g(j-1,\frac{n}{p_j})-\sum_{i=1}^{j-1}f(p_i)]=g(j-1,n)-p_j^k[g(j-1,\frac{n}{p_j})-\sum_{i=1}^{j-1}p_i^k]$$

当$p_j^2>n时，g(j,n)=g(j-1,n)$，所以只需要计算出$\sqrt{n}$以内的素数。

第二部分

预处理$id1[]，id2[]$和各种$\frac{n}{l}$的前缀和。

第三部分

设S(n,j)表示n以内所有最小质因子大于等于$p_j$的f(x)的和，即$S(n,j)=\sum_{i=1}^nF(i)[i最大质因子\ge p_j]$

状态转移：

$$S(n,j)=\sum{g_t(id[n])-sum_t[j-1]}+\sum_{k\ge y\;p_k^{e+1}\leq n}F(p_k^e)*S(\frac{n}{p_k^e},k+1)+F(p_k^{e+1})$$

$sum_t(j-1)=\sum_{i=1}^{j-1}p_i^t$为前j-1个质数的t次幂之和. $\sum{g_t(id[n])-sum_t[j-1]}$为多个多项式加减运算凑成的多项式的质数和$\sum_{p\epsilon primes}F(x)[p<P_j]$ 后面直接计算多项式函数$F(p_k^e),F(p_k^{e+1})$

递归完成即可。

例题

LOJ 6235. 区间素数个数

传送门：https://loj.ac/p/6235

求解1-n之间所有素数个数

$$g(j,n)=g(j-1,n)-f(p_j)[g(j-1,\frac{n}{p_j})-\sum_{i=1}^{j-1}f(p_i)]=g(j-1,n)-g(j-1,\frac{n}{p_j})-(j-1)$$

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
##include "bits/stdc++.h"
using namespace std;

typedef long long ll;

const int N = 1e6 + 10;

namespace Min_25 {
    int is_pr[N], pr[N], prs[N], cnt, m;
    ll n, T;
    ll id1[N], id2[N], a[N], g[N];

    inline int ID(ll x) {
        return x <= T ? id1[x] : id2[n / x];
    }

    inline ll calc(ll x) {
        return x - 1;
    }

    inline ll f(ll x) {
        return x;
    }

    void init() {
        cnt = m = 0;
        T = sqrt(n + 0.5);
        for(int i = 2;i <= T; i++) {
            if(!is_pr[i]) pr[++cnt] = i, prs[cnt] = prs[cnt - 1] + 1;
            for(int j = 1;j <= cnt && i * pr[j] <= T; j++) {
                is_pr[i * pr[j]] = true;
                if(i % pr[j] == 0) break;
            }
        }

        for(ll l = 1;l <= n; l = n / (n / l) + 1) {
            a[++m] = n / l;
            if(a[m] <= T) id1[a[m]] = m;
            else id2[n / a[m]] = m;
            g[m] = calc(a[m]);
        }

        for(int i = 1;i <= cnt; i++) {
            for(int j = 1;j <= m && (ll)pr[i] * pr[i] <= a[j]; j++) {
                g[j] = g[j] - (g[ID(a[j] / pr[i])] - prs[i - 1]);
            }
        }

    }

    ll Solve(ll x) {
        if(x <= 1) return x;
        return n = x, init(), g[ID(n)];
    }
}

void solve() {
    ll n; cin >> n;
    cout << Min_25::Solve(n) << endl;
}

signed main() {
    solve();
}

LOJ 6053. 简单的函数

传送门：https://loj.ac/p/6053

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 "bits/stdc++.h"
using namespace std;

typedef long long ll;

const ll mod = 1e9 + 7;
const int N = 1e6 + 10;

ll quick_pow(ll a, ll b) {
    ll ans = 1;
    while(b) {
        if(b & 1) ans = ans * a % mod;
         a = a * a % mod;
         b >>= 1;
    }
    return ans % mod;
}

namespace Min_25 {
    int is_pr[N], id1[N], id2[N], cnt;
    ll n, m, T;
    ll pr[N], prs[N], a[N], g1[N], g2[N];

    ll inv2;

    inline ll ID(ll x) {
        return x <= T ? id1[x] : id2[n / x];
    }

    inline ll calc1(ll x) {
        return (x - 1) % mod;
    }

    inline ll calc2(ll x) {
        return x % mod * (x + 1) % mod * inv2 % mod - 1;
    }

    void init() {
        cnt = m = 0;
        inv2 = quick_pow(2, mod - 2);
        T = sqrt(n + 0.5);
        for(int i = 2;i <= T; i++) {
            if(!is_pr[i]) pr[++cnt] = i, prs[cnt] = (prs[cnt - 1] + i) % mod;
            for(int j = 1;j <= cnt && i * pr[j] <= T; j++) {
                is_pr[i * pr[j]] = true;
                if(i % pr[j] == 0) break;
            }
        }

        for(ll l = 1;l <= n; l = n / (n / l) + 1) {
            a[++m] = n / l;
            if(a[m] <= T) id1[a[m]] = m;
            else id2[n / a[m]] = m;
            g1[m] = calc1(a[m]);
            g2[m] = calc2(a[m]);
        }

        for(int i = 1;i <= cnt; i++) {
            for(int j = 1;j <= m && pr[i] * pr[i] <= a[j]; j++) {
                g1[j] = (g1[j] - (g1[ID(a[j] / pr[i])] - (i - 1))) % mod;
                g2[j] = (g2[j] - pr[i] * (g2[ID(a[j] / pr[i])] - prs[i - 1]) % mod) % mod;
            }
        }
    }

    ll S(ll x, ll y) {
        if(x <= 1  x < pr[y]) return 0;
        ll ans = (g2[ID(x)] - prs[y - 1] - (g1[ID(x)] - (y - 1)) + (y == 1 ? 2 : 0) + mod) % mod;
        for(int i = y;i <= cnt && pr[i] * pr[i] <= x; i++) {
            ll pe = pr[i];
            for(int e = 1;pe * pr[i] <= x; e++, pe *= pr[i]) {
                ans = (ans + (pr[i] ^ e) * S(x / pe, i + 1) % mod + (pr[i] ^ e + 1) % mod) % mod;
            }
        }
        return ans % mod;
    }

    ll Solve(ll x) {
        return n = x, init(), S(n, 1) + 1;
    }
}

void solve() {
    ll n; cin >> n;
    ll ans = Min_25::Solve(n);
    cout << (ans % mod + mod) % mod << endl;
}

signed main() {
    solve();
}

本文作者：jujimeizuo
本文地址： https://blog.jujimeizuo.cn/2021/03/05/min-25/
本博客所有文章除特别声明外，均采用 CC BY-SA 3.0 协议。转载请注明出处！