Creating a polynomial with symbolic coefficients

[James Parson]

Dear sage-support group,

I am completely new to computer algebra systems and to computer
programming, and I hope you'll indulge the following beginner's
question. I was wondering if there is a simple way to create a
polynomial of degree d in x and y with symbolic coefficients in Sage.
Here is what I mean: if I were at the board in class, I might write
(in LaTeX transcription) something like

P(x,y) = \sum_{i+j\leq d} a_{ij} x^i y^j,

which I would view as an element of \mathbf{Z}[a_{ij},x,y]. I might
then impose some linear conditions on the a_{ij} by insisting that P
(x_t,y_t) = 0 for a list of points

(x_1,y_1), (x_2,y_2), ... .

Finally, I might solve the resulting system of linear equations.

How would you recommend that I set up something like the a_{ij} and P
(x,y) in Sage? In order to make the question more definite, I
illustrate it with an example that I took from a lecture of Doron
Zeilberger on experimental mathematics. He proposed the question of
finding a polynomial of degree d in x and y that vanishes when x and y
are specialized to consecutive Fibonacci numbers. The lines below are
my attempt at a Sage version of his suggested computer search for a
likely solution (originally written in Maple). The program should take
the degree d as an input and then provide a parameterized family of
polynomials of degree d that are likely candidates.

Here is what I came up with, after an enlightening afternoon of
studying computer manuals:

d = 4
e = d+1
L = []
M = []
for i in range(e):
for j in range(e-i):
L.append('a_%s_%s' %(i,j))
M.append([i,j])
V = var(' '.join(L))
P = sum(V[j]*x^(M[j][0])*y^(M[j][1]) for j in range(len(L)))
E = [P(x=fibonacci(n),y=fibonacci(n+1)) for n in range(1,len(V)+6)]
P.substitute(solve(E,V,solution_dict = True)[0])

I could not figure out how to create and reference the variables a_
{ij} conveniently, and so I ended up with the strange lists V and M
above. Even though I got my polynomial P and solved the original
problem to my satisfaction, I still don't think I know how I would
have Sage do something like sum the a_{i.i+1} for 2i+1<=d. Is there a
better way to do this sort of thing?


Thanks for your help and indulgence,

James Parson

[David Joyner]

I'm not sure if this helps, but you can create a polynomial
of the type you want a bit simpler:

sage: var("x,y")
(x, y)
sage: Inds = CartesianProduct(range(5), range(4))
sage: sum([var("a"+str(i)+str(j))*x^i*y^j for i,j in Inds])
a43*x^4*y^3 + a33*x^3*y^3 + a42*x^4*y^2 + a23*x^2*y^3 + a32*x^3*y^2 +
a41*x^4*y + a13*x*y^3 + a22*x^2*y^2 + a31*x^3*y + a40*x^4 + a03*y^3 +
a12*x*y^2 + a21*x^2*y + a30*x^3 + a02*y^2 + a11*x*y + a20*x^2 + a01*y
+ a10*x + a00

Now you can reference them on the fly like this:

sage: eval("a"+str(3)+str(2))
a32

[Robert Bradshaw]

Currently symbolic variables are un-indexable. What would people
think of having indexing create new subscripted variables?

sage: a = var('a')
sage: a[0]
a_0
sage: latex(a[1,2])
a_{1,2}

- Robert

[William Stein]

On Wed, Jun 3, 2009 at 11:45 AM, Robert Bradshaw
<robe...@math.washington.edu> wrote:
>
> Currently symbolic variables are un-indexable. What would people
> think of having indexing create new subscripted variables?
>
> sage: a = var('a')
> sage: a[0]
> a_0
> sage: latex(a[1,2])
> a_{1,2}

That's a pretty wild and crazy idea. Cool. Does any other math
software do that?
Are there any obvious gotcha's?

William

[Mike Hansen]

