Min_25筛代码及模板

最新推荐文章于 2020-10-07 15:30:24 发布

KXL5180

最新推荐文章于 2020-10-07 15:30:24 发布

阅读量577

点赞数 2

分类专栏：积性函数与莫比乌斯反演

本文链接：https://blog.csdn.net/KXL5180/article/details/104038247

版权

积性函数与莫比乌斯反演专栏收录该内容

11 篇文章 0 订阅

订阅专栏

LOJ简单的函数

以这道题为例子，讲一下min_25筛如何用代码实现。
首先min_25是一种亚线性筛，可以处理1e9以上的数据，老版min_25筛的复杂度为 $O(\frac{n^{0.75}}{logn})$ ，新版已经优化到 $O(n^{\frac{2}{3}})$ ，复杂度，我也不会证明。
据说常数比较大，确实是跑得快，但不过差距也不是特别大，我们还是说老版min_25筛吧。
min_25算法原理
上面这个博客写得比较简单，而且好理解，建议先大概看懂上面那个博客在看其他博客。
如果有错误，恳请各位指出来！
首先我先把min25的几部分摆出来吧，min25就是难在质数求和。

首先是运用的要求

首先对于积性函数 $f (x)$ ，当 $x\in prime$ 时， $f (x)$ 是一个关于 $x$ 的简单多项式，简单来说就是项数不能太多，幂次也不能太高，等会儿实现过程你就会发现。
$f(p^{k})，p \in prime$ 能够快速求得，这里快速求得，起码应该是log级别吧，等会儿代码实现过程你就会发现了。

我们发现其实对于一个积性函数来说，其实我们见过第一个条件很多积性函数都满足。
比如 $\mu(x)$ ，发现当为质数时，等于-1，而且第二个条件也很容易知道满足，所以莫比乌斯函数，可以用min25筛求他的前缀和。
同理，对于欧拉函数 $\varphi(x)$ ，当为质数时，等与 $p - 1$ ，第二个条件 $\varphi(p^{k})=p^{k-1}(p-1)$ ，同样也可以用min25筛。
再比如求 $\sum_{i=0}^{x-1} gcd(x,i)\mu(i)$ ，这个前缀和。总的来说很多积性函数的前缀和都可以秒杀。
我们看题 $f(1)=1,f(p^c)=p\bigoplus c,f(ab)=f(a)f(b)[gcd(a,b)==1]$

首先我们可以发现当 $p = 2$ 是 $f (2) = 3$ ，其他质数时 $f (p) = p - 1$ ，即 $f (p) = p - 1 + 2 [p = = 2]$ ，可以认为是多项式。同时第二个条件也满足。

接下讲一下第一部分，求g函数、

首先要处理关于质数的值，因为积性函数质数点的值是关于质数的多项式，我们接下来的假设是有意义的。

首先我们设一个函数 $g(n,j)=\sum_{i=1}^{n} i^k,(i\in prime) or(min(p)\geq p_{j},p|i,p\in prime)$ 换成中文就是， $i$ 为质数，或者 $i$ 的最小质因子大于 $p_{j}$ ， $p_j$ 表示第 $j$ 个质数。
然后我们发现上面的每一个 $n$ 和 $j$ 都对应了一个状态，我们考虑一些动态规划的思想，找一找这个式子的递推方程。我们考虑从 $g(n,j-1)\to g(n,j)$ ，当 $p_j^2>n$ 时，表面在 $1$ 到 $n$ 的范围内，不存在最小质因子为 $p_j$ 的数，因此对于总和没有贡献，即 $g(n,j)=g(n,j-1),p_j^2>n$ 时。
接下来考虑 $p_j^2\leq n$ ，这时从 $g(n,j-1)\to g(n,j)$ ，明显符合条件的数减少了，因此需要减去 $x$ 即 $g (n, j) = g (n, j - 1) - x$ ，我们现在考虑一下，这个 $x$ 显然应该是最小质因子是 $p_{j}$ 的一部分的值，我们减去最小质因子大于等于 $p_j$ 的那一部分值，我们不能从 $g (n, j - 1)$ 出发，这时，我们注意到， $i^k$ 是一个完成积性函数，因此我们可以这样状态转义，减去这一部分 $p_j^{k}g(\frac{n}{p_j},j-1)$ 表示所有含有 $p_j$ 的值，因为这里多减去了质数的部分的，所以要加上这一部分的值 $p_j^kg(p_{j-1},j-1)$ ，相当于一个容斥，如果不能很好的理解，可以手推模拟一下。
状态转移方程 $g(n,j)=\begin{cases}g(n,j-1), p_j^2>n \\g(n,j-1)-p_j^k[g(\frac{n}{p_j},j-1)-g(p_{j-1},j-1)], p_j^2\leq n\end{cases}$

