素数计数函数

来源：互联网发布：柏原崇知乎编辑：程序博客网时间：2024/06/07 06:09

一种计算π(n)的组合方法

这里π(n)指：不大于n的素数个数

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

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

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

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

$π (x) = O (x l n x)$

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

$T (n) = O (n \sum p \leq n 1 p)$

可以认为是线性的。

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

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

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

组合方法：

对于一个整数 q，q∈(n‾√,n]

q如果不能被[2,n‾√]的任何整数整除。

则：q是一个素。

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

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

$C (m, n) = ⌊ n ⌋ - \sum 1 \leq i \leq m ⌊ n P i ⌋ + \sum 1 \leq i < j \leq m ⌊ n P i P j ⌋ - \sum 1 \leq i < j < k \leq m ⌊ n P i P j P k ⌋ . . .$

Pi指第i个素数

则：

$C (π (n ‾ ‾ \sqrt), n) = π (n) - π (n ‾ ‾ \sqrt) + 1$

$C (m, n) = C (m - 1, n) - C (m - 1, ⌊ n P m ⌋)$

对于第二个式子。

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

可以被Pm整除的数字。但不能被前m−1个素数整除的数字

对于Pm的所有倍数：

$k P m, 其中 k \in [1, ⌊ n P m ⌋]$

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

$C (m - 1, ⌊ n P m ⌋)$

所以：

$C (m, n) = C (m - 1, n) - C (m - 1, ⌊ n P m ⌋)$

特别的：

$C (0, n) = ⌊ n ⌋$

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

对于 C(m,n)定义：

$D (m, d) = C (m, ⌊ n d ⌋)$

因为：

$⌊ ⌊ n a ⌋ b ⌋ = ⌊ n a b ⌋$

所以：

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

那么，通过D数组与C数组交替递推。并且选择一个恰当的分界B我们计算D[][1...B]和C[][1...⌊nB+1⌋]

有： $d P m \leq B : D [m] [d] = D [m - 1] [d] - D [m - 1] [d P m] d P m > B : D [m] [d] = D [m - 1] [d] - C [m - 1] [⌊ n d P m ⌋]$

$C [m] [j] = C [m - 1] [j] - C [m] [⌊ j P m ⌋]$

在B=⌊n‾‾√⌋时，这样计算的复杂度:O(nlogn)

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

优化1：

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

当Pm+1>n时：

$C (m, n) = 1$

则:P2m>n时：

$C (m, n) = C (m - 1, n) - C (m - 1, ⌊ n P m ⌋) = C (m - 1, n) - 1$

优化2:

当P2m≥⌊(nB+1)14⌋时：

$C (m, j) = π (m a x (P m, j)) - m + 1$

综上。我们分两段递推：

第一段，前：Pm≤⌊(nB+1)14⌋。正常计算

第二段。此时P2m>⌊(nB+1)14⌋

对于第二个优化。预处理出前⌊n‾√⌋个π[].

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

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

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

P2m>⌊nd⌋时,即：d>⌊nP2m⌋

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

可以肯定的是。不对d∈[⌊nP2m⌋+1,B]更新

仅仅维护d∈[1,⌊nP2m⌋]时。

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

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

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

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

所以第二阶段。我们仅仅计算了d∈[1,⌊nP2m⌋]

那么总多时间复杂度：

$T (n) = \sum P \leq (n B + 1) 1 / 4 (O (B) + O (n B)) + \sum (n B + 1) 1 / 4 < P \leq n 1 / 2 O (n p 2)$

在B=n‾‾√时：T(n)=O(n34logn)

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

使用数据结构的话是：T(n)=O(n34)

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

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

#include <algorithm>#include <stdio.h>#include <string.h>#include <cmath>#define MAXN 1111111using 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;}

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

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

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

那么有 $C (i, j) = \sum k = 0 \infty G k (i, j)$

那么： $C (π (j 1 3), j) = \sum k = 0 \infty G k (π (j 1 3), j) = \sum k = 0 2 G k (π (j 1 3), j) + \sum k = 3 \infty G k (π (j 1 3), j)$

明显1： $\sum k = 3 \infty G k (π (j 1 3), j) = 0$

明显2： $G 0 (i, j) = 1 G 1 (i, j) = π (m a x (j, P i)) - i$

对于G2(π(j13),j).它的可选素数是有P>j13

那么对于另一个素数P′>j13,且PP′≤j , P≤P′.

这样的素数P′个数为：

$π (j P) - π (P) + 1$

所以： $G 2 (π (j 1 3), j) = \sum j 1 / 3 < P \leq j (π (j P) - π (P) + 1)$

综上：

$π (j) = C (π (j 1 3), j) - G 0 (π (j 1 3), j) - G 2 (π (j 1 3), j) + π (j 1 3)$

这种方法。感觉前一部还是要打表出n1/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;}

阅读全文

1 0

素数计数函数

一种计算π(n)的组合方法

这里π(n)指：不大于n的素数个数

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

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

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

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

π(x)=O(xln x)

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

T(n)=O(n∑p≤n1p)

可以认为是线性的。

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

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

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

组合方法：

对于一个整数 q，q∈(n‾√,n]

q如果不能被[2,n‾√]的任何整数整除。

则：q是一个素。

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

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

C(m,n)=⌊n⌋−∑1≤i≤m⌊nPi⌋+∑1≤i<j≤m⌊nPiPj⌋−∑1≤i<j<k≤m⌊nPiPjPk⌋...

