fixed_discrete_log.hpp

SECTIONMath INCLUDEnoya/fixed_discrete_log.hpp

对固定底数和模数预处理，随后回答多次离散对数查询。

\[ \displaystyle a^x\equiv b\pmod p,\quad x=\min\{t\ge 0:a^t\equiv b\pmod p\} \]

Complexity: Time: O(p^(2/3) + sqrt(ppi(sqrt p)) log p) preprocessing and O(1) per logarithm query for prime p. Space: O(p^(2/3) + sqrt(ppi(sqrt p))).

AC 记录：discrete_logarithm_fixed_mod。

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

当前头文件，省略 include guard；依赖见 #include。

/// @complexity Time: O(p^(2/3) + sqrt(p*pi(sqrt p)) log p)
/// preprocessing and O(1) per logarithm query for prime p.
/// Space: O(p^(2/3) + sqrt(p*pi(sqrt p))).

#include "noya/factorize.hpp"

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <limits>
#include <utility>
#include <vector>

namespace noya {

/// @brief Preprocess discrete logarithms to one primitive root modulo a prime.
/// Farey approximation writes every nonzero x as x*q = +/-t (mod p) with
/// q <= p^(1/3) and |t| <= p^(2/3). Logs of the small integers are built from
/// a batched baby-step/giant-step pass for primes, multiplicativity for small
/// composites, and p = floor(p/i)*i + (p mod i) for the remaining interval.
/// A query then needs only the precomputed logs of q and |t|.
class fixed_discrete_log_table {
  using u32 = std::uint32_t;
  using u64 = std::uint64_t;
  using fraction = std::pair<u32, u32>;

public:
  fixed_discrete_log_table(u32 p, u32 pr) : p_(p), g_(pr), od_(p - 1) {
    assert(p >= 2 && is_prime(p));
    assert(g_ > 0 && g_ < p_);
    if (p_ == 2) {
      dl_ = {std::numeric_limits<u32>::max(), 0};
    } else if (p_ <= dl) {
      build_direct();
    } else {
      build_fast();
    }
  }

  u32 logarithm(u32 val) const {
    assert(val > 0 && val < p_);
    if (!dl_.empty()) {
      return dl_[val];
    }

    u32 idx = u32(u64(val) * fs / p_);
    auto [num, den] = pr_[idx];
    std::int64_t dif =
        std::int64_t(u64(val) * den) - std::int64_t(u64(p_) * num);
    if (std::uint64_t(std::abs(dif)) > fs) {
      num = nx_[idx].first;
      den = nx_[idx].second;
      dif = std::int64_t(u64(val) * den) - std::int64_t(u64(p_) * num);
    }
    assert(dif != 0 && std::uint64_t(std::abs(dif)) <= fs);
    u32 nl = sl_[std::size_t(std::abs(dif))];
    if (dif < 0) {
      nl = add_exponents(nl, od_ / 2);
    }
    return subtract_exponents(nl, sl_[den]);
  }

private:
  static constexpr u32 dl = 1'000'000;

  u32 p_;
  u32 g_;
  u32 od_;
  u32 fs = 0;
  std::vector<u32> dl_;
  std::vector<u32> sl_;
  std::vector<fraction> pr_;
  std::vector<fraction> nx_;

  u32 multiply(u32 a, u32 b) const { return u32(u64(a) * b % p_); }

  u32 power(u32 val, u64 exp) const {
    return u32(factorize_internal::power_mod(val, exp, p_));
  }

  u32 add_exponents(u32 a, u32 b) const {
    u32 res = a + b;
    return res >= od_ ? res - od_ : res;
  }

  u32 subtract_exponents(u32 a, u32 b) const {
    return a >= b ? a - b : a + od_ - b;
  }

  void build_direct() {
    dl_.assign(p_, std::numeric_limits<u32>::max());
    u32 val = 1;
    for (u32 exp = 0; exp < od_; exp++) {
      assert(dl_[val] == std::numeric_limits<u32>::max());
      dl_[val] = exp;
      val = multiply(val, g_);
    }
    assert(val == 1);
  }

