pari error at evaluating polynomial

[mmarco]

I get a strange error when i try to evaluate a polynomial in QQbar:

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)
---------------------------------------------------------------------------
PariError Traceback (most recent call
last)

/home/mmarco/<ipython console> in <module>()

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/
polynomial/multi_polynomial_element.pyc in __call__(self, *x,
**kwds)
154 y =
K(0)
155 for (m,c) in
self.element().dict().iteritems():
--> 156 y += c*misc.mul([ x[i]**m[i] for i in range(n) if
m[i] != 0])
157 return
y

158

/usr/local/sage/local/lib/python2.6/site-packages/sage/structure/
element.so in sage.structure.element.ModuleElement.__iadd__ (sage/
structure/element.c:8264)
()

/usr/local/sage/local/lib/python2.6/site-packages/sage/structure/
element.so in sage.structure.element.ModuleElement._add_ (sage/
structure/element.c:8117)
()

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/
polynomial/multi_polynomial_element.pyc in _add_(self,
right)
212 def _add_(self,
right):
213 #return self.parent()(self.__element +
right.__element)
--> 214 return self.__class__(self.parent(),self.__element +
right.__element)

215
216 def _sub_(self,
right):

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/
polynomial/polydict.so in
sage.rings.polynomial.polydict.PolyDict.__add__ (sage/rings/polynomial/
polydict.c:7503)()

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/
polynomial/polydict.so in
sage.rings.polynomial.polydict.PolyDict.__init__ (sage/rings/
polynomial/polydict.c:1781)()

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in __eq__(self, other)
3049 if self is other: return
True
3050 if other._descr.is_rational() and
other._descr.rational_value() == 0:
-> 3051 return not
self.__nonzero__()
3052 if self._descr.is_rational() and
self._descr.rational_value() == 0:
3053 return not
other.__nonzero__()

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in __nonzero__(self)
3081 if self._value.prec() <
128:
3082
self._more_precision()
-> 3083 return
self.__nonzero__()

3084
3085 #
Sigh...

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in __nonzero__(self)

3084
3085 #
Sigh...
-> 3086
self.exactify()
3087 return
self.__nonzero__()

3088

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in exactify(self)
2739 od =
self._descr
2740 if od.is_exact():
return
-> 2741
self._set_descr(self._descr.exactify())

2742
2743 def _set_descr(self,
new_descr):

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in exactify(self)
5891 left =
self._left
5892 right =
self._right
-> 5893
left.exactify()
5894
right.exactify()
5895 gen =
left._exact_field().union(right._exact_field())

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in exactify(self)
2739 od =
self._descr
2740 if od.is_exact():
return
-> 2741
self._set_descr(self._descr.exactify())

2742
2743 def _set_descr(self,
new_descr):

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in exactify(self)
5891 left =
self._left
5892 right =
self._right
-> 5893
left.exactify()
5894
right.exactify()
5895 gen =
left._exact_field().union(right._exact_field())

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in exactify(self)
2739 od =
self._descr
2740 if od.is_exact():
return
-> 2741
self._set_descr(self._descr.exactify())

2742
2743 def _set_descr(self,
new_descr):

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in exactify(self)
5891 left =
self._left
5892 right =
self._right
-> 5893
left.exactify()
5894
right.exactify()
5895 gen =
left._exact_field().union(right._exact_field())

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in exactify(self)
2739 od =
self._descr
2740 if od.is_exact():
return
-> 2741
self._set_descr(self._descr.exactify())

2742
2743 def _set_descr(self,
new_descr):

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in exactify(self)
5892 right =
self._right
5893
left.exactify()
-> 5894
right.exactify()
5895 gen =
left._exact_field().union(right._exact_field())
5896 left_value =
gen(left._exact_value())

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in exactify(self)
2739 od =
self._descr
2740 if od.is_exact():
return
-> 2741
self._set_descr(self._descr.exactify())

