min_25筛小记
引言
现在有一个积性函数 \(f(n)\) ,需要求出 \(f\) 的前缀和,即给定一个数 \(n\) ,求出:
显然,当 \(n\) 较小的时候,即可以用数组存下所有 \(f\) 的值时,可以使用线性筛,但是当 \(n\) 达到 \(1\times 10^{10}\) 时,线性筛显然不能胜任,需要使用 \(\texttt{min\_25}\) 筛 。
符号约定
- \(\mathbb{P}\) 表示全体素数,\(\mathbb{P}_i\) 表示第 \(i\) 小的素数,比如 \(\mathbb{P}_1=2, \mathbb{P}_3=5\) 。
- \(\mathrm{lpf}(n) = \min_{\mathbb{P}_i | n}\{\mathbb{P}_i\}\) ,表示 \(n\) 的最小素因子,特别的,\(n=1\) 的时候是 \(0\) 。
基本算法
算法分析(瞎吹)
这里的 \(f(i)\) 满足 \(f\left(\mathbb{P}_i^k\right) = \left(\mathbb{P_i^k}\right)^2 - \mathbb{P}_i^k\) 。
首先考虑将前缀和分成素数和合数两个部分,即:
注意,因为所有合数 \(a (a\not\in \mathbb{P})\)的最小质因子都是小于 \(\sqrt{a}\) 的,因此枚举最小质因子的时候只需要到 \(\sqrt{n}\) 就行了,这里用线性筛先将 \(\sqrt{n}\) 以内的所有素数筛出来,假设最大的素数是 \(\mathbb{P}_m\) 。
接着考虑对于 \(f(i)\) ,我们需要找到一个完全积性函数对其进行拟合,使其在 \(\mathbb{P}_i^k\) 处的值相等。通常我们可以找到一个低次多项式进行拟合,假设次数为 \(k\) ,次数为 \(k\) 的项的系数为 \(a_k\) ,比如对于这一题,选择函数 \(i^2-i\) 即可,此时 \(k=2, a_0=0,a_1=-1,a_2=1\) 。
设 \(g_k(n,x), h_k(n)\) ,这两个函数的定义为:
这里用递推求出 \(g_k(n,x)\) 。首先容易得到:
这里减一是为了将 \(1\) 去掉,接着转移,容易想到从 \(g_k(*,x-1)\) 转移到 \(g_k(*,x)\) ,这时候需要减去 \(\mathbb{P}_x\) 作为最小质因子的数。
注意,这里 \(g_k(\mathbb{P}_{x-1}, x-1)\) 是为了避免减去质数对 \(g_k(n,x)\) 的贡献,这里还可以用 \(h_k(x-1)\) 代替。
接下来考虑计算最终答案。有一个类似 \(\rm dp\) 的东西,设 \(S(n,x)\) 表示 \(1\) 到 \(n\) 中,最小质因子大于 \(\mathbb{P}_x\) 的所有 \(i\) 的 \(f(i)\) 的和,可以想到类似上面式子的转移:
特别的,令 \(S(n,x)=0 (\mathbb{P}_x>n\lor n=1)\) 。\(w(n,x)\) 表示 \(n\) 以内比 \(\mathbb{P}_x\) 大的素数的 \(f\) 值的和。在这道题中,有:
这里减去了小于等于 \(\mathbb{P}_x\) 的素数的贡献,提醒一下,这里的 \(m\) 指的是之前筛出的小于等于 \(\sqrt{n}\) 的最大素数的编号。
显然最终答案是 \(\mathrm{ANS} = S(n,0)\) 。因为玄学原因(我也不知道为什么),这个 \(S(n,0)\) 可以直接用递归计算,而且不用记忆化。
接下来考虑具体实现,关键在 \(g\) 的存储,可以发现形式都是 \(\left\lfloor\frac{n}{a}\right\rfloor\) ,同时,因此可以使用整除分块,同时因为小于 \(\sqrt{n}\) ,大于 \(\sqrt{n}\) 的位置的各自不超过 \(\sqrt{n}\) 个,因此可以用两个数组记一下下标。
整除性质 :
\[\left\lfloor\frac{\left\lfloor\frac{n}{a}\right\rfloor}{b}\right\rfloor=\left\lfloor\frac{n}{ab}\right\rfloor \]这个性质保证连续的整除,都可以表示整除另一个数。
同时,因为从 \(g_k(*,x-1)\) 转移到 \(g_k(*,x)\) ,因此可以用滚动数组,类似 \(\rm 01\) 背包的倒序枚举之类的进行优化,反正最后只需要使用 \(g_k(*, m)\) 。
代码
typedef Modint<1000000007> mint;
int m, pri[maxn], pri_cnt, ky_cnt, vis[maxn], id1[maxn], id2[maxn];
LL n, ky[maxn];
mint g[maxn][2], h[maxn][2];
inline int getid(LL x) { return x <= m ? id1[x] : id2[n / x]; };
void init() {
m = sqrt(n)/*抱歉,这里的m和讲解的m不是一个东西,讲解的m=pri_cnt*/;
mint inv2(500000004), inv6(166666668);
for (int x = 2; x <= m; x++) {
if (!vis[x]) {
pri[++pri_cnt] = x;
h[pri_cnt][0] = h[pri_cnt - 1][0] + mint(x);
h[pri_cnt][1] = h[pri_cnt - 1][1] + mint(x) * mint(x);
}
for (int j = 1; j <= pri_cnt && 1ll * x * pri[j] <= m; j++) {
vis[x * pri[j]] = 1;
if (x % pri[j] == 0) break;
}
}
//计算 g_k(*,0)
for (LL i = 1, j; i <= n; i = j + 1) {
j = n / (n / i), ky[++ky_cnt] = n / j;
mint t(ky[ky_cnt]);
g[ky_cnt][0] = t * (t + 1) * inv2 - 1;
g[ky_cnt][1] = t * (t + 1) * (mint(2) * t + 1) * inv6 - 1;
if (n / j <= m) id1[n / j] = ky_cnt; else id2[j] = ky_cnt;
}
//递推计算 g
for (int i = 1; i <= pri_cnt; i++) {
LL lim = 1ll * pri[i] * pri[i];
mint prii(pri[i]);
for (int j = 1; j <= ky_cnt && ky[j] >= lim; j++) {
int kyi = getid(ky[j] / pri[i]);
g[j][0] = g[j][0] - prii * (g[kyi][0] - h[i - 1][0]);
g[j][1] = g[j][1] - prii * prii * (g[kyi][1] - h[i - 1][1]);
}
}
}
mint S(LL x, LL y) {
if (pri[y] >= x) return mint(0);
LL id = getid(x);
mint ans = (g[id][1] - g[id][0]) - (h[y][1] - h[y][0]);
for (int i = y + 1; i <= pri_cnt && 1ll * pri[i] * pri[i] <= x; i++) {
for (LL c = 1, t = pri[i]; t <= x; c++, t *= pri[i]) {
mint mt(t);
ans = ans + mt * (mt - 1) * (S(x / t, i) + mint(c != 1));
}
}
return ans;
}
例题
求区间素数个数
就是 \(g_0(n,m)\) 。(草率)