cyh_toby 的洛谷博客

cyh_toby 的洛谷博客

The harder you roll, the more you get.

Min_25 筛

posted on 2022-01-12 12:33:49 | under 笔记 |

简介

对于积性函数 $f(x)$,如果能快速知道 $f(p^k)$ 的值(其中 $p$ 为质数,下同),且 $f(p)$ 可以表示为项数较少的多项式(或能表示成若干项使得每一项都是完全积性函数),那么 $\rm min\_25$ 筛可以在 $O(\frac{n^{\frac{3}{4}}}{\log n})$ 时间复杂度内求出 $\sum_{i=1}^nf(i)$。

由于其思想与埃氏筛类似,所以也被称作扩展埃氏筛。

内容

第一步

推导

目标:求 $g(n)=\sum_{p\le n} f(p)$ 。

不妨设 $f(p)$ 是完全积性函数,如果不是可以尝试拆成若干项完全积性函数,分别求然后相加。

首先要线性筛求出 $\sqrt n$ 以内的质数。

$g(n)$ 很难直接求解,考虑用 DP 计算。设 $g(n,j)=\sum_{i=1}^nf(i)[i\text{是质数或其最小质因子}>p_j]$,其中 $p_j$ 表示第 $j$ 个质数,那么我们要的就是 $g(n,k),k\text{为最小的满足}p_k\ge \sqrt n$。考虑从 $j-1$ 变到 $j$ ,那么最小质因子为 $p_j$ 的合数会被筛掉,那么它们的贡献要减去。则有转移$$ g(n,j)=g(n,j-1)-f(p_j)\left(g\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right)-g(p_{j-1},j-1)\right) $$ 系数 $f(p_j)$ 表示由于 $f(p)$ 是完全积性函数,所以可以把它从后面提出来。 $g\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right)$ 表示考虑所有 $p_j$ 的倍数,它们除以 $p_j$ 之后,最小质因子 $>p_{j-1}$ 的合数以及所有质数的贡献,应当减去。但是,这些质数中可能有 $\le p_{j-1}$ 的,它们在之前就被筛掉过了,所以要加回来,也就是 $g(p_{j-1},j-1)$。

由于有公式 $\lfloor\frac{\lfloor\frac{a}{b}\rfloor}{c}\rfloor=\lfloor\frac{a}{bc}\rfloor$ ,因此容易发现上述式子只会用到形如 $\lfloor\frac{n}{x}\rfloor,x\le n$ 的点处的 DP 值,即第一项的状态数是 $O(\sqrt n)$(实际实现的时候注意状态数是 $2\sqrt n$)。我们预处理出这 $O(\sqrt n)$ 个数,把他们离散化,顺带求出 $g(x,0)$,然后 DP 即可。

注意到上述转移式中还有一项 $g(p_{j-1},j-1)$,我们只处理了所有形如 $\lfloor\frac{n}{x}\rfloor$ 的数,是否会漏掉某些质数 $p$ 呢?其实不会漏。注意到,能用来转移的 $p$ 一定满足 $p \le \sqrt n$。我们只需要证明 $\forall k \le \sqrt n, \exists x,k=\lfloor\frac{n}{x}\rfloor$。

证明:

  1. 若 $k=\lfloor\sqrt n\rfloor$,设 $n=k^2+d$,则由于 $k^2\le n < (k+1)^2$,故 $d\in[0,2k]$。
    1. $d\in [0,k)$,那么 $\lfloor\frac{n}{k}\rfloor=k+\lfloor\frac{d}{k}\rfloor=k$。
    2. $d\in [k,2k]$,那么 $\lfloor\frac{n}{k+1}\rfloor=\lfloor\frac{k^2+k+d-k}{k+1}\rfloor=k+\lfloor\frac{d-k}{k+1}\rfloor=k$。
  2. 若 $k < \lfloor\sqrt n\rfloor$ 即 $k \le \sqrt n - 1$,假设存在 $i,\lfloor\frac{n}{i+1}\rfloor<k<\lfloor\frac{n}{i}\rfloor$,此时 $k$ 恰好夹在两个连续的 $\lfloor\frac{n}{x}\rfloor$ 之间,即不可被表出。则 $\frac{n}{i+1}<k$,故 $n<k(i+1)$,从而 $k<\lfloor\frac{n}{i}\rfloor<\lfloor\frac{k(i+1)}{i}\rfloor=k+\lfloor\frac{k}{i}\rfloor$。另一方面,$\frac{n}{i+1}<k<\sqrt n$,所以 $i+1>\sqrt n$,于是 $i>\sqrt n - 1 \ge k$,因此 $\lfloor\frac{k}{i}\rfloor=0$,于是得到 $k<k$,故假设不成立,原命题成立。

