Skip to content

prime.hpp

SECTIONMath INCLUDEnoya/prime.hpp

A linear sieve providing primes, smallest prime factors, Euler phi, and Mobius values.

\[ \displaystyle \varphi(n)=n\prod_{p\mid n}\left(1-\frac1p\right) \]

Implementation

View on GitHub

#ifndef NOYA_PRIME_HPP
#define NOYA_PRIME_HPP 1

/// @complexity Time: O(n) sieve build, O(log x) table factorization.
/// Space: O(n).

#include <algorithm>
#include <cassert>
#include <utility>
#include <vector>

namespace noya {

/// @brief Linear sieve with smallest prime factors, Euler phi, and Mobius
/// values.
struct prime_table {
  int n = 0;
  std::vector<int> primes;
  std::vector<int> min_factor;
  std::vector<int> phi;
  std::vector<int> mobius;

  prime_table() = default;
  explicit prime_table(int n_) { build(n_); }

  /// @brief Build tables for every integer in [0, n].
  void build(int n_) {
    assert(n_ >= 0);
    n = n_;
    primes.clear();
    min_factor.assign(n + 1, 0);
    phi.assign(n + 1, 0);
    mobius.assign(n + 1, 0);
    if (n >= 1) {
      phi[1] = 1;
      mobius[1] = 1;
    }
    for (int x = 2; x <= n; x++) {
      if (min_factor[x] == 0) {
        min_factor[x] = x;
        phi[x] = x - 1;
        mobius[x] = -1;
        primes.push_back(x);
      }
      for (int p : primes) {
        if (p > n / x) {
          break;
        }
        min_factor[x * p] = p;
        if (x % p == 0) {
          phi[x * p] = phi[x] * p;
          mobius[x * p] = 0;
          break;
        }
        phi[x * p] = phi[x] * (p - 1);
        mobius[x * p] = -mobius[x];
      }
    }
  }

  bool is_prime(int x) const {
    assert(0 <= x && x <= n);
    return x >= 2 && min_factor[x] == x;
  }

  /// @brief Factor x in O(log x); x must lie in [1, n].
  std::vector<std::pair<int, int>> factorize(int x) const {
    assert(1 <= x && x <= n);
    std::vector<std::pair<int, int>> result;
    while (x > 1) {
      int p = min_factor[x];
      int exponent = 0;
      do {
        x /= p;
        exponent++;
      } while (x % p == 0);
      result.emplace_back(p, exponent);
    }
    return result;
  }

  /// @brief Enumerate the positive divisors of x in increasing order.
  std::vector<int> divisors(int x) const {
    std::vector<int> result = {1};
    for (auto [p, exponent] : factorize(x)) {
      int old_size = int(result.size());
      int power = 1;
      for (int e = 1; e <= exponent; e++) {
        power *= p;
        for (int i = 0; i < old_size; i++) {
          result.push_back(result[i] * power);
        }
      }
    }
    std::sort(result.begin(), result.end());
    return result;
  }
};

} // namespace noya

#endif // NOYA_PRIME_HPP
#include <algorithm>
#include <cassert>
#include <utility>
#include <vector>

/// @complexity Time: O(n) sieve build, O(log x) table factorization.
/// Space: O(n).

namespace noya {

/// @brief Linear sieve with smallest prime factors, Euler phi, and Mobius
/// values.
struct prime_table {
  int n = 0;
  std::vector<int> primes;
  std::vector<int> min_factor;
  std::vector<int> phi;
  std::vector<int> mobius;

  prime_table() = default;
  explicit prime_table(int n_) { build(n_); }

  /// @brief Build tables for every integer in [0, n].
  void build(int n_) {
    assert(n_ >= 0);
    n = n_;
    primes.clear();
    min_factor.assign(n + 1, 0);
    phi.assign(n + 1, 0);
    mobius.assign(n + 1, 0);
    if (n >= 1) {
      phi[1] = 1;
      mobius[1] = 1;
    }
    for (int x = 2; x <= n; x++) {
      if (min_factor[x] == 0) {
        min_factor[x] = x;
        phi[x] = x - 1;
        mobius[x] = -1;
        primes.push_back(x);
      }
      for (int p : primes) {
        if (p > n / x) {
          break;
        }
        min_factor[x * p] = p;
        if (x % p == 0) {
          phi[x * p] = phi[x] * p;
          mobius[x * p] = 0;
          break;
        }
        phi[x * p] = phi[x] * (p - 1);
        mobius[x * p] = -mobius[x];
      }
    }
  }

  bool is_prime(int x) const {
    assert(0 <= x && x <= n);
    return x >= 2 && min_factor[x] == x;
  }

  /// @brief Factor x in O(log x); x must lie in [1, n].
  std::vector<std::pair<int, int>> factorize(int x) const {
    assert(1 <= x && x <= n);
    std::vector<std::pair<int, int>> result;
    while (x > 1) {
      int p = min_factor[x];
      int exponent = 0;
      do {
        x /= p;
        exponent++;
      } while (x % p == 0);
      result.emplace_back(p, exponent);
    }
    return result;
  }

  /// @brief Enumerate the positive divisors of x in increasing order.
  std::vector<int> divisors(int x) const {
    std::vector<int> result = {1};
    for (auto [p, exponent] : factorize(x)) {
      int old_size = int(result.size());
      int power = 1;
      for (int e = 1; e <= exponent; e++) {
        power *= p;
        for (int i = 0; i < old_size; i++) {
          result.push_back(result[i] * power);
        }
      }
    }
    std::sort(result.begin(), result.end());
    return result;
  }
};

} // namespace noya