素数计数函数

一种计算$\pi(n)$的组合方法

这里$\pi(n)$指：不大于$n$的素数个数

注：本文方法来自于维基百科对素数计数的一个组合方法：https://en.wikipedia.org/wiki/Prime-counting_function

由于最后的一些优化技巧还未掌握。文中的方法还有待优化。

与原文的方法略微不同。与洲阁筛的方法近似。。。

常识：一般来说，$x$以内的素数大约有：

$$\pi(x)=O(\frac{x}{ln\ x})$$

那么我们可以猜测埃式筛法的复杂度：

$$T(n)=O\Big(n\sum_{p\leq n}\frac{1}{p}\Big)$$

可以认为是线性的。

在$n$不是很大的时候。埃式筛法是个不错的选择。(不知道素数筛法的可以自行百度。)

埃式筛法在计算$\pi(n)$时。还会得到素数表。

像埃式筛法这样：通过得到所有不大于$n$的素数，来计算$\pi(n)$。算法复杂度的下界是$O(\frac{n}{log n})$。不会很快。

组合方法：

对于一个整数 $q$，$q\in(\sqrt{n},n]$

$q$如果不能被$[2,\sqrt{n}]$的任何整数整除。

则：$q$是一个素。

很容易得到一个基于容斥原理的式子：

$[1,n]$中，不能被前$m$个素数整除的数字个数为：

$$C(m,n)=\Big\lfloor n\Big\rfloor-\sum_{1\leq i\leq m}\Big\lfloor\frac{n}{P_i}\Big\rfloor+\sum_{1\leq i<j\leq m}\Big\lfloor\frac{n}{P_iP_j}\Big\rfloor-\sum_{1\leq i<j<k\leq m}\Big\lfloor\frac{n}{P_iP_jP_k}\Big\rfloor...$$

$P_i$指第$i$个素数

则：

$$C(\pi(\sqrt{n}),n)=\pi(n)-\pi(\sqrt{n})+1$$

$$C(m,n)=C(m-1,n)-C(m-1,\Big\lfloor\frac{n}{P_m}\Big\rfloor)$$

对于第二个式子。

$[1,n]$中不能被前$m-1$个素数整除的数字包含：

可以被$P_m$整除的数字。但不能被前$m-1$个素数整除的数字

对于$P_m$的所有倍数：

$$kP_m,\ \ 其中 \ k\in[1,\Big\lfloor\frac{n}{P_m}\Big\rfloor]$$

不能被前$m-1$个素数整除的$k$的数量为：

$$C(m-1,\Big\lfloor\frac{n}{P_m}\Big\rfloor)$$

所以：

$$C(m,n)=C(m-1,n)-C(m-1,\Big\lfloor\frac{n}{P_m}\Big\rfloor)$$

特别的：

$$C(0,n)=\lfloor n\rfloor$$

方便起见，下文所有数均为整数。

对于 $C(m,n)$定义：

$$D(m,d)=C(m,\Big\lfloor\frac{n}{d}\Big\rfloor)$$

因为：

$$\Big\lfloor\frac{\lfloor\frac{n}{a}\rfloor}{b}\Big\rfloor=\Big\lfloor\frac{n}{ab}\Big\rfloor$$

所以：

$$D(m,d)=D(m-1 , d)-D(m-1,dP_m)$$

那么，通过D数组与C数组交替递推。并且选择一个恰当的分界$B$我们计算$D[][1...B]$和$C[][1...\lfloor\frac{n}{B+1}\rfloor]$

有： $$dP_m\leq B:D[m][d]=D[m-1][d]-D[m-1][dP_m]\\dP_m>B:D[m][d]=D[m-1][d]-C[m-1][\lfloor\frac{n}{dP_m}\rfloor]$$

$$C[m][j]=C[m-1][j]-C[m][\lfloor\frac{j}{P_m}\rfloor]$$

在$B=\lfloor\sqrt{n}\rfloor$时，这样计算的复杂度:$O(\frac{n}{logn})$

并不比前面说到的方法快多少。

优化1：

通过上面的容斥原理很容易有：

当$P_{m+1}>n$时：

$$C(m,n)=1$$

则:$P_{m}^2>n$时：

$$C(m,n)=C(m-1,n)-C(m-1,\Big\lfloor\frac{n}{P_m}\Big\rfloor)\\=C(m-1,n)-1$$

优化2:

当$P_m^2\geq \lfloor\big(\frac{n}{B+1}\big)^{\frac{1}{4}}\rfloor$时：

$$C(m,j)=\pi\Big( max(P_m,j)\Big )-m+1$$

