finite_field_polynomial_factorization.hpp

SECTIONMath INCLUDEnoya/finite_field_polynomial_factorization.hpp

在任意素数有限域上分解多项式为不可约因子，并保留重数。

\[ \displaystyle f(x)=\prod_j (x-r_j)^{e_j} \]

Complexity: Time: Expected O(n^3 log p) with quadratic polynomial arithmetic. Space: O(n^2) across the factor queue and temporaries.

AC 记录：factorization_of_polynomials。

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

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

/// @complexity Time: Expected O(n^3 log p) with quadratic polynomial
/// arithmetic. Space: O(n^2) across the factor queue and temporaries.

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

namespace noya {

/// @brief One monic irreducible factor over F_p and its multiplicity.
struct finite_field_polynomial_factor {
  std::vector<std::uint64_t> pol;
  int cnt = 0;
};

namespace finite_field_polynomial_factorization_internal {

using polynomial = std::vector<std::uint64_t>;

inline void trim(polynomial &val) {
  while (!val.empty() && val.back() == 0) {
    val.pop_back();
  }
}

inline int degree(polynomial val) {
  trim(val);
  return int(val.size()) - 1;
}

inline std::uint64_t power(std::uint64_t bas, std::uint64_t exp,
                           std::uint64_t mod) {
  std::uint64_t res = 1 % mod;
  while (exp > 0) {
    if (exp & 1) {
      res = res * bas % mod;
    }
    bas = bas * bas % mod;
    exp >>= 1;
  }
  return res;
}

inline polynomial monic(polynomial val, std::uint64_t mod) {
  trim(val);
  if (val.empty()) {
    return val;
  }
  std::uint64_t scl = power(val.back(), mod - 2, mod);
  for (std::uint64_t &cf : val) {
    cf = cf * scl % mod;
  }
  return val;
}

inline polynomial add(polynomial l, const polynomial &r, std::uint64_t mod) {
  l.resize(std::max(l.size(), r.size()));
  for (int i = 0; i < int(r.size()); i++) {
    l[i] += r[i];
    if (l[i] >= mod) {
      l[i] -= mod;
    }
  }
  trim(l);
  return l;
}

inline polynomial subtract(polynomial l, const polynomial &r,
                           std::uint64_t mod) {
  l.resize(std::max(l.size(), r.size()));
  for (int i = 0; i < int(r.size()); i++) {
    if (l[i] >= r[i]) {
      l[i] -= r[i];
    } else {
      l[i] += mod - r[i];
    }
  }
  trim(l);
  return l;
}

inline polynomial multiply(const polynomial &l, const polynomial &r,
                           std::uint64_t mod) {
  if (l.empty() || r.empty()) {
    return {};
  }
  polynomial res(l.size() + r.size() - 1);
  for (int i = 0; i < int(l.size()); i++) {
    for (int j = 0; j < int(r.size()); j++) {
      res[i + j] = (res[i + j] + l[i] * r[j]) % mod;
    }
  }
  trim(res);
  return res;
}

inline std::pair<polynomial, polynomial> divmod(polynomial div, polynomial dvs,
                                                std::uint64_t mod) {
  trim(div);
  trim(dvs);
  assert(!dvs.empty());
  if (div.size() < dvs.size()) {
    return {{}, div};
  }
  polynomial quo(div.size() - dvs.size() + 1);
  std::uint64_t il = power(dvs.back(), mod - 2, mod);
  for (int pos = int(div.size() - dvs.size()); pos >= 0; pos--) {
    std::uint64_t cf = div[pos + dvs.size() - 1] * il % mod;
    quo[pos] = cf;
    for (int idx = 0; idx < int(dvs.size()); idx++) {
      std::uint64_t rem = cf * dvs[idx] % mod;
      std::uint64_t &tar = div[pos + idx];
      tar = tar >= rem ? tar - rem : tar + mod - rem;
    }
  }
  trim(quo);
  trim(div);
  return {quo, div};
}

inline polynomial remainder(const polynomial &val, const polynomial &mod,
                            std::uint64_t p) {
  return divmod(val, mod, p).second;
}

inline polynomial gcd(polynomial a, polynomial b, std::uint64_t mod) {
  trim(a);
  trim(b);
  while (!b.empty()) {
    polynomial nxt = remainder(a, b, mod);
    a = std::move(b);
    b = std::move(nxt);
  }
  return monic(std::move(a), mod);
}

inline polynomial multiply_mod(const polynomial &l, const polynomial &r,
                               const polynomial &mdf, std::uint64_t p) {
  return remainder(multiply(l, r, p), mdf, p);
}

inline polynomial power_mod(polynomial bas, std::uint64_t exp,
                            const polynomial &mdf, std::uint64_t p) {
  polynomial res = remainder({1}, mdf, p);
  bas = remainder(bas, mdf, p);
  while (exp > 0) {
    if (exp & 1) {
      res = multiply_mod(res, bas, mdf, p);
    }
    exp >>= 1;
    if (exp > 0) {
      bas = multiply_mod(bas, bas, mdf, p);
    }
  }
  return res;
}

inline std::uint64_t splitmix64(std::uint64_t &st) {
  std::uint64_t val = (st += 0x9e3779b97f4a7c15ULL);
  val = (val ^ (val >> 30)) * 0xbf58476d1ce4e5b9ULL;
  val = (val ^ (val >> 27)) * 0x94d049bb133111ebULL;
  return val ^ (val >> 31);
}

inline polynomial random_polynomial(int nc, std::uint64_t mod,
                                    std::uint64_t &st) {
  polynomial res(nc);
  for (std::uint64_t &cf : res) {
    cf = splitmix64(st) % mod;
  }
  trim(res);
  return res;
}

inline polynomial odd_character(const polynomial &val, int fd,
                                const polynomial &mdf, std::uint64_t p) {
  polynomial cjg = power_mod(val, (p - 1) / 2, mdf, p);
  polynomial res = remainder({1}, mdf, p);
  for (int idx = 0; idx < fd; idx++) {
    res = multiply_mod(res, cjg, mdf, p);
    if (idx + 1 < fd) {
      cjg = power_mod(cjg, p, mdf, p);
    }
  }
  return res;
}

inline polynomial binary_trace(const polynomial &val, int fd,
                               const polynomial &mdf) {
  polynomial cjg = remainder(val, mdf, 2);
  polynomial res;
  for (int idx = 0; idx < fd; idx++) {
    res = add(std::move(res), cjg, 2);
    if (idx + 1 < fd) {
      cjg = multiply_mod(cjg, cjg, mdf, 2);
    }
  }
  return res;
}

inline std::vector<polynomial> equal_degree_factorization(polynomial val,
                                                          int fd,
                                                          std::uint64_t p,
                                                          std::uint64_t &st) {
  val = monic(std::move(val), p);
  std::deque<polynomial> pen = {val};
  std::vector<polynomial> res;
  while (!pen.empty()) {
    polynomial cur = monic(std::move(pen.front()), p);
    pen.pop_front();
    int cd = degree(cur);
    if (cd == fd) {
      res.push_back(std::move(cur));
      continue;
    }
    while (true) {
      polynomial rng = random_polynomial(cd, p, st);
      polynomial sep = p == 2
                           ? binary_trace(rng, fd, cur)
                           : subtract(odd_character(rng, fd, cur, p), {1}, p);
      polynomial l = gcd(cur, sep, p);
      int ld = degree(l);
      if (ld <= 0 || ld == cd) {
        continue;
      }
      polynomial r = divmod(cur, l, p).first;
      pen.push_back(std::move(l));
      pen.push_back(std::move(r));
      break;
    }
  }
  return res;
}

} // namespace finite_field_polynomial_factorization_internal

/// @brief Factor a monic polynomial over F_p. Distinct-degree factorization
/// uses gcd(f,x^(p^d)-x); Cantor-Zassenhaus character or trace tests split each
/// equal-degree part, and exact repeated division recovers multiplicities.
inline std::vector<finite_field_polynomial_factor>
factor_finite_field_polynomial(std::vector<std::uint64_t> pol, std::uint64_t p,
                               std::uint64_t see = 0x13198a2e03707344ULL) {
  using namespace finite_field_polynomial_factorization_internal;
  assert(p >= 2);
  for (std::uint64_t &cf : pol) {
    cf %= p;
  }
  pol = monic(std::move(pol), p);
  std::vector<finite_field_polynomial_factor> res;
  if (degree(pol) <= 0) {
    return res;
  }

  const finite_field_polynomial_factorization_internal::polynomial x = {0, 1};
  auto frb = x;
  for (int fd = 1; degree(pol) > 0 && 2 * fd <= degree(pol); fd++) {
    frb = power_mod(frb, p, pol, p);
    auto sdp = gcd(pol, subtract(frb, x, p), p);
    if (degree(sdp) <= 0) {
      continue;
    }
    auto fs = equal_degree_factorization(std::move(sdp), fd, p, see);
    for (auto &fct : fs) {
      int cnt = 0;
      while (degree(pol) >= degree(fct)) {
        auto [quo, rv] = divmod(pol, fct, p);
        if (!rv.empty()) {
          break;
        }
        pol = std::move(quo);
        cnt++;
      }
      res.push_back({std::move(fct), cnt});
    }
  }
  if (degree(pol) > 0) {
    res.push_back({monic(std::move(pol), p), 1});
  }
  return res;
}

} // namespace noya