所以,上述担心就是多虑了。

我看到很多博客、题解都额外预处理了 $\sqrt n$ 内所有素数处的函数值的前缀和,但因为我“粗心”,代码中直接调用 $g$ 数组,却顺利通过了一堆题目,写笔记时意识到这一点,于是琢磨出了奥妙所在。其实,在看 zzt 集训队论文时注意到他说只需要用到 $\lfloor\frac{n}{x}\rfloor$ 处的值,就心生疑惑,现在终于明白了。因此,在我看来,预处理根号内素数处的函数前缀和的值是不必要的

时间复杂度

$$ \begin{aligned}&\sum_{i\le \sqrt n}O(\pi(\sqrt i))+\sum_{i\le \sqrt n}O(\pi(\sqrt\frac{n}{i}))\\ =&\sum_{i\le \sqrt n}O(\pi(\sqrt\frac{n}{i}))\\ =&\sum_{i\le \sqrt n}O\left(\frac{\sqrt\frac{n}{i}}{\log\sqrt\frac{n}{i}}\right)\\ =&O\left(\int_1^{\sqrt n} \frac{\sqrt\frac{n}{x}}{\log\sqrt\frac{n}{x}}\right)\\ =&O\left(\frac{n^\frac{3}{4}}{\log n}\right) \end{aligned} $$

第二步

推导

目标:求 $S(n)=\sum_{i\le n} f(i)$。与第一步类似,设 $S(n,j)=\sum_{i=1}^nf(i)[i\text{的最小质因子}>p_j]$。但此处 $f$ 不需要再拆分成单项式,直接是原函数即可(因为不需要依赖于完全积性,只需要积性即可)(但要能快速计算 $f(p^k)$ 的值)。

方法一

考虑把贡献拆成质数的和合数的,合数枚举最小质因子以及次数,于是有转移:$$ S(n,j)=g(n)-g(p_j)+\sum_{j<k,p_k\le \sqrt n,1\le e,p_k^e\le n} f(p_k^e)\left(S\left(\left\lfloor\frac{n}{p_k^e}\right\rfloor,k\right)+[e\neq 1]\right) $$ 最后一项 $[e\neq 1]$ 的意思是,对于 $e=1$ 的情况,$S$ 没有计算 $1$ 贡献,刚好,因为此时 $p_k\times 1$ 是质数,其贡献在之前计算过;对于 $e > 1$ 的情况,$p_k^e \times 1$ 是合数,贡献算漏了,要补上。直接暴力递归计算(并且不需要记忆化)。

方法二

也是把贡献拆成质数和合数,只是采用类似于第一步的递推方式:$$ S(n,j)=f(p_{j+1})+S(n,j+1)+\sum_{p_{j+1}\le \sqrt n,1\le e, p_{j+1}^e\le n}f(p_{j+1}^e)\left(S\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,j+1\right)+[e\neq 1]\right) $$ 但直接这样转移复杂度不对,我们要严格控制只有 $\le \sqrt n$ 的质数才发生转移。即,对于 $p_{j+1}>\sqrt n$ 的状态 $S(n,j)$,不显式地计算出来,用到的时候特判为 $g(n)-g(p_j)$。所以,我们要更新状态 $S(n,j)$ 时,一定有 $p_{j+1}\le \sqrt n$,此时需要用到 $S(n,j+1)$,即使要特判为 $g(n)-g(p_{j+1})$,由步骤一中的论述可知 $g(p_{j+1})$ 已经计算出,所以不会有问题。需要注意 $n=2,3$ 的状态不会被任何 $p$ 更新到,所以也要特判。一种较为简洁的方法是,像方法一的代码一样写一个函数 $S(x,y)$,只不过不用来递归,而是用来实现各种特判。

