Min_25 筛小结
Min_25 筛这个东西,完全理解花了我很长的时间,所以写点东西来记录一些自己的理解。
它能做什么
对于某个数论函数 \(f\),如果满足以下几个条件,那么它就可以用 Min_25 筛来快速求出这个函数的前缀和。
-
它是一个积性函数
-
对于一个质数 \(p\) ,\(f(p)\) 的表达式必须是一个项数比较小的多项式。即 \(\displaystyle f(p) = \sum a_ip^{b_i}\)。
-
对于一个质数 \(p\) ,\(f(p^k)\) 的表达式必须可以由 \(f(p)\) 快速得到。
看起来这些条件比较紧,其实仔细想想许多函数还是满足这些性质的。接下来讨论它的实现过程。
实现过程
Min_25 筛的核心思想就是对于 \(\le n\) 的数中,除掉质数之后,每个数的最小质因子大小不超过 \(\sqrt{n}\) 。即类似埃氏筛的方式,从每个质因子筛掉一部分的贡献。整个过程可以分成三部分,即质数的部分的贡献之和与合数的部分的贡献和,以及 \(1\) 的特判。
记号约定
这里用 \(P\) 表示质数集合,\(P_i\) 表示第 \(i\) 小的质数。
质数部分
首先我们不妨设 \(f(p) = p^k\) 。对于更加一般的项数小的多项式形式,我们可以分别记每一个项的贡献,然后分别相加。注意到这里 \(f(p)\) 是一个完全积性函数。对于 \(p\) 为合数的情况,为了方便,我们先姑且同样定义 \(f(p) = p^k\) ,到最后你会明白这样定义的原因。
我们再设 \(g(n,i)\) 表示 \(\le n\) 的所有数 \(x\) 中,\(x\) 的最小质因子 \(> P_i\) 或者 \(x\) 本身为质数的所有 \(f(x)\) 之和,这样 \(g(n,|P|)\) 就是所有质数的贡献之和了。而 \(g(n,0)\) 的答案也十分显然,即 \(\displaystyle \sum_{i = 2} ^{n} i^k\) 。我们只需要考虑,如何从 \(g(*,i - 1)\) 如何得到 \(g(*,i)\) 。
考虑我们从 \(i - 1\) 推到 \(i\) 的过程中,会减少哪个部分的贡献,容易发现最小质因子恰好为 \(P_i\) 的贡献被减掉了。如果 \(P_i^2 > n\) ,那么这一部分是没有贡献的,即 \(g(n,i) = g(n,i - 1)\)。
对于 \(P_i^2 \le n\) 的情况,我们需要把所有 \(P_i\) 的倍数,且满足其最小质因子除去 \(P_i\) 之外 \(> P_{i-1}\) 的贡献全部减掉。如果直接减去 \(f(P_i) \cdot g(\lfloor\frac{n}{P_i}\rfloor,i - 1)\) 你会发现贡献是减多了的,原因就是 \(g\) 函数的值还包含了所有质数的贡献之和,仔细分析一下,可以发现,对于 \(x = P_iP_j(j < i)\) 的贡献在 \(j\) 这个地方已经减过了,所以还要加回来,即 \(\displaystyle f(P_i) \cdot g(P_i - 1,i - 1) = \sum_{j = 1}^{i - 1} f(P_i)f(P_j)\) 。
式子汇总起来也就是
回头来,我们再来看最初我们提出的一些问题。
我们最开始为什么要定义,对于一个合数 \(p\) ,\(f(p) = p^k\) 呢?无非就是方便计算,方便利用完全积性函数的性质。你发现到最后,合数的贡献还是没有算进去,所以对于合数的定义是无所谓的。
我们到底在哪里利用了 \(f(p)\) 是一个完全积性函数的性质呢?这个性质是在上式中第二种情况中利用到的。注意到我们在减贡献的时候,并没有枚举 \(P_i\) 的次数,这也就意味着 \(g(\lfloor\frac{n}{P_i}\rfloor,i - 1)\) 中可能包含了某些数仍然是 \(P_i\) 的倍数。所以这一步需要有一个对于所有质数 \(p,f(p^2) = f^2(p)\) 的条件。由于之前我们假定了 \(f(p)\) 对于 \(p\) 为合数的情况也满足 \(f(p) = p^k\) ,因此这个条件是得到满足了的。
递归实现很显然。考虑一下非递归的实现。首先对于这个步骤,我们有用的 \(n\) 只有所有不同的 \(\lfloor\frac{n}{i}\rfloor,(1 \le i \le n)\),所以可以离散掉这些值,用一个大小为 \(2\sqrt{n}\) 的数组存下来。由于这个递推是一层一层推的,所以实现上可以不用开二维数组,可以直接离散之后用一个一维数组存下来,需要注意枚举顺序。
合数部分
考虑我们在求质数部分的时候,除了质数的答案之和之外,还求出了些什么东西。事实上,对于所有的 \(\lfloor\frac{n}{i}\rfloor,(1 \le i \le n)\) ,我们求出了 \(g(\lfloor\frac{n}{i}\rfloor,|P|)\) 的值。
我们设 \(S(n,i)\) 表示 \(\le n\) 的所有数 \(x\) 中,\(x\) 的最小质因子 \(\ge P_i\) 的 \(f(x)\) 之和。注意这里的 \(f(x)\) 在 \(x\) 为合数的时候不再沿用质数部分时的定义,也注意这里的定义与 \(g\) 的区别。
接下来我们先算上所有满足条件的质数贡献之和,即 \(\displaystyle g(n,|P|) - \sum_{j = 1}^{i - 1} f(P_j)\) 。对于合数,我们直接枚举其最小质因子以及质因子的个数,直接递归计算。注意在这里形如 \(p^k,(p \in P)\) 的贡献并没有算进去,所以还要单独加一下。式子的形式如下:
直接对着式子看,应该还是很好理解的。对这个式子,直接递归暴力计算即可,跑得会比第一部分快。
至此 Min_25 筛的全过程就讲完了。时间复杂度是 \(\displaystyle O(\frac{n^{\frac{3}{4}}}{\log n})\) 的,不大会证……
例题
首先可以用一个小问题来感受 Min_25 筛的威力。
SPOJ DIVCNT2
题意
求
\[\sum_{i = 1}^{n} \sigma_0(i^2) \]其中 \(\sigma_0(i)\) 表示 \(i\) 的约数个数。
\(n \le 10^{12},TimeLimit = 20s\) 。
题解
传统的杜教筛做法需要考虑将式子拆为每个质因子贡献相乘,然后反演 + 杜教筛。如果我们有 Min_25 筛呢?
我们发现 \(f(x) = \sigma_0(x^2)\) 满足 Min_25 筛的三个基本条件,\(f(p) = 3,f(p^k) = 2k + 1\) ,因此就是一道裸题了。
我们发现 SPOJ 还有两道题叫 DIVCNT3 / DIVCNTK。于是此题三倍经验……
代码
#include<bits/stdc++.h>
#define For(i,l,r) for(int i = (l),i##end = (r);i <= i##end;i++)
#define Fordown(i,r,l) for(int i = (r),i##end = (l);i >= i##end;i--)
#define debug(x) cout << #x << " = " << x << endl
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
template <typename T> inline bool chkmin(T &x,T y) { return y < x ? x = y,1 : 0; }
template <typename T> inline bool chkmax(T &x,T y) { return x < y ? x = y,1 : 0; }
const int INF = 0x3f3f3f3f;
const int N = 2e6 + 10;
int prime[N],vis[N],id1[N],id2[N];
ull B[N],Ans[N];
ll val[N],n;
int pcnt = 0;
inline ll read() {
ll x = 0,flag = 1;
char ch = getchar();
while(!isdigit(ch) && ch != '-') ch = getchar();
if(ch == '-') flag = -1,ch = getchar();
while(isdigit(ch)) x = (x << 3) + (x << 1) + (ch - '0'),ch = getchar();
return x * flag;
}
inline void init(int Lim) {
For(i,2,Lim) {
if(!vis[i]) prime[++pcnt] = i;
for(int j = 1;j <= pcnt && i * prime[j] <= Lim;j++) {
vis[i * prime[j]] = true;
if(i % prime[j] == 0) break;
}
}
}
#define id(x) (x <= d ? id1[x] : id2[n / (x)])
inline ull S(ll Lim,int cnt) {
if(Lim <= 1 || prime[cnt] > Lim) return 0;
int d = sqrt(n);
ull ans = 3 * (B[id(Lim)] - (cnt - 1));
for(int i = cnt;i <= pcnt && 1ll * prime[i] * prime[i] <= Lim;i++) {
ll prod = prime[i];
for(int times = 1;1ll * prod * prime[i] <= Lim;times++,prod = 1ll * prod * prime[i])
ans += (2 * times + 1) * S(Lim / prod,i + 1) + (2 * times + 3);
}
return ans;
}
inline void Calc(ll Lim) {
int d = sqrt(Lim),cnt = 0;
for(ll i = 1;i <= Lim;i = Lim / (Lim / i) + 1) {
val[++cnt] = Lim / i,id(val[cnt]) = cnt;
B[cnt] = val[cnt] - 1;
}
for(int i = 1;i <= pcnt && 1ll * prime[i] * prime[i] <= Lim;i++)
for(int j = 1;j <= cnt && 1ll * prime[i] * prime[i] <= val[j];j++)
B[j] = B[j] - B[id(val[j] / prime[i])] + (i - 1);
}
int main() {
#ifndef ONLINE_JUDGE
freopen("DIVCNT2.in","r",stdin);
freopen("DIVCNT2.out","w",stdout);
#endif
int Case = read();
init(2e6 + 5);
while(Case--) {
n = read();
Calc(n);
ull res = S(n,1) + 1;
printf("%llu\n",res);
}
return 0;
}
Min_25 筛在一些质因子相关的问题上,也有很强的适用性,比如下面的这些题。
UOJ188 Sanrd
题意
设 \(f(x)\) 为 \(x\) 的次大质因子,重复质因子算多个。求
\[\sum_{i = l} ^ {r} f(i) \]\(l,r \le 10^{11},TimeLimt = 5s\) 。
题解
这个函数显然不是积性函数了,但是这个次大质因子的形式,和 Min_25 筛枚举最小质因子的形式比较契合,因此可以想到在 Min_25 筛过程中统计答案。
考虑我们在递归计算 \(S\) 的过程中,假设当前递归到了 \(S(a,b)\) ,那么这些 \(\le a\) 且 \(> P_b\) 的质数,在最开始的时候,一定是 \(P_{b - 1}\) 的倍数(具体原因可以仔细看看 \(S\) 的递归计算的式子)。所以这些数的次大质因子都是 \(P_{b - 1}\) 。因此就可以在 Min_25 筛过程中直接计算了。
LOJ572 Misaka Network 与求和
题意
设 \(f(x)\) 为 \(x\) 的次大质因子,重复质因子算多个。求
\[\sum_{i = 1}^{n}\sum_{j = 1}^{n}f(\gcd(i,j))^k \bmod 2^{32} \]\(n,k \le 2 \times 10^9,TimeLimit = 4s\) 。
题解
前面经典的反演比较简单,式子可以转化为
设 \(F(x) = f(x)^k\) ,则所求即为
于是我们只要能够快速求出 \(F * \mu\) 的前缀和,就可以整除分块了。我们设 \(\displaystyle S(n) = \sum_{i = 1}^{n} (F * \mu)(n)\),可以写出一个杜教筛的式子。
由于 \(\mu * 1 = \epsilon\) ,容易想到令 \(g = 1\) 。那么 \((F * \mu * g)(i) = (F * (\mu * g))(i) = (F * \epsilon)(i) = F(i)\)
对于 \(F(i)\) 的前缀和,直接套用 Min_25 筛即可。