cpp/math/meissel-lehhmer.cpp

• View this file on GitHub
• Last update: 2024-01-24 16:14:34+09:00

Depends on

• cpp/data-structure/binary-indexed-tree.cpp
• cpp/math/exactsqrt.cpp
• cpp/math/sieve-of-eratosthenes.cpp

Required by

• cpp/geometry/two-points-distance.cpp

Verified with

• cpp/z_test/yosupo-counting_primes.cpp

Code

#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;
  }
};

#line 1 "cpp/data-structure/binary-indexed-tree.cpp"
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];
    return ret;
  }

  void add(int a,T w){
    for(int i = a; i <= n; i += i & -i)bit[i] += w;
  }

  T query(int r,int l){
    return sum(l-1)-sum(r-1);
  }


  int lower_bound(T x){
    if(x <= 0){
      return 0;
    }
    x--;
    int t = 1;
    while(t < n)t <<= 1;
    int st = 0;
    int base = 0;
    for(; t; t/=2){
      if(st+t <= n && base+bit[st+t] <= x){
        base += bit[st+t];
        st += t;
      }
    }
    return st+1;
  }
};
#line 1 "cpp/math/sieve-of-eratosthenes.cpp"
struct SieveEratos{
  int N;
  vector<int> minp;
  vector<bool> t;
  vector<int> primes;
  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++){
      if(t[i]){
        primes.emplace_back(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;
    }

  }
  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;}
};
#line 1 "cpp/math/exactsqrt.cpp"
// 整数型を取得することを想定
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;
}
#line 4 "cpp/math/meissel-lehhmer.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;
  }
};

Back to top page