Opened 5 years ago
Closed 5 years ago
#16885 closed defect (fixed)
random failure with factorization of polynomials over finite fields
Reported by:  vdelecroix  Owned by:  

Priority:  critical  Milestone:  sage6.4 
Component:  commutative algebra  Keywords:  random_fail 
Cc:  jpflori, vbraun  Merged in:  
Authors:  JeanPierre Flori  Reviewers:  Volker Braun 
Report Upstream:  N/A  Work issues:  
Branch:  7a1c115 (Commits)  Commit:  7a1c11575b725d703557c3c285ab0a3754a3f8bf 
Dependencies:  Stopgaps: 
Description (last modified by )
Ticket #16681 pops up a strange error when doing factorization of polynomials over finite fields. The problem occurs when we repeatedly create finite fields of different degree and ask for factorization of some random polynomial in that field as in the example below
sage: R = PolynomialRing(GF(3),'T') sage: for _ in xrange(5000): ....: p = R.random_element(degree=(2,8)) ....: d = p.degree() ....: if d <= 2: continue ....: q = p.change_ring(GF(3**d, 'z')) ....: assert q.factor().prod() == q, "problem with p={}".format(p) Traceback (most recent call last): ... AssertionError: problem with p=T^4 + T^2
You can also have a look at the file test_polynomials.py which contains more precise output for the same tests.
See also this sagedevel thread.
Attachments (1)
Change History (24)
Changed 5 years ago by
comment:1 Changed 5 years ago by
 Description modified (diff)
comment:2 Changed 5 years ago by
 Keywords random_fail added
comment:3 Changed 5 years ago by
 Description modified (diff)
comment:4 Changed 5 years ago by
comment:5 Changed 5 years ago by
 Description modified (diff)
comment:6 Changed 5 years ago by
 Description modified (diff)
