Program Listing for File fields.cpp

↰ Return to documentation for file (symengine/symengine/fields.cpp)

#include <symengine/fields.h>
#include <symengine/add.h>
#include <symengine/constants.h>
#include <symengine/mul.h>
#include <symengine/pow.h>
#include <symengine/symengine_exception.h>

namespace SymEngine
{
GaloisField::GaloisField(const RCP<const Basic> &var, GaloisFieldDict &&dict)
    : UIntPolyBase(var, std::move(dict))
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(get_poly()))
}

bool GaloisField::is_canonical(const GaloisFieldDict &dict) const
{
    // Check if dictionary contains terms with coeffienct 0
    if (dict.modulo_ <= integer_class(0))
        return false;
    if (not dict.empty())
        if (dict.dict_[dict.dict_.size() - 1] == integer_class(0))
            return false;
    return true;
}

hash_t GaloisField::__hash__() const
{
    hash_t seed = SYMENGINE_GALOISFIELD;

    seed += get_var()->hash();
    for (const auto &it : get_poly().dict_) {
        hash_t temp = SYMENGINE_GALOISFIELD;
        hash_combine<hash_t>(temp, mp_get_si(it));
        seed += temp;
    }
    return seed;
}

int GaloisField::compare(const Basic &o) const
{
    const GaloisField &s = down_cast<const GaloisField &>(o);

    if (get_poly().size() != s.get_poly().size())
        return (get_poly().size() < s.get_poly().size()) ? -1 : 1;

    int cmp = unified_compare(get_var(), s.get_var());
    if (cmp != 0)
        return cmp;

    cmp = unified_compare(get_poly().modulo_, s.get_poly().modulo_);
    if (cmp != 0)
        return cmp;

    return unified_compare(get_poly().dict_, s.get_poly().dict_);
}

RCP<const GaloisField> GaloisField::from_dict(const RCP<const Basic> &var,
                                              GaloisFieldDict &&d)
{
    return make_rcp<const GaloisField>(var, std::move(d));
}

RCP<const GaloisField>
GaloisField::from_vec(const RCP<const Basic> &var,
                      const std::vector<integer_class> &v,
                      const integer_class &modulo)
{
    return make_rcp<const GaloisField>(var,
                                       GaloisFieldDict::from_vec(v, modulo));
}

RCP<const GaloisField> GaloisField::from_uintpoly(const UIntPoly &a,
                                                  const integer_class &modulo)
{
    GaloisFieldDict wrapper(a.get_poly().get_dict(), modulo);
    return GaloisField::from_dict(a.get_var(), std::move(wrapper));
}

vec_basic GaloisField::get_args() const
{
    vec_basic args;
    if (get_poly().dict_.empty())
        args.push_back(zero);
    else {
        for (unsigned i = 0; i < get_poly().dict_.size(); i++) {
            if (get_poly().dict_[i] == integer_class(0))
                continue;
            if (i == 0) {
                args.push_back(integer(get_poly().dict_[i]));
            } else if (i == 1) {
                if (get_poly().dict_[i] == 1) {
                    args.push_back(get_var());
                } else {
                    args.push_back(Mul::from_dict(integer(get_poly().dict_[i]),
                                                  {{get_var(), one}}));
                }
            } else {
                if (get_poly().dict_[i] == 1) {
                    args.push_back(pow(get_var(), integer(i)));
                } else {
                    args.push_back(Mul::from_dict(integer(get_poly().dict_[i]),
                                                  {{get_var(), integer(i)}}));
                }
            }
        }
    }
    return args;
}

GaloisFieldDict::GaloisFieldDict(const int &i, const integer_class &mod)
    : modulo_(mod)
{
    integer_class temp;
    mp_fdiv_r(temp, integer_class(i), modulo_);
    if (temp != integer_class(0))
        dict_.insert(dict_.begin(), temp);
}

GaloisFieldDict::GaloisFieldDict(const map_uint_mpz &p,
                                 const integer_class &mod)
    : modulo_(mod)
{
    if (p.size() != 0) {
        dict_.resize(p.rbegin()->first + 1, integer_class(0));
        for (auto &iter : p) {
            integer_class temp;
            mp_fdiv_r(temp, iter.second, modulo_);
            dict_[iter.first] = temp;
        }
        gf_istrip();
    }
}