但是,这样代码不太好些,且实现出来常数比较大。我实现了一份:code。耗时几乎是别人的两倍。

方法二改进

既然采用了类似第一步的递推方式,干脆直接套用第一步的状态设计。设 $S(n,j)=\sum_{i=1}^nf(i)[i\text{是质数或其最小质因子}>p_j]$ ,那么参考第一步的方程,有转移:$$ S(n,j)=S(n,j+1)+\sum_{p_{j+1}\le \sqrt n,1\le e, p_{j+1}^e\le n}f(p_{j+1}^e)\left(S\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,j+1\right)-g(\min\{\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor p_{j+1}\})+[e\neq 1]\right) $$ 最后面的 $-g(\dots)$ 是因为,$S\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,j+1\right)$ 中包含了 $\le p_{j+1}$ 的质数的贡献,与转移式的意义不符,应当减去。但是 $g$ 里面有个 $\min$ 非常令人不爽,并且为了卡常,我们可以进一步简化上式:$$ S(n,j)=S(n,j+1)+\sum_{p_{j+1}\le \sqrt n,1\le e, p_{j+1}^{e+1}\le n}f(p_{j+1}^e)\left(S\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,j+1\right)-g(p_{j+1})\right)+f(p_{j+1}^{e+1}) $$ 这是因为当 $p_{j+1}^{e+1}>n$ 的时候 $s(\dots)-g(p_{j+1})$ 这一项不会产生贡献。(其实方法一也可以这样简化式子。)

所以只需要设置 $S(n,+\infty)$ 初值为 $g(n)$,即可进行和第一步几乎一样的 DP 过程了。

此外,注意这种方法不仅计算出了 $S(n)$,也顺带计算出了所有 $S(\lfloor n/x\rfloor)$,这是方法一做不到的。

时间复杂度

对于第一种方法,zzt 大佬在集训队论文中证明,一般情况下其复杂度将会是 $O(n^{1-\epsilon})$,其中 $\epsilon$ 代表一个无穷小量。但是,他也证明了,在 $n \le 10^{13}$ 时,复杂度是 $O\left(\frac{n^\frac{3}{4}}{\log n}\right)$!

对于第二种方法,其复杂度和第一步类似;虽然要枚举质数的指数,但随 $p$ 的增大,指数的枚举次数下降速率急剧趋缓,因而不会对复杂度产生太大影响,视为常数(?) 听说其密度确实仍然是 $O(\frac{1}{\log n})$ 级别。复杂度 $O\left(\frac{n^\frac{3}{4}}{\log n}\right)$。不过第一种方法常数更小,实际运行效可能会高那么一点点。

实现细节和代码

第一步中要使空间为 $O(\sqrt n)$,可以考虑根号分治表示下标。即,对于 $x< \sqrt n$,用 $x$ 映射到下标;对于 $x>\sqrt n$,用 $n/x$ 映射到下标。

P5325【模板】Min_25 筛

第一步把原函数拆成 $g_1(p) = p,g_2(p)=p^2$ 两个完全积性函数计算。

方法一

第二步中直接递归计算,但是用了化简过的递推式。

#include <bits/stdc++.h>

using namespace std;

const int N = 2e5 + 5, mod = 1e9 + 7, i2 = 5e8 + 4, i6 = 166666668;
typedef long long ll;

ll n;
int id1[N], id2[N], m;
int p[N], vis[N], cnt;
ll g1[N], g2[N], v[N];

inline int get(ll x) { return x < N ? id1[x] : id2[n/x]; }
inline ll S1(ll x) { return x %= mod, x * (x + 1) % mod * i2 % mod; }
inline ll S2(ll x) { return x %= mod, x * (x + 1) % mod * (2 * x + 1) % mod * i6 % mod; }
inline ll sq(ll x) { return x %= mod, x * x % mod; }
inline ll F(ll x) { return x %= mod, (sq(x) - x + mod) % mod; }