  std::vector<u32> smallest_prime_factors(u32 lim, std::vector<u32> &ps) const {
    std::vector<u32> mn(lim + 1);
    for (u32 val = 2; val <= lim; val++) {
      if (mn[val] == 0) {
        mn[val] = val;
        ps.push_back(val);
      }
      for (u32 p : ps) {
        if (p > mn[val] || u64(val) * p > lim) {
          break;
        }
        mn[val * p] = p;
      }
    }
    return mn;
  }

  std::vector<u32> batch_prime_logs(const std::vector<u32> &tar) const {
    if (tar.empty()) {
      return {};
    }
    u32 blk = u32(std::sqrt(static_cast<long double>(p_) / tar.size())) + 2;
    u32 ng = p_ / blk + 3;
    u32 gs = power(g_, blk);
    std::vector<std::pair<u32, u32>> gia;
    gia.reserve(ng);
    u32 val = gs;
    for (u32 x = 1; x <= ng; x++) {
      gia.emplace_back(val, x);
      val = multiply(val, gs);
    }
    std::sort(gia.begin(), gia.end());

    std::vector<u32> ans(tar.size(), std::numeric_limits<u32>::max());
    u32 bab = 1;
    for (u32 y = 0; y < blk; y++) {
      for (std::size_t idx = 0; idx < tar.size(); idx++) {
        u32 wan = multiply(tar[idx], bab);
        auto it = std::lower_bound(gia.begin(), gia.end(),
                                   std::pair<u32, u32>{wan, 0});
        if (it != gia.end() && it->first == wan) {
          u64 can = u64(it->second) * blk - y;
          if (can < ans[idx]) {
            ans[idx] = u32(can);
          }
        }
      }
      bab = multiply(bab, g_);
    }
    for (u32 exp : ans) {
      assert(exp < od_);
    }
    return ans;
  }

  void build_fast() {
    u32 fb = 1;
    while (u64(fb) * fb * fb <= p_) {
      fb *= 2;
    }
    fs = fb * fb;

    std::vector<fraction> exa(fs + 1);
    for (u32 num = 0; num <= fb; num++) {
      u32 fd = num == 1 ? 1 : num + 1;
      for (u32 den = fd; den <= fb; den++) {
        u32 idx = u32(u64(num) * fs / den);
        if (exa[idx].second == 0) {
          exa[idx] = {num, den};
        }
      }
    }
    pr_.resize(fs + 1);
    fraction cur{0, 1};
    for (u32 idx = 0; idx <= fs; idx++) {
      if (exa[idx].second != 0) {
        cur = exa[idx];
      }
      pr_[idx] = cur;
    }
    nx_.resize(fs + 1);
    cur = {1, 1};
    for (u32 idx = fs;; idx--) {
      if (exa[idx].second != 0) {
        cur = exa[idx];
      }
      nx_[idx] = cur;
      if (idx == 0) {
        break;
      }
    }

    u32 sqr = u32(std::sqrt(static_cast<long double>(p_)));
    while (u64(sqr) * sqr > p_) {
      sqr--;
    }
    while (u64(sqr + 1) * (sqr + 1) <= p_) {
      sqr++;
    }
    std::vector<u32> ps;
    std::vector<u32> mn = smallest_prime_factors(sqr, ps);
    std::vector<u32> pl = batch_prime_logs(ps);

    sl_.assign(fs + 1, 0);
    for (std::size_t idx = 0; idx < ps.size(); idx++) {
      sl_[ps[idx]] = pl[idx];
    }
    for (u32 val = 2; val <= sqr; val++) {
      if (mn[val] != val) {
        sl_[val] = add_exponents(sl_[mn[val]], sl_[val / mn[val]]);
      }
    }
    for (u32 val = sqr + 1; val <= fs; val++) {
      u32 quo = p_ / val;
      u32 rem = p_ % val;
      sl_[val] = subtract_exponents(add_exponents(od_ / 2, sl_[rem]), sl_[quo]);
    }
  }
};

} // namespace noya