On Wed, Jun 3, 2009 at 11:53 AM, William Stein <wst...@gmail.com> wrote:
>> Currently symbolic variables are un-indexable. What would people
>> think of having indexing create new subscripted variables?
>>
>> sage: a = var('a')
>> sage: a[0]
>> a_0
>> sage: latex(a[1,2])
>> a_{1,2}
>
> That's a pretty wild and crazy idea.  Cool.  Does any other math
> software do that?
> Are there any obvious gotcha's?

I think any system where things are left unevaluated can do something
like this. For example, I know people who use this in Maple to have
polynomials infinitely many variables. We have something like this in
Sage now:

sage: R.<x> = InfinitePolynomialRing(QQ)
sage: x[0] + x[100]
x100 + x0

I would still like to use __getitem__ to return the operands of a
symbolic expression. Since symbols have no operands, these two things
aren't totally incompatible.

--Mike

[Harald Schilly]

On Jun 3, 8:45 pm, Robert Bradshaw <rober...@math.washington.edu>
wrote:

> Currently symbolic variables are un-indexable. What would people  
> think of having indexing create new subscripted variables?
>
> sage: a = var('a')
> sage: a[0]
> a_0
> sage: latex(a[1,2])
> a_{1,2}
>

I like it, this idea could also be expanded to vectors, like
a[0:3] == vector([a_0, a_1, a_2, a_3]) # note, also a_3

and probably, the var constructor should also be more intelligent,
i.e. var('a[0:3]') creates all the a_i (and avoids the unindexed 'a'
you would need in your 2-step approach above)

h

[Jason Grout]

So as far as printing, a[0] would look the same as a0 would look the
same as a_0? Would a[0] actually be the variable a_0 or a0?

Do we ever want to make symbolic expressions indexable? If so, it would
be confusing to have:

(x+1)[0]

have totally different behavior than

(x)[0].

Jason

P.S. It seems like Maple did something like this---Maple experts can
comment on it, though.

--
Jason Grout

[Jason Grout]

Harald Schilly wrote:
> On Jun 3, 8:45 pm, Robert Bradshaw <rober...@math.washington.edu>
> wrote:
>> Currently symbolic variables are un-indexable. What would people
>> think of having indexing create new subscripted variables?
>>
>> sage: a = var('a')
>> sage: a[0]
>> a_0
>> sage: latex(a[1,2])
>> a_{1,2}
>>
>
> I like it, this idea could also be expanded to vectors, like
> a[0:3] == vector([a_0, a_1, a_2, a_3]) # note, also a_3

-1 to the a_3. That goes against all of the slicing conventions in Python.

It is kind of a cool idea, though. I'd like some time to try it out
before committing one way or the other.

What about having a "experiment mode" in Sage that turns on things like
this? Some variable in some module somewhere that people can set to
switch on some experimental behavior so they can test it out and give
feedback. In other words, setting sage.misc.misc.EXPERIMENT_MODE=True
turns it on.

Jason

--
Jason Grout

[Robert Bradshaw]

On Jun 3, 2009, at 12:11 PM, Jason Grout wrote:

> William Stein wrote:
>> On Wed, Jun 3, 2009 at 11:45 AM, Robert Bradshaw
>> <robe...@math.washington.edu> wrote:
>>> Currently symbolic variables are un-indexable. What would people
>>> think of having indexing create new subscripted variables?
>>>
>>> sage: a = var('a')
>>> sage: a[0]
>>> a_0
>>> sage: latex(a[1,2])
>>> a_{1,2}
>>
>> That's a pretty wild and crazy idea. Cool. Does any other math
>> software do that?
>> Are there any obvious gotcha's?
>
> So as far as printing, a[0] would look the same as a0 would look the
> same as a_0? Would a[0] actually be the variable a_0 or a0?

I'm not sure. My first intent was that a[0] would actually be a0, but
it's unclear how to format multiple indices (I'd want a[1, 23] != a
[12, 3]). Also, should we support a[i]? Then I'd want a[i].subs(i=5)
== a[5]. What about v[5].subs(v=vector(range(10)))?