#ifndef NOYA_FINITE_FIELD_POLYNOMIAL_FACTORIZATION_HPP
#define NOYA_FINITE_FIELD_POLYNOMIAL_FACTORIZATION_HPP 1

/// @complexity Time: Expected O(n^3 log p) with quadratic polynomial
/// arithmetic. Space: O(n^2) across the factor queue and temporaries.

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

namespace noya {

/// @brief One monic irreducible factor over F_p and its multiplicity.
struct finite_field_polynomial_factor {
  std::vector<std::uint64_t> pol;
  int cnt = 0;
};

namespace finite_field_polynomial_factorization_internal {

using polynomial = std::vector<std::uint64_t>;

inline void trim(polynomial &val) {
  while (!val.empty() && val.back() == 0) {
    val.pop_back();
  }
}

inline int degree(polynomial val) {
  trim(val);
  return int(val.size()) - 1;
}

inline std::uint64_t power(std::uint64_t bas, std::uint64_t exp,
                           std::uint64_t mod) {
  std::uint64_t res = 1 % mod;
  while (exp > 0) {
    if (exp & 1) {
      res = res * bas % mod;
    }
    bas = bas * bas % mod;
    exp >>= 1;
  }
  return res;
}

inline polynomial monic(polynomial val, std::uint64_t mod) {
  trim(val);
  if (val.empty()) {
    return val;
  }
  std::uint64_t scl = power(val.back(), mod - 2, mod);
  for (std::uint64_t &cf : val) {
    cf = cf * scl % mod;
  }
  return val;
}

inline polynomial add(polynomial l, const polynomial &r, std::uint64_t mod) {
  l.resize(std::max(l.size(), r.size()));
  for (int i = 0; i < int(r.size()); i++) {
    l[i] += r[i];
    if (l[i] >= mod) {
      l[i] -= mod;
    }
  }
  trim(l);
  return l;
}

inline polynomial subtract(polynomial l, const polynomial &r,
                           std::uint64_t mod) {
  l.resize(std::max(l.size(), r.size()));
  for (int i = 0; i < int(r.size()); i++) {
    if (l[i] >= r[i]) {
      l[i] -= r[i];
    } else {
      l[i] += mod - r[i];
    }
  }
  trim(l);
  return l;
}

inline polynomial multiply(const polynomial &l, const polynomial &r,
                           std::uint64_t mod) {
  if (l.empty() || r.empty()) {
    return {};
  }
  polynomial res(l.size() + r.size() - 1);
  for (int i = 0; i < int(l.size()); i++) {
    for (int j = 0; j < int(r.size()); j++) {
      res[i + j] = (res[i + j] + l[i] * r[j]) % mod;
    }
  }
  trim(res);
  return res;
}

inline std::pair<polynomial, polynomial> divmod(polynomial div, polynomial dvs,
                                                std::uint64_t mod) {
  trim(div);
  trim(dvs);
  assert(!dvs.empty());
  if (div.size() < dvs.size()) {
    return {{}, div};
  }
  polynomial quo(div.size() - dvs.size() + 1);
  std::uint64_t il = power(dvs.back(), mod - 2, mod);
  for (int pos = int(div.size() - dvs.size()); pos >= 0; pos--) {
    std::uint64_t cf = div[pos + dvs.size() - 1] * il % mod;
    quo[pos] = cf;
    for (int idx = 0; idx < int(dvs.size()); idx++) {
      std::uint64_t rem = cf * dvs[idx] % mod;
      std::uint64_t &tar = div[pos + idx];
      tar = tar >= rem ? tar - rem : tar + mod - rem;
    }
  }
  trim(quo);
  trim(div);
  return {quo, div};
}

inline polynomial remainder(const polynomial &val, const polynomial &mod,
                            std::uint64_t p) {
  return divmod(val, mod, p).second;
}

inline polynomial gcd(polynomial a, polynomial b, std::uint64_t mod) {
  trim(a);
  trim(b);
  while (!b.empty()) {
    polynomial nxt = remainder(a, b, mod);
    a = std::move(b);
    b = std::move(nxt);
  }
  return monic(std::move(a), mod);
}

inline polynomial multiply_mod(const polynomial &l, const polynomial &r,
                               const polynomial &mdf, std::uint64_t p) {
  return remainder(multiply(l, r, p), mdf, p);
}

inline polynomial power_mod(polynomial bas, std::uint64_t exp,
                            const polynomial &mdf, std::uint64_t p) {
  polynomial res = remainder({1}, mdf, p);
  bas = remainder(bas, mdf, p);
  while (exp > 0) {
    if (exp & 1) {
      res = multiply_mod(res, bas, mdf, p);
    }
    exp >>= 1;
    if (exp > 0) {
      bas = multiply_mod(bas, bas, mdf, p);
    }
  }
  return res;
}

inline std::uint64_t splitmix64(std::uint64_t &st) {
  std::uint64_t val = (st += 0x9e3779b97f4a7c15ULL);
  val = (val ^ (val >> 30)) * 0xbf58476d1ce4e5b9ULL;
  val = (val ^ (val >> 27)) * 0x94d049bb133111ebULL;
  return val ^ (val >> 31);
}

inline polynomial random_polynomial(int nc, std::uint64_t mod,
                                    std::uint64_t &st) {
  polynomial res(nc);
  for (std::uint64_t &cf : res) {
    cf = splitmix64(st) % mod;
  }
  trim(res);
  return res;
}

inline polynomial odd_character(const polynomial &val, int fd,
                                const polynomial &mdf, std::uint64_t p) {
  polynomial cjg = power_mod(val, (p - 1) / 2, mdf, p);
  polynomial res = remainder({1}, mdf, p);
  for (int idx = 0; idx < fd; idx++) {
    res = multiply_mod(res, cjg, mdf, p);
    if (idx + 1 < fd) {
      cjg = power_mod(cjg, p, mdf, p);
    }
  }
  return res;
}

inline polynomial binary_trace(const polynomial &val, int fd,
                               const polynomial &mdf) {
  polynomial cjg = remainder(val, mdf, 2);
  polynomial res;
  for (int idx = 0; idx < fd; idx++) {
    res = add(std::move(res), cjg, 2);
    if (idx + 1 < fd) {
      cjg = multiply_mod(cjg, cjg, mdf, 2);
    }
  }
  return res;
}

inline std::vector<polynomial> equal_degree_factorization(polynomial val,
                                                          int fd,
                                                          std::uint64_t p,
                                                          std::uint64_t &st) {
  val = monic(std::move(val), p);
  std::deque<polynomial> pen = {val};
  std::vector<polynomial> res;
  while (!pen.empty()) {
    polynomial cur = monic(std::move(pen.front()), p);
    pen.pop_front();
    int cd = degree(cur);
    if (cd == fd) {
      res.push_back(std::move(cur));
      continue;
    }
    while (true) {
      polynomial rng = random_polynomial(cd, p, st);
      polynomial sep = p == 2
                           ? binary_trace(rng, fd, cur)
                           : subtract(odd_character(rng, fd, cur, p), {1}, p);
      polynomial l = gcd(cur, sep, p);
      int ld = degree(l);
      if (ld <= 0 || ld == cd) {
        continue;
      }
      polynomial r = divmod(cur, l, p).first;
      pen.push_back(std::move(l));
      pen.push_back(std::move(r));
      break;
    }
  }
  return res;
}

} // namespace finite_field_polynomial_factorization_internal

/// @brief Factor a monic polynomial over F_p. Distinct-degree factorization
/// uses gcd(f,x^(p^d)-x); Cantor-Zassenhaus character or trace tests split each
/// equal-degree part, and exact repeated division recovers multiplicities.
inline std::vector<finite_field_polynomial_factor>
factor_finite_field_polynomial(std::vector<std::uint64_t> pol, std::uint64_t p,
                               std::uint64_t see = 0x13198a2e03707344ULL) {
  using namespace finite_field_polynomial_factorization_internal;
  assert(p >= 2);
  for (std::uint64_t &cf : pol) {
    cf %= p;
  }
  pol = monic(std::move(pol), p);
  std::vector<finite_field_polynomial_factor> res;
  if (degree(pol) <= 0) {
    return res;
  }

  const finite_field_polynomial_factorization_internal::polynomial x = {0, 1};
  auto frb = x;
  for (int fd = 1; degree(pol) > 0 && 2 * fd <= degree(pol); fd++) {
    frb = power_mod(frb, p, pol, p);
    auto sdp = gcd(pol, subtract(frb, x, p), p);
    if (degree(sdp) <= 0) {
      continue;
    }
    auto fs = equal_degree_factorization(std::move(sdp), fd, p, see);
    for (auto &fct : fs) {
      int cnt = 0;
      while (degree(pol) >= degree(fct)) {
        auto [quo, rv] = divmod(pol, fct, p);
        if (!rv.empty()) {
          break;
        }
        pol = std::move(quo);
        cnt++;
      }
      res.push_back({std::move(fct), cnt});
    }
  }
  if (degree(pol) > 0) {
    res.push_back({monic(std::move(pol), p), 1});
  }
  return res;
}

} // namespace noya

