finite_field_polynomial_factorization.hpp¶
在任意素数有限域上分解多项式为不可约因子，并保留重数。
\[
\displaystyle f(x)=\prod_j g_j(x)^{e_j}
\]
Complexity: Time: Expected \(O(n^{3} \log p)\) with quadratic polynomial arithmetic. Space: \(O(n)\) coefficients across the factor queue and temporaries.
AC 记录：factorization_of_polynomials。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: Expected O(n^3 log p) with quadratic polynomial
/// arithmetic. Space: O(n) coefficients 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(const polynomial &val) {
int pos = int(val.size()) - 1;
while (pos >= 0 && val[pos] == 0) --pos;
return pos;
}
inline std::uint64_t mul(std::uint64_t a, std::uint64_t b, std::uint64_t p) {
if (p <= UINT32_MAX) return a * b % p;
return static_cast<unsigned __int128>(a) * b % p;
}
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 = mul(res, bas, mod);
}
bas = mul(bas, bas, mod);
exp >>= 1;
}
return res;
}
inline polynomial monic(polynomial val, std::uint64_t mod) {
trim(val);
if (val.empty() || val.back() == 1) {
return val;
}
std::uint64_t scl = power(val.back(), mod - 2, mod);
for (std::uint64_t &cf : val) {
cf = mul(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] = l[i] >= mod - r[i] ? l[i] - (mod - r[i]) : l[i] + r[i];
}
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;
}
template <class Wide>
inline polynomial multiply_impl(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] = (Wide(res[i + j]) + Wide(l[i]) * r[j]) % mod;
}
}
trim(res);
return res;
}
inline polynomial multiply(const polynomial &l, const polynomial &r,
std::uint64_t mod) {
// (p-1)^2 + (p-1) fits 64 bits for every 32-bit prime.
if (mod <= UINT32_MAX) return multiply_impl<std::uint64_t>(l, r, mod);
return multiply_impl<unsigned __int128>(l, r, mod);
}
inline void divide(polynomial &div, const polynomial &dvs, std::uint64_t mod,
polynomial *quo = nullptr) {
trim(div);
int deg = degree(dvs);
assert(deg >= 0);
int len = std::max(0, int(div.size()) - deg);
if (quo) quo->assign(len, 0);
if (len == 0) return;
std::uint64_t il = dvs[deg] == 1 ? 1 : power(dvs[deg], mod - 2, mod);
for (int pos = len - 1; pos >= 0; pos--) {
std::uint64_t cf = il == 1 ? div[pos + deg] : mul(div[pos + deg], il, mod);
if (quo) (*quo)[pos] = cf;
div[pos + deg] = 0;
if (cf == 0) continue;
for (int idx = 0; idx < deg; idx++) {
std::uint64_t rem = mul(cf, dvs[idx], mod);
std::uint64_t &tar = div[pos + idx];
tar = tar >= rem ? tar - rem : tar + mod - rem;
}
}
trim(div);
}
inline std::pair<polynomial, polynomial> divmod(polynomial div,
const polynomial &dvs,
std::uint64_t mod) {
polynomial quo;
divide(div, dvs, mod, &quo);
trim(quo);
return {std::move(quo), std::move(div)};
}
inline polynomial remainder(polynomial val, const polynomial &mod,
std::uint64_t p) {
divide(val, mod, p);
return val;
}
inline polynomial gcd(polynomial a, polynomial b, std::uint64_t mod) {
trim(a);
trim(b);
while (!b.empty()) {
polynomial nxt = remainder(std::move(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(std::move(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;
pen.push_back(std::move(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.
/// p must be prime; coefficients are reduced modulo p. Nonzero inputs are
/// normalized to monic form. Constants and the zero polynomial return no factors.
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
// NOYA_FINITE_FIELD_POLYNOMIAL_FACTORIZATION_HPP
#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) coefficients 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(const polynomial &val) {
int pos = int(val.size()) - 1;
while (pos >= 0 && val[pos] == 0) --pos;
return pos;
}
inline std::uint64_t mul(std::uint64_t a, std::uint64_t b, std::uint64_t p) {
if (p <= UINT32_MAX) return a * b % p;
return static_cast<unsigned __int128>(a) * b % p;
}
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 = mul(res, bas, mod);
}
bas = mul(bas, bas, mod);
exp >>= 1;
}
return res;
}
inline polynomial monic(polynomial val, std::uint64_t mod) {
trim(val);
if (val.empty() || val.back() == 1) {
return val;
}
std::uint64_t scl = power(val.back(), mod - 2, mod);
for (std::uint64_t &cf : val) {
cf = mul(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] = l[i] >= mod - r[i] ? l[i] - (mod - r[i]) : l[i] + r[i];
}
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;
}
template <class Wide>
inline polynomial multiply_impl(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] = (Wide(res[i + j]) + Wide(l[i]) * r[j]) % mod;
}
}
trim(res);
return res;
}
inline polynomial multiply(const polynomial &l, const polynomial &r,
std::uint64_t mod) {
// (p-1)^2 + (p-1) fits 64 bits for every 32-bit prime.
if (mod <= UINT32_MAX) return multiply_impl<std::uint64_t>(l, r, mod);
return multiply_impl<unsigned __int128>(l, r, mod);
}
inline void divide(polynomial &div, const polynomial &dvs, std::uint64_t mod,
polynomial *quo = nullptr) {
trim(div);
int deg = degree(dvs);
assert(deg >= 0);
int len = std::max(0, int(div.size()) - deg);
if (quo) quo->assign(len, 0);
if (len == 0) return;
std::uint64_t il = dvs[deg] == 1 ? 1 : power(dvs[deg], mod - 2, mod);
for (int pos = len - 1; pos >= 0; pos--) {
std::uint64_t cf = il == 1 ? div[pos + deg] : mul(div[pos + deg], il, mod);
if (quo) (*quo)[pos] = cf;
div[pos + deg] = 0;
if (cf == 0) continue;
for (int idx = 0; idx < deg; idx++) {
std::uint64_t rem = mul(cf, dvs[idx], mod);
std::uint64_t &tar = div[pos + idx];
tar = tar >= rem ? tar - rem : tar + mod - rem;
}
}
trim(div);
}
inline std::pair<polynomial, polynomial> divmod(polynomial div,
const polynomial &dvs,
std::uint64_t mod) {
polynomial quo;
divide(div, dvs, mod, &quo);
trim(quo);
return {std::move(quo), std::move(div)};
}
inline polynomial remainder(polynomial val, const polynomial &mod,
std::uint64_t p) {
divide(val, mod, p);
return val;
}
inline polynomial gcd(polynomial a, polynomial b, std::uint64_t mod) {
trim(a);
trim(b);
while (!b.empty()) {
polynomial nxt = remainder(std::move(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(std::move(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;
pen.push_back(std::move(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.
/// p must be prime; coefficients are reduced modulo p. Nonzero inputs are
/// normalized to monic form. Constants and the zero polynomial return no factors.
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
/// @complexity Time: Expected O(n^3 log p) with quadratic polynomial
/// arithmetic. Space: O(n) coefficients 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(const polynomial &val) {
int pos = int(val.size()) - 1;
while (pos >= 0 && val[pos] == 0) --pos;
return pos;
}
inline std::uint64_t mul(std::uint64_t a, std::uint64_t b, std::uint64_t p) {
if (p <= UINT32_MAX) return a * b % p;
return static_cast<unsigned __int128>(a) * b % p;
}
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 = mul(res, bas, mod);
}
bas = mul(bas, bas, mod);
exp >>= 1;
}
return res;
}
inline polynomial monic(polynomial val, std::uint64_t mod) {
trim(val);
if (val.empty() || val.back() == 1) {
return val;
}
std::uint64_t scl = power(val.back(), mod - 2, mod);
for (std::uint64_t &cf : val) {
cf = mul(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] = l[i] >= mod - r[i] ? l[i] - (mod - r[i]) : l[i] + r[i];
}
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;
}
template <class Wide>
inline polynomial multiply_impl(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] = (Wide(res[i + j]) + Wide(l[i]) * r[j]) % mod;
}
}
trim(res);
return res;
}
inline polynomial multiply(const polynomial &l, const polynomial &r,
std::uint64_t mod) {
// (p-1)^2 + (p-1) fits 64 bits for every 32-bit prime.
if (mod <= UINT32_MAX) return multiply_impl<std::uint64_t>(l, r, mod);
return multiply_impl<unsigned __int128>(l, r, mod);
}
inline void divide(polynomial &div, const polynomial &dvs, std::uint64_t mod,
polynomial *quo = nullptr) {
trim(div);
int deg = degree(dvs);
assert(deg >= 0);
int len = std::max(0, int(div.size()) - deg);
if (quo) quo->assign(len, 0);
if (len == 0) return;
std::uint64_t il = dvs[deg] == 1 ? 1 : power(dvs[deg], mod - 2, mod);
for (int pos = len - 1; pos >= 0; pos--) {
std::uint64_t cf = il == 1 ? div[pos + deg] : mul(div[pos + deg], il, mod);
if (quo) (*quo)[pos] = cf;
div[pos + deg] = 0;
if (cf == 0) continue;
for (int idx = 0; idx < deg; idx++) {
std::uint64_t rem = mul(cf, dvs[idx], mod);
std::uint64_t &tar = div[pos + idx];
tar = tar >= rem ? tar - rem : tar + mod - rem;
}
}
trim(div);
}
inline std::pair<polynomial, polynomial> divmod(polynomial div,
const polynomial &dvs,
std::uint64_t mod) {
polynomial quo;
divide(div, dvs, mod, &quo);
trim(quo);
return {std::move(quo), std::move(div)};
}
inline polynomial remainder(polynomial val, const polynomial &mod,
std::uint64_t p) {
divide(val, mod, p);
return val;
}
inline polynomial gcd(polynomial a, polynomial b, std::uint64_t mod) {
trim(a);
trim(b);
while (!b.empty()) {
polynomial nxt = remainder(std::move(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(std::move(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;
pen.push_back(std::move(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.
/// p must be prime; coefficients are reduced modulo p. Nonzero inputs are
/// normalized to monic form. Constants and the zero polynomial return no factors.
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
// NOYA_FINITE_FIELD_POLYNOMIAL_FACTORIZATION_HPP