GaloisFieldDict::GaloisFieldDict(const integer_class &i,
                                 const integer_class &mod)
    : modulo_(mod)
{
    integer_class temp;
    mp_fdiv_r(temp, i, modulo_);
    if (temp != integer_class(0))
        dict_.insert(dict_.begin(), temp);
}

GaloisFieldDict GaloisFieldDict::from_vec(const std::vector<integer_class> &v,
                                          const integer_class &modulo)
{
    GaloisFieldDict x;
    x.modulo_ = modulo;
    x.dict_.resize(v.size());
    for (unsigned int i = 0; i < v.size(); ++i) {
        integer_class a;
        mp_fdiv_r(a, v[i], modulo);
        x.dict_[i] = a;
    }
    x.gf_istrip();
    return x;
}

GaloisFieldDict &GaloisFieldDict::negate()
{
    for (auto &a : dict_) {
        a *= -1;
        if (a != 0_z)
            a += modulo_;
    }
    return down_cast<GaloisFieldDict &>(*this);
}

void GaloisFieldDict::gf_istrip()
{
    for (auto i = dict_.size(); i-- != 0;) {
        if (dict_[i] == integer_class(0))
            dict_.pop_back();
        else
            break;
    }
}

GaloisFieldDict GaloisFieldDict::mul(const GaloisFieldDict &a,
                                     const GaloisFieldDict &b)
{
    if (a.modulo_ != b.modulo_)
        throw std::runtime_error("Error: field must be same.");
    if (a.get_dict().empty())
        return a;
    if (b.get_dict().empty())
        return b;

    GaloisFieldDict p;
    p.dict_.resize(a.degree() + b.degree() + 1, integer_class(0));
    p.modulo_ = a.modulo_;
    for (unsigned int i = 0; i <= a.degree(); i++)
        for (unsigned int j = 0; j <= b.degree(); j++) {
            auto temp = a.dict_[i];
            temp *= b.dict_[j];
            if (temp != integer_class(0)) {
                auto t = p.dict_[i + j];
                t += temp;
                mp_fdiv_r(t, t, a.modulo_);
                p.dict_[i + j] = t;
            }
        }
    p.gf_istrip();
    return p;
}

void GaloisFieldDict::gf_div(const GaloisFieldDict &o,
                             const Ptr<GaloisFieldDict> &quo,
                             const Ptr<GaloisFieldDict> &rem) const
{
    if (modulo_ != o.modulo_)
        throw SymEngineException("Error: field must be same.");
    if (o.dict_.empty())
        throw DivisionByZeroError("ZeroDivisionError");
    std::vector<integer_class> dict_out;
    if (dict_.empty()) {
        *quo = GaloisFieldDict::from_vec(dict_out, modulo_);
        *rem = GaloisFieldDict::from_vec(dict_, modulo_);
        return;
    }
    auto dict_divisor = o.dict_;
    auto deg_dividend = this->degree();
    auto deg_divisor = o.degree();
    if (deg_dividend < deg_divisor) {
        *quo = GaloisFieldDict::from_vec(dict_out, modulo_);
        *rem = GaloisFieldDict::from_vec(dict_, modulo_);
    } else {
        dict_out = dict_;
        integer_class inv;
        mp_invert(inv, *(dict_divisor.rbegin()), modulo_);
        integer_class coeff;
        for (auto it = deg_dividend + 1; it-- != 0;) {
            coeff = dict_out[it];
            auto lb = deg_divisor + it > deg_dividend
                          ? deg_divisor + it - deg_dividend
                          : 0;
            auto ub = std::min(it + 1, deg_divisor);
            for (size_t j = lb; j < ub; ++j) {
                mp_addmul(coeff, dict_out[it - j + deg_divisor],
                          -dict_divisor[j]);
            }
            if (it >= deg_divisor)
                coeff *= inv;
            mp_fdiv_r(coeff, coeff, modulo_);
            dict_out[it] = coeff;
        }
        std::vector<integer_class> dict_rem, dict_quo;
        dict_rem.resize(deg_divisor);
        dict_quo.resize(deg_dividend - deg_divisor + 1);
        for (unsigned it = 0; it < dict_out.size(); it++) {
            if (it < deg_divisor)
                dict_rem[it] = dict_out[it];
            else
                dict_quo[it - deg_divisor] = dict_out[it];
        }
        *quo = GaloisFieldDict::from_vec(dict_quo, modulo_);
        *rem = GaloisFieldDict::from_vec(dict_rem, modulo_);
    }
}

