make Expression.series return an element of SR[[]]

[rws]

Instead of duplicating power series functionality (existing in rings/) in Pynac, where there is momentarily only a skeleton implementation that does not even play nice with the rest of symbolics, the series method of Expression should return an element of PowerSeriesRing(SR).

It does not necessarily mean everything will be implemented in Python, Pynac could still compute the series as Python object.

[nbruin]

• What's the variable of SR[[]]? Does it support multivariate power series?
• How does f=(x^2+y^2+1) coerce into SR[['x']]? Note that f is already an element of SR, so by default it would be a "constant" in SR[['x']].

[rws]

Replying to nbruin:

• What's the variable of SR[[]]?

A symbolic variable is given as argument of series. The ring variable will have the same name and thus override it on new input from the user.

Does it support multivariate power series?

At the moment series takes one symbol.

• How does f=(x^2+y^2+1) coerce into SR[['x']]? Note that f is already an element of SR, so by default it would be a "constant" in SR[['x']].

As above the series argument determines the name of the series variable, so it also determines what is being done with the expression. Already now:

sage: (1+x+y).series(x,5)
(y + 1) + 1*x

(However, the coefficient functions don't work with "exact" series atm)

[nbruin]

Replying to rws:

• How does f=(x^2+y^2+1) coerce into SR[['x']]? Note that f is already an element of SR, so by default it would be a "constant" in SR[['x']].

As above the series argument determines the name of the series variable, so it also determines what is being done with the expression. Already now:

That's ambiguous and not what happens by default in our current coercion:

sage: R=PowerSeriesRing(SR,name='x',default_prec=10)
sage: F=1/(1-R.0)
sage: F-x
-x + 1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + O(x^10)
sage: F-R.0
1 + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + O(x^10)

You see that SR('x') doesn't coerce to R.0 in R. You'd have to override that somehow and that will be extremely messy: where does sin(x) coerce?

[rws]

Replying to nbruin:

That's ambiguous and not what happens by default in our current coercion:

sage: R=PowerSeriesRing(SR,name='x',default_prec=10)
sage: F=1/(1-R.0)
sage: F-x
-x + 1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + O(x^10)
sage: F-R.0
1 + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + O(x^10)

You see that SR('x') doesn't coerce to R.0 in R. You'd have to override that somehow and that will be extremely messy: where does sin(x) coerce?

To sin(R.0), where is the problem? Of course it would expand to

x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + 1/362880*x^9 - 1/39916800*x^11 + 1/6227020800*x^13 - 1/1307674368000*x^15 + 1/355687428096000*x^17 - 1/121645100408832000*x^19 + Order(x^20)

[nbruin]

Replying to rws:

where does sin(x) coerce? To sin(R.0), where is the problem? Of course it would expand to

x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + 1/362880*x^9 - 1/39916800*x^11 + 1/6227020800*x^13 - 1/1307674368000*x^15 + 1/355687428096000*x^17 - 1/121645100408832000*x^19 + Order(x^20)

That will be quite a bit of work to arrange, but that's not what I meant: What are the results of

SR[['x']]( sin(SR('x')) )
SR[['x']]( sin(SR('y')) )

(the main part: how are you going to arrive at those results? I'm not saying it's impossible to come up with something satisfactory. Just that there's some non-trivial work to be done)

[rws]

Replying to nbruin:

...how are you going to arrive at those results? I'm not saying it's impossible to come up with something satisfactory. Just that there's some non-trivial work to be done)

It's just another expression converter, isn't it? BTW, exp(x) already does the desired thing. Though it's not dependent the work would also benefit from #17759 which itself depends on #18255.

[rws]

Ah I see, sin(x) needs a second order DE. But Pari is able to do sin easily so there must be an algorithm.

[rws]

Power series via Pari would be #15601, but actually Pynac is more complete here. It just needs the conversion to PowerSeries.

[rws]

I think my new posting in sage-devel (​https://groups.google.com/forum/?hl=en#!topic/sage-devel/XhOhwjjeLTM) makes this ticket obsolete. Maybe we just need to change Pynac so that, if an expression contains a pseries it will be reexpanded (making the new order the min of all orders). What do you think?

Surprise.

[rws]

I think there is also a satisfying solution to the discussed interface problem. Instead of automatically creating a PowerSeriesRing, with all the problems of duplicate generators, the user will do this part and give the chosen generator as argument to series or simply use it as part of the expression.

[mmezzarobba]

Replying to rws:

I think my new posting in sage-devel (​https://groups.google.com/forum/?hl=en#!topic/sage-devel/XhOhwjjeLTM) makes this ticket obsolete. Maybe we just need to change Pynac so that, if an expression contains a pseries it will be reexpanded (making the new order the min of all orders). What do you think?

I'm not too convinced: calling series() is a bit like calling expand() or normal(), with symbolic expressions, you don't necessarily want that to happen automatically. And (FWIW) that's also what Maple does:

#-->s := series(sin(x),x);
s := series(1*x-1/6*x^3+1/120*x^5+O(x^7),x,7)
#-->t := series(exp(x),x);
t := series(1+1*x+1/2*x^2+1/6*x^3+1/24*x^4+1/120*x^5+O(x^6),x,6)
#-->s*t;
(series(1*x-1/6*x^3+1/120*x^5+O(x^7),x,7))*(series(1+1*x+1/2*x^2+1/6*x^3+1/24*x^4+1/120*x^5+O(x^6),x,6))
#-->series(s*t,x);
series(1*x+1*x^2+1/3*x^3-1/30*x^5+O(x^6),x,6)

But perhaps I misunderstood the question?

In any case, I agree that this ticket can be closed.

[rws]

Replying to mmezzarobba:

...And (FWIW) that's also what Maple does:

#-->s := series(sin(x),x);
s := series(1*x-1/6*x^3+1/120*x^5+O(x^7),x,7)
#-->t := series(exp(x),x);
t := series(1+1*x+1/2*x^2+1/6*x^3+1/24*x^4+1/120*x^5+O(x^6),x,6)
#-->s*t;
(series(1*x-1/6*x^3+1/120*x^5+O(x^7),x,7))*(series(1+1*x+1/2*x^2+1/6*x^3+1/24*x^4+1/120*x^5+O(x^6),x,6))
#-->series(s*t,x);
series(1*x+1*x^2+1/3*x^3-1/30*x^5+O(x^6),x,6)

But perhaps I misunderstood the question?

Interesting. But that can be considered a bug and see the sage-devel link for a bug report. I mean, what else do you expect as result when dealing with a series? But I'll open a new ticket for that.

[nbruin]

Replying to rws:

Interesting. But that can be considered a bug and see the sage-devel link for a bug report. I mean, what else do you expect as result when dealing with a series? But I'll open a new ticket for that.

The problem is: there are SR allows expressions of such generality that series expansions aren't always possible. Consider:

sage: function('f')
sage: s=series(sin(x))
sage: f(cos(x))+g

Should that be considered a series in x?

And then there are of course the cases with well-defined functions as ingredients where you really can't make a series expansion:

sage: s=series(sin(x))
sage: exp(1/x)+s

The general problem is: an element of SR (and the fundamental data type in maple) is just an expression tree. Mathematical meaning is only obtained by *what* you do with the element (it's like an untyped programming language). That means you can do algebra with objects in SR that represent a series in such a way that you get another object that can't possibly be interpreted as a series. maple (and apparently pynac by default) solves that problem by not even trying: they leave it to the user to do that, trusting that he/she will take care to only do that in situations where it is possible to interpret the result as a series.

[mmezzarobba]

Replying to nbruin:

The problem is: there are SR allows expressions of such generality that series expansions aren't always possible. Consider:

sage: function('f')
sage: s=series(sin(x))
sage: f(cos(x))+g

(I guess g should be s?)

Should that be considered a series in x?

Here Maple does (when one explicitly asks for a series expansion, of course):

#-->series(f(cos(x)) + sin(x), x);
series(f(1)+1*x+(-1/2*D(f)(1))*x^2-1/6*x^3+(1/24*D(f)(1)+1/8*`@@`(D,2)(f)(1))*x^4+1/120*x^5+O(x^6),x,6)

and I found that convenient on many occasions. But I seem to rememer that ginac can do something similar.

[rws]

Followup tickets: #18499, #18500