[GiNaC-devel] series() errors with divide by zero

[Richard B. Kreckel]

Hi!

On 06/26/2011 04:52 AM, J Sage wrote:
> There is an error in ginac when taking a series of a function at
> a singular point, for instance in ginsh,
>
>  > series(sqrt(x),x,3);
> power::eval(): division by zero
>
> from ginac-1.6 file power.cpp line 408.
>
> I instrumented power.cpp with some printf functions and found
> the offending dispatch to be in line 665:
>
>           coeff = deriv.subs(r, subs_options::no_pattern);
>
> while evaluating the third term in the expansion, presumably with
> a 1/x factor.  I'm new to ginac so its taking me a while to learn the
> code -- next step is to try and get it running in gdb -- so I'd appreciate
> any input from the list about how ginac is handling and/or is
> supposed to handle this kind of error.  Other software, like maxima,
> handles the problem fine:
>
> (%o6)                      maxima(sqrt(x), x, 0, 3)
> (%i7) taylor(sqrt(x),x,0,3);

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=y^q and compute the Laurent expansion in y. Setting 
q=2 in your case would give the desired result:

   series(sqrt(x),x,0,3)
     = series(sqrt(y^2),y,0,3*2)
     = y
     = 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.
-- 
Richard B. Kreckel
<http://www.ginac.de/~kreckel/>