[WIP] Hensel lift and zassenhaus

symengine/fields.cpp

@@ -316,10 +316,10 @@ GaloisFieldDict GaloisFieldDict::gf_pow(const unsigned int n) const
     }
 }
 
-void GaloisFieldDict::gf_monic(integer_class &res,
-                               const Ptr<GaloisFieldDict> &monic) const
+integer_class GaloisFieldDict::gf_monic(const Ptr<GaloisFieldDict> &monic) const
 {
     *monic = static_cast<GaloisFieldDict>(*this);
+    integer_class res;
     if (dict_.empty()) {
         res = integer_class(0);
     } else {
@@ -334,6 +334,7 @@ void GaloisFieldDict::gf_monic(integer_class &res,
             }
         }
     }
+    return res;
 }
 
 GaloisFieldDict GaloisFieldDict::gf_gcd(const GaloisFieldDict &o) const
@@ -347,11 +348,65 @@ GaloisFieldDict GaloisFieldDict::gf_gcd(const GaloisFieldDict &o) const
         f %= g; // f, g = f % g, g
         f.dict_.swap(g.dict_);
     }
-    integer_class temp_LC;
-    f.gf_monic(temp_LC, outArg(f));
+    f.gf_monic(outArg(f));
     return f;
 }
 
+void GaloisFieldDict::gf_gcdex(const GaloisFieldDict &f,
+                               const GaloisFieldDict &g,
+                               const Ptr<GaloisFieldDict> &s,
+                               const Ptr<GaloisFieldDict> &t,
+                               const Ptr<GaloisFieldDict> &h)
+{
+    if (f.modulo_ != g.modulo_)
+        throw std::runtime_error("Error: field must be same.");
+    s->dict_.clear();
+    t->dict_.clear();
+    h->dict_.clear();
+    s->modulo_ = t->modulo_ = h->modulo_ = f.modulo_;
+    if (f.dict_.empty() and g.dict_.empty()) {
+        s->dict_.push_back(1_z);
+        return;
+    }
+    integer_class p0, p1;
+    GaloisFieldDict r0;
+    p0 = f.gf_monic(outArg(r0));
+    p1 = g.gf_monic(outArg(*h));
+    if (f.dict_.empty()) {
+        mp_invert(p1, p1, f.modulo_);
+        t->dict_.push_back(p1);
+        return;
+    }
+    if (g.dict_.empty()) {
+        mp_invert(p0, p0, f.modulo_);
+        s->dict_.push_back(p0);
+        *h = r0;
+        return;
+    }
+    mp_invert(p0, p0, f.modulo_);
+    mp_invert(p1, p1, f.modulo_);
+    GaloisFieldDict s0, t0;
+    s0.modulo_ = (*s).modulo_ = t0.modulo_ = (*t).modulo_ = f.modulo_;
+    s0.dict_.push_back(p0);
+    (*t).dict_.push_back(p1);
+    while (1) {
+        GaloisFieldDict Q, R;
+        r0.gf_div((*h), outArg(Q), outArg(R));
+
+        if (R.dict_.empty())
+            break;
+        r0 = *h;
+        integer_class lc;
+        lc = R.gf_monic(outArg(*h));
+        mp_invert(lc, lc, f.modulo_);
+
+        s0 = (s0 - (*s) * Q) * lc;
+        std::swap(s0, (*s));
+        t0 = (t0 - (*t) * Q) * lc;
+        std::swap(t0, (*t));
+    }
+}
+
 GaloisFieldDict GaloisFieldDict::gf_lcm(const GaloisFieldDict &o) const
 {
     if (modulo_ != o.modulo_)
@@ -363,8 +418,7 @@ GaloisFieldDict GaloisFieldDict::gf_lcm(const GaloisFieldDict &o) const
     GaloisFieldDict out, temp_out;
     out = o * (*this);
     out /= gf_gcd(o);
-    integer_class temp_LC;
-    out.gf_monic(temp_LC, outArg(out));
+    out.gf_monic(outArg(out));
     return out;
 }
 
@@ -407,9 +461,8 @@ bool GaloisFieldDict::gf_is_sqf() const
 {
     if (dict_.empty())
         return true;
-    integer_class LC;
     GaloisFieldDict monic;
-    gf_monic(LC, outArg(monic));
+    gf_monic(outArg(monic));
     monic = monic.gf_gcd(monic.gf_diff());
     return monic.is_one();
 }
