【学习笔记】Min_25 筛
神仙 Min_25 发明的神奇筛子，用来筛积性函数前缀和。
对于要筛前缀和的积性函数 \(f\left(x\right)\)，其具体要求是：\(f\left(p\right)\) 是一个简单多项式，\(f\left(p^e\right)\) 可以快速计算。
复杂度是亚线性的，但我不会证。
Description
对于一个积性函数 \(f\left(x\right)\) 的前缀和，我们可以做如下变换：
第一行就是把 \(1\sim n\) 的数的贡献分为素数和合数分别进行计算。
第二行就是枚举每个合数的最小质因子，再计算剩下的，其中 \(\operatorname{lpf}\left(i\right)\) 代表 \(i\) 的最小质因子。
然后就可以分为两个部分来计算了！
Part 1 前半部分
已知 \(f\left(p\right)\) 是关于 \(p\) 的简单多项式，那么我们可以把这个多项式拆成若干个单项式来进行计算。也就是说，我们只需求出 \(\sum_{p\le n}p^k\) 即可，其中 \(p^k\) 是原多项式的一部分。
考虑进行一个神奇的 dp，定义 dp 数组 \(g\left(n,j\right)\) 为：
如果把算法过程想象成埃氏筛的话，\(g\left(n,j\right)\) 计算的其实就是筛去 \(p_j\) 及倍数后剩余的所有“素数”对答案的贡献。
考虑从 \(g\left(n,j-1\right)\) 转移过来，也就是说我们要计算筛掉 \(p_j\) 对答案的贡献。把 \(p_j\) 分成两类讨论：
- \(p_j>\sqrt n,p_j^2>n\)：可知最小质因数为 \(p_j\) 的最小合数为 \(p_j^2\)，因此筛去 \(p_j\) 对答案没有贡献，可得 \(g\left(n,j\right)=g\left(n,j-1\right)\)。
- \(p_j\le\sqrt n,p_j^2\le n\)：枚举除去最小质因子 \(p_j\) 后被筛掉的合数的值，可知筛去 \(p_j\) 对答案的贡献为 \(p_j^kg\left(\left\lfloor\frac n{p_j}\right\rfloor,j-1\right)\)。但是这样我们会多筛一次 \(p_1\sim p_{j-1}\)，因此答案要再加上 \(p_j^kg\left(p_j-1,j-1\right)\)。
于是我们得到了 \(g\left(n,j\right)\) 的递推式：
其中，\(g\left(p_j-1,j-1\right)\) 其实就是 \(\sum_{i=1}^{j-1}p_i^k\)，可以线性筛出。记 \(\sum_{i=1}^jp_i^k\) 为 \(\textit{Sp}_j\)。
注意到当 \(p_j>\sqrt n\) 时，\(g\left(n,j\right)\) 的值都是相同的，于是这样的 \(g\left(n,j\right)\) 的值可以不用筛。同时，对于每个 \(x\)，我们需要的只有所有 \(g\left(\left\lfloor\frac nd\right\rfloor,x\right)\) 的值，而这样的值只有 \(O\left(\sqrt n\right)\) 个，存储时映射一下下标即可。
例题：P5493 质数前缀统计
设 \(S\left(n\right)\) 为 \(n\) 以内质数的 \(k\) 次方和。
给定 \(N\)，求：\[\sum_{i=1}^{\left\lfloor\sqrt N\right\rfloor}i^2S\left(\left\lfloor\frac Ni\right\rfloor\right)\bmod p \]\(0 \le k \le 10\)，\(1 \le N \le 4\times {10}^{10}\)，\({10}^9 < p < 1.01 \times {10}^9\) 。
板子题，按照上面的过程实现即可。计算自然数 \(k\) 次方和可参考 CF622F。
code（用了 FastMod 优化取模）：
#include <cstdio>
#include <cmath>
const int maxn = 1e6 + 5;
typedef long long ll;
namespace FastModnmsp {
typedef unsigned int uint;
typedef __uint128_t uL;
typedef unsigned long long ull;
struct FastMod {
ull b,m;
FastMod() {}
FastMod(ull _b): b(_b),m(ull((uL(1) << 64) / _b)) {}
inline void init(ull _b) { b = _b,m = ull((uL(1) << 64) / _b); }
friend inline ull operator % (const ull& a,const FastMod& mod) {
ull r = a - (uL(mod.m) * a >> 64) * mod.b;
return r >= mod.b ? r - mod.b : r;
}
} Mod;
} using FastModnmsp::Mod;
ll n,sq,k,ans,P,finv[15],prsk[15];
ll cntP,pri[maxn],prik[maxn],sumk[maxn];
ll cntS,w[maxn],g[maxn];
int id1[maxn],id2[maxn];
bool isp[maxn];
inline ll read() {
#define gc c = getchar()
ll d = 0; int f = 0,gc;
for(;c < 48 || c > 57;gc) f |= (c == '-');
for(;c > 47 && c < 58;gc) d = (d << 1) + (d << 3) + (c ^ 48);
#undef gc
return f ? -d : d;
}
inline ll fpow(ll a,ll b) {
ll res = 1;
for(;b;a = a * a % Mod,b >>= 1) if(b & 1) res = res * a % Mod;
return res;
}
inline ll getPowSum(ll n) {
ll ans = 0,p = 1,pre = 0;
for(int i = 1;i <= k + 2;i ++) {
p = p * (n - i) % Mod,pre = (pre + prsk[i]) % Mod;
if(n == i) return pre;
ll inv = finv[i - 1] * finv[k + 2 - i] % Mod * fpow(n - i,P - 2) % Mod;
ll res = pre * inv % Mod;
if((k - i) & 1) res = P - res;
ans = (ans + res) % Mod;
}
return p * ans % Mod;
}
inline int getId(ll k) { return k <= sq ? id1[k] : id2[n / k]; }
inline void init() {
finv[0] = isp[1] = 1; sq = sqrt(n);
for(ll i = 1;i <= k + 2;i ++) finv[i] = finv[i - 1] * fpow(i,P - 2) % Mod,prsk[i] = fpow(i,k);
for(ll i = 2;i <= sq;i ++) {
if(!isp[i]) pri[++ cntP] = i,
prik[cntP] = fpow(i,k),
sumk[cntP] = (sumk[cntP - 1] + prik[cntP]) % Mod;
for(int j = 1;j <= cntP && i * pri[j] <= sq;j ++) {
isp[i * pri[j]] = true;
if(!(i % pri[j])) break;
}
}
for(ll r,l = 1;l <= n;l = r + 1) {
r = n / (n / l); w[++ cntS] = n / l;
g[cntS] = getPowSum(w[cntS] % Mod) - 1;
if(n / l <= sq) id1[n / l] = cntS;
else id2[r] = cntS;
}
for(int i = 1;i <= cntP;i ++)
for(int j = 1;j <= cntS && pri[i] * pri[i] <= w[j];j ++) {
ll h = getId(w[j] / pri[i]);
g[j] -= prik[i] * (g[h] - sumk[i - 1] + P) % Mod;
if(g[j] < 0) g[j] += P;
}
}
int main() {
n = read(),k = read(),P = read(); Mod.init(P);
init();
for(ll i = 1;i <= sq;i ++) ans = (ans + i * i % Mod * g[i] % Mod) % Mod;
printf("%lld\n",ans);
return 0;
}
Part 2 求解答案
还是 dp。考虑借鉴前半部分，定义 dp 数组 \(S\left(n,j\right)\)：
答案显然就是 \(S\left(n,0\right)+1\)。这里的 \(1\) 是 \(f\left(1\right)\) 的贡献，把它提出来计算。
接下来分两部分讨论贡献：
- 素数：所有素数的贡献为 \(g\left(n,j\right)\)，但是要减掉 \(p_1\sim p_{j}\) 的贡献，因此答案是 \(g\left(n,j\right)-\textit{Sp}_j\)。
- 合数：枚举最小质因数及其次数，可得贡献为 \(\sum_{\substack{k>j\\p_k^e\le n}}f\left(p_k^e\right)\left(S\left(\left\lfloor\frac n{p_k^e}\right\rfloor,k\right)+[e\not=1]\right)\)，其中 \([e\not=1]\) 计算的是 \(p_k^e\) 的贡献。
所以最终的式子就是：
其中 \(x\) 为 \(n\) 以内使得 \(p_x^2\le n\) 的最大的 \(x\)，即 \(x=\pi\left(\left\lfloor\sqrt n\right\rfloor\right)\)（因为 \(x\) 再大的话 \(g(n,x)\) 的值不会变化）。
很神奇的是，这玩意直接爆搜就行了，根本不需要记忆化，而且复杂度还是优秀的 \(O\left(\frac{n^{0.75}}{\log n}\right)\)。
例题：
令积性函数 \(f\left(p^k\right)=p^k\left(p^k-1\right)\)。给定 \(n\)，求：
\[\sum_{i=1}^nf\left(i\right)\bmod \left(10^9+7\right) \]\(1\le n\le 10^{10}\)。
在质数处，显然有 \(f\left(p\right)=p\left(p-1\right)=p^2-p\)，可以 Min_25 筛。
code：
#include <cstdio>
#include <cmath>
typedef long long ll;
const int maxn = 1e6 + 5;
const ll P = 1e9 + 7,inv3 = (P + 1) / 3;
ll n,sq,tot,cnt;
ll pri[maxn]; bool isp[maxn];
ll Sp1[maxn],Sp2[maxn];
ll g1[maxn],g2[maxn],w[maxn];
int id1[maxn],id2[maxn];
inline ll getLiSum(ll n) { return n * (n + 1) / 2 % P; }
inline ll getSqSum(ll n) { return n * (n + 1) / 2 % P * (2 * n + 1) % P * inv3 % P; }
inline ll Add(ll a,ll b) { a += b; return a >= P ? a - P : a; }
inline int getId(ll k) { return k <= sq ? id1[k] : id2[n / k]; }
inline void Init() {
isp[1] = true,sq = sqrt(n);
for(ll i = 2;i <= sq;i ++) {
if(!isp[i]) pri[++ cnt] = i,
Sp1[cnt] = Add(Sp1[cnt - 1],i),
Sp2[cnt] = Add(Sp2[cnt - 1],i * i % P);
for(int j = 1;j <= cnt && i * pri[j] <= sq;j ++) {
isp[i * pri[j]] = true;
if(!(i % pri[j])) break;
}
}
for(ll r,l = 1;l <= n;l = r + 1) {
r = n / (n / l);
w[++ tot] = n / l;
g1[tot] = getLiSum(w[tot] % P) - 1;
g2[tot] = getSqSum(w[tot] % P) - 1;
if(n / l <= sq) id1[n / l] = tot;
else id2[r] = tot;
}
for(int i = 1;i <= cnt;i ++)
for(int j = 1;j <= tot && pri[i] <= w[j] / pri[i];j ++) {
int k = getId(w[j] / pri[i]);
g1[j] -= pri[i] * (g1[k] - Sp1[i - 1] + P) % P;
g2[j] -= pri[i] * pri[i] % P * (g2[k] - Sp2[i - 1] + P) % P;
if(g1[j] <= 0) g1[j] += P; if(g2[j] <= 0) g2[j] += P;
}
}
inline ll S(ll n,int j) {
if(pri[j] >= n) return 0;
ll x = getId(n),ans = Add(((g2[x] - Sp2[j]) - (g1[x] - Sp1[j])) % P,P);
for(int k = j + 1;k <= cnt && pri[k] <= n / pri[k];k ++) {
ll p = pri[k];
for(int e = 1;p <= n;e ++,p *= pri[k]) {
ll tmp = p % P;
ans = Add(ans,tmp * (tmp - 1) % P * (S(n / p,k) + (e != 1)) % P);
}
}
return ans;
}
int main() {
scanf("%lld",&n); Init();
printf("%lld\n",Add(S(n,0),1));
return 0;
}
令 \(\sigma_0\left(n\right)\) 为 \(n\) 的约数个数。给定 \(n,k\)，求：
\[\sum_{i=1}^n\sigma_0\left(i^k\right)\bmod 2^{64} \]多测，\(1\le n,k\le10^{10}\)。
三倍经验题，写完这题就可以切了 DIVCNT2 和 DIVCNT3。
令 \(f\left(n\right)=\sigma_0\left(n^k\right)\)，则显然有：
可以 Min_25 筛出来。
code：
#include <cstring>
#include <cstdio>
#include <cmath>
typedef unsigned long long ull;
const int maxn = 1e6 + 5;
ull n,K,sq;
ull pri[maxn]; bool isp[maxn];
ull Sp[maxn],g[maxn],w[maxn];
int t,tot,cnt,id1[maxn],id2[maxn];
inline int getId(ull k) { return k <= sq ? id1[k] : id2[n / k]; }
inline void Init() {
isp[1] = true,sq = sqrt(n);
for(ull i = 2;i <= sq;i ++) {
if(!isp[i]) pri[++ cnt] = i,
Sp[cnt] = Sp[cnt - 1] + K + 1;
for(int j = 1;j <= cnt && i * pri[j] <= sq;j ++) {
isp[i * pri[j]] = true;
if(!(i % pri[j])) break;
}
}
for(ull r,l = 1;l <= n;l = r + 1) {
r = n / (n / l);
w[++ tot] = n / l;
g[tot] = w[tot] - 1;
if(n / l <= sq) id1[n / l] = tot;
else id2[r] = tot;
}
for(int i = 1;i <= cnt;i ++)
for(int j = 1;j <= tot && pri[i] <= w[j] / pri[i];j ++) {
int k = getId(w[j] / pri[i]);
g[j] -= g[k] - i + 1;
}
for(int i = 1;i <= tot;i ++) g[i] *= (K + 1);
}
inline ull S(ull n,int j) {
if(pri[j] >= n) return 0;
ull x = getId(n),ans = g[x] - Sp[j];
for(int k = j + 1;k <= cnt && pri[k] <= n / pri[k];k ++) {
ull p = pri[k];
for(ull e = 1;p <= n;e ++,p *= pri[k])
ans += (K * e + 1) * (S(n / p,k) + (e != 1));
}
return ans;
}
int main() {
scanf("%d",&t);
for(;t;t --) {
memset(isp,false,sizeof(isp)); tot = cnt = 0;
scanf("%llu%llu",&n,&K); Init();
printf("%llu\n",S(n,0) + 1);
}
return 0;
}
