min_25筛 学习笔记
min_25筛 学习笔记
简介
min_25筛是一种筛法,由min_25提出,能够 \(O(\dfrac{n^{\frac{3}{4}}}{\ln n})\) 的求积性函数 \(f\) 的前缀和。而 \(f\) 需要满足如下几个条件:
(设 \(p\) 为任意质数)
- \(f(p)\) 可以写成低次多项式的形式。设这个多项式为 \(g(x)\)。
- \(f(p^k)\) 很方便计算
约定
要求前缀和的积性函数为 \(f\),在质数位置可以写成多项式 \(g\)。
设第 \(i\) 个质数为 \(p_i\)。设 \([1,n]\) 中质数总共有 \(|p|\) 个。
\(n\) 的“基本和组”定义为集合 \(B(n)=\{x|x=n/k, k\in[1,n]\}\)。我们知道 \(|B(n)|\) 大约是 \(2\sqrt{n}\),根据整除分块。
求法
part-1 质数位置
首先我们求出 \([1,n]\) 中质数位置的 \(f\) 的和,也就是 \(g\) 的和。
\(g\) 是一个低次的多项式,因此我们可以把它的每一次单独拿出来做。那我们相当于要求 \([1,n]\) 中的质数的 \(k\) 次幂和。
考虑埃氏筛的过程。每次我们选一个质数 \(p\),并且把它的倍数筛去。假设我们只看 “有效筛” 的位置:即,如果这个位置已经被筛掉,就不去动它。如果 \(x\) 被“有效筛” 了,那 \(x\) 最小的质因子是 \(p\),且 \(x\) 为合数。那么 \(x\ge p^2\)。
也就是说,这个过程中,我们只需要 \(\le \sqrt{x}\) 的那些质数。这些质数是可以预先筛出来的。
我们记下每次筛完之后的结果。换句话说,记 \(G(n,i)\) 表示,\([1,n]\) 中,只留下素数与最小质因子 \(> p_i\) 的合数,这些数的 \(g\) 的和。
考虑 \(G(n,i-1)\to G(n,i)\),相当于要把 \(i\) 能 “有效筛” 的那些位置的贡献去掉。
设 \(x\) 能被“有效筛”,\(x\le n\)。那 \(x\) 一定是 \(p_i\) 的倍数,且 \(x/p_i\) 的最小质因子 \(\ge p_i\),也就是说,\(> p_{i-1}\)。
注意到 \(x\le n\),那 \(x/p_i\le n/p_i\)。且 \(x\) 的最小质因子 \(>p_{i-1}\)。满足条件的这些 \(x\) 大约就是 \(G(n/p_i,i-1)\)。但是,这个 \(G\) 还把质数算上去了,我们要抠掉比 \(p_i\) 小的质数。 (注意到,如果等于 \(p_i\) 是没问题的,这相当于 \(x=p_i^2\),满足的)
这相当于前缀若干个质数的 \(g\) 的和。可以在筛质数的时候顺便求个前缀和得到,记 \(sump(i)\) 表示前 \(i\) 个质数的 \(g\) 和。
那么我们得到
其中 \(\Delta\) 表示变化量,即 \(G(n,i-1)-G(n,i)\)。 也就是说我们令 \(G(n,i)=G(n,i-1)-\Delta\) 即可。
边界:\(G(n,0)=\sum\limits_{i=2}^{n} g(k)\),把 \(g\) 的每一项拆开,它就可以直接计算了。
备注:
我们在实现的时候,对 \(i\) 这一维滚动,且只保留基本和组位置的值。
由于后面只需要 \(G(*,tot)\) 位置的值,所以后面直接记
\[G(i)=\sum\limits_{p\le i,p\in prime} g(p) \]
可以用积分证明,这部分复杂度是 \(O(\dfrac{n^{\frac{3}{4}}}{\ln n})\) 的。
part-2 求答案
我们已经知道了质数位置的 \(g\) 和。
类似的考虑,设 \(S(n,i)\) 表示,最小质因子 \(\ge p_i\) 的那些数的 \(f\) 的和。
考虑一下边界:\(S(n,|p|+1)=0\) 。显然,因为 \([1,n]\) 中的数,最小质因子都 \(<p_{|p|}\)。
然后我们要求的答案是 \(S(n,1)+f(1)\)。这也显然,因为 \(1\) 以上所有数的最小质因子都 \(\ge 2\)。那 \(S(n,1)\) 就会覆盖到 \(2\) 以上所有数。
注意到我们是知道一个大值,然后要求的是一个小值。那考虑从大往小推。即,假设大的已知,做小的。
假设现在要做 \(S(n,i)\)。
首先加上满足条件质数的答案。这玩意就是 \([1,n]\) 中质数的 \(g\) 和,减去 \(<p_i\) 的质数的 \(g\) 和。即,\(G(n)-sump(i-1)\)。
然后是满足条件合数的答案。如何钦点合数呢?合数就至少有两个质因数。枚举它的最小质因数是谁,然后枚举它多少次,然后枚举后面选了啥数。
设枚举的最小质因数是 \(j\ (j\ge i)\),然后 \(p_j\) 有 \(e\) 次方。那它后面可能还跟了一个 \(p_j\),也可能跟的是比 \(p_j\) 大的质因数。
注意到 \(f(p_j^{e+1})\) 需要特判一下,其它的可以用 \(S\) 写出来。脑子转一下,得到它的贡献为
为啥后面直接乘起来是对的,因为 \(f\) 是积性函数。
需要满足:\(p_j^{e+1}\le n\)。这样枚举一遍,把贡献加起来即可。
最后 \(S(n,1)+1\) 就是答案力!!1
板子题
loj6053 简单的函数
分析一下,发现它在质数的时候可以写成 \(g(p)=p-1\) 的形式,特判 \(2\)。
然后就是个多项式,用 min-25 筛就可以做了。
代码中有注释,可以参考一下实现。
#include <bits/stdc++.h>
using namespace std;
namespace Flandre_Scarlet
{
#define N 200005
#define int long long
#define mod 1000000007
#define iv2 500000004
#define F(i,l,r) for(int i=l;i<=r;++i)
#define D(i,r,l) for(int i=r;i>=l;--i)
#define Tra(i,u) for(int i=G.h[u],v=G.to(i);~i;i=G.nx(i),v=G.to(i)) if (i>=0)
#define MEM(a,x) memset(a,x,sizeof(a))
#define FK(a) MEM(a,0)
#define sz(x) ((int)x.size())
#define all(x) x.begin(),x.end()
#define p_b push_back
#define pii pair<int,int>
#define fir first
#define sec second
int I() {char c=getchar(); int x=0; int f=1; while(c<'0' or c>'9') f=(c=='-')?-1:1,c=getchar(); while(c>='0' and c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar(); return ((f==1)?x:-x);}
template <typename T> void Rd(T& arg){arg=I();}
template <typename T,typename...Types> void Rd(T& arg,Types&...args){arg=I(); Rd(args...);}
void RA(int *p,int n) {F(i,1,n) *p=I(),++p;}
int pr[N]; bool notp[N];
int sp[N];
void Init() // 筛到sqrt(n)
{
int n=2e5;
int&c=pr[0];
notp[1]=1;
F(i,2,n)
{
if (!notp[i])
{
pr[++c]=i;
}
for(int j=1;j<=c and i*pr[j]<=n;++j)
{
int u=pr[j];
notp[i*u]=1;
if (i%u==0) break;
}
}
F(i,1,c) sp[i]=sp[i-1]+pr[i];
}
int n;
void Input()
{
n=I();
}
int g1[N],g0[N]; // 1次方, 0次方的和
int w[N],i1[N],i2[N];
// w: 存储基本和组位置; i1,i2: 给这些位置编号
// 若x<=sn,i1[x]保存x的编号; 否则, i2[n/x]保存x的编号
int sn,tot;
int s1(int x){x%=mod; return x*(x+1)%mod*iv2%mod;}
int S(int x,int p) // [1,x], min-p>=pr[p], f sum
{
if (x<=1 or pr[p]>x) return 0;
int k=(x<=sn)?(i1[x]):(i2[n/x]);
int ans=(g1[k]-g0[k])-(sp[p-1]-(p-1)); // 质数部分的答案
ans=(ans%mod+mod)%mod;
if (p==1) ans+=2; // 特判2: 我们把f(2)当成了1, 实际上它是3
for(int i=p;i<=pr[0] and pr[i]<=x/pr[i];++i)
{
int u=pr[i];
int e=1,ue=u;
// ue为u的e次方
while(ue*u<=x) // u的e+1次方<=x
{
ans+=((u^(e+1))+S(x/ue,i+1)*(u^e)%mod)%mod;
// u^e为 f(u的e次方), u^(e+1)同理
ans%=mod;
ue*=u; ++e;
}
}
return (ans%mod+mod)%mod;
}
void Sakuya()
{
sn=sqrt(n),tot=0;
for(int l=1,r;l<=n;l=r+1)
{
r=(n/(n/l));
w[++tot]=n/l;
if (n/l<=sn) i1[n/l]=tot;
else i2[n/(n/l)]=tot;
g0[tot]=((n/l)-1)%mod;
g1[tot]=(s1(n/l)-1)%mod;
// G(x)初始为g(2)+g(3)...+g(x)
// 这里拆成两个, 因此 g0=x-1, g1=2+3...+(x-1)
}
F(j,1,pr[0])
{
int u=pr[j];
for(int i=1;i<=tot and u*u<=w[i];++i)
{
int t=w[i]/u,k;
if (t<=sn) k=i1[t];
else k=i2[n/t];
g0[i]-=g0[k]-(j-1);
g1[i]-=(g1[k]-sp[j-1])%mod*u%mod;
g0[i]=(g0[i]%mod+mod)%mod;
g1[i]=(g1[i]%mod+mod)%mod;
}
}
// 计算G
int ans=S(n,1)+1; // f(1)=1, 加上
printf("%lld\n",ans);
}
void IsMyWife()
{
Init();
Input();
Sakuya();
}
}
#undef int //long long
int main()
{
Flandre_Scarlet::IsMyWife();
return 0;
}