[WIP] Hensel lift and zassenhaus #1066

Open
wants to merge 14 commits into
from

Projects

None yet

4 participants

@nishnik
Contributor
nishnik commented Aug 17, 2016

It factors a primitive polynomial in integer domain.
Implements extended GCD in finite fields.
And a small bug fix in UIntPoly::mul

@isuruf isuruf commented on an outdated diff Aug 17, 2016
symengine/polys/uintpoly.h
@@ -120,6 +124,19 @@ class UIntDict : public ODictWrapper<unsigned int, integer_class, UIntDict>
return curr;
}
+ void itrunc(const integer_class &mod);
+ UIntDict primitive() const;
+ integer_class l1_norm() const;
+ static void zz_hensel_step(const integer_class &m, const UIntDict &f,
+ const UIntDict &g, const UIntDict &h,
+ const UIntDict &s, const UIntDict &t,
+ const Ptr<UIntDict> &G, const Ptr<UIntDict> &H,
+ const Ptr<UIntDict> &S, const Ptr<UIntDict> &T);
+ std::vector<UIntDict> zz_hensel_lift(const integer_class &p,
+ const std::vector<UIntDict> &f_list,
+ unsigned int l) const;
+ static void zz_divide(const UIntDict &a, const UIntDict &b,
@isuruf
isuruf Aug 17, 2016 Member

Can you describe each of the methods above as a comment?

@certik certik commented on the diff Aug 17, 2016
symengine/tests/polynomial/test_uintpoly.cpp
+ REQUIRE(find_in_set(UIntPoly::from_vec(x, {1_z, -3_z}), out));
+ REQUIRE(find_in_set(UIntPoly::from_vec(x, {-1_z, -3_z}), out));
+ REQUIRE(out.size() == 2);
+
+ f = UIntPoly::from_vec(x, {-6_z, 11_z, -6_z, 1_z});
+ out = f->zz_zassenhaus();
+ REQUIRE(find_in_set(UIntPoly::from_vec(x, {-1_z, 1_z}), out));
+ REQUIRE(find_in_set(UIntPoly::from_vec(x, {-2_z, 1_z}), out));
+ REQUIRE(find_in_set(UIntPoly::from_vec(x, {-3_z, 1_z}), out));
+ REQUIRE(out.size() == 3);
+
+ f = UIntPoly::from_vec(x, {10_z, 13_z, 10_z, 3_z});
+ out = f->zz_zassenhaus();
+ REQUIRE(find_in_set(UIntPoly::from_vec(x, {2_z, 1_z}), out));
+ REQUIRE(find_in_set(UIntPoly::from_vec(x, {5_z, 4_z, 3_z}), out));
+ REQUIRE(out.size() == 2);
@certik
certik Aug 17, 2016 edited Contributor

Does the zz_zassenhaus give you a set of factors? E.g. the above 3 tests correspond to:

In [1]: factor(x**2-9)
Out[1]: (x - 3)⋅(x + 3)

In [2]: factor(-6*x**3+11*x**2-6*x+1)
Out[2]: -(x - 1)⋅(2⋅x - 1)⋅(3⋅x - 1)

In [3]: factor(10*x**3+13*x**2+10*x+3)
Out[3]: 
          ⎛   2          ⎞
(2⋅x + 1)⋅⎝5⋅x  + 4⋅x + 3⎠

It looks like all the factors are there, except in the first case, where you get the wrong overall sign, i.e. in your results
(x-3)*(-x-3) = -x^2+9, but you are factoring x^2-9.

@nishnik
nishnik Aug 17, 2016 Contributor

@certik This is not the actual factorization, zz_zassenhaus if given a square free polynomial with negative largest coefficient, will factor the negative of it (due to finite field algorithms).
So when we give -9x^2 + 1, it actually factors 9x^2 - 1.
Note that in SymEngine we are using a different convention that is {1_z, 0_z, -9_z} means -9x^2 + 1

@certik
Contributor
certik commented Aug 17, 2016

Otherwise it looks like it's working. Great job!

What needs to be done to finish this PR?

@nishnik
Contributor
nishnik commented Aug 17, 2016

@certik I will have to write the factor square free method then this can be merged.
After which in a different PR, I would like to implement the get square free part and factor method.

@isuruf isuruf commented on an outdated diff Aug 17, 2016
symengine/fields.cpp
+ if (R.dict_.empty())
+ break;
+ r0 = r1;
+ integer_class lc, inv;
+ R.gf_monic(lc, outArg(r1));
+ mp_invert(inv, lc, f.modulo_);
+
+ *s = s0 - s1 * Q;
+ *t = t0 - t1 * Q;
+
+ s0 = s1;
+ s1 = *s;
+ s1 *= inv;
+ t0 = t1;
+ t1 = *t;
+ t1 *= inv;
@isuruf
isuruf Aug 17, 2016 Member

This is making copies too much. Can you reduce copying here?

@isuruf isuruf commented on an outdated diff Aug 17, 2016
symengine/polys/uintpoly.cpp
+
+integer_class UIntDict::l1_norm() const
+{
+ integer_class out(0_z);
+ for (auto &a : dict_) {
+ out += mp_abs(a.second);
+ }
+ return out;
+}
+
+UIntDict UIntDict::primitive() const
+{
+ integer_class gcd(0_z);
+ UIntDict out(dict_);
+ for (auto &a : dict_) {
+ mp_gcd(gcd, gcd, mp_abs(a.second));
@isuruf
isuruf Aug 17, 2016 Member

You don't need mp_abs

@isuruf isuruf and 1 other commented on an outdated diff Aug 17, 2016
symengine/polys/uintpoly.cpp
+
+UIntDict UIntDict::primitive() const
+{
+ integer_class gcd(0_z);
+ UIntDict out(dict_);
+ for (auto &a : dict_) {
+ mp_gcd(gcd, gcd, mp_abs(a.second));
+ }
+ for (auto &a : out.dict_) {
+ a.second /= gcd;
+ }
+ return out;
+}
+
+std::vector<UIntDict>
+UIntDict::zz_hensel_lift(const integer_class &p,
@isuruf
isuruf Aug 17, 2016 Member

Isn't this supposed to be in GaloisField because it uses a mod ?

@nishnik
nishnik Aug 17, 2016 Contributor

Yes.
And I had read this algo from a source on internet, but after reading it from Modern Computer Algebra I am sure, I can make some more tweaks.

@isuruf isuruf and 2 others commented on an outdated diff Aug 17, 2016
symengine/polys/uintpoly.cpp
+}
+
+std::set<RCP<const UIntPoly>, RCPBasicKeyLess> UIntPoly::zz_zassenhaus() const
+{
+ std::set<RCP<const UIntPoly>, RCPBasicKeyLess> factors;
+ unsigned int n = get_degree();
+ UIntDict f(this->poly_);
+ if (n == 1) {
+ factors.insert(rcp_from_this_cast<const UIntPoly>());
+ return factors;
+ }
+ integer_class fc(this->poly_.dict_.begin()->second);
+ auto A = mp_get_ui(this->poly_.max_abs_coef());
+ auto b = get_lc();
+ auto B
+ = std::abs(int(std::sqrt(n + 1)) * std::pow(2, n) * A * mp_get_si(b));
@isuruf
isuruf Aug 17, 2016 Member

Are you sure mp_get_si won't overflow?
Also, when using doubles, you lose precision right?

@nishnik
nishnik Aug 18, 2016 Contributor

It might not fit, but we have a bottleneck here.
In the algorithm we have:

    auto l = std::ceil(std::log(2 * B + 1) / std::log(mp_get_ui(fsqf.first)));
    integer_class pl;
    mp_pow_ui(pl, fsqf.first, l);

as we have to use pl as modulo_ that means we need a similar sized dict_.
and pl is approximately equal to 2*B.
B = std::abs(int(std::sqrt(n + 1)) * std::pow(2, n) * A * mp_get_si(b));
and A is already greater than b, so I think we can use a signed int here, because if it is greater than a signed int, pl would also be and we have a bottleneck due to vector's size.

@isuruf
isuruf Aug 18, 2016 Member

Okay. then you should check that b fits a signed integer and if not throw an error. Also, what about the precision issue?

@nishnik
nishnik Aug 18, 2016 Contributor

We need the square root in integer domain, so thats not an issue.

@isuruf
isuruf Aug 18, 2016 Member

I meant using std::pow(2, n). This can overflow.

@isuruf
isuruf Aug 18, 2016 Member

Where is the dictionary of size pl made?

@nishnik
nishnik Aug 18, 2016 Contributor

When we do hensel lifting.
Here,
we convert the polynomial from Z to a polynomial in field p**l.
Then we do hensel lifting on all the factors of field p, and then here
we lift all the factors to field of p**l.
Though I see that, in the previous implementation, where we used polynomial
division in integer domain, this problem is solved. As we are using
UIntPoly representation, which is sparse there is no issue of space
(Though I am unsure of performance).
And here
we are implementing zassenhaus algorithm in field p, so a smaller
vector is required.
I personally think that we should move back to previous
implementation(SymPy does that, though its division is also in dense
polynomials).
I will look into the performance matter tonight and update you.
cc @certik

On Thu, Aug 18, 2016 at 11:16 PM, Isuru Fernando notifications@github.com
wrote:

In symengine/polys/uintpoly.cpp
#1066 (comment):

+}
+
+std::set<RCP, RCPBasicKeyLess> UIntPoly::zz_zassenhaus() const
+{

  • std::set<RCP, RCPBasicKeyLess> factors;
  • unsigned int n = get_degree();
  • UIntDict f(this->poly_);
  • if (n == 1) {
  •    factors.insert(rcp_from_this_cast<const UIntPoly>());
    
  •    return factors;
    
  • }
  • integer_class fc(this->poly_.dict_.begin()->second);
  • auto A = mp_get_ui(this->poly_.max_abs_coef());
  • auto b = get_lc();
  • auto B
  •    = std::abs(int(std::sqrt(n + 1)) \* std::pow(2, n) \* A \* mp_get_si(b));
    

Where is the dictionary of size pl made?


You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
https://github.com/symengine/symengine/pull/1066/files/e76c7a23a010d458384a14abb2973342023e1772#r75354742,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AI0uix0ExEpYj08AUdLIu1Y27wOEh5c6ks5qhJqTgaJpZM4Jl-E3
.

@abhinavagarwal07
abhinavagarwal07 Aug 19, 2016 Contributor

pow (2, n) can be replaced by (1u <<n) sincen is an unsigned int.
For n> 31 it will overflow

@isuruf
isuruf Aug 19, 2016 Member

Thanks for the links. Even though the field is p**l where p**l can be a large number, the vector is not going to be large right? For eg. (x**5 + 1) mod 10**100 needs only a vector of size 6.

@nishnik
nishnik Aug 19, 2016 Contributor

Sorry, I got confused.
You are right

@isuruf isuruf commented on the diff Aug 17, 2016
symengine/polys/uintpoly.cpp
+ integer_class q, r;
+ unsigned int a_deg, b_deg;
+ while (b_poly.degree() >= a.degree()) {
+ a_deg = a.degree();
+ b_deg = b_poly.degree();
+ q = b_poly.get_lc() / a.get_lc();
+ res[b_deg - a_deg] = q;
+ tmp = UIntDict(map_uint_mpz{{b_deg - a_deg, q}});
+ b_poly -= (a * tmp);
+ }
+ *quo = UIntDict(res);
+ *rem = b;
+ *rem -= (a * (*quo));
+}
+
+std::set<RCP<const UIntPoly>, RCPBasicKeyLess> UIntPoly::zz_zassenhaus() const
@isuruf
isuruf Aug 17, 2016 Member

This is a big method. Can you add a reference and also comment the code, so that I can understand what you are doing here?

@isuruf isuruf commented on an outdated diff Aug 17, 2016
symengine/polys/uintpoly.cpp
+ if (n == 1) {
+ factors.insert(rcp_from_this_cast<const UIntPoly>());
+ return factors;
+ }
+ integer_class fc(this->poly_.dict_.begin()->second);
+ auto A = mp_get_ui(this->poly_.max_abs_coef());
+ auto b = get_lc();
+ auto B
+ = std::abs(int(std::sqrt(n + 1)) * std::pow(2, n) * A * mp_get_si(b));
+ auto C = std::pow(n + 1, 2 * n) * std::pow(A, 2 * n - 1);
+ auto gamma = std::ceil(2 * std::log2(C));
+ auto bound = integer_class(int(2 * gamma * std::log(gamma)));
+ std::vector<std::pair<integer_class,
+ std::set<GaloisFieldDict, GaloisFieldDict::DictLess>>>
+ a;
+ for (integer_class i = 3_z; i <= bound; mp_nextprime(i, i)) {
@isuruf
isuruf Aug 17, 2016 Member

Use the prime sieve to generate primes.

@isuruf isuruf commented on an outdated diff Aug 17, 2016
symengine/polys/uintpoly.cpp
+ continue;
+ integer_class temp;
+ F.gf_monic(temp, outArg(F));
+ auto fsqfx = F.gf_zassenhaus();
+ a.push_back({i, fsqfx});
+ if (fsqfx.size() < 15 or a.size() > 4)
+ break;
+ }
+ auto fsqf = *(std::min_element(
+ a.begin(), a.end(),
+ [](const std::pair<integer_class,
+ std::set<GaloisFieldDict, GaloisFieldDict::DictLess>>
+ &lhs,
+ const std::pair<integer_class,
+ std::set<GaloisFieldDict, GaloisFieldDict::DictLess>>
+ &rhs) { return lhs.second.size() < rhs.second.size(); }));
@isuruf
isuruf Aug 17, 2016 Member

Sorry, didn't read this correctly before. Are you reusing a or just using it to get the min element? If so, you can get the min element, without even using the vector.

@isuruf isuruf commented on an outdated diff Aug 17, 2016
symengine/polys/uintpoly.cpp
+ } while (std::prev_permutation(bitmask.begin(), bitmask.end()));
+ return out;
+ };
+
+ auto _test_pl = [](const integer_class &fc, const integer_class &q,
+ const integer_class &pl) {
+ integer_class q_copy(q);
+ if (q_copy > pl / 2_z)
+ q_copy -= pl;
+ if (q_copy == 0_z)
+ return true;
+ return (fc % q_copy == 0_z);
+ };
+
+ while (2 * s <= T.size()) {
+ std::vector<std::vector<unsigned int>> subsets = _subsets(T, s);
@isuruf
isuruf Aug 17, 2016 Member

When s becomes large, this is going to be a huge vector.

@certik
Contributor
certik commented Aug 17, 2016
@nishnik
Contributor
nishnik commented Aug 18, 2016

@certik Yes sure! I will try to complete it by this week.

@isuruf isuruf commented on an outdated diff Aug 18, 2016
symengine/fields.cpp
+ while (1) {
+ GaloisFieldDict Q, R;
+ r0.gf_div(r1, outArg(Q), outArg(R));
+
+ if (R.dict_.empty())
+ break;
+ r0 = r1;
+ integer_class lc, inv;
+ R.gf_monic(lc, outArg(r1));
+ mp_invert(inv, lc, f.modulo_);
+
+ *s = s0 - s1 * Q;
+ *t = t0 - t1 * Q;
+
+ s0 = s1;
+ s1 = (*s) * inv;
@isuruf
isuruf Aug 18, 2016 Member

what you basically want to do here is s0, s1 = s1, (s0 - s1 *Q) * inv.

How about this?

s0 = (s0 - s1 *Q) * inv;
swap(s0, s1);
@isuruf isuruf commented on an outdated diff Aug 19, 2016
symengine/fields.cpp
+
+ if (R.dict_.empty())
+ break;
+ r0 = r1;
+ integer_class lc, inv;
+ lc = R.gf_monic(outArg(r1));
+ mp_invert(inv, lc, f.modulo_);
+
+ s0 = (s0 - s1 * Q) * inv;
+ std::swap(s0, s1);
+ t0 = (t0 - t1 * Q) * inv;
+ std::swap(t0, t1);
+ }
+ *s = s1;
+ *t = t1;
+ *h = r1;
@isuruf
isuruf Aug 19, 2016 Member

You can use s, t and h instead of new vectors s1, t1, r1

@isuruf isuruf commented on an outdated diff Aug 19, 2016
symengine/fields.cpp
+ 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);
+ GaloisFieldDict b(u);
@isuruf
isuruf Aug 19, 2016 Member

Why do you need an extra temporary b? You can use u instead of b

@isuruf isuruf commented on an outdated diff Aug 19, 2016
symengine/fields.cpp
+ u = t * b + 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()
+{
+ for (auto it = dict_.begin(); it != dict_.end(); ++it) {
+ if (*it > modulo_ / 2)
@isuruf
isuruf Aug 19, 2016 Member

You can take the computation, modulo_ / 2 outside the loop.

@isuruf isuruf commented on an outdated diff Aug 24, 2016
symengine/polys/uintpoly.cpp
+ } else {
+ G = UIntDict(b);
+ for (auto &i : S) {
+ G *= g[i];
+ }
+ G.itrunc(pl);
+ G.primitive(outArg(G));
+ auto q = G.dict_.begin()->second;
+ if (q != 0 and fc % q != 0)
+ continue;
+ }
+ UIntDict H = UIntDict(b);
+ std::vector<unsigned int> T_S(T.size());
+ auto it = std::set_difference(T.begin(), T.end(), S.begin(),
+ S.end(), T_S.begin());
+ T_S.resize(it - T_S.begin());
@isuruf
isuruf Aug 24, 2016 Member

What does this line do?

@isuruf isuruf and 1 other commented on an outdated diff Aug 24, 2016
symengine/polys/uintpoly.cpp
+ if (not _test_pl(fc, q, pl))
+ continue;
+ } else {
+ G = UIntDict(b);
+ for (auto &i : S) {
+ G *= g[i];
+ }
+ G.itrunc(pl);
+ G.primitive(outArg(G));
+ auto q = G.dict_.begin()->second;
+ if (q != 0 and fc % q != 0)
+ continue;
+ }
+ UIntDict H = UIntDict(b);
+ std::vector<unsigned int> T_S(T.size());
+ auto it = std::set_difference(T.begin(), T.end(), S.begin(),
@isuruf
isuruf Aug 24, 2016 Member

From what I gather, S is the set of numbers where bitmask is set to 1 and T_S is the set of numbers where bitmask is set to 0 right? If so, why not construct it like that, instead of this complicated way?

@nishnik
nishnik Aug 24, 2016 Contributor

Yes
That was needed in the previous implementation.
I will change that now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment