• HDU 5901 Count primes( Meisell-Lehmer算法模板 )



    求1-1e11 以内的素数。

    原理以后再搞,先把板子弄上:
     1 #include<cstdio>
     2 #include<cmath>
     3 
     4 using namespace std;
     5 
     6 #define LL long long
     7 const int N = 5e6 + 2;
     8 bool np[N];
     9 int prime[N], pi[N];
    10 
    11 int getprime(){
    12     int cnt = 0;
    13     np[0] = np[1] = true;
    14     pi[0] = pi[1] = 0;
    15     for(int i = 2; i < N; ++i){
    16     if(!np[i]) prime[++cnt] = i;
    17     pi[i] = cnt;
    18     for(int j = 1; j <= cnt && i * prime[j] < N; ++j){
    19         np[i * prime[j]] = true;
    20         if(i % prime[j] == 0)   break;
    21     }
    22     }
    23     return cnt;
    24 }
    25 
    26 const int M = 7;
    27 const int PM = 2 * 3 * 5 * 7 * 11 * 13 * 17;
    28 int phi[PM + 1][M + 1], sz[M + 1];
    29 
    30 void init(){
    31     getprime();
    32     sz[0] = 1;
    33     for(int i = 0; i <= PM; ++i)  phi[i][0] = i;
    34     for(int i = 1; i <= M; ++i){
    35     sz[i] = prime[i] * sz[i - 1];
    36     for(int j = 1; j <= PM; ++j) phi[j][i] = phi[j][i - 1] - phi[j / prime[i]][i - 1];
    37     }
    38 }
    39 
    40 int sqrt2(LL x){
    41     LL r = (LL)sqrt(x - 0.1);
    42     while(r * r <= x)   ++r;
    43     return int(r - 1);
    44 }
    45 
    46 int sqrt3(LL x){
    47     LL r = (LL)cbrt(x - 0.1);
    48     while(r * r * r <= x)   ++r;
    49     return int(r - 1);
    50 }
    51 
    52 LL getphi(LL x, int s){
    53     if(s == 0)  return x;
    54     if(s <= M)  return phi[x % sz[s]][s] + (x / sz[s]) * phi[sz[s]][s];
    55     if(x <= prime[s]*prime[s])   return pi[x] - s + 1;
    56     if(x <= prime[s]*prime[s]*prime[s] && x < N){
    57     int s2x = pi[sqrt2(x)];
    58     LL ans = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2;
    59     for(int i = s + 1; i <= s2x; ++i) ans += pi[x / prime[i]];
    60     return ans;
    61     }
    62     return getphi(x, s - 1) - getphi(x / prime[s], s - 1);
    63 }
    64 
    65 LL getpi(LL x){
    66     if(x < N)   return pi[x];
    67     LL ans = getphi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1;
    68     for(int i = pi[sqrt3(x)] + 1, ed = pi[sqrt2(x)]; i <= ed; ++i) ans -= getpi(x / prime[i]) - i + 1;
    69     return ans;
    70 }
    71 
    72 LL lehmer_pi(LL x){
    73     if(x < N)   return pi[x];
    74     int a = (int)lehmer_pi(sqrt2(sqrt2(x)));
    75     int b = (int)lehmer_pi(sqrt2(x));
    76     int c = (int)lehmer_pi(sqrt3(x));
    77     LL sum = getphi(x, a) +(LL)(b + a - 2) * (b - a + 1) / 2;
    78     for (int i = a + 1; i <= b; i++){
    79     LL w = x / prime[i];
    80     sum -= lehmer_pi(w);
    81     if (i > c) continue;
    82     LL lim = lehmer_pi(sqrt2(w));
    83     for (int j = i; j <= lim; j++) sum -= lehmer_pi(w / prime[j]) - (j - 1);
    84     }
    85     return sum;
    86 }
    87 
    88 int main(){
    89     init();
    90     LL n;
    91     while(~scanf("%lld",&n)){
    92     printf("%lld
    ",lehmer_pi(n));
    93     }
    94     return 0;
    95 }







  • 相关阅读:
    JavaEE学习之Spring Security3.x——模拟数据库实现用户,权限,资源的管理
    Couchbase入门——环境搭建以及HelloWorld
    VS2010如何使用Visual Studio Online在线服务管理团队资源(在线TFS)
    JSP报错:The superclass "javax.servlet.http.HttpServlet" was not found on the Java Build Path
    JavaEE学习之JAXB
    JavaEE学习之JPA中配置文件persistence.xml
    女朋友手速太慢,导致我无精打采。
    手机投屏之使命召唤
    智能灯改造计划
    女朋友老是埋怨我技术不行,于是我做了个辅助工具。
  • 原文地址:https://www.cnblogs.com/zhangbuang/p/10526120.html
Copyright © 2020-2023  润新知