-
Notifications
You must be signed in to change notification settings - Fork 0
/
findroots.hpp
73 lines (71 loc) · 1.75 KB
/
findroots.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
#pragma once
#include "Utility/random.hpp"
#include "FPS/halfgcd.hpp"
#include "FPS/prodofpolys.hpp"
namespace FindRoots {
template <typename T> Poly<T> powmod(Poly<T> &f, ll 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 /= 2;
}
return ret;
}
template <typename T> vector<Poly<T>> EDF(Poly<T> &f) {
if (f.deg() < 1)
return {};
if (f.deg() == 1)
return {f};
Poly<T> base(2);
base[0] = Random::get(T::get_mod() - 1);
base[1] = 1;
ll pw = (T::get_mod() - 1) / 2;
auto rem = powmod(base, pw, f);
if (rem.empty())
return EDF(f);
rem -= Poly<T>({1});
auto h = HalfGCD::gcd(rem, f);
auto ret = EDF(h);
auto fh = f / h;
auto add = EDF(fh);
ret.insert(ret.end(), ALL(add));
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);
g *= ProdOfPolys(base);
f /= g;
base.insert(base.begin(), f);
return base;
}
template <typename T> Poly<T> select(Poly<T> &f) {
Poly<T> X({0, 1});
auto xpmf = powmod(X, T::get_mod(), f);
xpmf -= X;
return HalfGCD::gcd(xpmf, f);
}
template <typename T> vector<T> run(Poly<T> &f) {
auto dec = SquarefreeDecomposition(f);
vector<T> ret;
rep(i, 0, SZ(dec)) {
auto g = select(dec[i]);
auto add = EDF(g);
for (auto &h : add) {
assert(h.deg() == 1);
rep(_, 0, i + 1) ret.push_back(-h[0]);
}
}
return ret;
}
}; // namespace FindRoots
/**
* @brief Find roots
*/