Are Laurent polynomials really this broken? (or is it me?)

[Andrew Mathas]

I have a Laurent polynomial ring in two variables: R = Z[u,u^-1,v,v^-1].

I want to implement the ring automorphism of R which sends u to u^-1 and v to v^-1.

If p(u,v) is in R then p.subsitutue(u,u^-1,v=v^-1) should work but sadly not:

sage: p=2*u**-1*v**-1+u*v
sage: p.substitute(u=u^-1,v=v^-1)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-39-eb90c4c3368e> in <module>()
----> 1 p.substitute(u=u**-Integer(1),v=v**-Integer(1))

/usr/local/src/sage/sage-5.11/local/lib/python2.7/site-packages/sage/structure/element.so in sage.structure.element.Element.substitute (sage/structure/element.c:6349)()

/usr/local/src/sage/sage-5.11/local/lib/python2.7/site-packages/sage/rings/polynomial/laurent_polynomial.so in sage.rings.polynomial.laurent_polynomial.LaurentPolynomial_mpair.subs (sage/rings/polynomial/laurent_polynomial.c:10753)()

/usr/local/src/sage/sage-5.11/local/lib/python2.7/site-packages/sage/structure/element.so in sage.structure.element.RingElement.__imul__ (sage/structure/element.c:14780)()

/usr/local/src/sage/sage-5.11/local/lib/python2.7/site-packages/sage/structure/coerce.so in sage.structure.coerce.CoercionModel_cache_maps.bin_op (sage/structure/coerce.c:8169)()

TypeError: unsupported operand parent(s) for '*': 'Multivariate Laurent Polynomial Ring in u, v over Rational Field' and 'Multivariate Laurent Polynomial Ring in u, v over Integer Ring'

I played with this for a little while and finally worked out where the prblem is. Consider:

sage: [ [c,u**-exp[0],v**-exp[1]] for (exp,c) in p.dict().iteritems()]
[[2, u, v], [1, u^-1, v^-1]]
sage: [ [c,u**-exp[0]*v**-exp[1]] for (exp,c) in p.dict().iteritems()]
[[2, u*v], [1, u^-1*v^-1]]
sage: [ [c*u**-exp[0]*v**-exp[1]] for (exp,c) in p.dict().iteritems()]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-36-e01dd53c18a2> in <module>()
----> 1 [ [c*u**-exp[Integer(0)]*v**-exp[Integer(1)]] for (exp,c) in p.dict().iteritems()]

/usr/local/src/sage/sage-5.11/local/lib/python2.7/site-packages/sage/structure/element.so in sage.structure.element.RingElement.__mul__ (sage/structure/element.c:14531)()

/usr/local/src/sage/sage-5.11/local/lib/python2.7/site-packages/sage/structure/coerce.so in sage.structure.coerce.CoercionModel_cache_maps.bin_op (sage/structure/coerce.c:8169)()

TypeError: unsupported operand parent(s) for '*': 'Multivariate Laurent Polynomial Ring in u, v over Rational Field' and 'Multivariate Laurent Polynomial Ring in u, v over Integer Ring'

So it seems that sage is happy to prduce the monomials but that when it comes to multiplying by the integer coefficeints it coerces itself into the field of fractions and blows up. On the plus side, at least I have a solution:

sage: sum( R(c)*u**-exp[0]*v**-exp[1] for (exp,c) in p.dict().iteritems() )
2*u*v + u^-1*v^-1

Am I missing something or is this the easiest hack to get this to work?

Andrew

[Nils Bruin]

I think the error is elsewhere:

sage: P.<u,v>=LaurentPolynomialRing(ZZ,2)
sage: p=2*u**-1*v**-1+u*v
sage: [ type(c) for (exp,c) in p.dict().iteritems()]
[sage.rings.rational.Rational, sage.rings.rational.Rational]

The coefficients are returned as rationals, not integers. Indeed:

sage: type(P(1).coefficients()[0])
<type 'sage.rings.integer.Integer'>
sage: type(P(1*u).coefficients()[0])
<type 'sage.rings.integer.Integer'>
sage: type(P(1*u^(-1)).coefficients()[0])
<type 'sage.rings.rational.Rational'>
sage: type(P(u^(-1)).coefficients()[0])
<type 'sage.rings.rational.Rational'>
sage: P(u^(-1)).coefficients()[0]*u+u
TypeError: unsupported operand parent(s) for '+': 'Multivariate Laurent Polynomial Ring in u, v over Rational Field' and 'Multivariate Laurent Polynomial Ring in u, v over Integer Ring'