2742
2743 def _set_descr(self,
new_descr):

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in exactify(self)
5891 left =
self._left
5892 right =
self._right
-> 5893
left.exactify()
5894
right.exactify()
5895 gen =
left._exact_field().union(right._exact_field())

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in exactify(self)
2739 od = self._descr
2740 if od.is_exact(): return
-> 2741 self._set_descr(self._descr.exactify())
2742
2743 def _set_descr(self, new_descr):

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in exactify(self)
5109 den, my_factor = clear_denominators(my_factor)
5110
-> 5111 red_elt, red_back, red_pol = do_polred(my_factor)
5112
5113 field = NumberField(red_pol, 'a', check=False)

/usr/local/sage/local/lib/python2.6/site-packages/sage/rings/qqbar.pyc
in do_polred(poly)
1206 assert(best is not None)
1207 parent = poly.parent()
-> 1208 rev = parent(best_elt.Mod(pari_poly).modreverse().lift())
1209 return parent(best_elt), rev, parent(best)
1210

/usr/local/sage/local/lib/python2.6/site-packages/sage/libs/pari/
gen.so in sage.libs.pari.gen._pari_trap (sage/libs/pari/gen.c:49380)()

PariError: (5)


Strange thing is that if i do only one of the two substitutions, it
works fine.

sage: f(x,y0)
(-0.4436467746868897? + 0.00818072177704596?*I)*x^4 + x^3 + x^2 -
0.006140793684108297? - 0.666096452500658?*I
sage: f(x+x0,y)
x^4*y^4 + 769/64*x^3*y^4 + 1774083/32768*x^2*y^4 +
454756609/4194304*x*y^4 + 349707832321/4294967296*y^4 + x^3 +
2563/256*x^2 - y^2 + 2167811/65536*x + 606145025/16777216



Someone has a clue about what the problem migh be?

[Nils Bruin]

Looks like a known problem:

http://trac.sagemath.org/sage_trac/ticket/13054

The routine in question tries to find a better defining polynomial for
a given field. With

R.<y>=QQ['y']
poly=y^4 - 4294967296*y^2 + 54265257667816538374400

it executes:

degree = poly.degree()
pari_poly = pari(poly)

red_table = pari_poly.polred(3)

best = None
best_discr = None

for i in range(red_table.nrows()):
red_poly = red_table[i,1]
if red_poly.poldegree() < degree:
continue
red_discr = red_poly.poldisc().abs()
if best_discr is None or red_discr < best_discr:
best = red_poly
best_discr = red_discr
best_elt = red_table[i,0]
assert(best is not None)

After which (this happened for i==3) we have best_discr == 0, so the
found element generates only a subfield and polred, quite
misleadingly, only reports the characteristic polynomial rather than
the minimal polynomial. We could avoid this problem by
testing "red_descr > 0 and ( best_discr is None or red_discr <
best_discr)". The pari documentation of polred suggests there is
always at least one element that passes that test, so the assertion
should still hold.

According to the documentation, this is an error in pari, since it
should be returning the minimal poly.

GP session illustrating the error: (GP/PARI 2.5.3 (development
git-6fd07f9))

? poly=x^4 - 4294967296*x^2 + 54265257667816538374400
%1 = x^4 - 4294967296*x^2 + 54265257667816538374400
? L=polred(poly,3)
%2 =
[1 x - 1]

[1/145522114880*x^3 - 67108864/2273783045*x x^4 - 11005853696*x^2 +
356327727107810941435025]

[1/4*x x^4 - 268435456*x^2 + 211973662764908353025]

[1/16*x^2 - 134217728 x^4 + 423911296732797742082*x^2 +
44925196874420524410175311836119348423681]

? elt=Mod(L[4,1],poly)
%3 = Mod(1/16*x^2 - 134217728, x^4 - 4294967296*x^2 +
54265257667816538374400)
? minpoly(elt)
%4 = x^2 + 211955648366398871041
? factor(L[4,2])
%5 =
[x^2 + 211955648366398871041 2]

As you can see L[4,2] is the square of the minimal polynomial.

[mmarco]

Thanks, i will try the workaround. Should this be reported upstream to
pari?

[Jeroen Demeyer]

On 2013-01-30 21:08, mmarco wrote:
> Should this be reported upstream to pari?

Just did so.