洲阁筛学习笔记

一、概述

在学习洲阁筛之前，请确保你对 \(\min25\) 筛有充分了解。

洲阁筛和 \(\min25\) 筛类似，都是对埃氏筛的优化。

不同之处在于，洲阁筛是块筛（对所有形如 \(\lfloor\frac ni\rfloor\) 的位置求前缀和），而 \(\min25\) 筛通常只求单点。

洲阁筛复杂度是严格 \(\mathcal O(\frac{n^{0.75}}{\log n})\) ，而 \(\min25\) 筛在 \(n\le 10^{13}\) 时可以被证明是 \(\mathcal O(\frac{n^{0.75}}{\log n})\) ，当 \(n\to\infty\) 时是 \(\mathcal O(n^{1-\epsilon})\) 。

数据范围较小时，\(\min25\) 筛因为常数较小所以更占优势，但是数据范围较大时更推荐用洲阁筛。

有空补一个效率对比图，先鸽着。

二、素数部分

洲阁筛和 \(\min25\) 筛的第一部分完全相同。

三、合数部分

温馨提示：下文沿用了 \(\min25\) 筛中的记号。

定义：

\[h(n,j)=\sum_{i=1}^n[i\in\text{prime}\vee\text{minp}(i)\ge j]f(i)\\ \]

即经过 \(j-1\) 轮埃氏筛后未被筛掉的所有数贡献之和。

初始值 \(h(n,cnt)=G(n)\) 。

转移枚举新加入的数含 \(p_j\) 的幂次 \(e\) ：

\[h(n,j)=h(n,j+1)+\sum_{e\ge 1,p_j^{e+1}\le n}\left(f(p_j^e)\cdot\big(h(\lfloor\frac n{p_j^e}\rfloor,j+1)-s_j\big)+f(p_j^{e+1})\right)\\ \]

for(int j=cnt;j>=1;j--)
    for(int i=1;p[j]*p[j]<=val[i];i++)
        for(int x=p[j],e=1;x*p[j]<=val[i];x*=p[j],e++)
            ///更新h[i]

最终 \(\sum_{i=1}^nf(i)=f(1)+h(n,1)\) 。

四、卡常技巧

都来学洲阁筛了，一定是碰到用 \(\min25\) 筛然后被卡常的题了吧。

• 整数除法改浮点数除法。

  整数除法比浮点数除法实在慢太多了！刚好整除分块里面又要大量使用除法，给这种优化方法留出用武之地。

  #define div _div///防止函数重名
  int div(int x,int j)
  {
      return (double)x/j;
  }
  

  分母大多是 \(\sqrt n\) 以内的数，理论上预处理 inv[j]=1.0/j 然后除法改乘法会更快，但是这样做掉精度非常严重。

• 更新 h[i] 前预处理 val[i]/p[j] 的值，从而省掉判断循环终止条件时的乘法。

  for(int j=cnt;j>=1;j--)
      for(int i=1;p[j]*p[j]<=val[i];i++)
          for(int x=p[j],e=1,lim=div(val[i],p[j]);x<=lim;x*=p[j],e++)
              ///更新h[i]
  

  这条优化还是比较有用的，因为如果 \(n^{1.5}\gt 2^{64}\) （比如 \(n=10^{13}\)）那么原版写法计算 x*p[j] 时需要强转 __int128 。

五、例题

放一道例题应该足够了。

例1、\(\texttt{LOJ6785 简单的函数 1e13版}\)

题目描述

\(f(x)\) 定义如下：

• \(f(1)=1\) 。
• \(f(p^c)=p\oplus c\) 。
• 若 \(\gcd(a,b)=1\) ，则 \(f(ab)=f(a)\cdot f(b)\) 。

令 \(S(m)=\big(\sum_{i=1}^mf(i)\big)\bmod 10^9+7\) 。

给定 \(n\) ，求 \(S(\lfloor\frac n1\rfloor),\cdots,S(\lfloor\frac nn\rfloor)\) 去重后的异或和。

数据范围

• \(1\le n\le 10^{13}\) 。

时间限制 \(\texttt{20s}\) ，空间限制 \(\texttt{512MB}\) 。

分析

令 \(f(2)=1\) （仅对素数部分有效），则 \(f(p)=p-1\) ，拆成完全积性函数 \(f_1(x)=x\) 和 \(f_2(x)=1\) 之差即可。

合数部分跑上述 \(\texttt{dp}\) ，别忘了当 \(n\ge 2\) 时最后要补上 \(f(2)\) 的贡献。