> Do we ever want to make symbolic expressions indexable? If so, it
> would
> be confusing to have:
>
> (x+1)[0]
>
> have totally different behavior than
>
> (x)[0].

They're not now. Having both would be confusing. I'd vote for the
latter--if (x+1)[0] worked, would it be the same as (1+x)[0], or (x+1-
x)[0]?

Answering your question about experimental mode, you're talking about
something easier and more permanent than applying a patch from trac?
Would people start depending on it, meaning we can't remove it
without sacrificing backwards compatibility (despite the name
"experimental")?

- Robert

[John H Palmieri]

On Jun 3, 12:16 pm, Jason Grout <jason-s...@creativetrax.com> wrote:
> What about having a "experiment mode" in Sage that turns on things like
> this?  Some variable in some module somewhere that people can set to
> switch on some experimental behavior so they can test it out and give
> feedback.  In other words, setting sage.misc.misc.EXPERIMENT_MODE=True
> turns it on.

It seems better to me to turn individual experimental features on or
off. So what about an argument to var, just like "ns"?

var(a, index=True)

for example. ("index" doesn't sound right to me, but I hope you get
the idea.)

John

[David Joyner]

+1

Is this possible?

sage: J = CartesianProduct(range(3),range(4))
sage var(a, index_set = J)


>
>  John
>
> >
>

[James Parson]

Thanks to David Joyner for his response to my original question. His
method worked nicely. Incidentally, here is the original Maple code
from the lecture of Doron Zeilberger that I was trying to translate
into Sage:

with(combinat): P:=(d,x,y)->add(add(a[i,j]*x**i*y**j,i=0..d-
j),j=0..d);
V:=d->fseq(seq(a[i,j],i=0..d-j),j=0..d)g;
E:=d->fseq(P(d,fibonacci(n),fibonacci(n+1)),n=1..nops(V(d))+5) g:
Q:=(d,x,y)->subs(solve(E(d),V(d)),P(d,x,y));

These lines feature the sort of indexed variables a[i,j] discussed
above.

(The full lecture from which I took these lines can be found at
http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/em.html.)

Here is a variant on the original question: suppose I wanted to write
a line that creates a polynomial ring whose variables are a_{ij} for i
+j<=d. How should I do it? I might want to set this up, for example,
so that I could tell Sage about an algebraic group action on the space
of polynomials of degree <=d. For a simpler variant: is there a
convenient way to construct QQ[x_{ij}] with 1\leq i,j\leq n? I am a
overwhelmed with the various ways to construct a polynomial ring, and
so I cannot tell if one of them would be appropriate for this purpose.
I can see how to make a polynomial ring in n^2 variables, but I do not
know how to name them x_{ij}.


Thanks again for your help,

James Parson

[David Joyner]

On Wed, Jun 3, 2009 at 5:54 PM, James Parson <par...@hood.edu> wrote:
>
> Thanks to David Joyner for his response to my original question. His
> method worked nicely. Incidentally, here is the original Maple code
> from the lecture of Doron Zeilberger that I was trying to translate
> into Sage:

BTW I think the more of Zeilberger's stuff that is translated into
Sage the better!
Please consider publishing your translation as a sagemath.org notebook
worksheet.

>
> with(combinat): P:=(d,x,y)->add(add(a[i,j]*x**i*y**j,i=0..d-
> j),j=0..d);
> V:=d->fseq(seq(a[i,j],i=0..d-j),j=0..d)g;
> E:=d->fseq(P(d,fibonacci(n),fibonacci(n+1)),n=1..nops(V(d))+5) g:
> Q:=(d,x,y)->subs(solve(E(d),V(d)),P(d,x,y));
>
> These lines feature the sort of indexed variables a[i,j] discussed
> above.
>
> (The full lecture from which I took these lines can be found at
> http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/em.html.)
>
> Here is a variant on the original question: suppose I wanted to write
> a line that creates a polynomial ring whose variables are a_{ij} for i
> +j<=d. How should I do it? I might want to set this up, for example,