@@ -423,9 +476,8 @@ GaloisFieldDict::gf_sqf_list() const
     unsigned n = 1;
     unsigned r = mp_get_ui(modulo_);
     bool sqf = false;
-    integer_class LC;
     GaloisFieldDict f;
-    gf_monic(LC, outArg(f));
+    gf_monic(outArg(f));
     while (true) {
         GaloisFieldDict F = f.gf_diff();
         if (not F.dict_.empty()) {
@@ -852,7 +904,7 @@ GaloisFieldDict::gf_factor() const
     integer_class lc;
     std::set<std::pair<GaloisFieldDict, unsigned>, DictLess> factors;
     GaloisFieldDict monic;
-    gf_monic(lc, outArg(monic));
+    lc = gf_monic(outArg(monic));
     if (monic.degree() < 1)
         return std::make_pair(lc, factors);
     std::vector<std::pair<GaloisFieldDict, unsigned>> sqf_list
@@ -864,4 +916,123 @@ GaloisFieldDict::gf_factor() const
     }
     return std::make_pair(lc, factors);
 }
+
+void GaloisFieldDict::zz_hensel_step(
+    const integer_class &m, const GaloisFieldDict &ff, GaloisFieldDict &g,
+    GaloisFieldDict &h, GaloisFieldDict &s, GaloisFieldDict &t,
+    const Ptr<GaloisFieldDict> &G, const Ptr<GaloisFieldDict> &H,
+    const Ptr<GaloisFieldDict> &S, const Ptr<GaloisFieldDict> &T)
+{
+    GaloisFieldDict q, r;
+    integer_class M;
+
+    mp_pow_ui(M, m, 2);
+    auto f = GaloisFieldDict::from_vec(ff.dict_, M);
+    g.modulo_ = M;
+    h.modulo_ = M;
+    s.modulo_ = M;
+    t.modulo_ = M;
+
+    auto e = f - (g * h);
+
+    (s * e).gf_div(h, outArg(q), outArg(r));
+    auto u = (t * e) + (q * g);
+    *G = (g + u);
+    *H = (h + r);
+
+    u = s * (*G) + t * (*H);
+    u.dict_[0] -= 1_z;
+    (s * u).gf_div(*H, outArg(q), outArg(r));
+    u = t * u + q * (*G);
+    *S = (s - r);
+    *T = (t - u);
+}
+
+integer_class GaloisFieldDict::get_lc() const
+{
+    if (dict_.empty())
+        return 0_z;
+    return *(dict_.rbegin());
+}
+
+void GaloisFieldDict::itrunc()
+{
+    integer_class modulo_2 = modulo_ / 2_z;
+    for (auto it = dict_.begin(); it != dict_.end(); ++it) {
+        if (*it > modulo_2)
+            *it -= modulo_;
+    }
+}
+
+std::vector<UIntDict>
+GaloisFieldDict::zz_hensel_lift(const UIntDict &f, const integer_class &p,
+                                const std::vector<UIntDict> &f_list,
+                                unsigned int l)
+{
+    integer_class pl;
+    mp_pow_ui(pl, p, l);
+    GaloisFieldDict ff(f.dict_, pl);
+    std::set<GaloisFieldDict, DictLess> s;
+    for (auto &a : f_list) {
+        s.insert(GaloisFieldDict(a.dict_, p));
+    }
+    return GaloisFieldDict::zz_hensel_lift(ff, p, s, l);
+}
+
+std::vector<UIntDict> GaloisFieldDict::zz_hensel_lift(
+    const GaloisFieldDict &f, const integer_class &p,
+    const std::set<GaloisFieldDict, GaloisFieldDict::DictLess> &f_list,
+    unsigned int l)
+{
+    std::vector<UIntDict> res;
+    size_t r = f_list.size();
+    auto lc = f.get_lc();
+
+    if (r == 1) {
+        integer_class g, r, s, pl;
+        mp_pow_ui(pl, p, l);
+        mp_gcdext(g, r, s, lc, pl);
+        UIntDict F = UIntDict::from_vec(f.dict_);
+        F *= r;
+        F.itrunc(pl);
+        res.push_back(F);
+        return res;
+    }
+
+    integer_class m(p);
+    size_t k(r / 2);
+    auto d = std::ceil(std::log2(l));
+    auto g = GaloisFieldDict::from_vec({lc}, p);
+    auto it_k = std::next(f_list.begin(), k);
+
+    for (auto iter = f_list.begin(); iter != it_k; ++iter) {
+        g *= (*iter);
+    }
+
+    auto h = GaloisFieldDict::from_vec({1_z}, p);
+
+    for (auto iter = it_k; iter != f_list.end(); ++iter) {
+        h *= (*iter);
+    }
+
+    std::set<GaloisFieldDict, GaloisFieldDict::DictLess> sub1(f_list.begin(),
+                                                              it_k),
+        sub2(it_k, f_list.end());
+
+    GaloisFieldDict temp2, s, t;
+    GaloisFieldDict::gf_gcdex(g, h, outArg(s), outArg(t), outArg(temp2));
+    g.itrunc();
+    h.itrunc();
+    s.itrunc();
+    t.itrunc();
+    for (unsigned int i = 1; i <= d; ++i) {
+        GaloisFieldDict::zz_hensel_step(m, f, g, h, s, t, outArg(g), outArg(h),
+                                        outArg(s), outArg(t));
+        mp_pow_ui(m, m, 2);
+    }
+    res = zz_hensel_lift(g, p, sub1, l);
+    auto temp = zz_hensel_lift(h, p, sub2, l);
+    res.insert(res.end(), temp.begin(), temp.end());
+    return res;
+}
 }

