min_25 筛学习
算法
min_25 筛是解决如下问题的：
设 \(f\) 为一个积性的数论函数，问求 \(\sum_{i=1}^n f(i)\) 。其中 \(f\) 满足若 \(i\) 为质数那么 \(f(i^k)\) 可以快速计算。
min_25 筛算法可以在 \(O\left(\frac{n^{\frac34}}{\log n}\right)\) （通常情况下）的时间复杂度内解决上述问题。
下面令 \(P\) 为质数集合。
我们记录 \(g\) 也是一个积性函数，且若 \(i\) 是质数，那么 \(f(i)=g(i)\) ，且 \(g(i)\) 需要满足如下两个性质：
- \(s(i)=\sum_{1\le j\le i}g(p_i)\) 可以快速计算。
- 若已知 \(\sum_{i\in S}g(i)\) ，那么对质数 \(j\)，\(T(\sum_{i\in S}g(i),j)=\sum_{i\in S}g(i\times j)\) 可以快速计算。
由于下式中在 \(2\) 到 \(n\) 中的合数的最小质因子一定 \(\le \sqrt n\)，于是设 \(p\) 是一个所有 \(\le \sqrt n\) 的质数的序列，长度为 \(m\) 。
令 \(G(i,j)\) 表示 \(2\) 到 \(i\) 中所有质数和所有最小质因子 \(>p_j\) 的合数的 \(g\) 的函数值的和。
令 \(s(i)=\sum_{1\le j\le i}g(p_i)\)。
边界条件为 \(G(i,0)=\sum_{k=2}^i g(i)\)。
对于 \(i>p_j^2\) 转移方程为 \(G(i,j)=G(i,j-1)-T\left(G\left(\left\lfloor\frac i{p_j}\right\rfloor,j-1\right),p_j\right)+s(j-1)\)。否则 \(G(i,j)=G(i,j-1)\)。
所有质数的 \(f\) 的和就是 \(G(n,m)\)，这里可以滚动数组优化。 \(i\) 的取值一定是 \(\left\lfloor\frac nd\right\rfloor\) 种类数量为是 \(O(\sqrt n)\) 的。可以 \(j\) 从小到大滚动数组。
令 \(F(i,j)\) 表示 \(2\) 到 \(i\) 中所有质数和所有最小质因子 \(>p_j\) 的合数的 \(f\) 的函数值的和。
初始时 \(F(n,m)=G(n,m)\)。
对于 \(i\le p_j^2\) 枚举 \(e\) 满足 \(p_j^e\le i\)，那么转移方程为 \(F(i,j)=G(i,m)-s(j)+\sum_{j<k\le m\wedge p_k^{e+1}\le i}f(p_k^e)G\left(\left\lfloor\frac i{p_k^e}\right\rfloor,k\right)+f(p_k^{e+1})\)。
同 \(G\)， \(i\) 的取值一定是 \(\left\lfloor\frac nd\right\rfloor\) 种类数量为是 \(O(\sqrt n)\) 的。可以记忆化搜索。
时间复杂度：
不会证明，但是根据论文所说，在 \(n\le 10^{13}\) 时时间复杂度 \(O\left(\frac{n^{3/4}}{\log n}\right)\)，在 \(n\) 足够大时，时间复杂度是 \(O(n^{1-\epsilon})\)，其中 \(\epsilon\) 为无穷小量。空间复杂度为 \(O(\sqrt n)\)。
例题
求 \(\sum_{i=1}^n\frac1{d(i)}\bmod p\) 。
题解
令 \(f(i)=\frac1{d(i)}\)，构造 \(g(i)=\frac12\) ，那么 \(s(i)=\frac i2\)，\(T(x,j)=x\)，\(f(p^k)=\frac1{k+1}\)，直接运用上述方法计算即可。
代码
#include<cstdio>
#include<algorithm>
#include<cmath>
typedef long long ll;
typedef __int128_t lll;
const ll N=2000005;
ll n,p,prime[N],inv[N],cnt,fs[N],fb[N],x[N],m;
bool vis[N];
inline void init(ll n){
cnt=0;
for(ll i=2;i*i<=n;i++){
if(!vis[i])prime[++cnt]=i;
for(ll j=1;j<=cnt&&i*i*prime[j]*prime[j]<=n;j++){
vis[i*prime[j]]=1;
if(i%prime[j]==0)break;
}
}
prime[cnt+1]=N;
}
inline ll&f(ll x){
return lll(x)*x<=n?fs[x]:fb[n/x];
}
inline ll fix(ll x){
return(x%p+p)%p;
}
inline ll s(ll x){
return fix(x*inv[2]);
}
inline void prework(){
init(n);
for(ll l=1,r;l<=n;l=r+1)r=n/(n/l),x[++m]=n/l,f(x[m])=fix((x[m]-1)*inv[2]);
for(ll i=1;i<=cnt;i++)for(ll j=1;x[j]>=prime[i]*prime[i];j++){
f(x[j])=fix(f(x[j])-f(x[j]/prime[i])+s(i-1));
}
}
ll F(ll n,ll i){
if(i<cnt&&prime[i+1]>n)return 0;
ll res=fix(f(n)-s(i));
for(ll j=i+1;prime[j]*prime[j]<=n;j++){
for(ll k=1,x=prime[j];x*prime[j]<=n;k++,x*=prime[j]){
res=fix(res+F(n/x,j)*inv[k+1]+inv[k+2]);
}
}
return res;
}
int main(){
scanf("%lld%lld",&n,&p),inv[1]=1;
for(ll i=2;i<=40;i++)inv[i]=inv[p%i]*(p-p/i)%p;
prework(),printf("%lld\n",fix(F(n,0)+1));
return fflush(stdout),fclose(stdin),fclose(stdout),0;
}