min25筛

        min25筛，要理解本文最好先阅读过杜教筛。

min25筛

        min25筛是一种可以在低于线性复杂度下求积性函数的前缀和的筛法。它要求对于质数，其积性函数为多项式的函数，并且质数方幂的积性函数要容易求。

$$\sum_{i=1}^nF(i)$$

        这里$F(i)$是一个积性函数。
        为了求这个前缀和，我们定义一个函数$S(n,j)$：

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

        这里$\min(x)$取$x$的最小质因子。换句话说，$S(n,j)$记录在$1$到$n$中，最小质因子大于$P_j$（$P_j$是第$j$个质数，从$1$开始）的积性函数之和。那么要求的前缀和就是$S(n,0)$。显然当$j>0$时，$F(1)$永远不会加入到答案中，我们完全可以最后再把$1$的答案加上，于是下面的过程均不考虑$1$。
        考虑把里面的数分类成质数和合数，如果我们可以求出$1$到$n$中大于$P_j$的所有质数对应的函数值之和（把它表示成$w(n,j)$），就可以表示成：

$$S(n,j)=w(n,j)+\sum_{k>j}\sum_{e}F(P_k^e)(S(\lfloor\frac {n}{P_k^e}\rfloor,k)+[e\not=1])$$

        前者是质数的和，后者是合数的和。对于后面的求和式，我们枚举最小质因子及其次数。当$e\not=1$时，它不是一个质数，还需要加上自己的一个一，这就是后面$[e\not=1]$的意义。

质数积性函数和求法

        现在问题转化为如何求$w(n,j)$。由于只是求质数的和，我们不妨降低一下标准，把积性函数$F(x)$换成完全积性函数$f(x)$，使它们在$x$是质数时相等。
        设一个函数$g(n,j)$：

$$g(n,j)=\sum_{i=1}^nf(i)[j\in P\lor\min(i)>P_j]$$

        即$g(n,j)$是$1$到$n$中所有的质数，以及所有的最小质因子大于$P_j$的合数对应的$f(i)$之和。考虑如何转移。
        当$P_j^2>n$时，此时不会有最小质因子大于$P_j$的合数加入到答案中，于是$g(n,j)=g(n,j-1)$。
        当$P_j^2\leq n$时，将从$g(n,j-1)$处转移过来，$g(n,j-1)$必然会比$g(n,j)$多一些东西，它们即是$1$到$n$中，
最小质因子恰好就是$P_j$的合数对应的$f(i)$之和，所以需要减去它。
        假设这样的合数形如$P_js$，可以发现$s$必然也满足$\min(s)\geq P_j$。于是可以表示成$g(\lfloor\frac {n}{P_j}\rfloor,j-1)$，但是$g(\lfloor\frac {n}{P_j}\rfloor,j-1)$还多包含$1$到$P_{j-1}$的所有质数，因此还需要减去这些质数对应的积性函数之和，即为$\sum_{i=1}^{j-1}f(P_i)$。
        这里可以保证$P_{j-1}\leq \lfloor\frac {n}{P_j}\rfloor$，根据$P_j^2\leq n$可得。
        综合来看，就有：

$$g(n,j)=\begin{cases} g(n,j-1)&P_j^2>n\\ g(n,j-1)-f(P_j)(g(\lfloor\frac{n}{P_j}\rfloor,j-1)-\sum_{i=1}^{j-1}f(P_i))&P_j^2\leq n \end{cases}$$

        那如何求$1$到$n$所有质数$f(i)$之和？很明显，它就是$g(n,x)$（这里$P_x^2\geq n$，简称$g(n)$）。由于$f(i)$在质数时与$F(i)$相同，所以它就是$F(i)$的和。