symengine/fields.h

@@ -65,8 +65,12 @@ class GaloisFieldDict
                    const Ptr<GaloisFieldDict> &rem) const;
     GaloisFieldDict gf_sqr() const;
     GaloisFieldDict gf_pow(const unsigned int n) const;
-    void gf_monic(integer_class &res, const Ptr<GaloisFieldDict> &monic) const;
+    integer_class gf_monic(const Ptr<GaloisFieldDict> &monic) const;
     GaloisFieldDict gf_gcd(const GaloisFieldDict &o) const;
+    static void gf_gcdex(const GaloisFieldDict &f, const GaloisFieldDict &g,
+                         const Ptr<GaloisFieldDict> &s,
+                         const Ptr<GaloisFieldDict> &t,
+                         const Ptr<GaloisFieldDict> &h);
     GaloisFieldDict gf_lcm(const GaloisFieldDict &o) const;
     GaloisFieldDict gf_diff() const;
     integer_class gf_eval(const integer_class &a) const;
@@ -146,6 +150,30 @@ class GaloisFieldDict
     std::pair<integer_class,
               std::set<std::pair<GaloisFieldDict, unsigned>, DictLess>>
     gf_factor() const;
+    // One step in Hensel lifting in `Z[x]`.
+    // References :
+    //     1.) J. von zur Gathen, J. Gerhard, Modern Computer Algebra, 1999,
+    //     page no: 445-446
+    static void
+    zz_hensel_step(const integer_class &m, const GaloisFieldDict &ff,
+                   GaloisFieldDict &g, GaloisFieldDict &h, GaloisFieldDict &s,
+                   GaloisFieldDict &t, const Ptr<GaloisFieldDict> &G,
+                   const Ptr<GaloisFieldDict> &H, const Ptr<GaloisFieldDict> &S,
+                   const Ptr<GaloisFieldDict> &T);
+    // Multifactor Hensel lifting in `Z[x]`.
+    // References :
+    //     1.) J. von zur Gathen, J. Gerhard, Modern Computer Algebra, 1999,
+    //     page no: 450-451
+    static std::vector<UIntDict> zz_hensel_lift(
+        const GaloisFieldDict &f, const integer_class &p,
+        const std::set<GaloisFieldDict, GaloisFieldDict::DictLess> &f_list,
+        unsigned int l);
+    std::vector<UIntDict> zz_hensel_lift(const UIntDict &f,
+                                         const integer_class &p,
+                                         const std::vector<UIntDict> &f_list,
+                                         unsigned int l);
+    integer_class get_lc() const;
+    void itrunc();
 
     GaloisFieldDict &operator=(GaloisFieldDict &&other) SYMENGINE_NOEXCEPT
     {
@@ -288,10 +316,12 @@ class GaloisFieldDict
     static GaloisFieldDict mul(const GaloisFieldDict &a,
                                const GaloisFieldDict &b);
 
-    friend GaloisFieldDict operator*(const GaloisFieldDict &a,
-                                     const GaloisFieldDict &b)
+    template <class T>
+    friend GaloisFieldDict operator*(const GaloisFieldDict &a, const T &b)
     {
-        return GaloisFieldDict::mul(a, b);
+        GaloisFieldDict c = a;
+        c *= b;
+        return c;
     }
 
     GaloisFieldDict &operator*=(const integer_class &other)