inline void init(int n) {
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) p[++cnt] = i;
        for (int j = 1; j <= cnt && p[j] <= n / i; j++) {
            vis[i*p[j]] = 1;
            if (i % p[j] == 0) break;
        }
    }
} 

ll S(ll x, int y) {
    if (p[y] >= x) return 0;
    ll res = (g2[get(x)] - g1[get(x)] - g2[get(p[y])] + g1[get(p[y])] + 2 * mod) % mod;
    for (int i = y + 1; i <= cnt && p[i] <= x / p[i]; i++) {
        ll w = p[i];
        for (int j = 1; w <= x / p[i]; j++, w = w * p[i])//注意后面有 x / w, 所以此处不能取模 
            res = (res + F(w) * S(x / w, i) % mod + F(w * p[i])) % mod;
    }
    return res;
}

int main()
{
    scanf("%lld", &n), init(sqrt(n) + 1);
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l), v[++m] = n / l;
        if (v[m] < N) id1[v[m]] = m;
        else id2[n/v[m]] = m;
        g1[m] = (S1(v[m]) - 1 + mod) % mod, g2[m] = (S2(v[m]) - 1 + mod) % mod;
    }
    for (int j = 1; j <= cnt; j++) {
        for (int i = 1; i <= m && p[j] <= v[i] / p[j]; i++) {
            g1[i] = (g1[i] - p[j] * (g1[get(v[i]/p[j])] - g1[get(p[j-1])]) % mod + mod) % mod;
            g2[i] = (g2[i] - sq(p[j]) * (g2[get(v[i]/p[j])] - g2[get(p[j-1])]) % mod + mod) % mod;
        }
    }
    printf("%lld\n", (S(n, 0) + 1) % mod);
    return 0;
}

评测记录

方法二

使用改进后的递推方式。

#include <bits/stdc++.h>

using namespace std;

const int N = 2e5 + 5, mod = 1e9 + 7, i2 = 5e8 + 4, i6 = 166666668;
typedef long long ll;

ll n;
int id1[N], id2[N], m;
int p[N], vis[N], cnt;
ll s[N], g1[N], g2[N], v[N];

inline int get(ll x) { return x < N ? id1[x] : id2[n/x]; }
inline ll S1(ll x) { return x %= mod, x * (x + 1) % mod * i2 % mod; }
inline ll S2(ll x) { return x %= mod, x * (x + 1) % mod * (2 * x + 1) % mod * i6 % mod; }
inline ll sq(ll x) { return x %= mod, x * x % mod; }
inline ll F(ll x) { return x %= mod, (sq(x) - x + mod) % mod; }

inline void init(int n) {
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) p[++cnt] = i;
        for (int j = 1; j <= cnt && p[j] <= n / i; j++) {
            vis[i*p[j]] = 1;
            if (i % p[j] == 0) break;
        }
    }
} 

int main()
{
    scanf("%lld", &n), init(sqrt(n) + 1);
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l), v[++m] = n / l;
        if (v[m] < N) id1[v[m]] = m;
        else id2[n/v[m]] = m;
        g1[m] = (S1(v[m]) - 1 + mod) % mod, g2[m] = (S2(v[m]) - 1 + mod) % mod;
    }
    for (int j = 1; j <= cnt; j++) {
        for (int i = 1; i <= m && p[j] <= v[i] / p[j]; i++) {
            g1[i] = (g1[i] - p[j] * (g1[get(v[i]/p[j])] - g1[get(p[j-1])]) % mod + mod) % mod;
            g2[i] = (g2[i] - sq(p[j]) * (g2[get(v[i]/p[j])] - g2[get(p[j-1])]) % mod + mod) % mod;
        }
    }
    for (int i = 1; i <= m; i++) s[i] = g1[i] = (g2[i] - g1[i] + mod) % mod;
    for (int j = cnt; j >= 1; j--) {
        for (int i = 1; i <= m && p[j] <= v[i] / p[j]; i++) {
            ll w = p[j];
            for (int k = 1; w <= v[i] / p[j]; k++, w *= p[j])
                s[i] = (s[i] + F(w) * (s[get(v[i]/w)] - g1[get(p[j])] + mod) % mod + F(w * p[j])) % mod;
        }
    }
    printf("%lld\n", (s[get(n)] + 1) % mod);
    return 0;
}

