素数の列挙 (src/NumberTheory/enumerate_primes.hpp)

View this file on GitHub
View document part on GitHub
Last update: 2024-02-19 15:31:52+09:00
Include: #include "src/NumberTheory/enumerate_primes.hpp"

エラトステネス篩による素数列挙の前処理を自動でおこなう．

関数

名前	概要	計算量
`enumerate_primes(N)`	$N$ 以下のすべての素数を返す. 返り値は `ListRange<int>`
`least_primes_factor(n)`	$n$ の最小素因数を返す.
`completely_multiplicative_table<T>(N, f)`	完全乗法的関数 $f$ の値が入ったテーブルを返す. $f$ は素冪のもの `f(p,e)`$=f(p^e)$ を渡す(完全だけど). サイズは $N+1$.	$\displaystyle O\left(N + \frac{N\alpha}{\log N}\right)$ ただし $f(p)$ の計算量を $O(\alpha)$ とおいた.
`pow_table<T>(N,k)`	$n^k$ の値が入ったテーブルを返す. サイズは $N+1$.	$\displaystyle O\left(N + \frac{N\log k}{\log N}\right)$

参考

https://37zigen.com/linear-sieve/

Depends on

CSR 表現を用いた二次元配列他 (src/Internal/ListRange.hpp)

Required by

Verified with

Code

#pragma once
#include <algorithm>
#include <cstdint>
#include "src/Internal/ListRange.hpp"
namespace nt_internal {
using namespace std;
vector<int> ps, lf;
void sieve(int N) {
 static int n= 2;
 if (n > N) return;
 if (lf.resize((N + 1) >> 1); n == 2) ps.push_back(n++);
 int M= (N - 1) / 2;
 for (int j= 1, e= ps.size(); j < e; ++j) {
  int p= ps[j];
  if (int64_t(p) * p > N) break;
  for (auto k= int64_t(p) * max(n / p / 2 * 2 + 1, p) / 2; k <= M; k+= p) lf[k]+= p * !lf[k];
 }
 for (; n <= N; n+= 2)
  if (!lf[n >> 1]) {
   ps.push_back(lf[n >> 1]= n);
   for (auto j= int64_t(n) * n / 2; j <= M; j+= n) lf[j]+= n * !lf[j];
  }
}
ConstListRange<int> enumerate_primes() { return {ps.cbegin(), ps.cend()}; }
ConstListRange<int> enumerate_primes(int N) {
 sieve(N);
 return {ps.cbegin(), upper_bound(ps.cbegin(), ps.cend(), N)};
}
int least_prime_factor(int n) { return n & 1 ? sieve(n), lf[(n >> 1)] : 2; }
// f(p,e) := f(p^e)
template <class T, class F> vector<T> completely_multiplicative_table(int N, const F &f) {
 vector<T> ret(N + 1);
 sieve(N);
 for (int n= 3, i= 1; n <= N; n+= 2, ++i) ret[n]= lf[i] == n ? f(n, 1) : ret[lf[i]] * ret[n / lf[i]];
 if (int n= 4; 2 <= N)
  for (T t= ret[2]= f(2, 1); n <= N; n+= 2) ret[n]= t * ret[n >> 1];
 return ret[1]= 1, ret;
}
}
using nt_internal::enumerate_primes, nt_internal::least_prime_factor, nt_internal::completely_multiplicative_table;
// O(N log k / log N + N)
template <class T> static std::vector<T> pow_table(int N, uint64_t k) {
 if (k == 0) return std::vector<T>(N + 1, 1);
 auto f= [k](int p, int) {
  T ret= 1, b= p;
  for (auto e= k;; b*= b) {
   if (e & 1) ret*= b;
   if (!(e>>= 1)) return ret;
  }
 };
 return completely_multiplicative_table<T>(N, f);
}