综上 。我们分两段递推：

第一段，前：$P_m\leq\lfloor\big(\frac{n}{B+1}\big)^{\frac{1}{4}}\rfloor$。正常计算

第二段。此时$P_m^2>\lfloor\big(\frac{n}{B+1}\big)^{\frac{1}{4}}\rfloor$

对于第二个优化。预处理出前$\lfloor\sqrt{n}\rfloor$个$\pi[]$.

此时不在更新$C[][]$，

如果需要C数组的信息。利用：优化2的式子得到。

对于D数组，应用优化1：

$P_m^2>\lfloor\frac{n}{d}\rfloor$时,即：$d>\lfloor\frac{n}{P_m^2}\rfloor$

$$D(m,d)=D(m-1,d)-1$$

可以肯定的是。不对$d\in[\lfloor\frac{n}{P_m^2}\rfloor+1,B]$更新

仅仅维护$d\in[1,\lfloor\frac{n}{P_m^2}\rfloor]$时。

我们记录最早开始不更新的哪个素数标号。并预处理前缀和。必要时刻查前缀和的表。即可。对于没有记录的区间有需要更新的区间。我们暴力更新。时间复杂度不会增加。（可以自行证明。很简单）

还有一种笨的方法。也不会很慢

笨的方法就是：使用区间数据结构来维护(总复杂度多了个 log？)。

建议使用数状数组。（常数小。也就慢了700ms。。。。

所以第二阶段。我们仅仅计算了$d\in[1,\lfloor\frac{n}{P_m^2}\rfloor]$

那么总多时间复杂度：

$$T(n)=\sum_{P\leq (\frac{n}{B+1})^{1/4}}\Big(O(B)+O(\frac{n}{B})\Big)+\sum_{ (\frac{n}{B+1})^{1/4}< P\leq n^{1/2}}O(\frac{n}{p^2})$$

在$B=\sqrt{n}$时：$T(n)=O(\frac{n^{\frac{3}{4}}}{log n})$

但是中间维护咱们不是 $O(1)$。

使用数据结构的话是：$T(n)=O(n^{\frac{3}{4}})$

注意：修改对C的定义。让其变为[2,n]上的数并不会优化计算。这是因为中间的递推依然会多出来一个余项。

下面是代码。。。(可以快速计算100亿以内的答案。)

#include <algorithm>
#include <stdio.h>
#include <string.h>
#include <cmath>
#define MAXN 1111111
using namespace std;
typedef long long LL;
struct arry
{
    int A[MAXN];
    int n;

    arry()
    {
        memset(A,0,sizeof A);
    }

    void clear(int m)
    {
        memset(A,0,(m+1)*sizeof(int));
        n=m;
    }

    int lowbit(int x)
    {
        return x&(-x);
    }

    void add(int x,int key)
    {
        while(x<=n)
        {
            A[x]+=key;
            x+=lowbit(x);
        }
    }

    int sum(int x)
    {
        int ans=0;
        while(x)
        {
            ans+=A[x];
            x-=lowbit(x);
        }
        return ans;
    }
}S;

int prim[MAXN],deep=1;
int pi[MAXN];
LL G[MAXN];
LL C[MAXN];

void init()
{
    for(int i=2;i<MAXN;i++)pi[i]=1;
    for(int i=2;i<MAXN;i++)
    {
        if(!pi[i])continue;
        prim[deep++]=i;
        for(int j=i<<1;j<MAXN;j+=i) pi[j]=0;
    }
    for(int i=2;i<MAXN;i++) pi[i]+=pi[i-1];
}

void clat_1(int m,int k,LL n)
{
    LL p=prim[k];
    for(int i=1;i<m;i++)
    {
        LL d=i*p;
        LL u=n/d;
        if(u<m)
            G[i]-=C[u];
        else
            G[i]-=G[d];
    }
    for(int i=m;i;i--)C[i]-=C[i/p];
}

LL slove(LL n)
{
    if(n<MAXN)return pi[n];
    int m=sqrt(n)+1.1;
    int n_4=pow(n,1.0/4.0)+1.1;
    for(int i=1;i<=m;i++)
    {
        C[i]=i;
        G[i]=n/i;
    }
    int k;
    for(k=1;prim[k]<n_4;k++)clat_1(m,k,n);
    S.clear(m);
    while(prim[k]<m)
    {
        LL p=(LL)prim[k]*prim[k];
        LL lim=n/p+1;
        for(int d=1;d<lim;d++)
        {
            LL u=(LL)d*prim[k];
            LL b=n/u;
            if(b<m)
            {
                if(b<=prim[k-1])
                    G[d]-=1;
                else
                    G[d]-=pi[b]-k+2;
            }
            else
                G[d]-=G[u]-S.sum((int)u);
        }
        S.add((int)lim,1);
        k++;
    }
    return k+G[1]-2;
}