#ifndef NOYA_FIXED_DISCRETE_LOG_HPP
#define NOYA_FIXED_DISCRETE_LOG_HPP 1

/// @complexity Time: O(p^(2/3) + sqrt(p*pi(sqrt p)) log p)
/// preprocessing and O(1) per logarithm query for prime p.
/// Space: O(p^(2/3) + sqrt(p*pi(sqrt p))).

#include "noya/factorize.hpp"

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <limits>
#include <utility>
#include <vector>

namespace noya {

/// @brief Preprocess discrete logarithms to one primitive root modulo a prime.
/// Farey approximation writes every nonzero x as x*q = +/-t (mod p) with
/// q <= p^(1/3) and |t| <= p^(2/3). Logs of the small integers are built from
/// a batched baby-step/giant-step pass for primes, multiplicativity for small
/// composites, and p = floor(p/i)*i + (p mod i) for the remaining interval.
/// A query then needs only the precomputed logs of q and |t|.
class fixed_discrete_log_table {
  using u32 = std::uint32_t;
  using u64 = std::uint64_t;
  using fraction = std::pair<u32, u32>;

public:
  fixed_discrete_log_table(u32 p, u32 pr) : p_(p), g_(pr), od_(p - 1) {
    assert(p >= 2 && is_prime(p));
    assert(g_ > 0 && g_ < p_);
    if (p_ == 2) {
      dl_ = {std::numeric_limits<u32>::max(), 0};
    } else if (p_ <= dl) {
      build_direct();
    } else {
      build_fast();
    }
  }

  u32 logarithm(u32 val) const {
    assert(val > 0 && val < p_);
    if (!dl_.empty()) {
      return dl_[val];
    }

    u32 idx = u32(u64(val) * fs / p_);
    auto [num, den] = pr_[idx];
    std::int64_t dif =
        std::int64_t(u64(val) * den) - std::int64_t(u64(p_) * num);
    if (std::uint64_t(std::abs(dif)) > fs) {
      num = nx_[idx].first;
      den = nx_[idx].second;
      dif = std::int64_t(u64(val) * den) - std::int64_t(u64(p_) * num);
    }
    assert(dif != 0 && std::uint64_t(std::abs(dif)) <= fs);
    u32 nl = sl_[std::size_t(std::abs(dif))];
    if (dif < 0) {
      nl = add_exponents(nl, od_ / 2);
    }
    return subtract_exponents(nl, sl_[den]);
  }

private:
  static constexpr u32 dl = 1'000'000;

  u32 p_;
  u32 g_;
  u32 od_;
  u32 fs = 0;
  std::vector<u32> dl_;
  std::vector<u32> sl_;
  std::vector<fraction> pr_;
  std::vector<fraction> nx_;

  u32 multiply(u32 a, u32 b) const { return u32(u64(a) * b % p_); }

  u32 power(u32 val, u64 exp) const {
    return u32(factorize_internal::power_mod(val, exp, p_));
  }

  u32 add_exponents(u32 a, u32 b) const {
    u32 res = a + b;
    return res >= od_ ? res - od_ : res;
  }

  u32 subtract_exponents(u32 a, u32 b) const {
    return a >= b ? a - b : a + od_ - b;
  }

  void build_direct() {
    dl_.assign(p_, std::numeric_limits<u32>::max());
    u32 val = 1;
    for (u32 exp = 0; exp < od_; exp++) {
      assert(dl_[val] == std::numeric_limits<u32>::max());
      dl_[val] = exp;
      val = multiply(val, g_);
    }
    assert(val == 1);
  }

