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: sage-6.4
Component: commutative algebra Keywords: random_fail
Cc: jpflori, vbraun Merged in:
Authors: Jean-Pierre Flori Reviewers: Volker Braun
Report Upstream: N/A Work issues:
Branch: 7a1c115 (Commits) Commit: 7a1c11575b725d703557c3c285ab0a3754a3f8bf
Dependencies: Stopgaps:

Description (last modified by vdelecroix)

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 sage-devel thread.

Attachments (1)

test_polynomials.py (3.0 KB) - added by vdelecroix 5 years ago.

Download all attachments as: .zip

Change History (24)

Changed 5 years ago by vdelecroix

comment:1 Changed 5 years ago by vdelecroix

  • Description modified (diff)

comment:2 Changed 5 years ago by vbraun

  • Keywords random_fail added

comment:3 Changed 5 years ago by jdemeyer

  • Description modified (diff)

comment:4 Changed 5 years ago by jdemeyer

Could you please describe the actual problem in the ticket description, preferably with a minimal example showing the bug.

Last edited 5 years ago by jdemeyer (previous) (diff)

comment:5 Changed 5 years ago by vdelecroix

  • Description modified (diff)

comment:6 Changed 5 years ago by vdelecroix

  • Description modified (diff)

comment:7 Changed 5 years ago by jpflori

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)*a4 + Mod(2, 3)*a3 + Mod(2, 3))*x + Mod(Mod(2, 3)*a3 +

Mod(2, 3)*a2, Mod(1, 3)*a4 + Mod(2, 3)*a3 + Mod(2, 3)) = T + 2*z3+z2 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 jpflori

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 jpflori

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 follow-up: Changed 5 years ago by jpflori

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 vdelecroix

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 pbruin

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 vdelecroix

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
Last edited 5 years ago by vdelecroix (previous) (diff)

comment:14 Changed 5 years ago by vdelecroix

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() with list(f._pari_()) in the code above then the RAM just fill in few seconds
  • if I replace f._pari_().list() with f.list() then I do not see error anymore

Vincent

comment:15 Changed 5 years ago by vdelecroix

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 jpflori

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 jpflori

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 jpflori

  • Authors set to Jean-Pierre Flori
  • 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:

7a1c115Perform additional contexts restoration for NTL.

comment:19 Changed 5 years ago by jpflori

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 jdemeyer

I don't think you need the <ntl_ZZ_pX> cast.

comment:21 Changed 5 years ago by jpflori

Yep, that's highly probable. I just got it from the previous snippet of code.

comment:22 Changed 5 years ago by vbraun

  • Reviewers set to Volker Braun
  • Status changed from needs_review to positive_review

comment:23 Changed 5 years ago by vbraun

  • Branch changed from u/jpflori/ticket/16885 to 7a1c11575b725d703557c3c285ab0a3754a3f8bf
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.