我们显然可以知道 $g(p_j,j)=\sum_{i=1}^{j}p_j^k$ ，这里我们用 $∣ p ∣$ 表示质数集的大小，因此 $g (n, ∣ p ∣)$ 就表示范围满足条件的质数的和，因此这里的质数我们，筛到 $\sqrt n$ 就可以了，同时 $g(n,0)=\sum_{i=2}^{n}i^k$ ，注意没有1哦。
其实这个一个过程就很像埃筛的思想，先筛出前 $j - 1$ 个质数的情况，在推出 $j$ 的情况。
这里我们上面要求 $f (p)$ 是关于 $p$ 的简单多项式，就是我们在求 $g$ 函数的时候方便求出来而已。我们来看看代码具体如何实现。

首先我们可以看到上面的递推，我们可以发现题目给出n的范围一般很大的，直接开数组是开不下的，但不过因为转移的时候只和整除有关，所以我们在代码实现时考虑数论分块后的值，这样的值大约有 $2\sqrt n$ 个。
我们首先要递推 $g$ 函数，因此要确定 $g$ 的边界，即 $g(n,0)=\sum_{i=2}^{n}f(i)=\sum_{i=2}^{n}\sum a_{x}i^{x}$ ，注意没有包括1，其实1包不包括不重要，最后加上就是了，记住这里的 $f (i)$ 是假的，就是把 $f (i)$ 当成全部是质数的情况来计算，才有后面的一部分，因为 $f (p)$ 是多项式，因此少不了自然数幂求和，以及拉格朗日插值等等，所以就和上面的对应起来了，如果复杂了，肯定不好求，而且效率堪忧，建议自然数的幂次不要超过100，否则插值，会花费大量时间。
接下来就是正常的类似动态规划的递推了，注意应该在上一步中处理好整除分块后的值。

int getid(ll x) {
    if (x <= sqr)
        return id1[x];
    return id2[n / x];
}//sqr=sqrt(n);
void calc1() {
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l), w[++m] = n / l;//假设全部是质数。
        h[m] = (w[m] - 1) % mod;//常数项
        g[m] = (w[m] % mod) * ((w[m] + 1) % mod) % mod;//一次项
        g[m] *= inv2, g[m]--;
        if (w[m] <= sqr)
            id1[w[m]] = m;
        else
            id2[r] = m;
    }
    for (int j = 1; j <= cnt; j++) {
        for (int i = 1; i <= m && pri[j] * pri[j] <= w[i]; i++) {
            int k = getid(w[i] / pri[j]);
            (g[i] -= 1ll * pri[j] * (g[k] - pre[j - 1]) % mod) %= mod;
            (h[i] -= h[k] - j + 1) %= mod;
        }
    }
}

时间复杂的的证明…我不会。。。。看这里吧时间复杂度

怎么求和

说到这里到底该怎么求和呢。。。

我们还是设一个函数 $s(n,j)=\sum_{i=2}^{n}f(i)[min(p) \geq p_j,p|i,p\in prime]$ ，中文就是 $i$ 的最小质因子大于等于 $p_j$ 注意这里和上面不同是大于等于，不是大于了，并且没有质数的必须计算的条件。
下面我们来思考，如何利用 $g$ 函数来辅助递推 $s$ 函数。我可以把 $s$ 分成两部分，一部分是质数的答案，一部分是合数的答案，关于质数部分的答案我们就可以了利用 $g$ 函数来递推了即 $\sum_{i=1}^{j-1}f(p_i)$ .这是质数部分的。
下面来考虑合数部分的值了。首先我们吧合数的最小质因子提出来,这些质因子都是大于等于 $p_j$ 的，可以得到 $\sum_{k=j}^{|p|}\sum_{e=1}^{p_{k}^{e+1}\leq n}[f(p_{k}^{e})s(\frac{n}{p_{k}^{e}},k+1)+f(p_{k}^{e+1})]$ 为什么这样就可以了？因为 $f (i)$ 是一个积性函数因此把最小质因子全部提出来然后相乘是肯定没有问题，因此我们可以枚举最小质因子以及幂次，就可以了，后面的那个一尾巴，是因为在处理过程种所以的 $f(p_{k}^{i})$ 都被筛掉了，因此要补充进来。
因此 $\sum_{i=1}^{j-1}f(p_i)+\sum_{k=j}^{|p|}\sum_{e=1}^{p_{k}^{e+1}\leq n}[f(p_{k}^{e})s(\frac{n}{p_{k}^{e}},k+1)+f(p_{k}^{e+1})]$ ，我们就可以利用递归进行状态转移了。同时从上面也可以看出来为什么要求条件二成立了。