时间复杂度 \(\mathcal O(\frac{n^{0.75}}{\log n})\) ，下面这份代码要跑 \(\texttt{19s}\) ：

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int maxn=7e6+5,mod=1e9+7;
int B,m,n,cnt,res;
int p[maxn>>3];
bitset<maxn> b;
int id1[maxn],id2[maxn],val[maxn];
int c[maxn],g1[maxn],g2[maxn],h[maxn],s[maxn];
void init(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!b[i]) p[++cnt]=i;
        for(int j=1;j<=cnt&&i*p[j]<=n;j++)
        {
            b[i*p[j]]=1;
            if(i%p[j]==0) break;
        }
    }
    for(int i=1;i<=cnt;i++) s[i]=(s[i-1]+p[i])%mod;
}
int id(int x)
{
    return x<=B?id1[x]:id2[n/x];
}
signed main()
{
    scanf("%lld",&n),init(B=sqrt(n));
    for(int l=1,r=0;l<=n;l=r+1)
    {
        int x=n/l;
        r=n/x,val[++m]=x,x<=B?id1[x]=m:id2[n/x]=m;
        x%=mod,g1[m]=x-1,g2[m]=(x*(x+1)/2-1)%mod;
    }
    for(int j=1;j<=cnt;j++)
        for(int i=1;p[j]*p[j]<=val[i];i++)
        {
            int k=id(val[i]/p[j]);
            g1[i]=(g1[i]-(g1[k]-(j-1)))%mod;
            g2[i]=(g2[i]-p[j]*(g2[k]-s[j-1]))%mod;
        }
    for(int i=1;i<=m;i++) h[i]=(g2[i]-g1[i])%mod;
    for(int j=cnt;j>=1;j--)
        for(int i=1;p[j]*p[j]<=val[i];i++)
            for(int x=p[j],e=1;(__int128)x*p[j]<=val[i];x*=p[j],e++)
                h[i]=(h[i]+(p[j]^e)*(h[id(val[i]/x)]-(s[j]-j))+(p[j]^(e+1)))%mod;
    for(int i=1;i<=m;i++) c[i]=(1+h[i]+2*(val[i]>=2)+mod)%mod;
    sort(c+1,c+m+1),m=unique(c+1,c+m+1)-c-1;
    for(int i=1;i<=m;i++) res^=c[i];
    printf("%lld\n",res);
    return 0;
}

优化后的代码（\(\texttt{17s}\)）：

#include<bits/stdc++.h>
#define int long long
#define div _div
using namespace std;
const int maxn=7e6+5,mod=1e9+7;
int B,m,n,cnt,res;
int p[maxn>>3];
bitset<maxn> b;
int id1[maxn],id2[maxn],val[maxn];
int c[maxn],g1[maxn],g2[maxn],h[maxn],s[maxn];
void init(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!b[i]) p[++cnt]=i;
        for(int j=1;j<=cnt&&i*p[j]<=n;j++)
        {
            b[i*p[j]]=1;
            if(i%p[j]==0) break;
        }
    }
    for(int i=1;i<=cnt;i++) s[i]=(s[i-1]+p[i])%mod;
}
int div(int x,int j)
{
    return (double)x/j;
}
int id(int x)
{
    return x<=B?id1[x]:id2[div(n,x)];
}
signed main()
{
    scanf("%lld",&n),init(B=sqrt(n));
    for(int l=1,r=0;l<=n;l=r+1)
    {
        int x=n/l;
        r=n/x,val[++m]=x,x<=B?id1[x]=m:id2[n/x]=m;
        x%=mod,g1[m]=x-1,g2[m]=(x*(x+1)/2-1)%mod;
    }
    for(int j=1;j<=cnt;j++)
        for(int i=1;p[j]*p[j]<=val[i];i++)
        {
            int k=id(div(val[i],p[j]));
            g1[i]=(g1[i]-(g1[k]-(j-1)))%mod;
            g2[i]=(g2[i]-p[j]*(g2[k]-s[j-1]))%mod;
        }
    for(int i=1;i<=m;i++) h[i]=(g2[i]-g1[i])%mod;
    for(int j=cnt;j>=1;j--)
        for(int i=1;p[j]*p[j]<=val[i];i++)
            for(int x=p[j],e=1,lim=div(val[i],p[j]);x<=lim;x*=p[j],e++)
                h[i]=(h[i]+(p[j]^e)*(h[id(val[i]/x)]-(s[j]-j))+(p[j]^(e+1)))%mod;
    for(int i=1;i<=m;i++) c[i]=(1+h[i]+2*(val[i]>=2)+mod)%mod;
    sort(c+1,c+m+1),m=unique(c+1,c+m+1)-c-1;
    for(int i=1;i<=m;i++) res^=c[i];
    printf("%lld\n",res);
    return 0;
}

注意第 \(52\) 行 val[i]/x 不要套 div 函数，因为 x 很大，大数相除反而会变慢。

不知为何本地特别吃这两个优化，但是 \(\texttt{loj}\) 上确实体现不明显，也许这就是评测机差异吧。