Merge pull request #74 from luckylat/Meissel-Lehmer/implement

feat(meissel-lehmer): implement

cpp/data-structure/binary-indexed-tree.cpp

@@ -2,8 +2,14 @@ template<typename T>
 struct BIT{//1_Indexed
   int n;
   vector<T> bit;
+  BIT(){}
   BIT(int n_):n(n_+1),bit(n,0){}
 
+  void set(int n_){
+    n = n_;
+    bit.assign(n,0);
+  }
+
   T sum(int a){
     T ret = 0;
     for(int i = a; i > 0; i -= i & -i) ret += bit[i];

cpp/math/exactsqrt.cpp

@@ -0,0 +1,41 @@
+// 整数型を取得することを想定
+template <typename T>
+T CeilExactSqrt(int c, T x){
+  T mn = 0;
+  T mx = x;
+  while(mx-mn > 1){
+    T ce = (mn+mx)/2;
+    T nw = 1;
+    for(int i = 0; c > i; i++){
+      if(x/ce < nw){
+        mx = ce;
+        break;
+      }
+      nw *= ce;
+      if(i+1 == c && x == nw){
+        mx = ce;
+      }
+    }
+    if(mx != ce)mn = ce;
+  }
+  return mx;
+}
+
+template <typename T>
+T FloorExactSqrt(int c, T x){
+  T mn = 0;
+  T mx = x;
+  while(mx-mn > 1){
+    T ce = (mn+mx)/2;
+    T nw = 1;
+    for(int i = 0; c > i; i++){
+      if(x/ce < nw){
+        mx = ce;
+        break;
+      }
+      nw *= ce;
+    }
+    if(mx != ce)mn = ce;
+  }
+  return mn;
+}

cpp/math/meissel-lehhmer.cpp

@@ -0,0 +1,84 @@
+#include "../data-structure/binary-indexed-tree.cpp"
+#include "sieve-of-eratosthenes.cpp"
+#include "exactsqrt.cpp"
+
+struct MeisselLehmer{
+  long long N;
+  long long alpha;
+  SieveEratos sieve;
+  BIT<int> varphiTable;
+  struct d{
+    long long m,b;
+    int cof;
+    bool operator<(const d x)const{
+      if(m == x.m) return b < x.b;
+      return m < x.m;
+    }
+  };
+  vector<d> varphiQueue; 
+  MeisselLehmer(long long n): N(n){
+    alpha = ceil(pow(N, 0.34));
+    sieve.set(N/alpha+10);
+    sieve.generate_invprime();
+    sieve.generate_minfactor();
+  }
+  long long varphiLoop(long long n, long long a, int cof = 1){
+    if(a == 0){
+      return n;
+    }else{
+      if(n/sieve.primes[a-1] <= N/alpha){
+        if(n/sieve.primes[a-1] >= 2){
+          varphiQueue.push_back({n/sieve.primes[a-1], a-1, cof*-1});
+          return varphiLoop(n, a-1, cof);
+        }else{
+          return varphiLoop(n, a-1, cof) - varphiLoop(n/sieve.primes[a-1], a-1, cof*-1);
+        } 
+      }else{
+        return varphiLoop(n, a-1, cof) - varphiLoop(n/sieve.primes[a-1], a-1, cof*-1);
+      }
+    }
+  }
+
+  long long varphi(long long n, long long a){
+    varphiTable.set(N/alpha+1);
+    long long val = varphiLoop(n,a);
+    sort(varphiQueue.begin(), varphiQueue.end());
+    int cur = 0;
+    for(int i = 2; N/alpha >= i; i++){
+      varphiTable.add(sieve.invprimes[sieve.minp[i]], 1);
+      while(cur < varphiQueue.size() && varphiQueue[cur].m == i){
+        val += varphiQueue[cur].cof * (varphiQueue[cur].m - varphiTable.query(1, varphiQueue[cur].b+1));
+        cur++;
+      }
+    }
+    varphiQueue.clear();
+    return val;
+  }
+
+  long long P2(long long n, long long a){
+    long long val = 0;
+    int cur = sieve.primes.size()-1;
+    for(int i = a; n/sieve.primes[a-1] > sieve.primes[i]; i++){
+      while((sieve.primes[cur] > n / sieve.primes[i] || (n%sieve.primes[i] == 0 && sieve.primes[cur] == n/sieve.primes[i]))){
+        cur--;
+      }
+      if(cur < i)break;
+      val += cur-i+1;
+    }
+    return val;
+  }
+
+
+  long long count(long long n = -1){
+    if(n == -1)n = N;
+    long long prevN = N;
+    N = n;
+    alpha = ceil(pow(N, 0.34));
+    if(N < 2)return 0;
+    else if(N < 3)return 1;
+    long long val = varphi(N, alpha) + alpha-1 - P2(N, alpha);
+    N = prevN;
+    return val;
+  }
+};
+

cpp/math/sieve-of-eratosthenes.cpp

@@ -1,17 +1,46 @@
 struct SieveEratos{
+  int N;
+  vector<int> minp;
   vector<bool> t;
   vector<int> primes;
-  SieveEratos(int n):t(n+1,true){
+  map<int,int> invprimes;
+  SieveEratos(){}
+  SieveEratos(int n):N(n+1){
+    generate();
+  }
+  void set(int n){
+    N = n+1;
+    generate();
+  }
+  void generate(){
+    t.assign(N, true);
     t[0] = t[1] = false;
-    for(int i = 2; n >= i; i++){
+    for(int i = 2; N > i; i++){
       if(t[i]){
         primes.emplace_back(i);
-        for(int j = i+i; n >= j; j+=i){
+        for(int j = i+i; N > j; j+=i){
           t[j] = false;
         }
       }
+    } 
+  }
+  void generate_invprime(){
+    for(size_t i = 0; primes.size() > i; i++){
+      invprimes[primes[i]] = i+1;
     }
+
   }
-  bool operator[](int x){return t[x];}
+  void generate_minfactor(){
+    minp.assign(N, N+2);
+    minp[0] = minp[1] = -1;
+    for(int i = 2; N > i; i++){
+      if(minp[i] == N+2){
+        minp[i] = i;
+        for(int j = i+i; N > j; j+=i){
+          minp[j] = min(i, minp[j]);
+        }
+      }
+    }
+  }
+  bool operator[](int x){return t[x] == x;}
 };
-

cpp/z_test/yosupo-counting_primes.cpp

@@ -0,0 +1,12 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/counting_primes"
+
+#include "../../cpp/template/template.cpp"
+
+#include "../../cpp/math/meissel-lehhmer.cpp"
+
+
+int main(){
+  long long n;;cin>>n;
+  MeisselLehmer P(n);
+  cout << P.count() << endl;
+}