  std::vector<u32> smallest_prime_factors(u32 lim, std::vector<u32> &ps) const {
    std::vector<u32> mn(lim + 1);
    for (u32 val = 2; val <= lim; val++) {
      if (mn[val] == 0) {
        mn[val] = val;
        ps.push_back(val);
      }
      for (u32 p : ps) {
        if (p > mn[val] || u64(val) * p > lim) {
          break;
        }
        mn[val * p] = p;
      }
    }
    return mn;
  }

  std::vector<u32> batch_prime_logs(const std::vector<u32> &tar) const {
    if (tar.empty()) {
      return {};
    }
    u32 blk = u32(std::sqrt(static_cast<long double>(p_) / tar.size())) + 2;
    u32 ng = p_ / blk + 3;
    u32 gs = power(g_, blk);
    std::vector<std::pair<u32, u32>> gia;
    gia.reserve(ng);
    u32 val = gs;
    for (u32 x = 1; x <= ng; x++) {
      gia.emplace_back(val, x);
      val = multiply(val, gs);
    }
    std::sort(gia.begin(), gia.end());

    std::vector<u32> ans(tar.size(), std::numeric_limits<u32>::max());
    u32 bab = 1;
    for (u32 y = 0; y < blk; y++) {
      for (std::size_t idx = 0; idx < tar.size(); idx++) {
        u32 wan = multiply(tar[idx], bab);
        auto it = std::lower_bound(gia.begin(), gia.end(),
                                   std::pair<u32, u32>{wan, 0});
        if (it != gia.end() && it->first == wan) {
          u64 can = u64(it->second) * blk - y;
          if (can < ans[idx]) {
            ans[idx] = u32(can);
          }
        }
      }
      bab = multiply(bab, g_);
    }
    for (u32 exp : ans) {
      assert(exp < od_);
    }
    return ans;
  }

  void build_fast() {
    u32 fb = 1;
    while (u64(fb) * fb * fb <= p_) {
      fb *= 2;
    }
    fs = fb * fb;

    std::vector<fraction> exa(fs + 1);
    for (u32 num = 0; num <= fb; num++) {
      u32 fd = num == 1 ? 1 : num + 1;
      for (u32 den = fd; den <= fb; den++) {
        u32 idx = u32(u64(num) * fs / den);
        if (exa[idx].second == 0) {
          exa[idx] = {num, den};
        }
      }
    }
    pr_.resize(fs + 1);
    fraction cur{0, 1};
    for (u32 idx = 0; idx <= fs; idx++) {
      if (exa[idx].second != 0) {
        cur = exa[idx];
      }
      pr_[idx] = cur;
    }
    nx_.resize(fs + 1);
    cur = {1, 1};
    for (u32 idx = fs;; idx--) {
      if (exa[idx].second != 0) {
        cur = exa[idx];
      }
      nx_[idx] = cur;
      if (idx == 0) {
        break;
      }
    }

    u32 sqr = u32(std::sqrt(static_cast<long double>(p_)));
    while (u64(sqr) * sqr > p_) {
      sqr--;
    }
    while (u64(sqr + 1) * (sqr + 1) <= p_) {
      sqr++;
    }
    std::vector<u32> ps;
    std::vector<u32> mn = smallest_prime_factors(sqr, ps);
    std::vector<u32> pl = batch_prime_logs(ps);

    sl_.assign(fs + 1, 0);
    for (std::size_t idx = 0; idx < ps.size(); idx++) {
      sl_[ps[idx]] = pl[idx];
    }
    for (u32 val = 2; val <= sqr; val++) {
      if (mn[val] != val) {
        sl_[val] = add_exponents(sl_[mn[val]], sl_[val / mn[val]]);
      }
    }
    for (u32 val = sqr + 1; val <= fs; val++) {
      u32 quo = p_ / val;
      u32 rem = p_ % val;
      sl_[val] = subtract_exponents(add_exponents(od_ / 2, sl_[rem]), sl_[quo]);
    }
  }
};

} // namespace noya

#endif // NOYA_FIXED_DISCRETE_LOG_HPP

