library
This documentation is automatically generated by online-judge-tools/verification-helper
Factorize Polynomial
(FPS/factorize.hpp)
- View this file on GitHub
- Last update: 2024-04-26 03:18:17+09:00
- Include:
#include "FPS/factorize.hpp"
Depends on
Random
(Utility/random.hpp)
Code
#pragma once
#include "Utility/random.hpp"
#include "FPS/halfgcd.hpp"
#include <boost/multiprecision/cpp_int.hpp>
namespace mp = boost::multiprecision;
using Bint = mp::cpp_int;
namespace FactorizePoly {
template <typename T> Poly<T> powmod(Poly<T> &f, Bint &n, Poly<T> &g) {
Poly<T> ret({1}), base = f;
while (n != 0) {
if (n & 1) {
ret *= base;
ret %= g;
}
base *= base;
base %= g;
n >>= 1;
}
return ret;
}
template <typename T> vector<Poly<T>> EDF(Poly<T> &f, int d) {
if (f.deg() < d)
return {};
if (f.deg() == d)
return {f};
Poly<T> base(SZ(f));
rep(i, 0, SZ(f)) base[i] = Random::get(T::get_mod() - 1);
auto g = HalfGCD::gcd(base, f);
if (g.deg() != 0) {
auto ret = EDF(g, d);
auto fg = f / g;
auto add = EDF(fg, d);
ret.insert(ret.end(), ALL(add));
return ret;
} else {
Bint pw = 1;
rep(_, 0, d) pw *= T::get_mod();
pw = (pw - 1) / 2;
auto rem = powmod(base, pw, f);
rem[0] -= 1;
g = HalfGCD::gcd(rem, f);
if (g.deg() != 0 and g != f) {
auto ret = EDF(g, d);
auto fg = f / g;
auto add = EDF(fg, d);
ret.insert(ret.end(), ALL(add));
return ret;
} else {
return EDF(f, d);
}
}
}
template <typename T> vector<Poly<T>> CantorZassenhaus(Poly<T> &f) {
auto base = Poly<T>({0, 1});
auto cur = f;
vector<Poly<T>> ret;
int d = 1;
Bint k = 1, md = T::get_mod();
for (;;) {
if (cur.deg() < d * 2) {
if (cur.deg())
ret.push_back(cur);
break;
}
k *= md;
auto rem = powmod(base, k, cur);
rem -= base;
auto g = HalfGCD::gcd(rem, cur);
auto add = EDF(g, d);
ret.insert(ret.end(), ALL(add));
cur /= g;
d++;
}
return ret;
}
template <typename T> vector<Poly<T>> SquarefreeDecomposition(Poly<T> f) {
if (f.deg() == 0)
return {f};
auto g = HalfGCD::gcd(f, f.diff());
auto base = SquarefreeDecomposition(g);
for (auto &b : base)
g *= b;
f /= g;
base.insert(base.begin(), f);
return base;
}
template <typename T> vector<Poly<T>> run(Poly<T> &f) { // f: monic
auto dec = SquarefreeDecomposition(f);
vector<Poly<T>> ret;
rep(i, 0, SZ(dec)) {
auto add = CantorZassenhaus(dec[i]);
rep(_, 0, i + 1) ret.insert(ret.end(), ALL(add));
}
return ret;
}
}; // namespace FactorizePoly
/**
* @brief Factorize Polynomial
*/#line 2 "Utility/random.hpp"
namespace Random {
mt19937_64 randgen(chrono::steady_clock::now().time_since_epoch().count());
using u64 = unsigned long long;
u64 get() {
return randgen();
}
template <typename T> T get(T L) { // [0,L]
return get() % (L + 1);
}
template <typename T> T get(T L, T R) { // [L,R]
return get(R - L) + L;
}
double uniform() {
return double(get(1000000000)) / 1000000000;
}
string str(int n) {
string ret;
rep(i, 0, n) ret += get('a', 'z');
return ret;
}
template <typename Iter> void shuffle(Iter first, Iter last) {
if (first == last)
return;
int len = 1;
for (auto it = first + 1; it != last; it++) {
len++;
int j = get(0, len - 1);
if (j != len - 1)
iter_swap(it, first + j);
}
}
template <typename T> vector<T> select(int n, T L, T R) { // [L,R]
if (n * 2 >= R - L + 1) {
vector<T> ret(R - L + 1);
iota(ALL(ret), L);
shuffle(ALL(ret));
ret.resize(n);
return ret;
} else {
unordered_set<T> used;
vector<T> ret;
while (SZ(used) < n) {
T x = get(L, R);
if (!used.count(x)) {
used.insert(x);
ret.push_back(x);
}
}
return ret;
}
}
void relabel(int n, vector<pair<int, int>> &es) {
shuffle(ALL(es));
vector<int> ord(n);
iota(ALL(ord), 0);
shuffle(ALL(ord));
for (auto &[u, v] : es)
u = ord[u], v = ord[v];
}
template <bool directed, bool simple> vector<pair<int, int>> genGraph(int n) {
vector<pair<int, int>> cand, es;
rep(u, 0, n) rep(v, 0, n) {
if (simple and u == v)
continue;
if (!directed and u > v)
continue;
cand.push_back({u, v});
}
int m = get(SZ(cand));
vector<int> ord;
if (simple)
ord = select(m, 0, SZ(cand) - 1);
else {
rep(_, 0, m) ord.push_back(get(SZ(cand) - 1));
}
for (auto &i : ord)
es.push_back(cand[i]);
relabel(n, es);
return es;
}
vector<pair<int, int>> genTree(int n) {
vector<pair<int, int>> es;
rep(i, 1, n) es.push_back({get(i - 1), i});
relabel(n, es);
return es;
}
}; // namespace Random
/**
* @brief Random
*/
#line 2 "FPS/halfgcd.hpp"
namespace HalfGCD{
template<typename T>using P=array<T,2>;
template<typename T>using Mat=array<T,4>;
template<typename T>P<T> operator*(const Mat<T>& a,const P<T>& b){
P<T> ret={a[0]*b[0]+a[1]*b[1],a[2]*b[0]+a[3]*b[1]};
rep(i,0,2)ret[i].shrink();
return ret;
}
template<typename T>Mat<T> operator*(const Mat<T>& a,const Mat<T>& b){
Mat<T> ret={a[0]*b[0]+a[1]*b[2],a[0]*b[1]+a[1]*b[3],
a[2]*b[0]+a[3]*b[2],a[2]*b[1]+a[3]*b[3]};
rep(i,0,4)ret[i].shrink();
return ret;
}
template<typename T>Mat<T> HGCD(P<T> a){
int m=(SZ(a[0])+1)>>1;
if(SZ(a[1])<=m){
Mat<T> ret;
ret[0]={1},ret[3]={1};
return ret;
}
auto R=HGCD(P<T>{a[0]>>m,a[1]>>m});
a=R*a;
if(SZ(a[1])<=m)return R;
Mat<T> Q;
Q[1]={1},Q[2]={1},Q[3]=-(a[0]/a[1]);
R=Q*R,a=Q*a;
if(SZ(a[1])<=m)return R;
int k=2*m+1-SZ(a[0]);
auto H=HGCD(P<T>{a[0]>>k,a[1]>>k});
return H*R;
}
template<typename T>Mat<T> InnerGCD(P<T> a){
if(SZ(a[0])<SZ(a[1])){
auto M=InnerGCD(P<T>{a[1],a[0]});
swap(M[0],M[1]);
swap(M[2],M[3]);
return M;
}
auto m0=HGCD(a);
a=m0*a;
if(a[1].empty())return m0;
Mat<T> Q;
Q[1]={1},Q[2]={1},Q[3]=-(a[0]/a[1]);
m0=Q*m0,a=Q*a;
if(a[1].empty())return m0;
return InnerGCD(a)*m0;
}
template<typename T>T gcd(const T& a,const T& b){
P<T> p({a,b});
auto M=InnerGCD(p);
p=M*p;
if(!p[0].empty()){
auto coeff=p[0].back().inv();
for(auto& x:p[0])x*=coeff;
}
return p[0];
}
template<typename T>pair<bool,T> PolyInv(const T& a,const T& b){
P<T> p({a,b});
auto M=InnerGCD(p);
T g=(M*p)[0];
if(g.size()!=1)return {false,{}};
P<T> x({T({1}),b});
auto ret=(M*x)[0]%b;
auto coeff=g[0].inv();
for(auto& x:ret)x*=coeff;
return {true,ret};
}
}
/**
* @brief Half GCD
*/
#line 4 "FPS/factorize.hpp"
#include <boost/multiprecision/cpp_int.hpp>
namespace mp = boost::multiprecision;
using Bint = mp::cpp_int;
namespace FactorizePoly {
template <typename T> Poly<T> powmod(Poly<T> &f, Bint &n, Poly<T> &g) {
Poly<T> ret({1}), base = f;
while (n != 0) {
if (n & 1) {
ret *= base;
ret %= g;
}
base *= base;
base %= g;
n >>= 1;
}
return ret;
}
template <typename T> vector<Poly<T>> EDF(Poly<T> &f, int d) {
if (f.deg() < d)
return {};
if (f.deg() == d)
return {f};
Poly<T> base(SZ(f));
rep(i, 0, SZ(f)) base[i] = Random::get(T::get_mod() - 1);
auto g = HalfGCD::gcd(base, f);
if (g.deg() != 0) {
auto ret = EDF(g, d);
auto fg = f / g;
auto add = EDF(fg, d);
ret.insert(ret.end(), ALL(add));
return ret;
} else {
Bint pw = 1;
rep(_, 0, d) pw *= T::get_mod();
pw = (pw - 1) / 2;
auto rem = powmod(base, pw, f);
rem[0] -= 1;
g = HalfGCD::gcd(rem, f);
if (g.deg() != 0 and g != f) {
auto ret = EDF(g, d);
auto fg = f / g;
auto add = EDF(fg, d);
ret.insert(ret.end(), ALL(add));
return ret;
} else {
return EDF(f, d);
}
}
}
template <typename T> vector<Poly<T>> CantorZassenhaus(Poly<T> &f) {
auto base = Poly<T>({0, 1});
auto cur = f;
vector<Poly<T>> ret;
int d = 1;
Bint k = 1, md = T::get_mod();
for (;;) {
if (cur.deg() < d * 2) {
if (cur.deg())
ret.push_back(cur);
break;
}
k *= md;
auto rem = powmod(base, k, cur);
rem -= base;
auto g = HalfGCD::gcd(rem, cur);
auto add = EDF(g, d);
ret.insert(ret.end(), ALL(add));
cur /= g;
d++;
}
return ret;
}
template <typename T> vector<Poly<T>> SquarefreeDecomposition(Poly<T> f) {
if (f.deg() == 0)
return {f};
auto g = HalfGCD::gcd(f, f.diff());
auto base = SquarefreeDecomposition(g);
for (auto &b : base)
g *= b;
f /= g;
base.insert(base.begin(), f);
return base;
}
template <typename T> vector<Poly<T>> run(Poly<T> &f) { // f: monic
auto dec = SquarefreeDecomposition(f);
vector<Poly<T>> ret;
rep(i, 0, SZ(dec)) {
auto add = CantorZassenhaus(dec[i]);
rep(_, 0, i + 1) ret.insert(ret.end(), ALL(add));
}
return ret;
}
}; // namespace FactorizePoly
/**
* @brief Factorize Polynomial
*/