PARI polred() bug

[rbeezer]

PARI's polred() returns reducible polynomials. GP session:

gp> pol = x^4 - 4294967296*x^2 + 54265257667816538374400;
gp> L = polred(pol);
gp> factor(L[4])
%14 =
[x^2 + 211955648366398871041 2]

Upstream:  http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=1395

This bug causes problem with arithmetic over QQbar:

Line 1563 (5.1.beta0) of sage/rings/qqbar.py:

rev = parent(best_elt.Mod(pari_poly).modreverse().lift())

Here, we are calling modreverse() on a element of a subfield. It would make a lot of exact linear algebra much more reliable if this was resolved.

Apply: 13054_polredbest.patch

spkg:  http://boxen.math.washington.edu/home/jdemeyer/spkg/pari-2.5.3.p3.spkg

[rbeezer]

This is not the same bug as above, but could be another symptom. Still hunting for an example to grab onto.

Run

n = 3
A = matrix(QQbar, n, n, [i/(i+1) + (i^2+3)*I for i in range(n^2)])
ev = A.eigenvalues()
(A-ev[0]).rref()

which should hangup and when killed with Ctrl-C yields

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_5.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("biA9IDMKQSA9IG1hdHJpeChRUWJhciwgbiwgbiwgW2kvKGkrMSkgKyAoaV4yKzMpKkkgZm9yIGkgaW4gcmFuZ2Uobl4yKV0pCmV2ID0gQS5laWdlbnZhbHVlcygpCihBLWV2WzBdKS5ycmVmKCk="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>

  File "/tmp/tmpiYGsKk/___code___.py", line 6, in <module>
    exec compile(u'(A-ev[_sage_const_0 ]).rref()
  File "", line 1, in <module>

  File "matrix2.pyx", line 5526, in sage.matrix.matrix2.Matrix.rref (sage/matrix/matrix2.c:28497)
  File "matrix2.pyx", line 5839, in sage.matrix.matrix2.Matrix.echelon_form (sage/matrix/matrix2.c:29882)
  File "matrix2.pyx", line 5745, in sage.matrix.matrix2.Matrix.echelonize (sage/matrix/matrix2.c:29429)
  File "matrix2.pyx", line 5897, in sage.matrix.matrix2.Matrix._echelon_in_place_classical (sage/matrix/matrix2.c:30399)
  File "/sage/sage-5.1.beta0/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3831, in __nonzero__
    return self.__nonzero__()
  File "/sage/sage-5.1.beta0/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3834, in __nonzero__
    self.exactify()
  File "/sage/sage-5.1.beta0/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3466, in exactify
    self._set_descr(self._descr.exactify())
  File "/sage/sage-5.1.beta0/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 7591, in exactify
    left.exactify()
  File "/sage/sage-5.1.beta0/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3466, in exactify
    self._set_descr(self._descr.exactify())
  File "/sage/sage-5.1.beta0/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 7591, in exactify
    left.exactify()
  File "/sage/sage-5.1.beta0/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3466, in exactify
    self._set_descr(self._descr.exactify())
  File "/sage/sage-5.1.beta0/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 7593, in exactify
    gen = left._exact_field().union(right._exact_field())
  File "/sage/sage-5.1.beta0/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 2251, in union
    pari_nf = self.pari_field()
  File "/sage/sage-5.1.beta0/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 2145, in pari_field
    self._pari_field = pari_pol.nfinit(1)
  File "gen.pyx", line 7470, in sage.libs.pari.gen.gen.nfinit (sage/libs/pari/gen.c:31383)
  File "gen.pyx", line 7484, in sage.libs.pari.gen.gen._nfinit_with_prec (sage/libs/pari/gen.c:31587)
KeyboardInterrupt
__SAGE__

[nbruin]

Found in  this sage-devel thread

sage: R.<x,y>=QQbar[]
sage: f=x^2+x^3-y^2+y^4*x^4
sage: RX=PolynomialRing(QQbar,'x')
sage: RY=PolynomialRing(QQbar,'y')
sage: x0=3+2^-8
sage: y0=RY(f((x0,y))).roots(multiplicities=False)[0]
sage: f(x+x0,y0)

sounds like it might be triggerering the same problem as reported in the ticket.

The slightest bit of %debug fun shows that the problem arises here > /usr/local/sage/5.6rc0/local/lib/python2.7/site-packages/sage/rings/qqbar.py(1563)do_polred():

   1562     parent = poly.parent()
-> 1563     rev = parent(best_elt.Mod(pari_poly).modreverse().lift())
   1564     return parent(best_elt), rev, parent(best)

with these inputs:

R.<y>=QQ['y']
best_elt = (1/16*y^2 - 134217728)._pari_()
pari_poly = y^4 - 4294967296*y^2 + 54265257667816538374400
rev = parent(best_elt.Mod(pari_poly).modreverse().lift())

confirming a degree mismatch problem.

EDIT: actually, from this example you don't see any degree mismatch problem immediately. The problem is that the minimal polynomial of best_elt.Mod(pari_poly) is quadratic rather than quartic. So modreverse is right in complaining. However, these polynomials were obtained from a routine that misreported a minimal poly somewhere.

[nbruin]

Workaround (it's really a bug in pari)

[jdemeyer]

Why is this a bug in PARI? It seems to me that PARI behaves as documented. The degree of the minimal polynomial is less than the degree of the field.

What would you expect PARI to do in this case?

sage: R.<y> = QQ['y']
sage: best_elt = (1/16*y^2 - 134217728)._pari_()
sage: pari_poly = y^4 - 4294967296*y^2 + 54265257667816538374400
sage: rev = parent(best_elt.Mod(pari_poly).modreverse().lift())

[nbruin]

GP instructions illustrating the problem:

poly=x^4 - 4294967296*x^2 + 54265257667816538374400
L=polred(poly,3)
elt=Mod(L[4,1],poly)
minpoly(elt)
factor(L[4,2])

The documentation of polred promises that L[4,2] is a minimal polynomial, so it should be square-free (and it's not ...)

[jdemeyer]

OK sorry, I misunderstood that there was a problem with modreverse(). But the bug description was very unclear...

[rbeezer]

Nils and Jeroen,

THANKS for making some progress on this one. David Roe and I had some "debug fun" with this last summer and I wasn't quite able to keep up and top of what was happening, so I'm glad we've got some better ideas about the source of the problem.

Hopefully this will make some matrix arithmetic over QQbar work better.

Rob

[jdemeyer]

The workaround looks good to me. You should add a comment though saying why we do this, refering to this ticket. And of course add the obligatory doctest.

[jdemeyer]

Diff for the PARI spkg, for review only

[jdemeyer]

The spkg is ready for review, but this also needs a Sage library patch.

[jdemeyer]

Any reviewers?...

[roed]

Doesn't the input to do_polred need to be irreducible (current documentation just says monic polynomial with integer coefficients)? Otherwise it looks good.