#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <limits>
#include <numeric>
#include <utility>
#include <vector>

/// @complexity Time: O(p^(2/3) + sqrt(p*pi(sqrt p)) log p)
/// preprocessing and O(1) per logarithm query for prime p.
/// Space: O(p^(2/3) + sqrt(p*pi(sqrt p))).

/// @complexity Time: O(log^3 n) primality testing; Pollard-rho factorization is expected about O(n^(1/4)).
/// Space: O(log n) recursion and factors.

namespace noya {
namespace factorize_internal {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

inline u64 multiply_mod(u64 a, u64 b, u64 mod) {
  return u64(u128(a) * b % mod);
}

inline u64 power_mod(u64 a, u64 exp, u64 mod) {
  u64 res = 1;
  while (exp > 0) {
    if (exp & 1) {
      res = multiply_mod(res, a, mod);
    }
    a = multiply_mod(a, a, mod);
    exp >>= 1;
  }
  return res;
}

inline bool miller_rabin(u64 n) {
  if (n < 2) {
    return false;
  }
  for (u64 p :
       std::array<u64, 12>{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
    if (n % p == 0) {
      return n == p;
    }
  }
  int shf = __builtin_ctzll(n - 1);
  u64 odd = (n - 1) >> shf;
  for (u64 bas :
       std::array<u64, 7>{2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
    if (bas % n == 0) {
      continue;
    }
    u64 val = power_mod(bas % n, odd, n);
    if (val == 1 || val == n - 1) {
      continue;
    }
    bool cmp = true;
    for (int i = 1; i < shf; i++) {
      val = multiply_mod(val, val, n);
      if (val == n - 1) {
        cmp = false;
        break;
      }
    }
    if (cmp) {
      return false;
    }
  }
  return true;
}

inline u64 splitmix64(u64 &st) {
  u64 z = (st += 0x9e3779b97f4a7c15ULL);
  z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
  z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
  return z ^ (z >> 31);
}

inline u64 pollard_rho(u64 n) {
  if (n % 2 == 0) {
    return 2;
  }
  if (n % 3 == 0) {
    return 3;
  }
  static u64 st = 0x123456789abcdef0ULL;
  while (true) {
    u64 y = splitmix64(st) % (n - 1) + 1;
    u64 c = splitmix64(st) % (n - 1) + 1;
    constexpr u64 blk = 128;
    u64 g = 1;
    u64 r = 1;
    u64 q = 1;
    u64 x = 0;
    u64 sy = 0;
    auto nxt = [&](u64 val) {
      return u64((u128(multiply_mod(val, val, n)) + c) % n);
    };
    while (g == 1) {
      x = y;
      for (u64 i = 0; i < r; i++) {
        y = nxt(y);
      }
      for (u64 off = 0; off < r && g == 1; off += blk) {
        sy = y;
        for (u64 i = 0; i < std::min(blk, r - off); i++) {
          y = nxt(y);
          u64 dif = x > y ? x - y : y - x;
          q = multiply_mod(q, dif, n);
        }
        g = std::gcd(q, n);
      }
      r <<= 1;
    }
    if (g == n) {
      do {
        sy = nxt(sy);
        u64 dif = x > sy ? x - sy : sy - x;
        g = std::gcd(dif, n);
      } while (g == 1);
    }
    if (g != n) {
      return g;
    }
  }
}

inline void collect_factors(u64 n, std::vector<u64> &res) {
  if (n == 1) {
    return;
  }
  if (miller_rabin(n)) {
    res.push_back(n);
    return;
  }
  u64 fct = pollard_rho(n);
  collect_factors(fct, res);
  collect_factors(n / fct, res);
}

} // namespace factorize_internal

/// @brief Deterministic Miller-Rabin primality test for unsigned 64-bit
/// integers.
inline bool is_prime(std::uint64_t n) {
  return factorize_internal::miller_rabin(n);
}

/// @brief Return the prime factors of n with multiplicity in increasing order.
inline std::vector<std::uint64_t> prime_factors(std::uint64_t n) {
  assert(n >= 1);
  std::vector<std::uint64_t> res;
  factorize_internal::collect_factors(n, res);
  std::sort(res.begin(), res.end());
  return res;
}

/// @brief Return the prime factorization of n as (prime, exponent) pairs.
inline std::vector<std::pair<std::uint64_t, int>> factorize(std::uint64_t n) {
  std::vector<std::pair<std::uint64_t, int>> res;
  for (std::uint64_t p : prime_factors(n)) {
    if (res.empty() || res.back().first != p) {
      res.emplace_back(p, 1);
    } else {
      res.back().second++;
    }
  }
  return res;
}

} // namespace noya

namespace noya {

/// @brief Preprocess discrete logarithms to one primitive root modulo a prime.
/// Farey approximation writes every nonzero x as x*q = +/-t (mod p) with
/// q <= p^(1/3) and |t| <= p^(2/3). Logs of the small integers are built from
/// a batched baby-step/giant-step pass for primes, multiplicativity for small
/// composites, and p = floor(p/i)*i + (p mod i) for the remaining interval.
/// A query then needs only the precomputed logs of q and |t|.
class fixed_discrete_log_table {
  using u32 = std::uint32_t;
  using u64 = std::uint64_t;
  using fraction = std::pair<u32, u32>;

public:
  fixed_discrete_log_table(u32 p, u32 pr) : p_(p), g_(pr), od_(p - 1) {
    assert(p >= 2 && is_prime(p));
    assert(g_ > 0 && g_ < p_);
    if (p_ == 2) {
      dl_ = {std::numeric_limits<u32>::max(), 0};
    } else if (p_ <= dl) {
      build_direct();
    } else {
      build_fast();
    }
  }

