快速求素数个数

最新推荐文章于 2018-11-09 22:24:09 发布

Drin_E

最新推荐文章于 2018-11-09 22:24:09 发布

阅读量2.2k

点赞数 1

分类专栏：数论

本文链接：https://blog.csdn.net/Drin_E/article/details/52015888

版权

数论专栏收录该内容

16 篇文章 0 订阅

订阅专栏

问题：求1~n出现了多少个素数。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
#define ll long long
using namespace std;
const int maxn=10000000;
int i,j;
ll n,m,c[maxn/10],f[maxn+5];
bool bz[maxn+5];
ll dfs(ll n,int m){
    if (m==0) return n;
    if (m==1) return n-n/2; 
    if (n<=c[m]) return 1;
     if(n<maxn){
        if(f[n]<=m) return 1;
        if(m>=f[(int)sqrt(n)]){
            return f[n]+1-m;
        }
    }
    return dfs(n,m-1)-dfs(n/c[m],m-1);
}
int main(){
    fo(i,2,maxn){
        f[i]=f[i-1];
        if (!bz[i]) c[++c[0]]=i,f[i]++;
        fo(j,1,c[0]){
            if (i>maxn/c[j]) break;
            bz[i*c[j]]=1;
            if (i%c[j]==0) break;
        }
    }
    scanf("%lld",&n);
    int m=sqrt(n);
    printf("%lld\n",dfs(n,f[m])+f[m]-1);
}

根据
Φ(m,n)=Φ(m,n-1)-Φ([m/p(i)],n-1)

给定一个自然数m,若n=Φ(m^(1/3))且μ=Φ(m^(1/2))-n,那么：

π(m)=Φ(m,n)+n(μ+1)+(μ^2-μ)/2-1-ΣΦ(m/p(n+k)) ,k={1,…,μ}.
我们可以得到一个类似的做法(并不会更高效)

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
#define ll long long
using namespace std;
const int maxn=5000000;
int i,j;
ll n,m,c[maxn/10],f[maxn+5];
bool bz[maxn+5];
ll b[maxn+5],h[maxn+5],l,r,t;
int hash(ll n){
    ll y=n;if (y>maxn) y=y%maxn;
    while (b[y]&&(b[y]!=n)) {
        y++;
        if (y>maxn) y=y%maxn;
    }
    return y;
}
ll dfs(ll n,int m){
    if (m==0) return n;
    if (m==1) return n-n/2; 
    if (n<=c[m]) return 1;
     if(n<maxn){
        if(f[n]<=m) return 1;
        if(m>=f[(int)sqrt(n)]){
            return f[n]+1-m;
        }
    }
    return dfs(n,m-1)-dfs(n/c[m],m-1);
}
ll pi(ll x){
    if (x<=maxn) return f[x];
    int i,j=hash(x);
    if (b[j]) return h[j];
    b[j]=x;
    ll n=pi(exp(log(x)/3)),mu=pi(sqrt(x))-n;
    ll s=dfs(x,n)+n*(mu+1)+(mu*mu-mu)/2-1;
    fo(i,1,mu) s=s-pi(x/c[n+i]);
    h[j]=s;
    return s;
}
int main(){
    fo(i,2,maxn){
        f[i]=f[i-1];
        if (!bz[i]) c[++c[0]]=i,f[i]++;
        fo(j,1,c[0]){
            if (i>maxn/c[j]) break;
            bz[i*c[j]]=1;
            if (i%c[j]==0) break;
        }
    }
    scanf("%lld",&n);
    printf("%lld\n",pi(n));
}