Lucy-Hedgehog 算法

Lucy-Hedgehog 算法用于求解 \(\pi(n)\) ，\(\pi(n)\) 表示 \(1 \sim n\) 的整数中质数的个数。

Lucy-Hedgehog 算法由 Project Euler 用户 Lucy-Hedgehog 提出。（Lucy-Hedgehog 原文章）（我的转载）

算法思想：

回顾埃氏筛的过程，对于每个质数 \(p\)，筛掉 \(p\) 的倍数。

定义函数 \(S(v,p)\) 表示用不超过 \(p\) 的质数筛 \(1 \sim v\) 后剩余的数的个数（除去 \(1\)）。

比如 \(S(10,3)\)，\(1 \sim 10\) 内，通过 \(2\) 筛去 \(4,6,8,10\)，通过 \(3\) 筛去 \(6,9\)，除去 \(1\) 剩余 \(2,3,5,7\)。因此 \(S(10,3)=4\).

对于 \(S(v,p)\)

• \(p\) 不是质数，\(S(v,p)=S(v,p-1)\)；

• \(p*p>v\)，\(p\) 不可能再筛任何数，\(S(v,p)=S(v,p-1)\)；

• \(p\) 是质数且 \(p*p \le v\)，假设 \(x=p \times k\) 能被筛掉，则 \(p \le k \le v/p\) 且 \(k\) 不被不超过 \(p-1\) 的质数整除，这样的 \(k\) 共有 \(S(v/p,p-1)-S(p-1,p-1)\) 个。

于是我们发现一个重要的式子:

\(S(v,p)=S(v,p-1)-(S(v/p,p-1)-S(p-1,p-1))\) （\(p\) 是质数且 \(p \le \sqrt{v}\)）

这个第二维可以滚动数组，于是可以用 DP 实现这个算法

最后的答案 \(\pi(n)=S(n,n)\)。

long long v[MAXN],s[MAXN];
long long Lucy_Hedgehog(long long x){
    int r=(int)sqrt(x)+1,m=x/r-1;
    int cnt=0;
    for(int i=1;i<=m;i++)
    {
        v[++cnt]=i;
        s[v[cnt]]=v[cnt]-1;
    }
    for(int i=r;i>=1;i--)
    {
        v[++cnt]=x/i;
        s[v[cnt]]=v[cnt]-1;
    }
    for(int i=2;i<=r;i++)
    {
        if(s[i]>s[i-1]){
            long long sp=s[i-1];
            long long lim=1ll*i*i;
            for(int j=cnt;v[j]>=lim&&j;j--)
            {
                s[v[j]]-=(s[v[j]/i]-sp);
            }
        }
    }
    return s[x];
}

又有数论分块可以把数组空间开到 \(O(\sqrt{n})\)

long long v[MAXN],s[MAXN];
long long Lucy_Hedgehog(long long x){
    int r=(int)sqrt(x)+1,m=x/r-1;
    auto id=[&](long long i){
        if(i<=m){
            return (int)i;
        }
        else{
            return lower_bound(v+1,v+1+cnt,i)-v;
        }
    };
    int cnt=0;
    for(int i=1;i<=m;i++)
    {
        v[++cnt]=i;
        s[cnt]=v[cnt]-1;
    }
    for(int i=r;i>=1;i--)
    {
        v[++cnt]=x/i;
        s[cnt]=v[cnt]-1;
    }
    for(int i=2;i<=r;i++)
    {
        if(s[i]>s[i-1]){
            long long sp=s[i-1];
            long long lim=1ll*i*i;
            for(int j=cnt;v[j]>=lim&&j;j--)
            {
                s[j]-=(s[id(v[j]/i)]-sp);
            }
        }
    }
    return s[id(x)];
}

这个算法还可以求 \(1 \sim n\) 的整数中质数的和。

只需要把式子改为 \(S(v,p)=S(v,p-1)-p \times (S(v/p,p-1)-S(p-1,p-1))\) 就行了。

long long v[MAXN];
int128 s[MAXN];
int128 Lucy_Hedgehog(long long x){
    if(x<=1){
        return 0;
    }
    int r=(int)sqrt(x)+1,m=x/r-1;
    auto id=[&](long long i){
        if(i<=m){
            return (int)i;
        }
        else{
            return lower_bound(v+1,v+1+cnt,i)-v;
        }
    };
    int cnt=0;
    for(int i=1;i<=m;i++)
    {
        v[++cnt]=i;
        s[cnt]=(int128)v[cnt]*(v[cnt]+1)/2-1;
    }
    for(int i=r;i>=1;i--)
    {
        v[++cnt]=x/i;
        s[cnt]=(int128)v[cnt]*(v[cnt]+1)/2-1;
    }
    for(int i=2;i<=r;i++)
    {
        if(s[i]>s[i-1]){
            int128 sp=s[i-1];
            long long lim=1ll*i*i;
            for(int j=cnt;v[j]>=lim&&j;j--)
            {
                s[j]-=1ll*i*(s[id(v[j]/i)]-sp);
            }
        }
    }
    return s[id(x)];
}

总体时间复杂度 \(O(x^{\frac{3}{4}})\)，空间复杂度 \(O(\sqrt{x})\)。

这篇文章 写了如何将时间复杂度优化至 \(O(x^{\frac{2}{3}})\)