代码实现

        在实现方面，我们先求$g$，再求$S$。
        预处理出$\sqrt n$范围内的所有质数，并求出相应的积性函数$f(i)$前缀和。
        注意到$g(n,j)$的第一维一定是$\lfloor\frac{n}{x}\rfloor$的形式，它最多有$2\sqrt n$个。并且可以发现只有满足$P_j^2\leq n$的质数$P_j$才有意义，于是这里的$n$和$j$其实都是$O(\sqrt n)$的。
        首先我们对$n$，进行一步整数分块，对于所有可能的$\lfloor\frac {n}{x}\rfloor$，进行离散化。离散化不建议使用$map$或排序二分，这样会引入$logn$的额外复杂度。一个优雅的方法是直接放到数组里，记录$\lfloor\frac{n}{x}\rfloor$对应的数组下标。但是这样数可能过大，数组开不下，于是可以在$\lfloor\frac{n}{x}\rfloor\leq\sqrt n$时才放入数组，否则将$\frac {n}{\lfloor\frac{n}{x}\rfloor}$放入数组。这样只需要记两个大小为$2\sqrt n$的数组就可以离散化。
        然后根据上文给出的公式，递推出$g$数组。求出$g$数组后，就可以表示出一开始提到的$w(n,j)$。
        由于$w(n,j)$中不能含有小于等于$P_j$的质数对应的$F(i)$和，于是需要减去，即：

$$S(n,j)=g(n)-\sum_{i=1}^jF(P_i)+\sum_{k>j}\sum_{e}F(P_k^e)(S(\lfloor\frac {n}{P_k^e}\rfloor,k)+[e\not=1])$$

        再根据以上式子递推即可。min25筛时间复杂度$O(\frac {n^{\frac {3}{4}}}{logn})$。
        模板题。这里先构造一个完全积性函数$f(i)$。考虑到$f(p)=p(p-1)$不是积性函数，但是可以拆成$p^2-p$，这样前者和后者都是积性函数，分两次去求。

#include <bits/stdc++.h>
using namespace std;
const int N = 200000 + 5;
const int mod = 1e9 + 7;
typedef long long ll;
int idx1[N], idx2[N], SQRT, mk[N], pr[N], tot, fs1[N], fs2[N], num;
ll to[N], g1[N], g2[N], n;
inline int F(int x)
{
    return 1ll * x * (x - 1) % mod;
}

inline int ID(ll x)
{
    return x > SQRT ? idx2[n / x] : idx1[x];
}
inline void pre(int x)
{
    for (int i = 2; i <= x; i++) {
        if (!mk[i]) pr[++tot] = i;
        for (int j = 1; j <= tot; j++) {
            if (1ll * pr[j] * i > x) break;
            mk[pr[j] * i] = 1;
            if (i % pr[j] == 0) break;
        }
    }
    for (int i = 1; i <= tot; i++) {
        fs1[i] = (fs1[i - 1] + pr[i]) % mod;
        fs2[i] = (fs2[i - 1] + 1ll * pr[i] * pr[i] % mod) % mod;
    }
}
int S(ll x, int y)
{
    if (pr[y] >= x) return 0;
    int ans = ((g2[ID(x)] - g1[ID(x)]) % mod - (fs2[y] - fs1[y]) % mod) % mod;
    for (int i = y + 1; i <= tot && 1ll * pr[i] * pr[i] <= x; i++) {
        ll pe = pr[i];
        for (int e = 1; pe <= x; e++, pe = pe * pr[i]) {
            ll xx = pe % mod;
            ans = (ans + 1ll * F(xx) * (S(x / pe, i) + (e != 1))) % mod;
        }
    }
    return ans;
}
int main()
{
    // freopen("out.txt", "w", stdout);
    scanf("%lld", &n), pre(SQRT = round(sqrt(n)));
    for (ll l = 1, r, x; l <= n; l = r + 1) {
        r = n / (n / l), x = n / l;
        if (x > SQRT)
            idx2[n / x] = ++num;
        else
            idx1[x] = ++num;
        to[num] = x, x %= mod;
        g1[num] = 1ll * x * (x + 1) % mod * 500000004ll % mod - 1;
        g2[num] = 1ll * x * (x + 1) % mod * (2ll * x + 1) % mod * 166666668ll % mod - 1;
    }
    for (int i = 1; i <= tot; i++) {
        for (int j = 1; j <= num && 1ll * pr[i] * pr[i] <= to[j]; j++) {
            g1[j] -= 1ll * pr[i] * ((g1[ID(to[j] / pr[i])] - fs1[i - 1]) % mod) % mod;
            g1[j] %= mod;
            g2[j] -= 1ll * pr[i] * pr[i] % mod * ((g2[ID(to[j] / pr[i])] - fs2[i - 1]) % mod) % mod;
            g2[j] %= mod;
        }
    }

    int ans = (S(n, 0) + 1) % mod;
    printf("%d", (ans + mod) % mod);
    return 0;
}

