Sage: Ticket #9555: Series expansions at singularities don't work
https://trac.sagemath.org/ticket/9555
<p>
Calling the series() method on a symbolic expression at a singularity (algebraic, logarithmic, essential) returns nonsense, inconsistent results, or raises an exception:
</p>
<p>
Examples:
</p>
<pre class="wiki">sage: sqrt(x).series(x,0)
Order(1)
sage: sqrt(x).series(x,1)
Order(x)
sage: sqrt(x).series(x,2)
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
/home/fredrik/sage/<ipython console> in <module>()
/home/fredrik/sage/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.series (sage/symbolic/expression.cpp:12811)()
RuntimeError: power::eval(): division by zero
</pre><pre class="wiki">sage: (log(x) + x).series(x,0)
Order(1)
sage: (log(x) + x).series(x,1)
(log(x)) + Order(x)
sage: (log(x) + x).series(x,2)
(log(x)) + 1*x
</pre><pre class="wiki">sage: exp(1/x).series(x,0)
Order(1)
sage: exp(1/x).series(x,1)
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
/home/fredrik/sage/<ipython console> in <module>()
/home/fredrik/sage/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.series (sage/symbolic/expression.cpp:12811)()
RuntimeError: power::eval(): division by zero
sage: exp(1/x).series(x,2)
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
/home/fredrik/sage/<ipython console> in <module>()
/home/fredrik/sage/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.series (sage/symbolic/expression.cpp:12811)()
RuntimeError: power::eval(): division by zero
</pre>en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/9555
Trac 1.1.6fredrik.johanssonTue, 20 Jul 2010 12:17:11 GMT
https://trac.sagemath.org/ticket/9555#comment:1
https://trac.sagemath.org/ticket/9555#comment:1
<p>
Another example:
</p>
<pre class="wiki">sage: (x^n).series(x,0)
Order(1)
sage: (x^n).series(x,1)
(0^n) + Order(x)
sage: (x^n).series(x,2)
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
/home/fredrik/sage/<ipython console> in <module>()
/home/fredrik/sage/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.series (sage/symbolic/expression.cpp:12811)()
RuntimeError: power::eval(): division by zero
</pre>
TicketjlhsageSun, 26 Jun 2011 20:57:57 GMT
https://trac.sagemath.org/ticket/9555#comment:2
https://trac.sagemath.org/ticket/9555#comment:2
<p>
This problem shows up in pure ginac-1.60, as you can check in their ginsh shell. I posted about this error on the ginac-devel list, and their reply is here:
</p>
<p>
<a class="ext-link" href="http://www.cebix.net/pipermail/ginac-devel/2011-June/001946.html"><span class="icon"></span>http://www.cebix.net/pipermail/ginac-devel/2011-June/001946.html</a>
</p>
<p>
Essentially, they don't claim to work with fractional power series. Here's the reply:
</p>
<blockquote class="citation">
<p>
GiNaC can only compute Taylor and Laurent series. Yours is a Puiseux
series: a series not in integer powers of x but in rational powers of x.
</p>
<p>
The Puiseux expansion of sqrt(x) is, well, x<sup>(1/2)</sup>.
</p>
<p>
You may try to set x=y<sup>q</sup> and compute the Laurent expansion in y. Setting
q=2 in your case would give the desired result:
</p>
<pre class="wiki">series(sqrt(x), x, 0, 3)
= series(sqrt(y^2), y, 0, 3*2)
= y
= x^(1/2).
</pre><p>
Note that the member functions degree() and ldegree() currently return
int, so this would have to be generalized somehow, when implementing
Puiseux series directly in GiNaC.
</p>
<p>
Bye
</p>
<blockquote>
<p>
-richy.
</p>
</blockquote>
</blockquote>
<p>
Note that Richy's suggestion about using <code>x=y^2</code> doesn't work either in ginsh.
Does anybody know how maxima resolves this issue?
</p>
TicketkcrismanTue, 28 Jun 2011 19:09:44 GMTcc set
https://trac.sagemath.org/ticket/9555#comment:3
https://trac.sagemath.org/ticket/9555#comment:3
<ul>
<li><strong>cc</strong>
<em>kcrisman</em> added
</li>
</ul>
<p>
Maxima:
</p>
<pre class="wiki">(%i1) taylor(sqrt(x),x,0,3);
(%o1)/T/ sqrt(x) + . . .
(%i3) taylor(%e^(-x)/x,x,0,5);
2 3 4 5
1 x x x x x
(%o3)/T/ - - 1 + - - -- + -- - --- + --- + . . .
x 2 6 24 120 720
</pre><p>
The first one in particular seems odd. Maybe that should be considered a bug? I've submitted <a class="ext-link" href="https://sourceforge.net/tracker/?func=detail&aid=3341693&group_id=4933&atid=104933"><span class="icon"></span>this Maxima bug</a>.
</p>
<hr />
<p>
Anyway, perhaps that's not relevant. There's a pretty relevant followup to the discussion referenced above:
</p>
<pre class="wiki">
> Actually, ginac cannot do a series on x^n or on sqrt(x^2). Here's the
> ginsh output:
>
> > series(x^n,x,1);
> (0^n)+Order(x)
> > series(x^n,x,2);
> power::eval(): division by zero
> > series(sqrt(x^2),x,3);
> power::eval(): division by zero
Oh, you're so right! Actually, that's because the result of
series(sqrt(x^2),x==0,3) is not necessarly x: it could just as well be
-x. The presence of the branch point makes it more tricky.
Sorry for not being able to help.
</pre><p>
That doesn't resolve things here. My guess is that the current behavior in Ginac could actually be considered correct, so that we should just add some documentation saying so?
</p>
TicketjlhsageWed, 29 Jun 2011 01:20:44 GMT
https://trac.sagemath.org/ticket/9555#comment:4
https://trac.sagemath.org/ticket/9555#comment:4
<p>
We should probably fix pynac/ginac to not give an error. From the referenced ginac-devel thread, the correct response to series(x<sup>{1/2}} at 0 is a Puiseux series and:
</sup></p>
<blockquote class="citation">
<p>
The Puiseux expansion of sqrt(x) is, well, x<sup>(1/2).
</sup></p>
</blockquote>
<p>
which ginac is not currently able to do, but which maxima in fact returns correctly as x<sup>(1/2).
Also, here's maxima on the series of x</sup>2 at 0:
</p>
<p>
(%i1) taylor(x<sup>2,x,0,3);
</sup></p>
<blockquote>
<p>
2
</p>
</blockquote>
<p>
(%o1)/T/ x + . . .
</p>
<p>
which again looks to be correct. So far there is not much response to the referenced thread on ginac-devel so we may have to fix this ourselves, and we might start by understanding how maxima recognizes Puiseux conditions. Towards that end I've started reading maxima's series.lisp code but, being a newbie to maxima development, pointers from the maxima experts would be very helpful.
</p>
<p>
John Hoebing
</p>
TicketkcrismanWed, 29 Jun 2011 01:59:33 GMT
https://trac.sagemath.org/ticket/9555#comment:5
https://trac.sagemath.org/ticket/9555#comment:5
<p>
Well, you can certainly implement Puiseux series. See <a class="closed ticket" href="https://trac.sagemath.org/ticket/4618" title="enhancement: Puiseux series (closed: fixed)">#4618</a>, where this is already requested.
</p>
<p>
But having a command called "taylor" return things that are not Taylor series is ... problematic. And I don't think that Maxima is recognizing this necessarily.
</p>
<pre class="wiki">>], ...)
`taylor (<expr>, <x>, <a>, <n>)' expands the expression <expr> in
a truncated Taylor or Laurent series in the variable <x> around
the point <a>, containing terms through `(<x> - <a>)^<n>'.
</pre><p>
nothing about Puiseux, nor do any examples.
</p>
<p>
Now, our <code>.series()</code> is not named quite that way, so I suppose it's possible we could have that. I still think it might be a little confusing, so I would at the very least recommend that implementing P. series be at <a class="closed ticket" href="https://trac.sagemath.org/ticket/4618" title="enhancement: Puiseux series (closed: fixed)">#4618</a> and then adding a keyword or something (even default, but maybe so that it could be disabled) here for the <code>.series()</code> method would be what we would do here.
</p>
TicketburcinWed, 29 Jun 2011 13:28:52 GMT
https://trac.sagemath.org/ticket/9555#comment:6
https://trac.sagemath.org/ticket/9555#comment:6
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9555#comment:5" title="Comment 5">kcrisman</a>:
</p>
<blockquote class="citation">
<p>
But having a command called "taylor" return things that are not Taylor series is ... problematic. And I don't think that Maxima is recognizing this necessarily.
</p>
</blockquote>
<pre class="wiki">>], ...)
`taylor (<expr>, <x>, <a>, <n>)' expands the expression <expr> in
a truncated Taylor or Laurent series in the variable <x> around
the point <a>, containing terms through `(<x> - <a>)^<n>'.
</pre><blockquote class="citation">
<p>
nothing about Puiseux, nor do any examples.
</p>
</blockquote>
<p>
+1
</p>
<blockquote class="citation">
<p>
Now, our <code>.series()</code> is not named quite that way, so I suppose it's possible we could have that.
</p>
</blockquote>
<p>
The <code>taylor()</code> method is cruft left over from pre-pynac symbolics. We should depreciate it in favor of the <code>series()</code> method. It's perfectly acceptable to give Puiseux series as a result of a call to <code>.series()</code>. I expect this to be done in <a class="new ticket" href="https://trac.sagemath.org/ticket/6119" title="enhancement: deprecate taylor() in favor of series() (new)">#6119</a>, where we add an <code>algorithm=</code> option to <code>.series()</code>. The default behavior can be to call pynac and fall back to maxima if that returns an error.
</p>
<p>
This ticket can then be about fixing pynac/ginac to allow these series expansions, in which case, it should be moved to the issue tracker on bitbucket:
</p>
<p>
<a class="ext-link" href="https://bitbucket.org/burcin/pynac/issues"><span class="icon"></span>https://bitbucket.org/burcin/pynac/issues</a>
</p>
TicketkcrismanWed, 29 Jun 2011 13:38:31 GMT
https://trac.sagemath.org/ticket/9555#comment:7
https://trac.sagemath.org/ticket/9555#comment:7
<blockquote class="citation">
<blockquote class="citation">
<p>
Now, our <code>.series()</code> is not named quite that way, so I suppose it's possible we could have that.
</p>
</blockquote>
<p>
The <code>taylor()</code> method is cruft left over from pre-pynac symbolics. We should depreciate it in favor of the <code>series()</code> method.
</p>
</blockquote>
<p>
I don't know if I'd go that far. I think that especially the top-level <code>taylor()</code> (which calls <code>SR.taylor()</code>) is very useful; having the top command be <code>series()</code> is more confusing, as there are lots of series out there.
</p>
<p>
That said, <code>SR.taylor()</code> should probably be calling Pynac at this point, probably just the <code>.series()</code> method (without Puiseaux). I agree that the code itself is cruft, just not the name :)
</p>
TicketburcinWed, 29 Jun 2011 13:48:21 GMT
https://trac.sagemath.org/ticket/9555#comment:8
https://trac.sagemath.org/ticket/9555#comment:8
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9555#comment:7" title="Comment 7">kcrisman</a>:
</p>
<blockquote class="citation">
<blockquote class="citation">
<p>
The <code>taylor()</code> method is cruft left over from pre-pynac symbolics. We should depreciate it in favor of the <code>series()</code> method.
</p>
</blockquote>
<p>
I don't know if I'd go that far. I think that especially the top-level <code>taylor()</code> (which calls <code>SR.taylor()</code>) is very useful; having the top command be <code>series()</code> is more confusing, as there are lots of series out there.
</p>
</blockquote>
<p>
OK. I can see the educational appeal of a <code>taylor()</code> method.
</p>
<blockquote class="citation">
<p>
That said, <code>SR.taylor()</code> should probably be calling Pynac at this point, probably just the <code>.series()</code> method (without Puiseaux). I agree that the code itself is cruft, just not the name :)
</p>
</blockquote>
<p>
This would work, given the current behavior of Pynac/GiNaC's <code>series()</code> function. It would cause problems in the future if that ever changes. We should also add some doctests to make sure it only returns Taylor expansions.
</p>
<p>
This can also be done as a part of <a class="new ticket" href="https://trac.sagemath.org/ticket/6119" title="enhancement: deprecate taylor() in favor of series() (new)">#6119</a>. Then, instead of a deprecation message we give an error pointing to the <code>.series()</code> method.
</p>
TicketkcrismanWed, 29 Jun 2011 13:50:07 GMT
https://trac.sagemath.org/ticket/9555#comment:9
https://trac.sagemath.org/ticket/9555#comment:9
<blockquote class="citation">
<blockquote class="citation">
<p>
That said, <code>SR.taylor()</code> should probably be calling Pynac at this point, probably just the <code>.series()</code> method (without Puiseaux). I agree that the code itself is cruft, just not the name :)
</p>
</blockquote>
<p>
This would work, given the current behavior of Pynac/GiNaC's <code>series()</code> function. It would cause problems in the future if that ever changes. We should also add some doctests to make sure it only returns Taylor expansions.
</p>
</blockquote>
<p>
Well, I would envision it would have a special keyword that would require returning only Taylor series. Presumably <code>.series()</code> should also have this, don't you think? Even if that weren't the default behavior.
</p>
<blockquote class="citation">
<p>
This can also be done as a part of <a class="new ticket" href="https://trac.sagemath.org/ticket/6119" title="enhancement: deprecate taylor() in favor of series() (new)">#6119</a>. Then, instead of a deprecation message we give an error pointing to the <code>.series()</code> method.
</p>
</blockquote>
<p>
Sounds good.
</p>
TicketkcrismanSat, 06 Dec 2014 20:27:29 GMTcc, upstream changed
https://trac.sagemath.org/ticket/9555#comment:10
https://trac.sagemath.org/ticket/9555#comment:10
<ul>
<li><strong>cc</strong>
<em>rws</em> added
</li>
<li><strong>upstream</strong>
changed from <em>N/A</em> to <em>Reported upstream. Developers acknowledge bug.</em>
</li>
</ul>
<p>
I completely forgot about this ticket when I reported <a class="ext-link" href="http://sourceforge.net/p/maxima/bugs/2850/"><span class="icon"></span>http://sourceforge.net/p/maxima/bugs/2850/</a> Maxima (and we) need better documentation about this stuff. That said, I may have changed my mind about this. Maybe we can make <code>taylor</code> clearly only Taylor or something... Incidentally, someone said that another name for these series at singular points are Frobenius series.
</p>
TicketjakobkroekerWed, 06 May 2015 13:34:02 GMTcc changed
https://trac.sagemath.org/ticket/9555#comment:11
https://trac.sagemath.org/ticket/9555#comment:11
<ul>
<li><strong>cc</strong>
<em>jakobkroeker</em> added
</li>
</ul>
Ticket