GaloisFieldDict GaloisFieldDict::gf_lshift(const integer_class n) const
{
    std::vector<integer_class> dict_out;
    auto to_ret = GaloisFieldDict::from_vec(dict_out, modulo_);
    if (!dict_.empty()) {
        auto n_val = mp_get_ui(n);
        to_ret.dict_.resize(n_val, integer_class(0));
        to_ret.dict_.insert(to_ret.dict_.end(), dict_.begin(), dict_.end());
    }
    return to_ret;
}

void GaloisFieldDict::gf_rshift(const integer_class n,
                                const Ptr<GaloisFieldDict> &quo,
                                const Ptr<GaloisFieldDict> &rem) const
{
    std::vector<integer_class> dict_quo;
    *quo = GaloisFieldDict::from_vec(dict_quo, modulo_);
    auto n_val = mp_get_ui(n);
    if (n_val < dict_.size()) {
        quo->dict_.insert(quo->dict_.end(), dict_.begin() + n_val, dict_.end());
        std::vector<integer_class> dict_rem(dict_.begin(),
                                            dict_.begin() + n_val);
        *rem = GaloisFieldDict::from_vec(dict_rem, modulo_);
    } else {
        *rem = down_cast<const GaloisFieldDict &>(*this);
    }
}

GaloisFieldDict GaloisFieldDict::gf_sqr() const
{
    return (*this * *this);
}

GaloisFieldDict GaloisFieldDict::gf_pow(const unsigned long n) const
{
    if (n == 0) {
        return GaloisFieldDict({integer_class(1)}, modulo_);
    }
    if (n == 1)
        return down_cast<const GaloisFieldDict &>(*this);
    if (n == 2)
        return gf_sqr();
    auto num = n;
    GaloisFieldDict to_sq = down_cast<const GaloisFieldDict &>(*this);
    GaloisFieldDict to_ret = GaloisFieldDict({integer_class(1)}, modulo_);
    while (1) {
        if (num & 1) {
            to_ret *= to_sq;
        }
        num >>= 1;
        if (num == 0)
            return to_ret;
        to_sq = to_sq.gf_sqr();
    }
}

void GaloisFieldDict::gf_monic(integer_class &res,
                               const Ptr<GaloisFieldDict> &monic) const
{
    *monic = down_cast<const GaloisFieldDict &>(*this);
    if (dict_.empty()) {
        res = integer_class(0);
    } else {
        res = *dict_.rbegin();
        if (res != integer_class(1)) {
            integer_class inv, temp;
            mp_invert(inv, res, modulo_);
            for (auto &iter : monic->dict_) {
                temp = inv;
                temp *= iter;
                mp_fdiv_r(iter, temp, modulo_);
            }
        }
    }
}

GaloisFieldDict GaloisFieldDict::gf_gcd(const GaloisFieldDict &o) const
{
    if (modulo_ != o.modulo_)
        throw SymEngineException("Error: field must be same.");
    GaloisFieldDict f = down_cast<const GaloisFieldDict &>(*this);
    GaloisFieldDict g = o;
    GaloisFieldDict temp_out;
    while (not g.dict_.empty()) {
        f %= g; // f, g = f % g, g
        f.dict_.swap(g.dict_);
    }
    integer_class temp_LC;
    f.gf_monic(temp_LC, outArg(f));
    return f;
}

GaloisFieldDict GaloisFieldDict::gf_lcm(const GaloisFieldDict &o) const
{
    if (modulo_ != o.modulo_)
        throw SymEngineException("Error: field must be same.");
    if (dict_.empty())
        return down_cast<const GaloisFieldDict &>(*this);
    if (o.dict_.empty())
        return o;
    GaloisFieldDict out, temp_out;
    out = o * (*this);
    out /= gf_gcd(o);
    integer_class temp_LC;
    out.gf_monic(temp_LC, outArg(out));
    return out;
}