Pi指第i个素数

则：

C(π(n‾‾√),n)=π(n)−π(n‾‾√)+1

C(m,n)=C(m−1,n)−C(m−1,⌊nPm⌋)

对于第二个式子。

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

可以被Pm整除的数字。但不能被前m−1个素数整除的数字

对于Pm的所有倍数：

kPm, 其中 k∈[1,⌊nPm⌋]

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

C(m−1,⌊nPm⌋)

所以：

C(m,n)=C(m−1,n)−C(m−1,⌊nPm⌋)

特别的：

C(0,n)=⌊n⌋

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

对于 C(m,n)定义：

D(m,d)=C(m,⌊nd⌋)

因为：

⌊⌊na⌋b⌋=⌊nab⌋

所以：

D(m,d)=D(m−1,d)−D(m−1,dPm)

那么，通过D数组与C数组交替递推。并且选择一个恰当的分界B我们计算D[][1...B]和C[][1...⌊nB+1⌋]

有：dPm≤B:D[m][d]=D[m−1][d]−D[m−1][dPm]dPm>B:D[m][d]=D[m−1][d]−C[m−1][⌊ndPm⌋]

C[m][j]=C[m−1][j]−C[m][⌊jPm⌋]

在B=⌊n‾‾√⌋时，这样计算的复杂度:O(nlogn)

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

优化1：

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

当Pm+1>n时：

C(m,n)=1

则:P2m>n时：

C(m,n)=C(m−1,n)−C(m−1,⌊nPm⌋)=C(m−1,n)−1

优化2:

当P2m≥⌊(nB+1)14⌋时：

C(m,j)=π(max(Pm,j))−m+1

综上 。我们分两段递推：

第一段，前：Pm≤⌊(nB+1)14⌋。正常计算

第二段。此时P2m>⌊(nB+1)14⌋

对于第二个优化。预处理出前⌊n‾√⌋个π[].

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

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

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

P2m>⌊nd⌋时,即：d>⌊nP2m⌋

D(m,d)=D(m−1,d)−1

可以肯定的是。不对d∈[⌊nP2m⌋+1,B]更新

仅仅维护d∈[1,⌊nP2m⌋]时。

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

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

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

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

所以第二阶段。我们仅仅计算了d∈[1,⌊nP2m⌋]

那么总多时间复杂度：

T(n)=∑P≤(nB+1)1/4(O(B)+O(nB))+∑(nB+1)1/4<P≤n1/2O(np2)

在B=n‾‾√时：T(n)=O(n34logn)

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

使用数据结构的话是：T(n)=O(n34)

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

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

$π (x) = O (x l n x)$

$T (n) = O (n \sum p \leq n 1 p)$

$C (m, n) = ⌊ n ⌋ - \sum 1 \leq i \leq m ⌊ n P i ⌋ + \sum 1 \leq i < j \leq m ⌊ n P i P j ⌋ - \sum 1 \leq i < j < k \leq m ⌊ n P i P j P k ⌋ . . .$

$C (π (n ‾ ‾ \sqrt), n) = π (n) - π (n ‾ ‾ \sqrt) + 1$

$C (m, n) = C (m - 1, n) - C (m - 1, ⌊ n P m ⌋)$

$k P m, 其中 k \in [1, ⌊ n P m ⌋]$

$C (m - 1, ⌊ n P m ⌋)$

$C (m, n) = C (m - 1, n) - C (m - 1, ⌊ n P m ⌋)$

$C (0, n) = ⌊ n ⌋$

$D (m, d) = C (m, ⌊ n d ⌋)$

$⌊ ⌊ n a ⌋ b ⌋ = ⌊ n a b ⌋$

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

有： $d P m \leq B : D [m] [d] = D [m - 1] [d] - D [m - 1] [d P m] d P m > B : D [m] [d] = D [m - 1] [d] - C [m - 1] [⌊ n d P m ⌋]$

$C [m] [j] = C [m - 1] [j] - C [m] [⌊ j P m ⌋]$

$C (m, n) = 1$

$C (m, n) = C (m - 1, n) - C (m - 1, ⌊ n P m ⌋) = C (m - 1, n) - 1$

$C (m, j) = π (m a x (P m, j)) - m + 1$

综上。我们分两段递推：

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

$T (n) = \sum P \leq (n B + 1) 1 / 4 (O (B) + O (n B)) + \sum (n B + 1) 1 / 4 < P \leq n 1 / 2 O (n p 2)$

那么有 $C (i, j) = \sum k = 0 \infty G k (i, j)$

那么： $C (π (j 1 3), j) = \sum k = 0 \infty G k (π (j 1 3), j) = \sum k = 0 2 G k (π (j 1 3), j) + \sum k = 3 \infty G k (π (j 1 3), j)$

明显1： $\sum k = 3 \infty G k (π (j 1 3), j) = 0$

明显2： $G 0 (i, j) = 1 G 1 (i, j) = π (m a x (j, P i)) - i$

$π (j P) - π (P) + 1$

所以： $G 2 (π (j 1 3), j) = \sum j 1 / 3 < P \leq j (π (j P) - π (P) + 1)$

$π (j) = C (π (j 1 3), j) - G 0 (π (j 1 3), j) - G 2 (π (j 1 3), j) + π (j 1 3)$

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