评测记录

可以看出方法二实际效率确实劣于方法一。

例题

求质数个数,联想到 min_25 筛的第一步。相当于求 $f(p)=1$ 这个函数在质数处取值的前缀和。

评测记录

容易发现两题本质相同,以第二道为例。

求合法的质数个数,容易联想到 min_25 筛的第一步。可以设 $f(n,j,k)$ 表示前 $n$ 个数中,最小质因子 $>p_j$,且 $\bmod m=k$ 的数的个数。转移的时候,设 $p_j\bmod m=c$,那么枚举 $k\in[0,m-1]$,把 $k$ 的状态贡献到 $kc\bmod m$ 的状态即可。

评测记录

考虑函数质数处的取值。容易发现,$f(p)=\begin{cases}p+1 &p=2\\p-1 &p\neq 2\end{cases}$。

第一步中,应当直接把 $f(p)$ 当成 $p-1$,然后拆成 $g_1(p)=p,g_0(p)=1$ 计算,最后整合起来,并特判含 $2$ 的前缀和。注意,我们应当保证拆成的若干个函数是完全积性函数,因此 $g_0(p)=-1$ 这样的拆分是不可以的。

第二步没有特别之处,套用模板即可。

评测记录

观察发现,当 $\mu(n)=0$ 时,$f(n)=0$。当 $n=1$ 时 $f(n)=1$。接下来考虑 $\mu(n)\neq 0$ 且 $n>1$ 的情况。

根据 $\mu$ 的定义知,此时 $n$ 可以表示成 $p_1p_2\dots p_m$。而 $f(n)$ 取 $+1$ 还是 $-1$ 取决于 $\sum_{d|n}[\mu(d)=-1]$ 的奇偶性。我们发现,$\mu(d)=-1$ 当且仅当 $d$ 中包含奇数个质因子,因此上式等价于 $\sum_{i=2k+1,1\le i \le m} {m\choose i}$。由简单组合知识可知这个式子等于 $2^{m-1}$。于是,当 $m=1$ 即 $n$ 为质数时,$f(n)=-1$,否则 $f(n)=1$。

我们惊奇地发现 $f(n)=\begin{cases}\mu^2(n) & n\text{为合数}\\ -1 &n\text{为质数}\end{cases}$。所以要求 $\sum_{i=1}^nf(i)$,就是求 $\sum_{i=1}^n\mu^2(i)-2\cdot\sum_{i=1}^n[i\text{是质数}]$。

前一项是积性函数前缀和,可以 min_25 筛。后一项显然套用 min_25 筛的第一步即可。

评测记录

容易知道就是求 $\sum_{i=1}^n{d(i)\choose 2}$,其中 $d(n)$ 表示 $n$ 的因数个数。展开变成 $\frac{1}{2}\left(\sum_{i=1}^nd^2(i)-\sum_{i=1}^nd(i)\right)$。

对于第一项,$d^2(n)$ 是积性函数,且 $d^2(p)=4,d^2(p^k)=(k+1)^2$,考虑 min_25 筛。第一步要求函数是完全积性函数,所以提一个 $4$ 出来,$d'(p)=1$,最后再乘回去。

对于第二项,枚举因数:$\sum_{d=1}^n\lfloor\frac{n}{d}\rfloor$。数论分块即可。

评测记录

就是说要知道所有 $S(\lfloor\frac{n}{x}\rfloor)$ 的值。刚好 min_25 筛的第二种写法可以求出这些东西。

评测记录

不过 $10^{13}$ 的版本似乎不是为 min_25 筛准备的(或者我的写法常数太大?),反正我觉得很卡长,尽管有大佬的代码跑得飞快。最后卡了一晚上+一早上,终于离 TL $500ms$ 左右通过了,但代码几乎惨不忍睹(。评测记录

有没有大佬教教我正解的啊(悲)

小结

好像重要的点在简介里都提到了。

参考资料