GaloisFieldDict GaloisFieldDict::gf_diff() const
{
    auto df = degree();
    GaloisFieldDict out = GaloisFieldDict({}, modulo_);
    out.dict_.resize(df, integer_class(0));
    for (unsigned i = 1; i <= df; i++) {
        if (dict_[i] != integer_class(0)) {
            out.dict_[i - 1] = i * dict_[i];
            mp_fdiv_r(out.dict_[i - 1], out.dict_[i - 1], modulo_);
        }
    }
    out.gf_istrip();
    return out;
}

integer_class GaloisFieldDict::gf_eval(const integer_class &a) const
{
    integer_class res = 0_z;
    for (auto rit = dict_.rbegin(); rit != dict_.rend(); ++rit) {
        res *= a;
        res += (*rit);
        res %= modulo_;
    }
    return res;
}

vec_integer_class
GaloisFieldDict::gf_multi_eval(const vec_integer_class &v) const
{
    vec_integer_class res(v.size());
    for (unsigned int i = 0; i < v.size(); ++i)
        res[i] = gf_eval(v[i]);
    return res;
}

bool GaloisFieldDict::gf_is_sqf() const
{
    if (dict_.empty())
        return true;
    integer_class LC;
    GaloisFieldDict monic;
    gf_monic(LC, outArg(monic));
    monic = monic.gf_gcd(monic.gf_diff());
    return monic.is_one();
}

std::vector<std::pair<GaloisFieldDict, unsigned>>
GaloisFieldDict::gf_sqf_list() const
{
    std::vector<std::pair<GaloisFieldDict, unsigned>> vec_out;
    if (degree() < 1)
        return vec_out;
    unsigned n = 1;
    // This cast is okay, because the multiplicities are unsigned
    unsigned r = numeric_cast<unsigned>(mp_get_ui(modulo_));
    bool sqf = false;
    integer_class LC;
    GaloisFieldDict f;
    gf_monic(LC, outArg(f));
    while (true) {
        GaloisFieldDict F = f.gf_diff();
        if (not F.dict_.empty()) {
            GaloisFieldDict g = f.gf_gcd(F);
            GaloisFieldDict h = f / g;

            unsigned i = 1;

            while (not h.is_one()) {
                GaloisFieldDict G = h.gf_gcd(g);
                GaloisFieldDict H = h / G;

                if (H.degree() > 0)
                    vec_out.push_back({H, i * n});

                ++i;
                g /= G;
                h = G;
            }
            if (g.is_one())
                sqf = true;
            else
                f = g;
        }
        if (not sqf) {
            auto deg = f.degree();
            auto d = deg / r;
            GaloisFieldDict temp = f;
            for (unsigned int i = 0; i <= d; ++i) {
                f.dict_[d - i] = temp.dict_[deg - i * r];
            }
            n *= r;
            f.dict_.resize(d + 1);
            f.gf_istrip();
        } else
            break;
    }
    return vec_out;
}

GaloisFieldDict GaloisFieldDict::gf_sqf_part() const
{
    auto sqf = gf_sqf_list();
    GaloisFieldDict g = GaloisFieldDict::from_vec({1_z}, modulo_);

    for (auto &f : sqf)
        g *= f.first;

    return g;
}

GaloisFieldDict GaloisFieldDict::gf_compose_mod(const GaloisFieldDict &g,
                                                const GaloisFieldDict &h) const
{
    if (g.modulo_ != h.modulo_)
        throw SymEngineException("Error: field must be same.");
    if (g.modulo_ != modulo_)
        throw SymEngineException("Error: field must be same.");
    if (g.dict_.size() == 0)
        return g;
    GaloisFieldDict out
        = GaloisFieldDict::from_vec({*(g.dict_.rbegin())}, modulo_);
    if (g.dict_.size() >= 2) {
        for (auto i = g.dict_.size() - 2;; --i) {
            out *= h;
            out += g.dict_[i];
            out %= (*this);
            if (i == 0)
                break;
        }
    }
    return out;
}

