series(sqrt(x),x,3);
Hi, There is an error in ginac when taking a series of a function at a singular point, for instance in ginsh, 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); Thanks in advance, John * *
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/>
...
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/>
series(x^n,x,1); (0^n)+Order(x) series(x^n,x,2);
series(sqrt(x^2),x,3);
Actually, ginac cannot do a series on x^n or on sqrt(x^2). Here's the ginsh output: power::eval(): division by zero power::eval(): division by zero
John
Hi! On 06/27/2011 07:09 AM, J Sage wrote:
... 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/ <http://www.ginac.de/%7Ekreckel/>>
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. Bye -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
participants (2)
-
J Sage
-
Richard B. Kreckel