comment:7 Changed 5 years ago by
Here is a little more detailed failure output:
<type 'list'> [0, Mod(Mod(1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3)), Mod(Mod(2, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3)), 0, Mod(Mod(2, 3), Mod(1, 3)* a^4 + Mod(2, 3)*a^3 + Mod(2, 3))] <type 'sage.libs.pari.gen.gen'> Mod(Mod(2, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x^4 + Mod(Mod(2, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x^2 + Mod( Mod(1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x <type 'sage.libs.pari.gen.gen'> [Mod(Mod(1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x + Mod(Mod(0, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3)), 1; Mod(Mod (1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x + Mod(Mod(1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3)), 1; Mod(Mod(1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x + Mod(Mod(1, 3)*a^3 + Mod(1, 3)*a^2 + Mod(2, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3)), 1; Mod(Mod(1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3 ))*x + Mod(Mod(2, 3)*a^3 + Mod(2, 3)*a^2, Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3)), 1] <type 'list'> [[Mod(Mod(1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x + Mod(Mod(0, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3)), Mod(Mod(1, 3), Mod(1, 3)*a^ 4 + Mod(2, 3)*a^3 + Mod(2, 3))*x + Mod(Mod(1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3)), Mod(Mod(1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x + Mod(Mod( 1, 3)*a^3 + Mod(1, 3)*a^2 + Mod(2, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3)), Mod(Mod(1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x + Mod(Mod(2, 3)*a^3 + Mod(2, 3)*a^2, Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))]~, [1, 1, 1, 1]~] <type 'sage.libs.pari.gen.gen'> Mod(Mod(1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x^4 + Mod(Mod(0, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x^3 + Mod( Mod(1, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x^2 + Mod(Mod(2, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a^3 + Mod(2, 3))*x + Mod(Mod(0, 3), Mod(1, 3)*a^4 + Mod(2, 3)*a ^3 + Mod(2, 3)) ERROR (after 761 iterations): polynomial : 2*T^4 + 2*T^2 + T polynomial parent : <class 'sage.rings.finite_rings.finite_field_givaro.FiniteField_givaro_with_category'>, x^4 + 2*x^3 + 2 factorization: <class 'sage.structure.factorization.Factorization'>, (2) * T * (T + 1) * (T + 2*z^2 + 2*z + 1) * (T + z^3 + z^2 + 2) product : <type 'sage.rings.polynomial.polynomial_zz_pex.Polynomial_ZZ_pEX'>, 2*T^4 + (2*z^3 + z + 2)*T^3 + (2*z^3 + z^2 + z + 1)*T^2 + (z^2 + 1)*T product parent : <class 'sage.rings.finite_rings.finite_field_givaro.FiniteField_givaro_with_category'>, x^4 + 2*x^3 + 2 product : 2*T^4 + (2*z^3 + z + 2)*T^3 + (2*z^3 + z^2 + z + 1)*T^2 + (z^2 + 1)*T
basically printing what happens when PARI is called before the slightly modified Vincent output.
As you might see, the factorization at the end of the call to factor() in polynomial element is correct. More precisely, the first two factors in the PARI factorization are the same as what we get back in Vincent's file in the Factorization object, but the third one Mod(Mod(1, 3), Mod(1, 3)*a^{4 + Mod(2, 3)*a}3 + Mod(2, 3))*x + Mod(Mod(2, 3)*a^{3 + }
Mod(2, 3)*a^{2, Mod(1, 3)*a}4 + Mod(2, 3)*a^{3 + Mod(2, 3)) = T + 2*z}3+z^{2 does not appear there and the fourth one becomes the third one. }
The fourth one (T + 2*z2 + 2*z + 1) in Sage's Factorizaiton object is mysterious (and different from the one which disappeared from PARI fac object.
comment:8 Changed 5 years ago by
Something wrong happens in:
F = [(R(f), int(e)) for f, e in zip(pols, exps)]
With further debugging, it seems doing that without R and int is fine, and putting R and int breaks some values in the list.
comment:9 Changed 5 years ago by
Here is one quite smallish example:
[(Mod(Mod(1, 3), Mod(1, 3)*a^5 + Mod(2, 3)*a + Mod(1, 3))*x + Mod(Mod(1, 3)*a^3 + Mod(2, 3)*a^2 + Mod(1, 3), Mod(1, 3)*a^5 + Mod(2, 3)*a + Mod(1, 3)), 1), (Mod(Mod(1, 3), Mod(1, 3)*a^5 + Mod(2, 3)*a + Mod(1, 3))*x + Mod(Mod(1, 3)*a^4 + Mod(2, 3), Mod(1, 3)*a^5 + Mod(2, 3)*a + Mod(1, 3)), 1), (Mod(Mod(1, 3), Mod(1, 3)*a^5 + Mod(2, 3)*a + Mod(1, 3))*x + Mod(Mod(1, 3)*a^4 + Mod(1, 3)*a^3 + Mod(1, 3)*a^2 + Mod(2, 3), Mod(1, 3)*a^5 + Mod(2, 3)*a + Mod(1, 3)), 1), (Mod(Mod(1, 3), Mod(1, 3)*a^5 + Mod(2, 3)*a + Mod(1, 3))*x + Mod(Mod(2, 3)*a^4 + Mod(2, 3)*a^2 + Mod(2, 3)*a, Mod(1, 3)*a^5 + Mod(2, 3)*a + Mod(1, 3)), 1), (Mod(Mod(1, 3), Mod(1, 3)*a^5 + Mod(2, 3)*a + Mod(1, 3))*x + Mod(Mod(2, 3)*a^4 + Mod(1, 3)*a^3 + Mod(1, 3)*a^2 + Mod(1, 3)*a, Mod(1, 3)*a^5 + Mod(2, 3)*a + Mod(1, 3)), 1)] <type 'list'> [(T + z^3 + 2*z^2 + 1, 1), (T + z^3, 1), (T + z^4 + z^3 + z^2 + 2, 1), (T + 2*z^4 + 2*z^2 + 2*z, 1), (T + 2*z^4 + z^3 + z^2 + z, 1)]
and the ring/pol were:
polynomial : T^5 + 2*T^4 + 1 polynomial parent : <class 'sage.rings.finite_rings.finite_field_givaro.FiniteField_givaro_with_category'>, x^5 + 2*x + 1 factorization: <class 'sage.structure.factorization.Factorization'>, (T + z^3) * (T + z^3 + 2*z^2 + 1) * (T + z^4 + z^3 + z^2 + 2) * (T + 2*z^4 + 2*z^2 + 2*z) * (T + 2*z^4 + z^3 + z^2 + z)
}}}
comment:10 followup: ↓ 11 Changed 5 years ago by
Simpler test without factor stuff:
while True: d = ZZ.random_element(2,10) K = GF(3**d, 'a') R = PolynomialRing(K, 'x') for i in xrange(10**4): f = R.random_element() if f != pari([c._pari_("a") for c in R(pari(f))]).Polrev(): print f, pari(f), R(pari(f)), i raise ValueError
Does someone confirm failures with that one?
comment:11 in reply to: ↑ 10 Changed 5 years ago by
Replying to jpflori:
Does someone confirm failures with that one?
Yes, got a failure as well. You might update the description of the ticket.
Vincent
comment:12 Changed 5 years ago by
I am also getting failures with the following variant:
while True: d = ZZ.random_element(2, 8) K = GF(3**d, 'a') R = PolynomialRing(K, 'x') f = R.random_element() pf = pari(f) pfl = pf.list() Kpfl = [K(c) for c in pfl] Rf = R(Kpfl) if f != Rf: print f, pf, pfl, Kpfl, Rf raise ValueError
The values in one particularly simple case are as follows (applied lift(lift(...))
to PARI objects to make them readable):
f = a^3*x^2 + (a^3 + 2*a^2 + a)*x pf = a^3*x^2 + (a^3 + 2*a^2 + a)*x pfl = [0, a^3 + 2*a^2 + a, a^3] Kpfl = [0, a^3 + 2*a^2 + a, a^3] Rf = a^3*x^2 + (2*a^2 + 2*a + 2)*x
Everything is consistent up to the conversion of Kpfl
(the list of coefficients) to a polynomial, where the coefficient a^3 + 2*a^2 + a
changes into 2*a^2 + 2*a + 2
.
comment:13 Changed 5 years ago by
Hi Peter,
Just to mention, using your code, we also get the fact that it works the second time
sage: while True: ....: d = ZZ.random_element(2, 8) ....: K = GF(3**d, 'a') ....: R = PolynomialRing(K, 'x') ....: f = R.random_element() ....: pf = pari(f) ....: pfl = pf.list() ....: Kpfl = [K(c) for c in pfl] ....: Rf = R(Kpfl) ....: if f != Rf: ....: print f, pf, pfl, Kpfl, Rf ....: raise ValueError <BOOM> sage: R(Kpfl) == f True
comment:14 Changed 5 years ago by
The error occurs more precisely in the initialization of Polynomial_ZZ_pEX
. It can be tracked with the following modifications
 a/src/sage/rings/polynomial/polynomial_zz_pex.pyx +++ b/src/sage/rings/polynomial/polynomial_zz_pex.pyx @@ 132,6 +132,7 @@ cdef class Polynomial_ZZ_pEX(Polynomial_template): x = x.list() if PY_TYPE_CHECK(x, list) or PY_TYPE_CHECK(x, tuple): + y = x[:] Polynomial.__init__(self, parent, is_gen=is_gen) (<Polynomial_template>self)._cparent = get_cparent(parent) celement_construct(&self.x, (<Polynomial_template>self)._cparent) @@ 143,6 +144,19 @@ cdef class Polynomial_ZZ_pEX(Polynomial_template): e = K(e) d = parent._modulus.ZZ_pE(list(e.polynomial())) ZZ_pEX_SetCoeff(self.x, i, d.x) + while y and y[1] == 0: + y.pop(1) + if list(self) != y: + print "ERROR in Polynomial_ZZ_pEX constructor" + print "call with args:" + print " x = {}".format(x) + print " check = {}".format(check) + print " is_gen = {}".format(is_gen) + print " construct = {}".format(construct) + print "but got the result {}".format(self) + print "variables" + print " K = {}".format(K) + raise RuntimeError return
and running
while True: d = ZZ.random_element(2, 8) K = GF(3**d, 'a') R = PolynomialRing(K, 'x') f = R.random_element() l = [K(c) for c in f._pari_().list()] g = R(l)
I got the following kind of error
ERROR in Polynomial_ZZ_pEX constructor call with args: x = [a^5 + a^4 + 2*a^3 + 2*a^2 + 2*a + 1, a^6 + 2*a^5 + a^4 + a^3 + 2*a^2 + 2*a + 1, a^6 + 2*a^5 + 2*a^4 + 2*a^3 + a^2 + 2] check = True is_gen = False construct = False but got the result (a^6 + 2*a^5 + 2*a^4 + 2*a^3 + a^2 + 2)*x^2 + (2*a^5 + 2*a^4 + a^3 + a^2 + 2)*x + a^5 + a^4 + 2*a^3 + 2*a^2 + 2*a + 1 variables K = Finite Field in a of size 3^7
Might be related or not:
 if I replace
f._pari_().list()
withlist(f._pari_())
in the code above then the RAM just fill in few seconds  if I replace
f._pari_().list()
withf.list()
then I do not see error anymore
Vincent
comment:15 Changed 5 years ago by
I go a bit further and seems that the trouble comes from the initialization in the Sage NTL wrapper. More precisely, in the constructor of Polynomial_ZZ_pEX
there is a problem in the following loop
for i,e in enumerate(x): # self(x) is supposed to be a conversion, # not necessarily a coercion. So, we must # not do K.coerce(e) but K(e). e = K(e) d = parent._modulus.ZZ_pE(list(e.polynomial())) ZZ_pEX_SetCoeff(self.x, i, d.x)
When an error occurs the coefficient list(e.polynomial())
is right but d
is wrong. The mysterious parent._modulus
is of class ntl_ZZ_pEContext_class
in libs/ntl/ntl_ZZ_pEContext.pyx
. The file is poorly documented and I do not know NTL enough to dig further.
Vincent
comment:16 Changed 5 years ago by
I suspected something with NTL and contexts wrongly restored at some point but did not investigate far enough and thought "well it doesn't seem NTL is used that much in the problematic code, just PARI and factorization". But with your info it seems it might indeed be some problem with NTL contexts restoration (at the C++ level). I'll try to track it down.
comment:17 Changed 5 years ago by
Might be something like #9524. There we found out that the NTL context stuff for the polynomial defining the finite field was restored, but the context for the integer giving the characteristic was not.
comment:18 Changed 5 years ago by
 Branch set to u/jpflori/ticket/16885
 Commit set to 7a1c11575b725d703557c3c285ab0a3754a3f8bf
 Status changed from new to needs_review
It seems adding a little contexts restoration helps solves this.
New commits:
7a1c115  Perform additional contexts restoration for NTL.

comment:19 Changed 5 years ago by
Maybe we should add a test, but it implies finding an easily and quickly reproducible test...
Fixing the typos I did would be nice as well :)
comment:20 Changed 5 years ago by
I don't think you need the <ntl_ZZ_pX>
cast.
comment:21 Changed 5 years ago by
Yep, that's highly probable. I just got it from the previous snippet of code.
comment:22 Changed 5 years ago by
 Reviewers set to Volker Braun
 Status changed from needs_review to positive_review
comment:23 Changed 5 years ago by
 Branch changed from u/jpflori/ticket/16885 to 7a1c11575b725d703557c3c285ab0a3754a3f8bf
 Resolution set to fixed
 Status changed from positive_review to closed
Could you please describe the actual problem in the ticket description, preferably with a minimal example showing the bug.