sage: Inds = CartesianProduct(range(5), range(5))
sage: vars = ["a"+str(i)+str(j) for i,j in Inds]
sage: PolynomialRing(QQ,25,vars)
Multivariate Polynomial Ring in a00, a01, a02, a03, a04, a10, a11,
a12, a13, a14, a20, a21, a22, a23, a24, a30, a31, a32, a33, a34, a40,
a41, a42, a43, a44 over Rational Field

[kcrisman]

>
> > Here is a variant on the original question: suppose I wanted to write
> > a line that creates a polynomial ring whose variables are a_{ij} for i
> > +j<=d. How should I do it? I might want to set this up, for example,
>
> sage: Inds = CartesianProduct(range(5), range(5))
> sage: vars = ["a"+str(i)+str(j) for i,j in Inds]

I think it would actually be something like

sage: vars = ["a"+str(i)+str(j) for i,j in Inds if i+j<5]

Of course you'd have to change 25 below. Actually, I think the number
of variables is optional, though if you include it you need to
calculate the correct number of them :)

> sage: PolynomialRing(QQ,25,vars)
> Multivariate Polynomial Ring in a00, a01, a02, a03, a04, a10, a11,
> a12, a13, a14, a20, a21, a22, a23, a24, a30, a31, a32, a33, a34, a40,
> a41, a42, a43, a44 over Rational Field

- kcrisman

[Bill Page]

On Wed, Jun 3, 2009 at 5:54 PM, James Parson wrote:
>
> Thanks to David Joyner for his response to my original question. His
> method worked nicely. Incidentally, here is the original Maple code
> from the lecture of Doron Zeilberger that I was trying to translate
> into Sage:
>
> with(combinat): P:=(d,x,y)->add(add(a[i,j]*x**i*y**j,i=0..d-
> j),j=0..d);
> V:=d->fseq(seq(a[i,j],i=0..d-j),j=0..d)g;
> E:=d->fseq(P(d,fibonacci(n),fibonacci(n+1)),n=1..nops(V(d))+5) g:
> Q:=(d,x,y)->subs(solve(E(d),V(d)),P(d,x,y));
>
> These lines feature the sort of indexed variables a[i,j] discussed
> above.
>
> (The full lecture from which I took these lines can be found at
> http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/em.html.)
>

There is some pdf cut-and-paste corruption of the example above. The
correct Maple code is:

with(combinat):
P:=(d,x,y)->add(add(a[i,j]*x**i*y**j,i=0..d-j),j=0..d);

V:=d->{seq(seq(a[i,j],i=0..d-j),j=0..d)};
E:=d->{seq(P(d,fibonacci(n),fibonacci(n+1)),n=1..nops(V(d))+5) };
Q:=(d,x,y)->subs(solve(E(d),V(d)),P(d,x,y));

A more direct translation to Sage might be something like this:

sage: P=lambda d,x,y: sum([ sum([ var("a"+str(i)+str(j))*x^i*y^j for i
in [0..d-j]]) for j in [0..d]])
sage: V=lambda d:sum([[var("a"+str(i)+str(j)) for i in [0..d-j]] for j
in [0..d]],[])
sage: E=lambda d: [ P(d,fibonacci(n),fibonacci(n+1)) for n in [1..len(V(d))+5] ]
sage: Q=lambda d,x,y:P(d,x,y).subs_expr(solve(E(d),V(d),solution_dict=True)[0])

sage: Q(1,x,y)
0
sage: Q(2,x,y)
0
sage: Q(3,x,y)
0
sage: Q(4,x,y).factor()
r21*(-y^2 + x*y + x^2 - 1)*(-y^2 + x*y + x^2 + 1)

Regards,
Bill Page.

[James Parson]