#line 2 "src/NumberTheory/enumerate_primes.hpp"
#include <algorithm>
#include <cstdint>
#line 2 "src/Internal/ListRange.hpp"
#include <vector>
#include <iostream>
#include <iterator>
#include <type_traits>
#define _LR(name, IT, CT) \
 template <class T> struct name { \
  using Iterator= typename std::vector<T>::IT; \
  Iterator bg, ed; \
  Iterator begin() const { return bg; } \
  Iterator end() const { return ed; } \
  size_t size() const { return std::distance(bg, ed); } \
  CT &operator[](int i) const { return bg[i]; } \
 }
_LR(ListRange, iterator, T);
_LR(ConstListRange, const_iterator, const T);
#undef _LR
template <class T> struct CSRArray {
 std::vector<T> dat;
 std::vector<int> p;
 size_t size() const { return p.size() - 1; }
 ListRange<T> operator[](int i) { return {dat.begin() + p[i], dat.begin() + p[i + 1]}; }
 ConstListRange<T> operator[](int i) const { return {dat.cbegin() + p[i], dat.cbegin() + p[i + 1]}; }
};
template <template <class> class F, class T> std::enable_if_t<std::disjunction_v<std::is_same<F<T>, ListRange<T>>, std::is_same<F<T>, ConstListRange<T>>, std::is_same<F<T>, CSRArray<T>>>, std::ostream &> operator<<(std::ostream &os, const F<T> &r) {
 os << '[';
 for (int _= 0, __= r.size(); _ < __; ++_) os << (_ ? ", " : "") << r[_];
 return os << ']';
}
#line 5 "src/NumberTheory/enumerate_primes.hpp"
namespace nt_internal {
using namespace std;
vector<int> ps, lf;
void sieve(int N) {
 static int n= 2;
 if (n > N) return;
 if (lf.resize((N + 1) >> 1); n == 2) ps.push_back(n++);
 int M= (N - 1) / 2;
 for (int j= 1, e= ps.size(); j < e; ++j) {
  int p= ps[j];
  if (int64_t(p) * p > N) break;
  for (auto k= int64_t(p) * max(n / p / 2 * 2 + 1, p) / 2; k <= M; k+= p) lf[k]+= p * !lf[k];
 }
 for (; n <= N; n+= 2)
  if (!lf[n >> 1]) {
   ps.push_back(lf[n >> 1]= n);
   for (auto j= int64_t(n) * n / 2; j <= M; j+= n) lf[j]+= n * !lf[j];
  }
}
ConstListRange<int> enumerate_primes() { return {ps.cbegin(), ps.cend()}; }
ConstListRange<int> enumerate_primes(int N) {
 sieve(N);
 return {ps.cbegin(), upper_bound(ps.cbegin(), ps.cend(), N)};
}
int least_prime_factor(int n) { return n & 1 ? sieve(n), lf[(n >> 1)] : 2; }
// f(p,e) := f(p^e)
template <class T, class F> vector<T> completely_multiplicative_table(int N, const F &f) {
 vector<T> ret(N + 1);
 sieve(N);
 for (int n= 3, i= 1; n <= N; n+= 2, ++i) ret[n]= lf[i] == n ? f(n, 1) : ret[lf[i]] * ret[n / lf[i]];
 if (int n= 4; 2 <= N)
  for (T t= ret[2]= f(2, 1); n <= N; n+= 2) ret[n]= t * ret[n >> 1];
 return ret[1]= 1, ret;
}
}
using nt_internal::enumerate_primes, nt_internal::least_prime_factor, nt_internal::completely_multiplicative_table;
// O(N log k / log N + N)
template <class T> static std::vector<T> pow_table(int N, uint64_t k) {
 if (k == 0) return std::vector<T>(N + 1, 1);
 auto f= [k](int p, int) {
  T ret= 1, b= p;
  for (auto e= k;; b*= b) {
   if (e & 1) ret*= b;
   if (!(e>>= 1)) return ret;
  }
 };
 return completely_multiplicative_table<T>(N, f);
}