random subs() surprise that has just kept me guessing for over an hour

[Kevin Buzzard]

This took me a fair while to track down.

sage: m=matrix(GF(3),[[0,1,0,0,0,0],[0,0,1,0,0,0],[0,0,0,1,0,0],[0,0,0,0,1,0],[0,0,0,0,0,1],[2,2,0,1,0,0]]) # random 6x6 matrix
sage: g=m.charpoly('x')
sage: g
x^6 + x^3 + 2*x + 2

[I've just build a degree 6 poly. Now let's build a degree 12 one]

sage: h=expand((g.subs(x+2/x))*x^6)

[easy check shows that this will have degree 12]

sage: h
2*x^5 + x^3 + 1

[gaargh]

Am I some sort of a victim of some secret property of the letter 'x'? I wanted to use x rather than another letter because that would save me from having to define the ring GF(3)[x,1/x] which would have involved some thinking on my part ;-)

PS change x+2/x to x+GF(3)(2)/x and it works fine ;-)

Kevin

[Simon King]

Hi Kevin,

On 2014-07-27, Kevin Buzzard <kevin.m...@gmail.com> wrote:
> [I've just build a degree 6 poly. Now let's build a degree 12 one]
>
> sage: h=expand((g.subs(x+2/x))*x^6)

Let's work without the x^6 factor:

sage: g
x^6 + 2*x^3 + x + 1
sage: g.subs(x+2/x).expand()
2/x + 1/x^3 + 1/x^6

No idea what is going wrong here. Let's try different ways to
substitute:

sage: g(x=x+2/x).expand()
2/x + 1/x^3 + 1/x^6

Same problem. Whithout expanding, it seems to be fine:

sage: g.subs(x+2/x)
(((x + 2/x)^3 + 2)*(x + 2/x)^2 + 1)*(x + 2/x) + 1

I guess the above form comes from Horner's scheme.

Now, really strange: If one copy-and-pastes the above expression and
expands it, then all is fine!

sage: ((((x + 2/x)^3 + 2)*(x + 2/x)^2 + 1)*(x + 2/x) + 1).expand()
x^6 + 12*x^4 + 2*x^3 + 60*x^2 + 13*x + 26/x + 240/x^2 + 16/x^3 + 192/x^4 + 64/x^6 + 161

That is all very puzzling to me---and is yet another reason for my creed
that one should use *proper* polynomials and not symbolic expressions
(see below for the underlying types).

> Am I some sort of a victim of some secret property of the letter 'x'?

Not of the letter 'x', but of the fact that the variable x is
pre-defined in Sage as a symbolic expression, which often enough causes
trouble to those who expect to work with polynomials. Note that x is *not*
the generator of a polynomial ring:

sage: type(x)
<type 'sage.symbolic.expression.Expression'>

whereas your polynomial g is a proper polynomial (there are various underlying
implementations of polynomials in Sage, for different purposes, and this
one relies on FLINT):

sage: type(g)
<type 'sage.rings.polynomial.polynomial_zmod_flint.Polynomial_zmod_flint'>

> I
> wanted to use x rather than another letter because that would save me from
> having to define the ring GF(3)[x,1/x] which would have involved some
> thinking on my part ;-)

Well, the polynomial ring already is there:

sage: g.parent()
Univariate Polynomial Ring in x over Finite Field of size 3

and the fraction field will automatically be created when needed. So,
you could do:

sage: g.parent().inject_variables()
Defining x

The preceding line works if you are in an interactive session.
Alternatively (or if you write a program), you could explicitly define
it:

sage: x = g.parent().gen()
sage: type(x)
<type 'sage.rings.polynomial.polynomial_zmod_flint.Polynomial_zmod_flint'>

Now, the substitution works out of the box, and "expansion" of symbolic
expressions is not needed, since we work with proper polynomials and
their quotients.

sage: g.subs(x+2/x)
(x^12 + 2*x^9 + x^7 + 2*x^6 + 2*x^5 + x^3 + 1)/x^6

My general recommendation is: Unless you want to do calculus, get used
to define x as the generator of an appropriate polynomial ring, taking
care of the underlying range of coefficients. This is totally easy in an
interactive session:

sage: P.<x> = QQ[] # or GF(3)['x'] in your case

The above line both defines P and x, as polynomial ring over the
rationals, and its generator. In a program, you should write

P = QQ['x'] # or GF(3)['x'] in your case
x = P.gen()

instead.

Best regards,
Simon

[John Cremona]

Minor remark: if you have defined x to be a polynomial ring generator
then Sage will happily form rational functions without you having to
define new rings:

sage: R.<x> = GF(3)[]
sage: R

Univariate Polynomial Ring in x over Finite Field of size 3

sage: f = x+1/x
sage: f.parent()
Fraction Field of Univariate Polynomial Ring in x over Finite Field of size 3

> --
> You received this message because you are subscribed to the Google Groups "sage-support" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.
> To post to this group, send email to sage-s...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sage-support.
> For more options, visit https://groups.google.com/d/optout.

[William Stein]

On Sun, Jul 27, 2014 at 12:52 PM, John Cremona <john.c...@gmail.com> wrote:
> Minor remark: if you have defined x to be a polynomial ring generator
> then Sage will happily form rational functions without you having to
> define new rings:
>
> sage: R.<x> = GF(3)[]
> sage: R
> Univariate Polynomial Ring in x over Finite Field of size 3
> sage: f = x+1/x
> sage: f.parent()
> Fraction Field of Univariate Polynomial Ring in x over Finite Field of size 3

Or, as David Kohel used to always remind me, making "/" something that
always constructs something in the fraction field was a good design
decision in Magma, which he made sure I maintained in Sage.

>>> R<x> := PolynomialRing(GF(3));
>>> Parent(1/x);
Univariate rational function field over GF(3) Variables: x
>>> Parent(2/1);
Rational Field

--
William Stein
Professor of Mathematics
University of Washington
http://wstein.org

[Nils Bruin]

I don't think we'll ever get SR to operate properly in positive characteristics; especially because it would allow completely arbitrary characteristic combinations in the first place, but perhaps the cases below help in tracking down if we can so something to improve the situation a bit:

sage: ((1)*((1)*x+(1)/x)^3+(1)).expand()
x^3 + 3*x + 3/x + 1/x^3 + 1
sage: ((k(1)*x+k(1)/x)^3+k(1)).expand()
x^3 + 3*x + 3/x + 1/x^3 + 1
sage: (k(1)*(k(1)*x+k(1)/x)^3).expand()
x^3 + 1/x^3
sage: (k(1)*(k(1)*x+k(1)/x)^3+k(1)).expand()
0/x + 1/x^3 + 1
sage: ((k(1)*x+k(1)/x)^3+k(1)).expand()
x^3 + 3*x + 3/x + 1/x^3 + 1

so it's the combination of multiplying the third power by k(1) and adding k(1) to it that leads "expand" to producing an expression that is really invalid.

My guess is that it finds a nonzero factor 3 somewhere, which it assumes is invertible, and ends up multiplying both numerator and denominator of some fraction somewhere. The numerator ends up being combined in GF(3), so that one gets 0 and SR takes it happily from there, removing the factor 3 from the denominator again.

[Kevin Buzzard]

Thanks to everyone who commented. I am a fool. I got tangled up in some similar sort of situation a few months ago and David Loeffler emailed me to basically warn me off sage's x and to use anything else instead. I didn't follow his advice this time because I couldn't get an R.<t> constructor to work in a script because I was using execfile instead of %runfile. We live and learn.

Kevin