> > Here is a variant on the original question: suppose I wanted to write
> > a line that creates a polynomial ring whose variables are a_{ij} for i
> > +j<=d. How should I do it? I might want to set this up, for example,
>
> sage: Inds = CartesianProduct(range(5), range(5))
> sage: vars = ["a"+str(i)+str(j) for i,j in Inds]
> sage: PolynomialRing(QQ,25,vars)
> Multivariate Polynomial Ring in a00, a01, a02, a03, a04, a10, a11,
> a12, a13, a14, a20, a21, a22, a23, a24, a30, a31, a32, a33, a34, a40,
> a41, a42, a43, a44 over Rational Field

Thanks again for the suggestions. I have one more foolish question
about this sort of construction: if I type those lines into Sage and
then type something like

a00 + a11,

I get an error

NameError: name 'a00' is not defined.

I read about this sort of thing in the Sage Tutorial, but I couldn't
understand it well enough to figure out how to name the variables what
I wanted. Is there any easy way to do this?


Regards,

James Parson

[David Joyner]

Sorry, I'm stuck here too. Can you just write R("a00")+R("a11") instead?



>
>
> Regards,
>
> James Parson
> >
>

[William Stein]

Two options:

(1) Just type inject_on() and then aij will be defined:

sage: inject_on()
Redefining: FiniteField Frac FractionField FreeMonoid GF
LaurentSeriesRing NumberField PolynomialRing quo quotient

sage: Inds = CartesianProduct(range(5), range(5))
sage: vars = ["a"+str(i)+str(j) for i,j in Inds]
sage: PolynomialRing(QQ,25,vars)

Defining a00, a01, a02, a03, a04, a10, a11, a12, a13, a14, a20, a21,

a22, a23, a24, a30, a31, a32, a33, a34, a40, a41, a42, a43, a44
Multivariate Polynomial Ring in a00, a01, a02, a03, a04, a10, a11,
a12, a13, a14, a20, a21, a22, a23, a24, a30, a31, a32, a33, a34, a40,
a41, a42, a43, a44 over Rational Field

sage: a00 + a01
a00 + a01
sage: (a00 + a01)^3
a00^3 + 3*a00^2*a01 + 3*a00*a01^2 + a01^3

(2) Edit the globals dictionary:

sage: Inds = CartesianProduct(range(5), range(5))
sage: vars = ["a"+str(i)+str(j) for i,j in Inds]

sage: R = PolynomialRing(QQ,25,vars)
sage: for v in R.gens(): globals()[str(v)] = v
....:
sage: (a00 + a01)^3
a00^3 + 3*a00^2*a01 + 3*a00*a01^2 + a01^3

[Robert Dodier]

William Stein wrote:

> On Wed, Jun 3, 2009 at 11:45 AM, Robert Bradshaw
> <robe...@math.washington.edu> wrote:
> >
> > Currently symbolic variables are un-indexable. What would people
> > think of having indexing create new subscripted variables?

> That's a pretty wild and crazy idea. Cool. Does any other math
> software do that?

For the record, Maxima treats subscripted variables somewhat the
same as simple variables. (They should be the same but Maxima
is less than entirely consistent ....) A subscripted variable x[0] is
distinct from x_0 and x0.

Subscripted variables are the same as unevaluated function-like
expressions except that they have an extra flag which shows that
it's a subscript instead of a function argument.
I think it's useful to consider subscripted variables as a subset
of functional expressions; after all a subscripted variable is
function
which maps its set of indices to whatever.

> Are there any obvious gotcha's?

I think of two general shortcomings of Maxima's scheme.
One is that x[m] and x[n] have no known relation when m and n
are different; there isn't any way to apply some property across
all of the subscripted variables x[foo].
The other is that x[foo] could be any of several things ---
could be a list element, a matrix row, an array element, a
hash table element, as well as a subscripted variable.
This multitude of interpretations of x[foo] can lead to confusion.

FWIW

Robert Dodier