代码实现如下，我可以利用dfs实现状态转移，剩下的就是按照公式来就是了。那么我们要求的就是 $s (n, 1)$ 最后加上1对应的函数值就可以了。

ll solve(ll x, int y) {
    if (x <= 1 || pri[y] > x)
        return 0;//边界条件
    int k = getid(x);//同样注意分块
    ll ret = (g[k] - pre[y - 1] - h[k] + y - 1) % mod;
    if (y == 1)
        ret += 2;//注意2的特殊性
    for (int i = y; i <= cnt && 1ll * pri[i] * pri[i] <= x; i++) {
        ll t1 = pri[i], t2 = 1ll * pri[i] * pri[i];
        for (int e = 1; t2 <= x; ++e, t1 = t2, t2 *= pri[i]) {
            (ret += ((1ll * solve(x / t1, i + 1) * (pri[i] ^ e) % mod + (pri[i] ^ (e + 1)) % mod))) %= mod;
        }
    }
    return ret;
}

好了LOJ6053已经讲解完了。
总结一下

首先看函数满不满足条件。
其次找 $f (p)$ ,然后分开项计算 $g$ 函数。
然后计算 $s$ 函数

这应该是标准模板吧。
然后min25还可以处理特殊质因子，最大质因子，次大质因子，最小质因子等问题。。。
感觉迟早要完。。。
下面来一道，洛谷上面的模板题
https://www.luogu.com.cn/problem/P5325

#include "bits/stdc++.h"

using namespace std;
inline int read() {
    int x = 0;
    bool f = 1;
    char c = getchar();
    for (; !isdigit(c); c = getchar()) if (c == '-') f = 0;
    for (; isdigit(c); c = getchar()) x = (x << 3) + (x << 1) + c - '0';
    if (f) return x;
    return 0 - x;
}
typedef long long ll;
const int maxn = 100010;
const ll mod = 1000000000 + 7;
ll quick(ll a, ll n, ll p) {
    ll ans = 1;
    for (; n; n >>= 1, a = a * a % p)
        if (n & 1) ans = ans * a % p;
    return ans;
}
int vis[maxn], sqr, cnt = 0;
ll suf1[maxn], suf2[maxn], inv2, inv6, n, pri[maxn];
void init(int lim) {
    vis[1] = 1;
    for (int i = 2; i <= lim; i++) {
        if (!vis[i]) {
            pri[++cnt] = i;
            suf1[cnt] = suf1[cnt - 1] + i;
            suf2[cnt] = (suf2[cnt - 1] + 1ll * i * i % mod) % mod;
        }
        for (int j = 1; j <= cnt && i * pri[j] <= lim; j++) {
            vis[i * pri[j]] = 1;
            if (i % pri[j] == 0) break;
        }
    }
}
int id1[maxn], id2[maxn], m = 0;
ll w[maxn << 1], g[maxn << 1], h[maxn << 1];
int getid(ll x) {
    if (x <= sqr) return id1[x];
    return id2[n / x];
}
ll sum1(ll x) {
    x %= mod;
    return (x + 1) * x % mod * inv2 % mod;
}
ll sum2(ll x) {
    x %= mod;
    return x * (x + 1) % mod * (2ll * x % mod + 1) % mod * inv6 % mod;
}
void calc1() {
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        w[++m] = n / l;
        g[m] = sum2(w[m]) - 1;
        h[m] = sum1(w[m]) - 1;
        if (w[m] <= sqr) id1[w[m]] = m;
        else id2[r] = m;
    }
    for (int j = 1; j <= cnt; j++) {
        for (int i = 1; i <= m && pri[j] * pri[j] <= w[i]; i++) {
            int k = getid(w[i] / pri[j]);
            (g[i] -= pri[j] * pri[j] % mod * (g[k] - suf2[j - 1]) % mod) %= mod;
            (h[i] -= pri[j] * (h[k] - suf1[j - 1]) % mod) %= mod;
        }
    }
}
ll solve(ll x, int y) {
    if (x <= 1 || pri[y] > x) return 0;
    int k = getid(x);
    ll ret = (g[k] - suf2[y - 1] - h[k] + suf1[y - 1]) % mod;
    for (int i = y; i <= cnt && pri[i] * pri[i] <= x; i++) {
        ll t1 = pri[i], t2 = pri[i] * pri[i];
        for (int e = 1; t2 <= x; ++e, t1 = t2, t2 *= pri[i]) {
            ll pk = quick(pri[i], e, mod);
            (ret += pk * (pk - 1) % mod * solve(x / t1, i + 1) % mod) %= mod;
            pk = pk * pri[i] % mod;
            (ret += pk * (pk - 1) % mod) %= mod;
        }
    }
    return ret;
}
int main() {
    inv2 = quick(2, mod - 2, mod);
    inv6 = quick(6, mod - 2, mod);
    cin >> n;
    sqr = sqrt(n);
    init(sqr);
    calc1();
    ll ans = solve(n, 1) + 1;
    cout << (ans + mod) % mod << endl;
    return 0;
}

