51nod 1184 第n个素数

这题大概有两种解法：

1.预先打表放到数组里，大概100万间隔的素数的值$p_n$，然后用区间筛法。可以分的再细一点。

2.先估计，再用区间筛法。你需要一个第n个素数的估计值，一个素数统计函数和区间筛法，根据wiki上的关于第n个素数的公式，通过打表，我发现他有时大于实际值，有时小于实际值，我觉得让估计值保持小于实际值比较方便，这样区间筛的时候，只要向上筛就好了。经过修正我预估的第n个素数

$ p_n=n (n \log +\log  (n \log )-1)+\frac{n (\log  (n \log )-2)}{\log }-6*n/1000 $

我加入了一个向左的移动，通过打表验证，它一定会比实际的$p_n$小，最大误差也在3e6左右。求出预估值以内有多少素数$\pi \left(p_n\right)$,记下与目标的差值，然后向上区间筛。

素数计数函数$\pi \left(n\right)$的算法来源于project euler problem 10的Forum中第5页Lucy_Hedgehog的python代码，原理基本与wiki百科中的方法类似，虽然problem10是求素数的和，但简单修改2处之后就会变成素数统计函数，下文代码中会有注释。区间筛事先预处理出16万以内的素数就够了。

#include<iostream>
#include<cmath>
#include<assert.h>
using namespace std;
typedef long long LINT;
LINT a,b,goal,n;
int mark[160000],prime[160000],e,bl[10000005];
LINT pn(int n)
{
    LINT s=LINT(n*(log(n)+log(log(n))-1)+n*(log(log(n))-2)/log(n)-6.0*n/1000.0);
    return s<=1?1:s;
}

inline LINT V2IDX(LINT v, LINT N, LINT Ndr, LINT nv) {
    return v >= Ndr ? (N/v - 1) : (nv - v);
}

LINT primesum(LINT N) {
    LINT *S;
    LINT *V;

    LINT r = (LINT)sqrt(N);
    LINT Ndr = N/r;

    assert(r*r <= N and (r+1)*(r+1) > N);

    LINT nv = r + Ndr - 1;

    V = new LINT[nv];
    S = new LINT[nv];

    for (LINT i=0; i<r; i++) {
        V[i] = N/(i+1);
    }
    for (LINT i=r; i<nv; i++) {
        V[i] = V[i-1] - 1;

    }

    for (LINT i=0; i<nv; i++) {
      //S[i] = V[i] * (V[i] + 1) / 2 - 1;若求素数和，使用此处
        S[i]=V[i] - 1;
      //若求素数个数使用此处
    }

    for (LINT p=2; p<=r; p++) {
        if (S[nv-p] > S[nv-p+1]) {
            LINT sp = S[nv-p+1];
            LINT p2 = p*p;
            for (LINT i=0; i<nv; i++) {
                if (V[i] >= p2) {
                //S[i] -= p* (S[V2IDX(V[i]/p, N, Ndr, nv)] - sp);若求素数和，使用此处
                    S[i] -= 1* (S[V2IDX(V[i]/p, N, Ndr, nv)] - sp);
                  //若求素数个数，使用此处
                } else {
                    break;
                }
            }
        }
    }

    return S[0];
}

int main()
{
    cin>>n;
    a=pn(n);
    if(a%2&&a>1)a=a-1;//防止预估值本身就是素数的情况，刚开始被这里坑了3个test
    b=a+7000000;
    goal=n-primesum(a);
    for(int i=2;i<=160000;i++)
    {
        if(!mark[i])
        {
            prime[++e]=i;
            for(LINT j=(LINT)i*i;j<=160000;j+=i)
            {
                mark[j]=1;
            }
        }
    }
    LINT xxx,c;
    for(int i=1;i<=e;i++)
    {
        xxx=(LINT)ceil(1.0*a/prime[i]);
        if(xxx==1) xxx++;
        for(LINT j=xxx;(c=j*prime[i])<b;j++)
        {
            bl[c-a]=1;
        }
    }
    int ans=0;
    c=b-a;
    if(a==1) ans--;
    for(int i=0;i<c;i++)
    {
        if(!bl[i]) ans++;
        if(ans==goal)
        {
            cout<<i+a<<endl;
            break;
        }
    }
    return 0;
}

View Code

 

 

转载于:https://www.cnblogs.com/czp001/p/6338045.html