GaloisFieldDict GaloisFieldDict::gf_pow_mod(const GaloisFieldDict &f,
                                            const unsigned long &n) const
{
    if (modulo_ != f.modulo_)
        throw SymEngineException("Error: field must be same.");
    if (n == 0)
        return GaloisFieldDict::from_vec({1_z}, modulo_);
    GaloisFieldDict in = f;
    if (n == 1) {
        return f % (*this);
    }
    if (n == 2) {
        return f.gf_sqr() % (*this);
    }
    GaloisFieldDict h = GaloisFieldDict::from_vec({1_z}, modulo_);
    auto mod = n;
    while (true) {
        if (mod & 1) {
            h *= in;
            h %= *this;
        }
        mod >>= 1;

        if (mod == 0)
            break;

        in = in.gf_sqr() % *this;
    }
    return h;
}

std::vector<GaloisFieldDict> GaloisFieldDict::gf_frobenius_monomial_base() const
{
    auto n = degree();
    std::vector<GaloisFieldDict> b;
    if (n == 0)
        return b;
    b.resize(n);
    b[0] = GaloisFieldDict::from_vec({1_z}, modulo_);
    GaloisFieldDict temp_out;
    if (mp_get_ui(modulo_) < n) {
        for (unsigned i = 1; i < n; ++i) {
            b[i] = b[i - 1].gf_lshift(modulo_);
            b[i] %= (*this);
        }
    } else if (n > 1) {
        b[1] = gf_pow_mod(GaloisFieldDict::from_vec({0_z, 1_z}, modulo_),
                          mp_get_ui(modulo_));
        for (unsigned i = 2; i < n; ++i) {
            b[i] = b[i - 1] * b[1];
            b[i] %= (*this);
        }
    }
    return b;
}

GaloisFieldDict
GaloisFieldDict::gf_frobenius_map(const GaloisFieldDict &g,
                                  const std::vector<GaloisFieldDict> &b) const
{
    if (modulo_ != g.modulo_)
        throw SymEngineException("Error: field must be same.");
    auto m = g.degree();
    GaloisFieldDict temp_out(*this), out;
    if (this->degree() >= m) {
        temp_out %= g;
    }
    if (temp_out.empty()) {
        return temp_out;
    }
    m = temp_out.degree();
    out = GaloisFieldDict::from_vec({temp_out.dict_[0]}, modulo_);
    for (unsigned i = 1; i <= m; ++i) {
        auto v = b[i];
        v *= temp_out.dict_[i];
        out += v;
    }
    out.gf_istrip();
    return out;
}

std::pair<GaloisFieldDict, GaloisFieldDict> GaloisFieldDict::gf_trace_map(
    const GaloisFieldDict &a, const GaloisFieldDict &b,
    const GaloisFieldDict &c, const unsigned long &n) const
{
    unsigned long n_val(n);
    auto u = this->gf_compose_mod(a, b);
    GaloisFieldDict v(b), U, V;
    if (n_val & 1) {
        U = a + u;
        V = b;
    } else {
        U = a;
        V = c;
    }
    n_val >>= 1;
    while (n_val) {
        u += this->gf_compose_mod(u, v);
        v = gf_compose_mod(v, v);

        if (n_val & 1) {
            auto temp = gf_compose_mod(u, V);
            U += temp;
            V = gf_compose_mod(v, V);
        }
        n_val >>= 1;
    }
    return std::make_pair(gf_compose_mod(a, V), U);
}

GaloisFieldDict
GaloisFieldDict::_gf_trace_map(const GaloisFieldDict &f, const unsigned long &n,
                               const std::vector<GaloisFieldDict> &b) const
{
    GaloisFieldDict x = f % (*this);
    auto h = f;
    auto r = f;
    for (unsigned i = 1; i < n; ++i) {
        h = gf_frobenius_map(h, b);
        r += h;
        r %= (*this);
    }
    return r;
}

