洲阁筛 学习笔记
洲阁筛
对具有如下性质的积性函数 \(F\)
-
\(p\in prime\) 时，\(F(p)=G(p)\)
-
$p\in prime $ 时，\(F(p^c)=T(p^c),c>1\)(如无特殊说明，和式里面的 \(p\) 同样指代质数)
-
其中，\(G,T\) 为关于 \(p,c\) 的低阶多项式或幂函数
洲阁筛的作用是在 \(O(\frac{n^{\frac{3}{4}}}{\log n})\) 的复杂度内求出 \(\sum\limits_{i\le n} F(i)\)的一种筛法
将所有数分成两个集合
- \(A\) 集合： \(i\) 没有大于 \(\sqrt n\) 的质因数
- \(B\) 集合： \(i\) 有大于 \(\sqrt n\) 的质因数
我们前半部分预处理 \(F\) ，对 \(G\) 整除分块，后面的 \(F\) 独立求。
考虑对每个 \(y=\lfloor \frac{n}{x}\rfloor\) 的情况求出
设 \(p_1,p_2,\dots,p_m\) 表示$\le \sqrt n $ 的从小到大的 \(m\) 个质数
设 \(g_{i,j}\) 表示在 \([1,j]\) 中所有与前 \(i\) 个质数互质的数 \(x\) 的低阶多项式的幂之和 \(\sum x^c\) (例如 \(G(p)=p\)，则需要统计互质的数 \(x\) 的 \(\sum x^1\) )
这样我们可以用 \(g_{m,y}\) 得到 \(\sqrt n <p\le y\) 中质数的相关信息(低阶多项式幂之和)
考虑转移
考虑到所有与前 \(i-1\) 个数互质且不与前 \(i\) 个数互质的数 \(x\) 一定可以被拆成 \(x=p\times res\)，我们减去这部分就可以了
考虑这个朴素转移的的复杂度，我们有
- 质数密度 \(\pi(n)=\frac{n}{\ln n}\)
- \(\lfloor \frac{n}{xp_i}\rfloor\) 的总个数不超过 \(\sqrt n\) 个（这里讨论了每个 \(y\) 的情况）
因此状态数是 \(O(\frac{n}{\log n})\) 的，转移是 \(O(1)\) 的，因此复杂度也是这个。
这里的第二维使用 \(Hash\) 表进行存储
参考Code:
const int mod=1145141;
int head[mod],Next[N<<1],tim[N<<1],cnt;
ll val[N<<1],g[N<<1],f[N<<1];
void ins(ll v)
{
int u=v%mod;
val[++cnt]=v,Next[cnt]=head[u],head[u]=cnt;
}
int qry(ll v)
{
int u=v%mod;
for(int i=head[u];i;i=Next[i])
if(val[i]==v)
return i;
return 0;
}
void Insert()
{
for(ll i=1;i<=n;i=(n/(n/i))+1) head[(n/i)%mod]=0;
cnt=0;
for(ll i=1;i<=n;i=(n/(n/i))+1)
{
ins(n/i);
//something else
}
}
- 这份代码是对 \(\lfloor \frac{n}{xp_i}\rfloor\) 从大到小插入的
- 长度要开 \(2\sqrt n\)
考虑到当 \(p_{i-1}^2\le j < p_i^2\) 时
其中 \(g_{i-1,\lfloor\frac{j}{p_i}\rfloor}=1\)，因为 \(\lfloor\frac{j}{p_i}\rfloor <p_i\)，则 \([1,\lfloor\frac{j}{p_i}\rfloor]\) 中没有与前 \(i-1\) 个质数互质的数字
则
这样的话我们考虑当 \(j<p_i^2\) 时先不计算这一部分的式子，当需要用它转移的时候我们再一下子把所有没计算的一起加上就可以了。
这里还是蛮抽象的，因此给出 \(c=1\) 时的代码
参考Code:
void getg()
{
for(ll i=1;i<=n;i=(n/(n/i))+1)
{
g[cnt]=n/i;
tim[cnt]=0;
}
for(int i=1;i<=m;i++)
for(int j=1;j<=cnt&&1ll*pri[i]*pri[i]<=val[j];j++)
{
int v=qry(val[j]/pri[i]);
g[j]-=g[v]-(i-1-tim[v]);
tim[j]=i;
}
}
- 采用了滚动数组
- \(tim_j\) 表示 \(j\) 最后一次被转移到的时间
复杂度？
这样我们就解决了第一个部分
考虑计算
还是将 \(\sqrt n\) 的质数从小到大排列成 \(m\) 个
令 \(f_{i,j}\) 代表在 \([1,j]\) 中分解质因数后只包含 后 \(i\) 个质数的数的 \(F(x)\)
则有
朴素转移
复杂度分析和 \(g\) 一样
考虑到当 \(p_{i-1}^2\le j < p_i^2\) 时
和上面一样，这一部分是不需要计算的，我们只需要预处理后在需要的时候把 \(\sum F(p_i)\) 累加进去就可以了
复杂度分析同上，最后的复杂度也是 \(O(\frac{n^{\frac{3}{4}}}{\log n})\) 的
以 \(T(p^c)=3c+1\) 为例
参考Code:
void getf()
{
for(int i=m;i;i--)
{
for(int j=1;j<=cnt&&1ll*pri[i]*pri[i]<=val[j];j++)
{
if(val[j]<pri[i+1])
f[j]=1;
else if(val[j]<1ll*pri[i+1]*pri[i+1])
f[j]=(pris[min(N-1ll,val[j])]-pris[pri[i+1]-1])*4ll+1;
for(ll c=1,pr=pri[i];pr<=val[j];pr=pr*pri[i],c++)
{
ll k=val[j]/pr,v;
if(k<pri[i+1])
v=1;
else if(k<1ll*pri[i+1]*pri[i+1])
v=(pris[min(N-1ll,k)]-pris[pri[i+1]-1])*4ll+1;
else
v=f[qry(k)];
f[j]+=v*(3*c+1);
}
}
}
}
- \(pris_i\) 表示 \([1,i]\) 中素数数量
最后的合并就比较简单，直接按照整除分块差不多的方法就可以了
以为DIVCNT3 - Counting Divisors (cube)为例
参考Code:
#include <cstdio>
#include <cctype>
#include <algorithm>
#define ll long long
using std::min;
const int SIZE=1<<21;
char ibuf[SIZE],*iS,*iT;
//#define gc() (iS==iT?(iT=(iS=ibuf)+fread(ibuf,1,SIZE,stdin),iS==iT?EOF:*iS++):*iS++)
#define gc() getchar()
template <class T>
void read(T &x)
{
x=0;char c=gc();
while(!isdigit(c)) c=gc();
while(isdigit(c)) x=x*10+c-'0',c=gc();
}
const int N=316241;
int ispri[N],pri[N],pmi[N],c[N],pris[N],m;
ll F[N],SF[N],n;
void init()
{
ispri[1]=F[1]=1;
for(int i=2;i<N;i++)
{
pris[i]=pris[i-1]+(!ispri[i]);
if(!ispri[i])
{
pri[++m]=i;
pmi[i]=i;
c[i]=1;
F[i]=4;
}
for(int j=1;j<=m&&pri[j]*i<N;j++)
{
int x=pri[j]*i;
ispri[x]=1;
if(i%pri[j])
{
c[x]=1;
pmi[x]=pri[j];
F[x]=F[i]*4;
}
else
{
c[x]=c[i]+1;
pmi[x]=pmi[i]*pri[j];
F[x]=F[i/pmi[i]]*(3*c[x]+1);
break;
}
}
}
pri[m+1]=N;//写这个少一个特判
for(int i=1;i<N;i++) SF[i]=SF[i-1]+F[i];
}
//g_{i,j} 在[1,j]中与前i个质数(从小到大)互质的数p的\sum p^0(即个数)
//这里标准写法应该看做0阶多项式
//g_{i,j}=g_{i-1,j}-p^0*g_{i-1,j/p}
//g_{0,j}=j
const int mod=1145141;
int head[mod],Next[N<<1],tim[N<<1],cnt;
ll val[N<<1],g[N<<1],f[N<<1];
void ins(ll v)
{
int u=v%mod;
val[++cnt]=v,Next[cnt]=head[u],head[u]=cnt;
}
int qry(ll v)
{
int u=v%mod;
for(int i=head[u];i;i=Next[i])
if(val[i]==v)
return i;
return 0;
}
void getg()
{
for(ll i=1;i<=n;i=(n/(n/i))+1) head[(n/i)%mod]=0;
cnt=0;
for(ll i=1;i<=n;i=(n/(n/i))+1)//将需要的值 从 大 到 小 的编号塞到hash表去
{
ins(n/i);
g[cnt]=n/i;
tim[cnt]=0;
}
//p<=j<p^2时
//g_{i,j}=g_{i-1,j}-1,记录tim[j]从哪个i开始，j没有处理
//询问到它的时候删除即可
for(int i=1;i<=m;i++)
for(int j=1;j<=cnt&&1ll*pri[i]*pri[i]<=val[j];j++)//p^2>j时，懒惰处理
{
int v=qry(val[j]/pri[i]);
g[j]-=g[v]-(i-1-tim[v]);
tim[j]=i;
}
}
//f_{i,j} 在[1,j]中分解后仅为前i个质数(从大到小)的数x的\sum F(x)
//f_{i,j}=f_{i-1,j}+\sum_{c} f_{i-1,j/(p_i^c)}
//f_{0,j}=1,即F(1)
void getf()
{
//p<=j<p^2时
//f_{i,j}=f_{i+1,j}+F(p_i)
for(int i=m;i;i--)
{
for(int j=1;j<=cnt&&1ll*pri[i]*pri[i]<=val[j];j++)
{
//判断是不是从这一步开始才可以被转移
if(val[j]<pri[i+1])
f[j]=1;//F(1)
else if(val[j]<1ll*pri[i+1]*pri[i+1])
f[j]=(pris[min(N-1ll,val[j])]-pris[pri[i+1]-1])*4ll+1;//\sum_{p} F(p)+F(1)
for(ll c=1,pr=pri[i];pr<=val[j];pr=pr*pri[i],c++)
{
ll k=val[j]/pr,v;
if(k<pri[i+1])
v=1;
else if(k<1ll*pri[i+1]*pri[i+1])
v=(pris[min(N-1ll,k)]-pris[pri[i+1]-1])*4ll+1;
else
v=f[qry(k)];
f[j]+=v*(3*c+1);
}
}
}
}
int main()
{
init();
int T;read(T);
while(T--)
{
read(n);
if(n<N)
{
printf("%lld\n",SF[n]);
continue;
}
getg(),getf();
ll ans=0;
for(ll l=1,r,t;l<N;l=r+1)
{
int v=qry(n/l);
r=min(N-1ll,n/(n/l));
if(n/l<N) t=1;
else t=g[v]-(m-tim[v]);
ans+=(SF[r]-SF[l-1])*(t-1)*4;
}
printf("%lld\n",ans+f[1]);
}
return 0;
}
2019.6.12
