Opened 11 years ago

Last modified 15 months ago

# Series expansions at singularities don't work

Reported by: Owned by: fredrik.johansson burcin major symbolics kcrisman, rws, jakobkroeker Reported upstream. Developers acknowledge bug.

### Description

Calling the series() method on a symbolic expression at a singularity (algebraic, logarithmic, essential) returns nonsense, inconsistent results, or raises an exception:

Examples:

```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
```
```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
```
```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
```

### comment:1 Changed 11 years ago by fredrik.johansson

Another example:

```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
```

### comment:2 Changed 10 years ago by jlhsage

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:

Essentially, they don't claim to work with fractional power series. Here's the reply:

GiNaC can only compute Taylor and Laurent series. Your's is a Puiseux series: a series not in integer powers of x but in rational powers of x.

The Puiseux expansion of sqrt(x) is, well, x(1/2).

You may try to set x=yq and compute the Laurent expansion in y. Setting q=2 in your case would give the desired result:

series(sqrt(x),x,0,3)

# x(1/2).

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.

Bye

-richy.

Note that Richy's suggestion about using 'x=y2' doesn't work either in ginsh. Does anybody know how maxima resolves this issue?

Version 0, edited 10 years ago by jlhsage (next)

### comment:3 Changed 10 years ago by kcrisman

• Cc kcrisman added

Maxima:

```(%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
```

The first one in particular seems odd. Maybe that should be considered a bug? I've submitted this Maxima bug.

Anyway, perhaps that's not relevant. There's a pretty relevant followup to the discussion referenced above:

```> 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.
```

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?

### comment:4 Changed 10 years ago by jlhsage

We should probably fix pynac/ginac to not give an error. From the referenced ginac-devel thread, the correct response to series(x{1/2}} at 0 is a Puiseux series and:

The Puiseux expansion of sqrt(x) is, well, x(1/2).

which ginac is not currently able to do, but which maxima in fact returns correctly as x(1/2). Also, here's maxima on the series of x2 at 0:

(%i1) taylor(x2,x,0,3);

2

(%o1)/T/ x + . . .

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.

John Hoebing

### comment:5 follow-up: ↓ 6 Changed 10 years ago by kcrisman

Well, you can certainly implement Puiseux series. See #4618, where this is already requested.

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.

```>], ...)
`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>'.
```

nothing about Puiseux, nor do any examples.

Now, our `.series()` 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 #4618 and then adding a keyword or something (even default, but maybe so that it could be disabled) here for the `.series()` method would be what we would do here.

### comment:6 in reply to: ↑ 5 ; follow-up: ↓ 7 Changed 10 years ago by burcin

Replying to kcrisman:

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.

```>], ...)
`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>'.
```

nothing about Puiseux, nor do any examples.

+1

Now, our `.series()` is not named quite that way, so I suppose it's possible we could have that.

The `taylor()` method is cruft left over from pre-pynac symbolics. We should depreciate it in favor of the `series()` method. It's perfectly acceptable to give Puiseux series as a result of a call to `.series()`. I expect this to be done in #6119, where we add an `algorithm=` option to `.series()`. The default behavior can be to call pynac and fall back to maxima if that returns an error.

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:

### comment:7 in reply to: ↑ 6 ; follow-up: ↓ 8 Changed 10 years ago by kcrisman

Now, our `.series()` is not named quite that way, so I suppose it's possible we could have that.

The `taylor()` method is cruft left over from pre-pynac symbolics. We should depreciate it in favor of the `series()` method.

I don't know if I'd go that far. I think that especially the top-level `taylor()` (which calls `SR.taylor()`) is very useful; having the top command be `series()` is more confusing, as there are lots of series out there.

That said, `SR.taylor()` should probably be calling Pynac at this point, probably just the `.series()` method (without Puiseaux). I agree that the code itself is cruft, just not the name :)

### comment:8 in reply to: ↑ 7 ; follow-up: ↓ 9 Changed 10 years ago by burcin

Replying to kcrisman:

The `taylor()` method is cruft left over from pre-pynac symbolics. We should depreciate it in favor of the `series()` method.

I don't know if I'd go that far. I think that especially the top-level `taylor()` (which calls `SR.taylor()`) is very useful; having the top command be `series()` is more confusing, as there are lots of series out there.

OK. I can see the educational appeal of a `taylor()` method.

That said, `SR.taylor()` should probably be calling Pynac at this point, probably just the `.series()` method (without Puiseaux). I agree that the code itself is cruft, just not the name :)

This would work, given the current behavior of Pynac/GiNaC's `series()` 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.

This can also be done as a part of #6119. Then, instead of a deprecation message we give an error pointing to the `.series()` method.

### comment:9 in reply to: ↑ 8 Changed 10 years ago by kcrisman

That said, `SR.taylor()` should probably be calling Pynac at this point, probably just the `.series()` method (without Puiseaux). I agree that the code itself is cruft, just not the name :)

This would work, given the current behavior of Pynac/GiNaC's `series()` 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.

Well, I would envision it would have a special keyword that would require returning only Taylor series. Presumably `.series()` should also have this, don't you think? Even if that weren't the default behavior.

This can also be done as a part of #6119. Then, instead of a deprecation message we give an error pointing to the `.series()` method.

Sounds good.

### comment:10 Changed 7 years ago by kcrisman

• Cc rws added
• Report Upstream changed from N/A to Reported upstream. Developers acknowledge bug.

I completely forgot about this ticket when I reported http://sourceforge.net/p/maxima/bugs/2850/ Maxima (and we) need better documentation about this stuff. That said, I may have changed my mind about this. Maybe we can make `taylor` clearly only Taylor or something... Incidentally, someone said that another name for these series at singular points are Frobenius series.

### comment:11 Changed 6 years ago by jakobkroeker

• Cc jakobkroeker added
Note: See TracTickets for help on using tickets.