区間篩
(number-theory/range-sieve.hpp)

View this file on GitHub
Last update: 2026-07-28 20:47:27+09:00
Include: #include "number-theory/range-sieve.hpp"

区間篩

$L$ 以上 $R$ 以下の整数を区間篩により処理する．$1\leq L\leq R$ とする．

RangeSieve::lpf(L, R)：lpf[x - L] が $x$ の最小素因数となる長さ $R-L+1$ の列を返す．ただし，$x=1$ のときは $1$ とする．
RangeSieve::table(L, R)：$L$ 以上 $R$ 以下の素数を昇順に並べた列を返す．
RangeSieve::factorize(L, R)：factors[x - L] が $x$ の素因数と指数の組を昇順に並べた列となる，長さ $R-L+1$ の列を返す．$1$ に対応する列は空である．

素数表の前計算は static に保持し，上限を $2$ 倍ずつ拡張する．一連の呼び出しで必要な $\lfloor\sqrt R\rfloor$ の最大値を $X$ とすると，前計算の合計は $O(X\log\log X)$ 時間，$O(X)$ 空間である．各呼び出しの区間処理に $O((R-L+1)\log\log R)$ 時間を要する．

Depends on

Verified with

verify/number-theory/UNIT_range_sieve_table.test.cpp

Code

#pragma once

#include "math/util.hpp"
#include "number-theory/prime-sieve.hpp"

namespace RangeSieve {
using ll = long long;

namespace internal {
const vector<int>& primes(int limit) {
  static int n = 1;
  static vector<int> ps;
  if (n < limit) {
    while (n < limit) n *= 2;
    ps = PrimeSieve::table(n);
  }
  return ps;
}
};  // namespace internal

// lpf of [l,r]
vector<ll> lpf(ll l, ll r) {
  assert(1 <= l && l <= r);
  int limit = Math::isqrt(r);
  vector<ll> ret(r - l + 1);
  for (size_t i = 0; i < ret.size(); i++) ret[i] = l + (ll)i;
  for (ll p : internal::primes(limit)) {
    if (p > limit) break;
    for (ll x = Math::ceil(l, p) * p; x <= r; x += p)
      if (ret[x - l] > p) ret[x - l] = p;
  }
  return ret;
}

vector<ll> table(ll l, ll r) {
  assert(1 <= l && l <= r);
  int limit = Math::isqrt(r);
  vector<bool> composite(r - l + 1, false);
  for (ll p : internal::primes(limit)) {
    if (p > limit) break;
    for (ll x = max(Math::ceil(l, p), p) * p; x <= r; x += p)
      composite[x - l] = true;
  }
  vector<ll> ps;
  for (size_t i = 0; i < composite.size(); i++) {
    ll x = l + (ll)i;
    if (x >= 2 && !composite[i]) ps.push_back(x);
  }
  return ps;
}

vector<vector<pair<ll, int>>> factorize(ll l, ll r) {
  assert(1 <= l && l <= r);
  int limit = Math::isqrt(r);
  vector<vector<pair<ll, int>>> factors(r - l + 1);
  vector<ll> rem(r - l + 1);
  for (size_t i = 0; i < rem.size(); i++) rem[i] = l + (ll)i;
  for (ll p : internal::primes(limit)) {
    if (p > limit) break;
    for (ll x = Math::ceil(l, p) * p; x <= r; x += p) {
      int e = 0;
      while (rem[x - l] % p == 0) rem[x - l] /= p, e++;
      factors[x - l].emplace_back(p, e);
    }
  }
  for (size_t i = 0; i < rem.size(); i++)
    if (rem[i] > 1) factors[i].emplace_back(rem[i], 1);
  return factors;
}
};  // namespace RangeSieve

/**
 * @brief 区間篩
 * @docs docs/number-theory/range-sieve.md
 */

#line 2 "number-theory/range-sieve.hpp"

#line 2 "math/util.hpp"