int main ()
{
    init();
    LL n;
    while(scanf("%lld",&n)==1)  printf("%lld\n",slove(n));
    return 0;
}

增加内容(与维基百科上的一致):

上面的递推有一个简化和推广。（其实是针对上面递推变形出来的一个）

记：$G_k(i,j)$表示：$[1,j]$上，由$k$个大于$P_i$的素数组成的数的数量。

那么有 $$C(i,j)=\sum_{k=0}^{∞}G_k(i,j)$$

那么： $$C(\pi(j^{\frac{1}{3}}),j)=\sum_{k=0}^{∞}G_k(\pi(j^{\frac{1}{3}}),j)\\=\sum_{k=0}^{2}G_k(\pi(j^{\frac{1}{3}}),j)+\sum_{k=3}^{∞}G_k(\pi(j^{\frac{1}{3}}),j)$$

明显1： $$\sum_{k=3}^{∞}G_k(\pi(j^{\frac{1}{3}}),j)=0$$

明显2： $$G_0(i,j)=1\\G_1(i,j)=\pi\big(max(j,P_i)\big)-i$$

对于$G_2(\pi(j^{\frac{1}{3}}),j)$.它的可选素数是有$P>j^{\frac{1}{3}}$

那么对于另一个素数$P'>j^{\frac{1}{3}}$,且$PP'\leq j\ ,\ P\leq P'$.

这样的素数$P'$个数为：

$$\pi(\frac{j}{P})-\pi(P)+1$$

所以： $$G_2(\pi(j^{\frac{1}{3}}),j)=\sum_{j^{1/3}<P\leq j}\Big(\pi(\frac{j}{P})-\pi(P)+1\Big)$$

综上：

$$\pi(j)=C(\pi(j^{\frac{1}{3}}),j)-G_0(\pi(j^{\frac{1}{3}}),j)-G_2(\pi(j^{\frac{1}{3}}),j)+\pi(j^{\frac{1}{3}})$$

这种方法 。感觉前一部还是要打表出$n^{1/4}$比较好。我比较弱。还是不知道优化的优雅方法。

这种方法空间复杂度比较高。有些人貌似是部分记忆话。效率还很高呢。

#include<stdio.h>
#include<algorithm>
#include<string.h>
#include<math.h>
using namespace std;
typedef long long ll;
#define INF 30000000000
#define chk(pos, i) (((pos[i/64])&(1<<((i>>1)&31))))
#define set(pos, i) (((pos[i/64])|=(1<<((i>>1)&31))))
#define check(x) ((x&&(x&1)&&(!chk(pos, x)))||(x==2))
const int N = 101;
const int M = 49500;
const int P = 700000;
const int UP = 5000000;
ll tmp[N][M];
unsigned int pos[157000] = {0};
int len = 0, pm[P], cnt[UP];
void init(){
    set(pos, 0), set(pos, 1);
    for(ll i = 3; i*i < UP; i += 2){
        if(!chk(pos, i)){
            ll k = i<<1;
            for (ll j = (i * i); j < UP; j += k) set(pos, j);
        }
    }
    for(int i = 1; i < UP; ++i){
        cnt[i] = cnt[i-1];
        if(check(i)) pm[len++] = i, cnt[i]++;
    }
    for(ll n = 0; n < N; ++n){
        for(int m = 0; m < M; ++m){
            if(!n){ tmp[n][m] = m; continue; }
            tmp[n][m] = tmp[n-1][m]-tmp[n-1][m/pm[n - 1]];
        }
    }
}
ll euler(ll m, int n){
    if(n == 0) return m;
    if(pm[n - 1] >= m) return 1;
    if(m < M && n < N) return tmp[n][m];

    return euler(m, n - 1) - euler(m / pm[n - 1], n - 1);
}
ll solve(ll m){
    if(m < UP) return cnt[m];
    int s = sqrt(m+0.9);
    int y = pow(m+0.9,1.0/3.0);
    ll res = euler(m, cnt[y])+cnt[y]-1;
    for (int i = cnt[y]; pm[i] <= s; i++){
        res += - solve(m/pm[i]) + solve(pm[i]) - 1;
    }
    return res;
}
int main() {
    init();
    ll n;
    while(scanf("%lld",&n)==1)printf("%lld\n",solve(n));
    return 0;
}