过年在家无事(太堕落了。。)无事于是想了想和杜教筛比一比速度。。。。
就拿洛谷上面那一道杜教筛模板题来对比一下

上面min25下面杜教筛。。
看来差距不是很大，杜教筛稍优

#include "bits/stdc++.h"

using namespace std;
inline int read() {
    int x = 0;
    bool f = 1;
    char c = getchar();
    for (; !isdigit(c); c = getchar()) if (c == '-') f = 0;
    for (; isdigit(c); c = getchar()) x = (x << 3) + (x << 1) + c - '0';
    if (f) return x;
    return 0 - x;
}
typedef long long ll;
const int maxn = 100010;
const ll mod = 1000000000 + 7;
ll quick(ll a, ll n) {
    ll ans = 1;
    for (; n; n >>= 1, a = a * a)
        if (n & 1) ans = ans * a;
    return ans;
}
int vis[maxn], sqr, cnt = 0;
ll suf1[maxn], suf2[maxn], inv2, inv6, n, pri[maxn];
void init(int lim) {
    vis[1] = 1;
    for (int i = 2; i <= lim; i++) {
        if (!vis[i]) {
            pri[++cnt] = i;
            suf1[cnt] = suf1[cnt - 1] + i;
        }
        for (int j = 1; j <= cnt && i * pri[j] <= lim; j++) {
            vis[i * pri[j]] = 1;
            if (i % pri[j] == 0) break;
        }
    }
}
int id1[maxn], id2[maxn], m = 0;
ll w[maxn << 1], g[maxn << 1], h[maxn << 1];
int getid(ll x) {
    if (x <= sqr) return id1[x];
    return id2[n / x];
}

ll sum2(ll x) {
    return x * (x + 1) / 2;
}

void calc1() {
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        w[++m] = n / l;
        g[m] = sum2(w[m]) - 1;
        h[m] = w[m] - 1;
        if (w[m] <= sqr) id1[w[m]] = m;
        else id2[r] = m;
    }
    for (int j = 1; j <= cnt; j++) {
        for (int i = 1; i <= m && pri[j] * pri[j] <= w[i]; i++) {
            int k = getid(w[i] / pri[j]);
            g[i] -= 1ll * pri[j] * (g[k] - suf1[j - 1]);
            h[i] -= h[k] - j + 1;
        }
    }
}

ll phi(ll p, ll k) {
    return quick(p, k - 1) * (p - 1);
}

ll solve(ll x, int y) {
    if (x <= 1 || pri[y] > x) return 0;
    int k = getid(x);
    ll ret = g[k] - suf1[y - 1] - h[k] + y - 1;
    for (int i = y; i <= cnt && 1ll * pri[i] * pri[i] <= x; i++) {
        ll t1 = pri[i], t2 = 1ll * pri[i] * pri[i];
        for (int e = 1; t2 <= x; ++e, t1 = t2, t2 *= pri[i]) {
            ret += 1ll * solve(x / t1, i + 1) * phi(pri[i], e) + phi(pri[i], e + 1);
        }
    }
    return ret;
}
ll mu(ll k) {
    if (k == 1) return -1;
    return 0;
}

ll solve1(ll x, int y) {
    if (x <= 1 || pri[y] > x) return 0;
    int k = getid(x);
    ll ret = -h[k] + y - 1;
    for (int i = y; i <= cnt && 1ll * pri[i] * pri[i] <= x; i++) {
        ll t1 = pri[i], t2 = 1ll * pri[i] * pri[i];
        for (int e = 1; t2 <= x; ++e, t1 = t2, t2 *= pri[i]) {
            ret += 1ll * solve1(x / t1, i + 1) * mu(e) + mu(e + 1);
        }
    }
    return ret;
}

int main() {
    int T = read();
    init(100000);
    while (T--) {
        cin >> n;
        sqr = sqrt(n);
        m = 0;
        calc1();
        cout << solve(n, 1) + 1 << " ";
        cout << solve1(n, 1) + 1 << "\n";
    }
    return 0;
}

KXL5180

关注

2
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
Min_25筛代码及模板

LOJ简单的函数以这道题为例子，讲一下min_25筛如何用代码实现。首先min_25是一种亚线性筛，可以处理1e9以上的数据，老版min_25筛的复杂度为O(n0.75logn)O(\frac{n^{0.75}}{logn})O(lognn0.75)，新版已经优化到O(n23)O(n^{\frac{2}{3}})O(n32)，复杂度，我也不会证明。据说常数比较大，确实是跑得快，但不过差距...
复制链接

扫一扫