namespace Math {
template <class T>
T safe_mod(T a, T b) {
  assert(b != 0);
  if (b < 0) a = -a, b = -b;
  a %= b;
  return a >= 0 ? a : a + b;
}
template <class T>
T floor(T a, T b) {
  assert(b != 0);
  if (b < 0) a = -a, b = -b;
  return a >= 0 ? a / b : (a + 1) / b - 1;
}
template <class T>
T ceil(T a, T b) {
  assert(b != 0);
  if (b < 0) a = -a, b = -b;
  return a > 0 ? (a - 1) / b + 1 : a / b;
}
long long isqrt(long long n) {
  if (n <= 0) return 0;
  long long x = sqrt(n);
  while ((__int128)(x + 1) * (x + 1) <= n) x++;
  while ((__int128)x * x > n) x--;
  return x;
}
long long floor_root(long long n, int k) {
  assert(n >= 0);
  if (n == 0) return 0;
  assert(k >= 1);
  if (k == 1) return n;
  if (k > 64) return 1;
  long long x = round(pow((long double)n, 1.0L / k));
  auto check = [&](long long a) {
    if (a <= 0) return true;
    __int128_t p = 1;
    for (int i = 0; i < k; ++i)
      if ((p *= a) > n) return false;
    return true;
  };
  while (check(x + 1)) x++;
  while (!check(x)) x--;
  return x;
}
unsigned long long floor_root_unsigned(unsigned long long n, int k) {
  assert(k >= 1);
  if (n <= 1 || k == 1) return n;
  if (k >= 64) return 1;
  int bits = (64 + k - 1) / k;
  unsigned long long ok = 1, ng = min(n, 1ULL << bits);
  auto check = [&](unsigned long long a) {
    __uint128_t p = 1;
    for (int i = 0; i < k; i++) {
      p *= a;
      if (p > n) return false;
    }
    return true;
  };
  while (ok + 1 < ng) {
    unsigned long long mid = ok + (ng - ok) / 2;
    (check(mid) ? ok : ng) = mid;
  }
  return ok;
}
// return g=gcd(a,b)
// a*x+b*y=g
// - b!=0 -> 0<=x<|b|/g
// - b=0  -> ax=g
template <class T>
T ext_gcd(T a, T b, T& x, T& y) {
  T a0 = a, b0 = b;
  bool sgn_a = a < 0, sgn_b = b < 0;
  if (sgn_a) a = -a;
  if (sgn_b) b = -b;
  if (b == 0) {
    x = sgn_a ? -1 : 1;
    y = 0;
    return a;
  }
  T x00 = 1, x01 = 0, x10 = 0, x11 = 1;
  while (b != 0) {
    T q = a / b, r = a - b * q;
    x00 -= q * x01;
    x10 -= q * x11;
    swap(x00, x01);
    swap(x10, x11);
    a = b, b = r;
  }
  x = x00, y = x10;
  if (sgn_a) x = -x;
  if (sgn_b) y = -y;
  if (b0 != 0) {
    a0 /= a, b0 /= a;
    if (b0 < 0) a0 = -a0, b0 = -b0;
    T q = x >= 0 ? x / b0 : (x + 1) / b0 - 1;
    x -= b0 * q;
    y += a0 * q;
  }
  return a;
}
constexpr long long inv_mod(long long x, long long m) {
  x %= m;
  if (x < 0) x += m;
  long long a = m, b = x;
  long long y0 = 0, y1 = 1;
  while (b > 0) {
    long long q = a / b;
    swap(a -= q * b, b);
    swap(y0 -= q * y1, y1);
  }
  if (y0 < 0) y0 += m / a;
  return y0;
}
long long pow_mod(long long x, long long n, long long m) {
  if (m == 1) return 0;
  x = (x % m + m) % m;
  long long y = 1;
  while (n) {
    if (n & 1) y = y * x % m;
    x = x * x % m;
    n >>= 1;
  }
  return y;
}
constexpr long long pow_mod_constexpr(long long x, long long n, int m) {
  if (m == 1) return 0;
  unsigned int _m = (unsigned int)(m);
  unsigned long long r = 1;
  unsigned long long y = x % m;
  if (y >= m) y += m;
  while (n) {
    if (n & 1) r = (r * y) % _m;
    y = (y * y) % _m;
    n >>= 1;
  }
  return r;
}
constexpr bool is_prime_constexpr(int n) {
  if (n <= 1) return false;
  if (n == 2 || n == 7 || n == 61) return true;
  if (n % 2 == 0) return false;
  long long d = n - 1;
  while (d % 2 == 0) d /= 2;
  constexpr long long bases[3] = {2, 7, 61};
  for (long long a : bases) {
    long long t = d;
    long long y = pow_mod_constexpr(a, t, n);
    while (t != n - 1 && y != 1 && y != n - 1) {
      y = y * y % n;
      t <<= 1;
    }
    if (y != n - 1 && t % 2 == 0) {
      return false;
    }
  }
  return true;
}
template <int n>
constexpr bool is_prime = is_prime_constexpr(n);
};  // namespace Math
#line 2 "number-theory/prime-sieve.hpp"