  u32 logarithm(u32 val) const {
    assert(val > 0 && val < p_);
    if (!dl_.empty()) {
      return dl_[val];
    }

    u32 idx = u32(u64(val) * fs / p_);
    auto [num, den] = pr_[idx];
    std::int64_t dif =
        std::int64_t(u64(val) * den) - std::int64_t(u64(p_) * num);
    if (std::uint64_t(std::abs(dif)) > fs) {
      num = nx_[idx].first;
      den = nx_[idx].second;
      dif = std::int64_t(u64(val) * den) - std::int64_t(u64(p_) * num);
    }
    assert(dif != 0 && std::uint64_t(std::abs(dif)) <= fs);
    u32 nl = sl_[std::size_t(std::abs(dif))];
    if (dif < 0) {
      nl = add_exponents(nl, od_ / 2);
    }
    return subtract_exponents(nl, sl_[den]);
  }

private:
  static constexpr u32 dl = 1'000'000;

  u32 p_;
  u32 g_;
  u32 od_;
  u32 fs = 0;
  std::vector<u32> dl_;
  std::vector<u32> sl_;
  std::vector<fraction> pr_;
  std::vector<fraction> nx_;

  u32 multiply(u32 a, u32 b) const { return u32(u64(a) * b % p_); }

  u32 power(u32 val, u64 exp) const {
    return u32(factorize_internal::power_mod(val, exp, p_));
  }

  u32 add_exponents(u32 a, u32 b) const {
    u32 res = a + b;
    return res >= od_ ? res - od_ : res;
  }

  u32 subtract_exponents(u32 a, u32 b) const {
    return a >= b ? a - b : a + od_ - b;
  }

  void build_direct() {
    dl_.assign(p_, std::numeric_limits<u32>::max());
    u32 val = 1;
    for (u32 exp = 0; exp < od_; exp++) {
      assert(dl_[val] == std::numeric_limits<u32>::max());
      dl_[val] = exp;
      val = multiply(val, g_);
    }
    assert(val == 1);
  }

