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(lfloorfrac{n}{P_i} floor,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(lfloorfrac{n}{P_i} floor,i - 1)) 中可能包含了某些数仍然是 (P_i) 的倍数。所以这一步需要有一个对于所有质数 (p,f(p^2) = f^2(p)) 的条件。由于之前我们假定了 (f(p)) 对于 (p) 为合数的情况也满足 (f(p) = p^k) ,因此这个条件是得到满足了的。
递归实现很显然。考虑一下非递归的实现。首先对于这个步骤,我们有用的 (n) 只有所有不同的 (lfloorfrac{n}{i} floor,(1 le i le n)),所以可以离散掉这些值,用一个大小为 (2sqrt{n}) 的数组存下来。由于这个递推是一层一层推的,所以实现上可以不用开二维数组,可以直接离散之后用一个一维数组存下来,需要注意枚举顺序。
合数部分
考虑我们在求质数部分的时候,除了质数的答案之和之外,还求出了些什么东西。事实上,对于所有的 (lfloorfrac{n}{i} floor,(1 le i le n)) ,我们求出了 (g(lfloorfrac{n}{i} floor,|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
",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 mod 2^{32} ](n,k le 2 imes 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 筛即可。