namespace PrimeSieve {
using ll = long long;

vector<int> lpf(int n) {
  assert(n >= 0);
  vector<int> ret(n + 1);
  for (size_t i = 0; i < ret.size(); i++) ret[i] = (int)i;
  for (int p = 2; (ll)p * p <= n; p++) {
    if (ret[p] != p) continue;
    for (ll x = (ll)p * p;; x += p) {
      if (ret[x] == x) ret[x] = p;
      if (n - x < p) break;
    }
  }
  return ret;
}

vector<int> table(int n) {
  assert(n >= 0);
  vector<bool> composite(n + 1, false);
  for (int p = 2; (ll)p * p <= n; p += (p & 1) + 1) {
    if (composite[p]) continue;
    for (ll x = (ll)p * p;; x += p) {
      composite[x] = true;
      if (n - x < p) break;
    }
  }
  vector<int> ps;
  for (int p = 2; p <= n;) {
    if (!composite[p]) ps.push_back(p);
    int step = (p & 1) + 1;
    if (n - p < step) break;
    p += step;
  }
  return ps;
}

vector<vector<pair<ll, int>>> factorize(int n) {
  assert(n >= 0);
  vector<vector<pair<ll, int>>> factors(n + 1);
  auto lp = lpf(n);
  for (int x = 2; x <= n;) {
    int y = x;
    while (y > 1) {
      int p = lp[y], e = 0;
      while (y % p == 0) y /= p, e++;
      factors[x].emplace_back(p, e);
    }
    if (x == n) break;
    x++;
  }
  return factors;
}
};  // namespace PrimeSieve

/**
 * @brief 素数篩
 * @docs docs/number-theory/prime-sieve.md
 */
#line 5 "number-theory/range-sieve.hpp"

namespace RangeSieve {
using ll = long long;

namespace internal {
const vector<int>& primes(int limit) {
  static int n = 1;
  static vector<int> ps;
  if (n < limit) {
    while (n < limit) n *= 2;
    ps = PrimeSieve::table(n);
  }
  return ps;
}
};  // namespace internal

// lpf of [l,r]
vector<ll> lpf(ll l, ll r) {
  assert(1 <= l && l <= r);
  int limit = Math::isqrt(r);
  vector<ll> ret(r - l + 1);
  for (size_t i = 0; i < ret.size(); i++) ret[i] = l + (ll)i;
  for (ll p : internal::primes(limit)) {
    if (p > limit) break;
    for (ll x = Math::ceil(l, p) * p; x <= r; x += p)
      if (ret[x - l] > p) ret[x - l] = p;
  }
  return ret;
}

vector<ll> table(ll l, ll r) {
  assert(1 <= l && l <= r);
  int limit = Math::isqrt(r);
  vector<bool> composite(r - l + 1, false);
  for (ll p : internal::primes(limit)) {
    if (p > limit) break;
    for (ll x = max(Math::ceil(l, p), p) * p; x <= r; x += p)
      composite[x - l] = true;
  }
  vector<ll> ps;
  for (size_t i = 0; i < composite.size(); i++) {
    ll x = l + (ll)i;
    if (x >= 2 && !composite[i]) ps.push_back(x);
  }
  return ps;
}

vector<vector<pair<ll, int>>> factorize(ll l, ll r) {
  assert(1 <= l && l <= r);
  int limit = Math::isqrt(r);
  vector<vector<pair<ll, int>>> factors(r - l + 1);
  vector<ll> rem(r - l + 1);
  for (size_t i = 0; i < rem.size(); i++) rem[i] = l + (ll)i;
  for (ll p : internal::primes(limit)) {
    if (p > limit) break;
    for (ll x = Math::ceil(l, p) * p; x <= r; x += p) {
      int e = 0;
      while (rem[x - l] % p == 0) rem[x - l] /= p, e++;
      factors[x - l].emplace_back(p, e);
    }
  }
  for (size_t i = 0; i < rem.size(); i++)
    if (rem[i] > 1) factors[i].emplace_back(rem[i], 1);
  return factors;
}
};  // namespace RangeSieve

/**
 * @brief 区間篩
 * @docs docs/number-theory/range-sieve.md
 */

区間篩 (number-theory/range-sieve.hpp)

区間篩

Depends on

Verified with

Code

区間篩
(number-theory/range-sieve.hpp)