莫比乌斯函数与欧拉函数

        既然是一个积性函数筛，必然脱不开求两个经典函数——莫比乌斯函数与欧拉函数的前缀和。这里只提一下思路。
        对于莫比乌斯函数，它在$x$为质数时有$\mu(x)=-1$，这可以很容易地把$f(x)=1$当做处理$g$数组的完全积性函数，而且在求$S$时，$\mu(p^e)=-1[e=1]$，也十分容易求。并且发现在$e>1$时，$\mu(p^e)=0$，这样甚至可以省略第二个枚举次数的for循环。
        对于欧拉函数，其在质数时有$\phi(p)=p-1$，这处理成$f_1(x)=x$以及$f_2(x)=1$即可。求$S$时则可以使用如下公式：

$$\phi(p^e)=p^e(1-\frac{1}{p})$$

素数和以及素数数量

        根据$g(n)$函数的性质，我们可以得到一个亚线性复杂度下，求素数和以及素数数量的算法。
        规定$f(p)=1$，显然这是一个积性函数，而且它在素数中的前缀和显然就是素数的数量。修改筛法模板，可以得到如下代码：

//计算1到n的素数个数
#include <bits/stdc++.h>
using namespace std;
const int N = 5000000 + 5;
typedef long long ll;
int idx1[N], idx2[N], SQRT, mk[N], pr[N], tot, num;
ll to[N], g[N], n, fs[N];

inline int ID(ll x)
{
    return x > SQRT ? idx2[n / x] : idx1[x];
}
inline void pre(int x)
{
    for (int i = 2; i <= x; i++) {
        if (!mk[i]) pr[++tot] = i;
        for (int j = 1; j <= tot; j++) {
            if (1ll * pr[j] * i > x) break;
            mk[pr[j] * i] = 1;
            if (i % pr[j] == 0) break;
        }
    }
    for (int i = 1; i <= tot; i++) fs[i] = fs[i - 1] + 1;
}
int main()
{
    scanf("%lld", &n), pre(SQRT = round(sqrt(n)));
    for (ll l = 1, r, x; l <= n; l = r + 1) {
        r = n / (n / l), x = n / l;
        if (x > SQRT)
            idx2[n / x] = ++num;
        else
            idx1[x] = ++num;
        to[num] = x, g[num] = x - 1;
    }
    for (int i = 1; i <= tot; i++) {
        for (int j = 1; j <= num && 1ll * pr[i] * pr[i] <= to[j]; j++) {
            g[j] -= (g[ID(to[j] / pr[i])] - fs[i - 1]);
        }
    }
    printf("%lld", g[ID(n)]);
    return 0;
}

        同理，我们可以修改$f(x)=x$，进而得到素数和：

//求1到n的素数和
#include <bits/stdc++.h>
using namespace std;
const int N = 5000000 + 5;
typedef long long ll;
int idx1[N], idx2[N], SQRT, mk[N], pr[N], tot, num;
ll to[N], g[N], n, fs[N];

inline int ID(ll x)
{
    return x > SQRT ? idx2[n / x] : idx1[x];
}
inline void pre(int x)
{
    for (int i = 2; i <= x; i++) {
        if (!mk[i]) pr[++tot] = i;
        for (int j = 1; j <= tot; j++) {
            if (1ll * pr[j] * i > x) break;
            mk[pr[j] * i] = 1;
            if (i % pr[j] == 0) break;
        }
    }
    for (int i = 1; i <= tot; i++) fs[i] = fs[i - 1] + pr[i];
}
int main()
{
    scanf("%lld", &n), pre(SQRT = round(sqrt(n)));
    for (ll l = 1, r, x; l <= n; l = r + 1) {
        r = n / (n / l), x = n / l;
        if (x > SQRT)
            idx2[n / x] = ++num;
        else
            idx1[x] = ++num;
        to[num] = x, g[num] = x * (x + 1) / 2 - 1;
    }
    for (int i = 1; i <= tot; i++) {
        for (int j = 1; j <= num && 1ll * pr[i] * pr[i] <= to[j]; j++) {
            g[j] -= 1ll * pr[i] * (g[ID(to[j] / pr[i])] - fs[i - 1]);
        }
    }
    printf("%lld", g[ID(n)]);
    return 0;
}

例题

