Find roots
(FPS/findroots.hpp)

• View this file on GitHub
• Last update: 2024-04-26 03:18:17+09:00
• Include: #include "FPS/findroots.hpp"

Depends on

• Half GCD (FPS/halfgcd.hpp)
• Product of Polynomials (FPS/prodofpolys.hpp)
• Random (Utility/random.hpp)

Code

#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
 */

#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 2 "FPS/prodofpolys.hpp"

template<typename T>Poly<T> ProdOfPolys(vector<Poly<T>>& fs){
    if(fs.empty())return Poly<T>({T(1)});
    sort(ALL(fs),[&](Poly<T>& a,Poly<T>& b){return a.size()<b.size();});
    deque<Poly<T>> deq;
    for(auto& f:fs)deq.push_back(f);
    while(deq.size()>1){
        deq.push_back(deq[0]*deq[1]);
        deq.pop_front();
        deq.pop_front();
    }
    return deq[0];
}

/**
 * @brief Product of Polynomials
*/
#line 5 "FPS/findroots.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
 */

Back to top page