Opened 9 years ago

Closed 8 years ago

#13054 closed defect (fixed)

PARI polred() bug

Reported by: rbeezer Owned by: davidloeffler
Priority: major Milestone: sage-5.8
Component: number fields Keywords: sd40.5
Cc: roed Merged in: sage-5.8.beta2
Authors: Jeroen Demeyer Reviewers: David Roe
Report Upstream: Fixed upstream, but not in a stable release. Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description (last modified by jdemeyer)

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

Attachments (3)

trac_13054-workaround.patch (717 bytes) - added by nbruin 8 years ago.
Workaround (it's really a bug in pari)
pari-2.5.3.p3.diff (25.9 KB) - added by jdemeyer 8 years ago.
Diff for the PARI spkg, for review only
13054_polredbest.patch (14.3 KB) - added by jdemeyer 8 years ago.

Download all attachments as: .zip

Change History (23)

comment:1 Changed 9 years ago by 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__

comment:2 Changed 8 years ago by 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.

Last edited 8 years ago by nbruin (previous) (diff)

Changed 8 years ago by nbruin

Workaround (it's really a bug in pari)

comment:3 Changed 8 years ago by 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())

comment:4 Changed 8 years ago by 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 ...)

comment:5 Changed 8 years ago by jdemeyer

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

comment:6 Changed 8 years ago by jdemeyer

  • Description modified (diff)

comment:7 Changed 8 years ago by jdemeyer

  • Report Upstream changed from N/A to Reported upstream. No feedback yet.

comment:8 Changed 8 years ago by 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

comment:9 Changed 8 years ago by jdemeyer

  • Description modified (diff)

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.

comment:10 Changed 8 years ago by jdemeyer

  • Authors set to Jeroen Demeyer
  • Description modified (diff)

Changed 8 years ago by jdemeyer

Diff for the PARI spkg, for review only

comment:11 Changed 8 years ago by jdemeyer

  • Report Upstream changed from Reported upstream. No feedback yet. to Fixed upstream, but not in a stable release.

comment:12 Changed 8 years ago by jdemeyer

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

comment:13 Changed 8 years ago by jdemeyer

  • Description modified (diff)
  • Status changed from new to needs_review

comment:14 Changed 8 years ago by jdemeyer

  • Summary changed from Pernicious bug for algebraic numbers to PARI polred() bug

comment:15 Changed 8 years ago by jdemeyer

Any reviewers?...

comment:16 Changed 8 years ago by roed

  • Status changed from needs_review to needs_info

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

Changed 8 years ago by jdemeyer

comment:17 Changed 8 years ago by jdemeyer

  • Status changed from needs_info to needs_review

comment:18 Changed 8 years ago by roed

  • Status changed from needs_review to positive_review

comment:19 Changed 8 years ago by roed

  • Reviewers set to David Roe

comment:20 Changed 8 years ago by jdemeyer

  • Merged in set to sage-5.8.beta2
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.