Loj6053简单的函数

        定义完全积性函数$f(p)=p\oplus 1$，这样对于$p>2$时的质数其实就是$p-1$，拆成$p$和$1$两个积性函数就和模板题一样了。处理完$g$数组后，再修改其中的值，把$2$缺少的部分补回来。然后再求$S$。

#include <bits/stdc++.h>
using namespace std;
const int N = 200000 + 5;
const int mod = 1e9 + 7;
typedef long long ll;
int idx1[N], idx2[N], SQRT, mk[N], pr[N], tot, fs1[N], fs2[N], num;
ll to[N], g1[N], g2[N], n;
inline int F(int x, int y) { return x ^ y; }

inline int ID(ll x) { return x > SQRT ? idx2[n / x] : idx1[x]; }
inline void pre(int x) {
    for (int i = 2; i <= x; i++) {
        if (!mk[i])
            pr[++tot] = i;
        for (int j = 1; j <= tot; j++) {
            if (1ll * pr[j] * i > x)
                break;
            mk[pr[j] * i] = 1;
            if (i % pr[j] == 0)
                break;
        }
    }
    for (int i = 1; i <= tot; i++) {
        fs1[i] = (fs1[i - 1] + pr[i]) % mod;
        fs2[i] = (fs2[i - 1] + 1) % mod;
    }
}
int S(ll x, int y) {
    if (pr[y] >= x)
        return 0;
    int ans = ((g1[ID(x)] - g2[ID(x)]) % mod - (fs1[y] - fs2[y]) % mod) % mod;
    if (x >= 2 && y < 1)
        ans += 2;
    for (int i = y + 1; i <= tot && 1ll * pr[i] * pr[i] <= x; i++) {
        ll pe = pr[i];
        for (int e = 1; pe <= x; e++, pe = pe * pr[i]) {
            ans = (ans + 1ll * F(pr[i], e) * (S(x / pe, i) + (e != 1))) % mod;
        }
    }
    return ans;
}
int main() {
    scanf("%lld", &n), pre(SQRT = round(sqrt(n)));
    for (ll l = 1, r, x; l <= n; l = r + 1) {
        r = n / (n / l), x = n / l;
        if (x > SQRT)
            idx2[n / x] = ++num;
        else
            idx1[x] = ++num;
        to[num] = x, x %= mod;
        g1[num] = 1ll * x * (x + 1) % mod * 500000004ll % mod - 1;
        g2[num] = (x - 1) % mod;
    }
    for (int i = 1; i <= tot; i++) {
        for (int j = 1; j <= num && 1ll * pr[i] * pr[i] <= to[j]; j++) {
            g1[j] -= 1ll * pr[i] * ((g1[ID(to[j] / pr[i])] - fs1[i - 1]) % mod) % mod;
            g1[j] %= mod;
            g2[j] -= ((g2[ID(to[j] / pr[i])] - fs2[i - 1]) % mod) % mod;
            g2[j] %= mod;
        }
    }
    int ans = (S(n, 0) + 1) % mod;
    printf("%d", (ans + mod) % mod);
    return 0;
}

HDU6834-Yukikaze and Smooth numbers

        问题问的很清楚，求$1$到$n$中，质因子不大于$k$的数的数量。
        考虑大于$k$但是不大于$n$的质数有$s$个，它们分别是$P_1,\cdots,P_s$。需要删去的数必然就是它们倍数的并集，这个并集的大小可以使用容斥原理去求。
        对于多个质数$P_i$来说，$1$到$n$中它们的倍数并集的数量就是$\lfloor\frac{n}{p}\rfloor$，其中$p=\prod_iP_i$。
        这样根据容斥原理得到：

$$\sum_{i=1}^s\lfloor\frac {n}{P_i}\rfloor-\sum_{1\leq l< r\leq n}\lfloor\frac {n}{P_lP_r}\rfloor\cdots+(-1)^{s+1}\lfloor\frac{n}{\prod_jP_j}\rfloor$$

        注意到这里每一项前的$(-1)^i$可以用莫比乌斯函数来代替，因为对于形如若干互异质数的乘积的数，它的莫比乌斯函数恰好就是$(-1)^s$，其中$s$是质因子的数量。但是这样也需要枚举所有的乘积组合，没有多少优化的余地。考虑到对于含有两个及以上的相同质因子的数，莫比乌斯函数为$0$，这样可以把其它的数也加进来，乘上它们对应的莫比乌斯函数即可。
        综上我们需要枚举最小质因子大于$k$，不小于$n$的所有数，则整个题的答案就可以表示成：