  std::vector<u32> smallest_prime_factors(u32 lim, std::vector<u32> &ps) const {
    std::vector<u32> mn(lim + 1);
    for (u32 val = 2; val <= lim; val++) {
      if (mn[val] == 0) {
        mn[val] = val;
        ps.push_back(val);
      }
      for (u32 p : ps) {
        if (p > mn[val] || u64(val) * p > lim) {
          break;
        }
        mn[val * p] = p;
      }
    }
    return mn;
  }

  std::vector<u32> batch_prime_logs(const std::vector<u32> &tar) const {
    if (tar.empty()) {
      return {};
    }
    u32 blk = u32(std::sqrt(static_cast<long double>(p_) / tar.size())) + 2;
    u32 ng = p_ / blk + 3;
    u32 gs = power(g_, blk);
    std::vector<std::pair<u32, u32>> gia;
    gia.reserve(ng);
    u32 val = gs;
    for (u32 x = 1; x <= ng; x++) {
      gia.emplace_back(val, x);
      val = multiply(val, gs);
    }
    std::sort(gia.begin(), gia.end());

    std::vector<u32> ans(tar.size(), std::numeric_limits<u32>::max());
    u32 bab = 1;
    for (u32 y = 0; y < blk; y++) {
      for (std::size_t idx = 0; idx < tar.size(); idx++) {
        u32 wan = multiply(tar[idx], bab);
        auto it = std::lower_bound(gia.begin(), gia.end(),
                                   std::pair<u32, u32>{wan, 0});
        if (it != gia.end() && it->first == wan) {
          u64 can = u64(it->second) * blk - y;
          if (can < ans[idx]) {
            ans[idx] = u32(can);
          }
        }
      }
      bab = multiply(bab, g_);
    }
    for (u32 exp : ans) {
      assert(exp < od_);
    }
    return ans;
  }

  void build_fast() {
    u32 fb = 1;
    while (u64(fb) * fb * fb <= p_) {
      fb *= 2;
    }
    fs = fb * fb;

    std::vector<fraction> exa(fs + 1);
    for (u32 num = 0; num <= fb; num++) {
      u32 fd = num == 1 ? 1 : num + 1;
      for (u32 den = fd; den <= fb; den++) {
        u32 idx = u32(u64(num) * fs / den);
        if (exa[idx].second == 0) {
          exa[idx] = {num, den};
        }
      }
    }
    pr_.resize(fs + 1);
    fraction cur{0, 1};
    for (u32 idx = 0; idx <= fs; idx++) {
      if (exa[idx].second != 0) {
        cur = exa[idx];
      }
      pr_[idx] = cur;
    }
    nx_.resize(fs + 1);
    cur = {1, 1};
    for (u32 idx = fs;; idx--) {
      if (exa[idx].second != 0) {
        cur = exa[idx];
      }
      nx_[idx] = cur;
      if (idx == 0) {
        break;
      }
    }

    u32 sqr = u32(std::sqrt(static_cast<long double>(p_)));
    while (u64(sqr) * sqr > p_) {
      sqr--;
    }
    while (u64(sqr + 1) * (sqr + 1) <= p_) {
      sqr++;
    }
    std::vector<u32> ps;
    std::vector<u32> mn = smallest_prime_factors(sqr, ps);
    std::vector<u32> pl = batch_prime_logs(ps);

    sl_.assign(fs + 1, 0);
    for (std::size_t idx = 0; idx < ps.size(); idx++) {
      sl_[ps[idx]] = pl[idx];
    }
    for (u32 val = 2; val <= sqr; val++) {
      if (mn[val] != val) {
        sl_[val] = add_exponents(sl_[mn[val]], sl_[val / mn[val]]);
      }
    }
    for (u32 val = sqr + 1; val <= fs; val++) {
      u32 quo = p_ / val;
      u32 rem = p_ % val;
      sl_[val] = subtract_exponents(add_exponents(od_ / 2, sl_[rem]), sl_[quo]);
    }
  }
};

} // namespace noya