std::vector<std::pair<GaloisFieldDict, unsigned>>
GaloisFieldDict::gf_ddf_zassenhaus() const
{
    unsigned i = 1;
    GaloisFieldDict f(*this);
    GaloisFieldDict g = GaloisFieldDict::from_vec({0_z, 1_z}, modulo_);
    GaloisFieldDict to_sub(g);
    std::vector<std::pair<GaloisFieldDict, unsigned>> factors;

    auto b = f.gf_frobenius_monomial_base();
    while (2 * i <= f.degree()) {
        g = g.gf_frobenius_map(f, b);

        GaloisFieldDict h = f.gf_gcd(g - to_sub);

        if (not h.is_one()) {
            factors.push_back({h, i});
            f /= h;
            g %= f;
            b = f.gf_frobenius_monomial_base();
        }
        ++i;
    }
    if (not(f.is_one() || f.empty())) {
        factors.push_back({f, f.degree()});
    }
    return factors;
}

GaloisFieldDict
GaloisFieldDict::_gf_pow_pnm1d2(const GaloisFieldDict &f, const unsigned &n,
                                const std::vector<GaloisFieldDict> &b) const
{
    GaloisFieldDict f_in(f);
    f_in %= *this;
    GaloisFieldDict h, r;
    h = r = f_in;
    for (unsigned i = 1; i < n; ++i) {
        h = h.gf_frobenius_map(*this, b);
        r *= h;
        r %= *this;
    }
    auto res = gf_pow_mod(r, (mp_get_ui(modulo_) - 1) / 2);
    return res;
}

GaloisFieldDict GaloisFieldDict::gf_random(const unsigned int &n_val,
                                           mp_randstate &state) const
{
    std::vector<integer_class> v(n_val + 1);
    for (unsigned i = 0; i < n_val; ++i) {
        state.urandomint(v[i], modulo_);
    }
    v[n_val] = 1_z;
    return GaloisFieldDict::from_vec(v, modulo_);
}

std::set<GaloisFieldDict, GaloisFieldDict::DictLess>
GaloisFieldDict::gf_edf_zassenhaus(const unsigned &n) const
{
    std::set<GaloisFieldDict, DictLess> factors;
    factors.insert(*this);
    if (this->degree() <= n)
        return factors;

    auto N = this->degree() / n;

    std::vector<GaloisFieldDict> b;
    if (modulo_ != 2_z)
        b = this->gf_frobenius_monomial_base();
    mp_randstate state;
    while (factors.size() < N) {
        auto r = gf_random(2 * n - 1, state);
        GaloisFieldDict g;
        if (modulo_ == 2_z) {
            GaloisFieldDict h = r;
            unsigned ub = 1 << (n * N - 1);
            for (unsigned i = 0; i < ub; ++i) {
                r = gf_pow_mod(r, 2);
                h += r;
            }
            g = this->gf_gcd(h);
        } else {
            GaloisFieldDict h = _gf_pow_pnm1d2(r, n, b);
            h -= 1_z;
            g = this->gf_gcd(h);
        }

        if (!g.is_one() and g != (*this)) {
            factors = g.gf_edf_zassenhaus(n);
            auto to_add = ((*this) / g).gf_edf_zassenhaus(n);
            if (not to_add.empty())
                factors.insert(to_add.begin(), to_add.end());
        }
    }
    return factors;
}

std::vector<std::pair<GaloisFieldDict, unsigned>>
GaloisFieldDict::gf_ddf_shoup() const
{
    std::vector<std::pair<GaloisFieldDict, unsigned>> factors;
    if (dict_.empty())
        return factors;
    GaloisFieldDict f(*this);
    auto n = this->degree();
    auto k = static_cast<unsigned>(std::ceil(std::sqrt(n / 2)));
    auto b = gf_frobenius_monomial_base();
    auto x = GaloisFieldDict::from_vec({0_z, 1_z}, modulo_);
    auto h = x.gf_frobenius_map(f, b);

    std::vector<GaloisFieldDict> U;
    U.push_back(x);
    U.push_back(h);
    U.resize(k + 1);
    for (unsigned i = 2; i <= k; ++i)
        U[i] = U[i - 1].gf_frobenius_map(*this, b);
    h = U[k];
    U.resize(k);
    std::vector<GaloisFieldDict> V;
    V.push_back(h);
    V.resize(k);
    for (unsigned i = 1; i + 1 <= k; ++i)
        V[i] = this->gf_compose_mod(V[i - 1], h);
    for (unsigned i = 0; i < V.size(); i++) {
        h = GaloisFieldDict::from_vec({1_z}, modulo_);
        auto j = k - 1;
        GaloisFieldDict g;
        for (auto &u : U) {
            g = V[i] - u;
            h *= g;
            h %= f;
        }
        g = f.gf_gcd(h);
        f /= g;
        for (auto rit = U.rbegin(); rit != U.rend(); ++rit) {
            h = V[i] - (*rit);
            auto F = g.gf_gcd(h);
            if (not F.is_one()) {
                unsigned temp = k * (i + 1) - j;
                factors.push_back({F, temp});
            }
            g /= F;
            --j;
        }
    }
    if (not f.is_one())
        factors.push_back({f, f.degree()});
    return factors;
}

