51nod 1184 第n个素数

Posted 2020-08-28

这题大概有两种解法：

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 $

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

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

  1 #include<iostream>
  2 #include<cmath>
  3 #include<assert.h>
  4 using namespace std;
  5 typedef long long LINT;
  6 LINT a,b,goal,n;
  7 int mark[160000],prime[160000],e,bl[10000005];
  8 LINT pn(int n)
  9 {
 10     LINT s=LINT(n*(log(n)+log(log(n))-1)+n*(log(log(n))-2)/log(n)-6.0*n/1000.0);
 11     return s<=1?1:s;
 12 }
 13 
 14 inline LINT V2IDX(LINT v, LINT N, LINT Ndr, LINT nv) {
 15     return v >= Ndr ? (N/v - 1) : (nv - v);
 16 }
 17 
 18 LINT primesum(LINT N) {
 19     LINT *S;
 20     LINT *V;
 21 
 22     LINT r = (LINT)sqrt(N);
 23     LINT Ndr = N/r;
 24 
 25     assert(r*r <= N and (r+1)*(r+1) > N);
 26 
 27     LINT nv = r + Ndr - 1;
 28 
 29     V = new LINT[nv];
 30     S = new LINT[nv];
 31 
 32     for (LINT i=0; i<r; i++) {
 33         V[i] = N/(i+1);
 34     }
 35     for (LINT i=r; i<nv; i++) {
 36         V[i] = V[i-1] - 1;
 37 
 38     }
 39 
 40     for (LINT i=0; i<nv; i++) {
 41       //S[i] = V[i] * (V[i] + 1) / 2 - 1;若求素数和，使用此处
 42         S[i]=V[i] - 1;
 43       //若求素数个数使用此处
 44     }
 45 
 46     for (LINT p=2; p<=r; p++) {
 47         if (S[nv-p] > S[nv-p+1]) {
 48             LINT sp = S[nv-p+1];
 49             LINT p2 = p*p;
 50             for (LINT i=0; i<nv; i++) {
 51                 if (V[i] >= p2) {
 52                 //S[i] -= p* (S[V2IDX(V[i]/p, N, Ndr, nv)] - sp);若求素数和，使用此处
 53                     S[i] -= 1* (S[V2IDX(V[i]/p, N, Ndr, nv)] - sp);
 54                   //若求素数个数，使用此处
 55                 } else {
 56                     break;
 57                 }
 58             }
 59         }
 60     }
 61 
 62     return S[0];
 63 }
 64 
 65 int main()
 66 {
 67     cin>>n;
 68     a=pn(n);
 69     if(a%2&&a>1)a=a-1;//防止预估值本身就是素数的情况，刚开始被这里坑了3个test
 70     b=a+7000000;
 71     goal=n-primesum(a);
 72     for(int i=2;i<=160000;i++)
 73     {
 74         if(!mark[i])
 75         {
 76             prime[++e]=i;
 77             for(LINT j=(LINT)i*i;j<=160000;j+=i)
 78             {
 79                 mark[j]=1;
 80             }
 81         }
 82     }
 83     LINT xxx,c;
 84     for(int i=1;i<=e;i++)
 85     {
 86         xxx=(LINT)ceil(1.0*a/prime[i]);
 87         if(xxx==1) xxx++;
 88         for(LINT j=xxx;(c=j*prime[i])<b;j++)
 89         {
 90             bl[c-a]=1;
 91         }
 92     }
 93     int ans=0;
 94     c=b-a;
 95     if(a==1) ans--;
 96     for(int i=0;i<c;i++)
 97     {
 98         if(!bl[i]) ans++;
 99         if(ans==goal)
100         {
101             cout<<i+a<<endl;
102             break;
103         }
104     }
105     return 0;
106 }

View Code