#endif // NOYA_FINITE_FIELD_POLYNOMIAL_FACTORIZATION_HPP

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

/// @complexity Time: Expected O(n^3 log p) with quadratic polynomial
/// arithmetic. Space: O(n^2) across the factor queue and temporaries.

namespace noya {

/// @brief One monic irreducible factor over F_p and its multiplicity.
struct finite_field_polynomial_factor {
  std::vector<std::uint64_t> pol;
  int cnt = 0;
};

namespace finite_field_polynomial_factorization_internal {

using polynomial = std::vector<std::uint64_t>;

inline void trim(polynomial &val) {
  while (!val.empty() && val.back() == 0) {
    val.pop_back();
  }
}

inline int degree(polynomial val) {
  trim(val);
  return int(val.size()) - 1;
}

inline std::uint64_t power(std::uint64_t bas, std::uint64_t exp,
                           std::uint64_t mod) {
  std::uint64_t res = 1 % mod;
  while (exp > 0) {
    if (exp & 1) {
      res = res * bas % mod;
    }
    bas = bas * bas % mod;
    exp >>= 1;
  }
  return res;
}

inline polynomial monic(polynomial val, std::uint64_t mod) {
  trim(val);
  if (val.empty()) {
    return val;
  }
  std::uint64_t scl = power(val.back(), mod - 2, mod);
  for (std::uint64_t &cf : val) {
    cf = cf * scl % mod;
  }
  return val;
}

inline polynomial add(polynomial l, const polynomial &r, std::uint64_t mod) {
  l.resize(std::max(l.size(), r.size()));
  for (int i = 0; i < int(r.size()); i++) {
    l[i] += r[i];
    if (l[i] >= mod) {
      l[i] -= mod;
    }
  }
  trim(l);
  return l;
}

inline polynomial subtract(polynomial l, const polynomial &r,
                           std::uint64_t mod) {
  l.resize(std::max(l.size(), r.size()));
  for (int i = 0; i < int(r.size()); i++) {
    if (l[i] >= r[i]) {
      l[i] -= r[i];
    } else {
      l[i] += mod - r[i];
    }
  }
  trim(l);
  return l;
}

inline polynomial multiply(const polynomial &l, const polynomial &r,
                           std::uint64_t mod) {
  if (l.empty() || r.empty()) {
    return {};
  }
  polynomial res(l.size() + r.size() - 1);
  for (int i = 0; i < int(l.size()); i++) {
    for (int j = 0; j < int(r.size()); j++) {
      res[i + j] = (res[i + j] + l[i] * r[j]) % mod;
    }
  }
  trim(res);
  return res;
}

inline std::pair<polynomial, polynomial> divmod(polynomial div, polynomial dvs,
                                                std::uint64_t mod) {
  trim(div);
  trim(dvs);
  assert(!dvs.empty());
  if (div.size() < dvs.size()) {
    return {{}, div};
  }
  polynomial quo(div.size() - dvs.size() + 1);
  std::uint64_t il = power(dvs.back(), mod - 2, mod);
  for (int pos = int(div.size() - dvs.size()); pos >= 0; pos--) {
    std::uint64_t cf = div[pos + dvs.size() - 1] * il % mod;
    quo[pos] = cf;
    for (int idx = 0; idx < int(dvs.size()); idx++) {
      std::uint64_t rem = cf * dvs[idx] % mod;
      std::uint64_t &tar = div[pos + idx];
      tar = tar >= rem ? tar - rem : tar + mod - rem;
    }
  }
  trim(quo);
  trim(div);
  return {quo, div};
}

inline polynomial remainder(const polynomial &val, const polynomial &mod,
                            std::uint64_t p) {
  return divmod(val, mod, p).second;
}

inline polynomial gcd(polynomial a, polynomial b, std::uint64_t mod) {
  trim(a);
  trim(b);
  while (!b.empty()) {
    polynomial nxt = remainder(a, b, mod);
    a = std::move(b);
    b = std::move(nxt);
  }
  return monic(std::move(a), mod);
}

inline polynomial multiply_mod(const polynomial &l, const polynomial &r,
                               const polynomial &mdf, std::uint64_t p) {
  return remainder(multiply(l, r, p), mdf, p);
}

inline polynomial power_mod(polynomial bas, std::uint64_t exp,
                            const polynomial &mdf, std::uint64_t p) {
  polynomial res = remainder({1}, mdf, p);
  bas = remainder(bas, mdf, p);
  while (exp > 0) {
    if (exp & 1) {
      res = multiply_mod(res, bas, mdf, p);
    }
    exp >>= 1;
    if (exp > 0) {
      bas = multiply_mod(bas, bas, mdf, p);
    }
  }
  return res;
}

inline std::uint64_t splitmix64(std::uint64_t &st) {
  std::uint64_t val = (st += 0x9e3779b97f4a7c15ULL);
  val = (val ^ (val >> 30)) * 0xbf58476d1ce4e5b9ULL;
  val = (val ^ (val >> 27)) * 0x94d049bb133111ebULL;
  return val ^ (val >> 31);
}

inline polynomial random_polynomial(int nc, std::uint64_t mod,
                                    std::uint64_t &st) {
  polynomial res(nc);
  for (std::uint64_t &cf : res) {
    cf = splitmix64(st) % mod;
  }
  trim(res);
  return res;
}

inline polynomial odd_character(const polynomial &val, int fd,
                                const polynomial &mdf, std::uint64_t p) {
  polynomial cjg = power_mod(val, (p - 1) / 2, mdf, p);
  polynomial res = remainder({1}, mdf, p);
  for (int idx = 0; idx < fd; idx++) {
    res = multiply_mod(res, cjg, mdf, p);
    if (idx + 1 < fd) {
      cjg = power_mod(cjg, p, mdf, p);
    }
  }
  return res;
}

inline polynomial binary_trace(const polynomial &val, int fd,
                               const polynomial &mdf) {
  polynomial cjg = remainder(val, mdf, 2);
  polynomial res;
  for (int idx = 0; idx < fd; idx++) {
    res = add(std::move(res), cjg, 2);
    if (idx + 1 < fd) {
      cjg = multiply_mod(cjg, cjg, mdf, 2);
    }
  }
  return res;
}

inline std::vector<polynomial> equal_degree_factorization(polynomial val,
                                                          int fd,
                                                          std::uint64_t p,
                                                          std::uint64_t &st) {
  val = monic(std::move(val), p);
  std::deque<polynomial> pen = {val};
  std::vector<polynomial> res;
  while (!pen.empty()) {
    polynomial cur = monic(std::move(pen.front()), p);
    pen.pop_front();
    int cd = degree(cur);
    if (cd == fd) {
      res.push_back(std::move(cur));
      continue;
    }
    while (true) {
      polynomial rng = random_polynomial(cd, p, st);
      polynomial sep = p == 2
                           ? binary_trace(rng, fd, cur)
                           : subtract(odd_character(rng, fd, cur, p), {1}, p);
      polynomial l = gcd(cur, sep, p);
      int ld = degree(l);
      if (ld <= 0 || ld == cd) {
        continue;
      }
      polynomial r = divmod(cur, l, p).first;
      pen.push_back(std::move(l));
      pen.push_back(std::move(r));
      break;
    }
  }
  return res;
}

} // namespace finite_field_polynomial_factorization_internal

/// @brief Factor a monic polynomial over F_p. Distinct-degree factorization
/// uses gcd(f,x^(p^d)-x); Cantor-Zassenhaus character or trace tests split each
/// equal-degree part, and exact repeated division recovers multiplicities.
inline std::vector<finite_field_polynomial_factor>
factor_finite_field_polynomial(std::vector<std::uint64_t> pol, std::uint64_t p,
                               std::uint64_t see = 0x13198a2e03707344ULL) {
  using namespace finite_field_polynomial_factorization_internal;
  assert(p >= 2);
  for (std::uint64_t &cf : pol) {
    cf %= p;
  }
  pol = monic(std::move(pol), p);
  std::vector<finite_field_polynomial_factor> res;
  if (degree(pol) <= 0) {
    return res;
  }

  const finite_field_polynomial_factorization_internal::polynomial x = {0, 1};
  auto frb = x;
  for (int fd = 1; degree(pol) > 0 && 2 * fd <= degree(pol); fd++) {
    frb = power_mod(frb, p, pol, p);
    auto sdp = gcd(pol, subtract(frb, x, p), p);
    if (degree(sdp) <= 0) {
      continue;
    }
    auto fs = equal_degree_factorization(std::move(sdp), fd, p, see);
    for (auto &fct : fs) {
      int cnt = 0;
      while (degree(pol) >= degree(fct)) {
        auto [quo, rv] = divmod(pol, fct, p);
        if (!rv.empty()) {
          break;
        }
        pol = std::move(quo);
        cnt++;
      }
      res.push_back({std::move(fct), cnt});
    }
  }
  if (degree(pol) > 0) {
    res.push_back({monic(std::move(pol), p), 1});
  }
  return res;
}

} // namespace noya