std::set<GaloisFieldDict, GaloisFieldDict::DictLess>
GaloisFieldDict::gf_edf_shoup(const unsigned &n) const
{
    auto N = this->degree();
    std::set<GaloisFieldDict, DictLess> factors;
    if (N <= n) {
        if (N != 0)
            factors.insert(*this);
        return factors;
    }
    auto x = GaloisFieldDict::from_vec({0_z, 1_z}, modulo_);
    mp_randstate state;
    auto r = gf_random(N - 1, state);
    if (modulo_ == 2_z) {
        auto h = gf_pow_mod(x, mp_get_ui(modulo_));
        auto H = gf_trace_map(r, h, x, n - 1).second;
        auto h1 = gf_gcd(H);
        auto h2 = (*this) / h1;
        factors = h1.gf_edf_shoup(n);
        auto temp = h2.gf_edf_shoup(n);
        factors.insert(temp.begin(), temp.end());
    } else {
        auto b = gf_frobenius_monomial_base();
        auto H = _gf_trace_map(r, n, b);
        auto h = gf_pow_mod(H, (mp_get_ui(modulo_) - 1) / 2);
        auto h1 = gf_gcd(h);
        auto h2 = gf_gcd(h - 1_z);
        auto h3 = (*this) / (h1 * h2);
        factors = h1.gf_edf_shoup(n);
        auto temp = h2.gf_edf_shoup(n);
        factors.insert(temp.begin(), temp.end());
        temp = h3.gf_edf_shoup(n);
        factors.insert(temp.begin(), temp.end());
    }
    return factors;
}

std::set<GaloisFieldDict, GaloisFieldDict::DictLess>
GaloisFieldDict::gf_zassenhaus() const
{
    std::set<GaloisFieldDict, DictLess> factors;
    auto temp1 = gf_ddf_zassenhaus();
    for (auto &f : temp1) {
        auto temp2 = f.first.gf_edf_zassenhaus(f.second);
        factors.insert(temp2.begin(), temp2.end());
    }
    return factors;
}

std::set<GaloisFieldDict, GaloisFieldDict::DictLess>
GaloisFieldDict::gf_shoup() const
{
    std::set<GaloisFieldDict, DictLess> factors;
    auto temp1 = gf_ddf_shoup();
    for (auto &f : temp1) {
        auto temp2 = f.first.gf_edf_shoup(f.second);
        factors.insert(temp2.begin(), temp2.end());
    }
    return factors;
}

std::pair<integer_class, std::set<std::pair<GaloisFieldDict, unsigned>,
                                  GaloisFieldDict::DictLess>>
GaloisFieldDict::gf_factor() const
{
    integer_class lc;
    std::set<std::pair<GaloisFieldDict, unsigned>, DictLess> factors;
    GaloisFieldDict monic;
    gf_monic(lc, outArg(monic));
    if (monic.degree() < 1)
        return std::make_pair(lc, factors);
    std::vector<std::pair<GaloisFieldDict, unsigned>> sqf_list
        = monic.gf_sqf_list();
    for (auto a : sqf_list) {
        auto temp = (a.first).gf_zassenhaus();
        for (auto f : temp)
            factors.insert({f, a.second});
    }
    return std::make_pair(lc, factors);
}
}