$$n+\sum_{\min(i)>k}^{i\leq n}\mu(i)\lfloor\frac{n}{i}\rfloor$$

        这样就可以进行整除分块，然后配合min25筛解决这个问题。
        具体来说，先考虑$k^2>n$的情况，此时最小质因子大于$k$的必然都是质数，这样只需要求质数对应的前缀和，根本用不到$S$（见上文min25筛原理），只用$g$即可。将求和的范围设置成$[k+1,r]$，在这个区间整除分块，然后乘上对应的莫比乌斯函数和。注意到min25筛中需要给定$n$，并且所有的端点都需要是$\lfloor\frac{n}{x}\rfloor$的形式，这里的整除分块恰好满足这个要求（考虑到整除分块的右端点是$\frac{n}{n/l}$）。唯一可能不满足的就是$k+1$时，这里额外算一步即可。
        当$k^2\leq n$时，情况就比较复杂了。此时既可能有质数，又可能有合数，需要用到$S$来求，但原理和上面相同。我们需要先找出最后一个不大于$k$的质数，然后把求和范围设置成$[1,n]$，整除分块然后利用$S$求最小质因子大于$k$的对应莫比乌斯函数和。

#include <bits/stdc++.h>
using namespace std;
const int N = 500000 + 5;
typedef long long ll;
int idx1[N], idx2[N], SQRT, mk[N], pr[N], tot, fs[N], num, to[N], g[N], nn;
inline int F(int x, int y)
{
    return y > 1 ? 0 : -1;
}

inline int ID(int x)
{
    return x > SQRT ? idx2[nn / x] : idx1[x];
}
inline void pre(int x)
{
    for (int i = 2; i <= x; i++) mk[i] = 0;
    for (int i = 2; i <= x; i++) {
        if (!mk[i]) pr[++tot] = i;
        for (int j = 1; j <= tot; j++) {
            if (1ll * pr[j] * i > x) break;
            mk[pr[j] * i] = 1;
            if (i % pr[j] == 0) break;
        }
    }
    for (int i = 1; i <= tot; i++) fs[i] = fs[i - 1] + 1;
}
int S(int x, int y)
{
    if (x == 0) return 0;
    if (pr[y] >= x) return 0;
    int ans = -g[ID(x)] + fs[y];
    for (int i = y + 1; i <= tot && 1ll * pr[i] * pr[i] <= x; i++) ans -= S(x / pr[i], i);
    return ans;
}
inline void solve_g(int x)
{
    nn = x, num = tot = 0, pre(SQRT = sqrt(x));
    for (int l = 1, r, x; l <= nn; l = r + 1) {
        r = nn / (nn / l), x = nn / l;
        if (x > SQRT)
            idx2[nn / x] = ++num;
        else
            idx1[x] = ++num;
        to[num] = x, g[num] = x - 1;
    }
    for (int i = 1; i <= tot; i++) {
        for (int j = 1; j <= num && 1ll * pr[i] * pr[i] <= to[j]; j++) {
            g[j] -= g[ID(to[j] / pr[i])] - fs[i - 1];
        }
    }
}
int main()
{
    // DEBUG;
    int t;
    scanf("%lld", &t);
    while (t--) {
        int n, k, ans = 0, cnt = 0;
        scanf("%lld%lld", &n, &k);
        if (k >= n) {
            ans = n;
        } else if (1ll * k * k > n) {
            int l = k + 1, r = n / (n / l);
            solve_g(k), cnt -= g[ID(k)];
            solve_g(n), cnt += g[ID(r)];
            ans -= (n / l) * cnt, l = r + 1, ans += n;
            while (l <= n) r = n / (n / l), ans -= (g[ID(r)] - g[ID(l - 1)]) * (n / l), l = r + 1;
        } else {
            solve_g(n);
            int t = 1, l = 1, r;
            while (t <= tot && pr[t] <= k) ++t;
            --t, ans = n;
            while (l <= n) {
                r = n / (n / l);
                ans += (S(r, t) - S(l - 1, t)) * (n / l);
                l = r + 1;
            }
        }
        printf("%d\n", ans);
    }
    return 0;
}