As you can see, even the coefficient of the monomial u^(-1) is already a rational number. I can see how that would confuse the coercion system. I suspect that's because we're just hitting generic inversion code when we do u^(-1), i.e., the same code that can compute (3*u)^(-1). One can construct the monomial properly:

sage: m=P({(-1,0):1})
sage: type(m.coefficients()[0]
<type 'sage.rings.integer.Integer'>

Using generic division for these operations is a little problematic:

sage: parent((u)^(-1))

Multivariate Laurent Polynomial Ring in u, v over Integer Ring

sage: parent((1+u)^(-1))
Fraction Field of Multivariate Polynomial Ring in u, v over Integer Ring
sage: parent(1^(-1))
Rational Field

As you can see, 1^(-1) is returned as a rational number even though it would fit in ZZ again. Laurent polynomials behave differently, but as it turns out, their coefficients don't.

So yes, I'd say laurent polynomials are broken, especially if defined over non-fields. But it's not easy to fix, because Laurent polynomials want u^(-1) to behave differently from (1+u)^(-1) and sage isn't set up for that.

[Andrew Mathas]

On Thursday, 19 September 2013 08:51:02 UTC+2, Nils Bruin wrote:

I think the error is elsewhere:

sage: P.<u,v>=LaurentPolynomialRing(ZZ,2)
sage: p=2*u**-1*v**-1+u*v
sage: [ type(c) for (exp,c) in p.dict().iteritems()]
[sage.rings.rational.Rational, sage.rings.rational.Rational]

The coefficients are returned as rationals, not integers.

Thanks for the clarification Nils. So Laurent polynomials are not playing well with their coefficients, but I wonder if the real issue is that even polynomials over rings have not been implemented carefully eough. Consider:

sage: A.<x>=PolynomialRing(ZZ)
sage: x in A, x^-1 in A
(True, False)
sage: x^-1
1/x

Everything is correct except, I think, the last line: x^-1 it NOT defined so this should be returning an error rather than coercing into the field of fractions. If I want x^-1 to be defined then I should define it properly. I think do not think that sage should try and second guess what I mean and put the result into a different ring. There are scenarios I can think of where by not triggering an error code further down the track will produce wrong answers. As a result, to compensate for this it will sometimes be necessary litter our code with checks to make sure that our coefficients belong to the ring that we think they do.

This is essentially the same problem that LaurentPolynomials is having:

     
sage: A.<x>=LaurentPolynomialRing(ZZ)
sage: x^-1 in A
True
sage: (2*x)^-1
1/2*x^-1

As you say, for x^-1 and (2*x)^-1 sage is returning what it thinks is are elements of LaurentPolynomialRing(QQ) rather then elements of LaurentPolynomialRing (ZZ) as requested. I think this is wrong.

Andrew

Andrew

[Simon King]

Hi Andrew,

On 2013-09-19, Andrew Mathas <andrew...@gmail.com> wrote:
> Thanks for the clarification Nils. So Laurent polynomials are not playing
> well with their coefficients, but I wonder if the real issue is that even
> polynomials over rings have not been implemented carefully eough. Consider:
>
> sage: A.<x>=PolynomialRing(ZZ)
> sage: x in A, x^-1 in A
> (True, False)
> sage: x^-1
> 1/x
>
> Everything is correct except, I think, the last line: x^-1 it NOT defined
> so this should be returning an error rather than coercing into the field of
> fractions.

The polynomial ring is an integral domain, thus, it HAS a fraction
field. So, I don't see why one should not pass to the fraction field
upon division. It is like 2^-1 yielding a rational number, even though 2
is an integer.

> If I want x^-1 to be defined then I should define it properly.

Do you really want that the user defines 2^-1 before allowing him/her to
do division of integers???

> I
> think do not think that sage should try and second guess what I mean and
> put the result into a different ring.

Why not? You asked for it!

Really: I think we agree that division of integers should work without
an extra definition. So, why do you want to forbid divisions for other
integral domains? That would be plainly inconsistent.

Best regards,
Simon

[Nils Bruin]

On Thursday, September 19, 2013 1:34:02 AM UTC-7, Andrew Mathas wrote:

 

Everything is correct except, I think, the last line: x^-1 it NOT defined so this should be returning an error rather than coercing into the field of fractions.

If you're mainly interested in laurent polynomials then perhaps, but if you're also interested in rational functions (as Simon pointed out, the field of fractions of a polynomial ring) then probably not.

There is a "special" division, "//", which is normally "euclidean division", but one could (and there may be examples of it already) where it's "divide if the quotient is in the ring and throw an error otherwise". Then you'd be able to write

1+u+v+2//(u*v)

probably also not something you're waiting for ...

In the mean time, you can do

sage: def norm_coeffs(p):
....:     return parent(p)({k:ZZ(c) for k,c in p.dict().iteritems()})
sage: P.<u,v>=LaurentPolynomialRing(ZZ,2)
sage: p=norm_coeffs(2*u**-1*v**-1+u*v)
sage: ui=norm_coeffs(u^(-1))
sage: vi=norm_coeffs(v^(-1))
sage: p(ui,vi)
2*u*v + u^-1*v^-1
sage: q= u+v+ui
sage: q(ui,vi)

which is not great, but not as heavy-handed as your other solution. Plus, ui is a couple of characters shorter than u^(-1).

By the way, in laurent polynomial rings, u^(-1) and 1/u do not follow the same code path:

sage: parent(u^(-1))

Multivariate Laurent Polynomial Ring in u, v over Integer Ring

sage: parent(1/u)
Fraction Field of Multivariate Laurent Polynomial Ring in u, v over Integer Ring

so wherever the exponentiation code for laurent polynomials lives, it seems to special-case inversion of one-term polynomials. If we just add to that special case the check if the coefficient is invertible in the base ring, we can probably avoid ending up with rational numbers all the time.

In fact, in laurent_polynomial, line 202 we have

        if n < 0:
            E = self._poly.exponents()
            if len(E) == 0:
                raise ZeroDivisionError
            elif len(E) > 1:
                return self.parent()._R.fraction_field()(self)**n
            else:
                ans = self._new_c()
                ans._poly = self._poly.parent().change_ring(self.parent().base_ring().fraction_field())(self._poly.coefficients()[0]**n)
                ans._mon = self._mon.eadd(E[0]).emul(n)

so the coercion of the coefficient happens rather explicitly (and wrongly! If you're going to change the ring of a coefficient you should change the polynomial ring as well).

[Travis Scrimshaw]

Hey everyone,

<shameless needs review plug>

   At least some of this for univariate Laurent rings is fixed in #11726 (although this currently depends on #14261, which is where this problem comes from):

sage: R.<t> = LaurentPolynomialRing(ZZ)

sage: 1 / t

t^-1

</end shameless needs review plug>

Let me also make explicit the question: Given ZZ[t, t^-1], do we want 1 / (2*t) to be in the fraction field or in QQ[t, t^-1]?

I'm leaning towards the latter since it is a "smaller/more simple" ring (and I realize now that #11726 currently does the former).

   However I think there is a bigger question: if we have a / b and b divides a, do we want, by default, to return an object in the original ring and not its fraction field? I think yes and since the fraction field would run quo_rem() anyways, but we could run it in the _div_ method and return an object in the original ring if the remainder is 0. Thus there should not be a slow-down. The fact that exact division ends up in the fraction field can irritate me when using polynomials that I know results in a polynomial, because methods I expect are not there, and I have forgotten about the //.

Best,

Travis

[Nils Bruin]

On Thursday, September 19, 2013 10:10:26 AM UTC-7, Travis Scrimshaw wrote:

Let me also make explicit the question: Given ZZ[t, t^-1], do we want 1 / (2*t) to be in the fraction field or in QQ[t, t^-1]?

I'm leaning towards the latter since it is a "smaller/more simple" ring (and I realize now that #11726 currently does the former).

It currently lies in the fraction field and that is consistent with other choices in sage: The parent of the result of a/b depends on the parents of a,b; not on the particular values of a,b. If you don't do that kind of thing, you'll soon find your code littered with tests for where the result you just got actually lives.

 

For LaurentPolynomials to be consistent with the rest of sage, a/b must lie in the fraction field, even though by the time you've constructed laurent polynomials, you're obviously not interested in the fration field. Really ZZ[x,x^(-1)] doesn't have division: it just happens to have an extra unit x with inverse "x^(-1)".

There is precedent for a^(-1) being a little more varied. In particular, for Laurent polynomials to be usable we must have that x^(-n) is short-hand for "the n-th power of the inverse of x". But I think redefining "/" to be at odds with the rest of sage is going too far [and bringing the rest of sage in line with LaurentPolynomials is not going to happen due to backwards compatibility/inertia, even if it were desirable]

Note that a^(-1) for matrices seems to prefer extending scalars in a regular way rather than a minimal way:

sage: (matrix([[1,1],[0,1]]).is_invertible())
True
sage: (matrix([[2,1],[0,1]]).is_invertible()) #determinant not a unit in ZZ
False
sage: parent(matrix([[1,1],[0,1]])^(-1))
Full MatrixSpace of 2 by 2 dense matrices over Rational Field
sage: parent(matrix([[2,1],[0,1]])^(-1))
Full MatrixSpace of 2 by 2 dense matrices over Rational Field

so it looks like "inverse over a non-field" isn't really implemented for matrices.