X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns_trans.cpp;h=dbbb122c87e8cbe0c5279a7a21a64782f4fe3765;hp=390f84f1a0472d25d82116f95a48d46e92e12921;hb=d9927cf4bc27d542d5f6487d31f5d43dc2e39af7;hpb=3ad4e97e949298e5d1733843130e399f4e171d9d diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp index 390f84f1..dbbb122c 100644 --- a/ginac/inifcns_trans.cpp +++ b/ginac/inifcns_trans.cpp @@ -34,9 +34,7 @@ #include "pseries.h" #include "utils.h" -#ifndef NO_NAMESPACE_GINAC namespace GiNaC { -#endif // ndef NO_NAMESPACE_GINAC ////////// // exponential function @@ -158,7 +156,7 @@ static ex log_series(const ex &arg, } catch (pole_error) { must_expand_arg = true; } - // or we are at the branch cut anyways + // or we are at the branch point anyways if (arg_pt.is_zero()) must_expand_arg = true; @@ -176,23 +174,31 @@ static ex log_series(const ex &arg, const ex point = rel.rhs(); const int n = argser.ldegree(*s); epvector seq; - seq.push_back(expair(n*log(*s-point), _ex0())); + // construct what we carelessly called the n*log(x) term above + ex coeff = argser.coeff(*s, n); + // expand the log, but only if coeff is real and > 0, since otherwise + // it would make the branch cut run into the wrong direction + if (coeff.info(info_flags::positive)) + seq.push_back(expair(n*log(*s-point)+log(coeff), _ex0())); + else + seq.push_back(expair(log(coeff*pow(*s-point, n)), _ex0())); if (!argser.is_terminating() || argser.nops()!=1) { // in this case n more terms are needed - ex newarg = ex_to_pseries(arg.series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true); + // (sadly, to generate them, we have to start from the beginning) + ex newarg = ex_to_pseries((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true); return pseries(rel, seq).add_series(ex_to_pseries(log(newarg).series(rel, order, options))); } else // it was a monomial return pseries(rel, seq); } if (!(options & series_options::suppress_branchcut) && - arg_pt.info(info_flags::negative)) { + arg_pt.info(info_flags::negative)) { // method: // This is the branch cut: assemble the primitive series manually and // then add the corresponding complex step function. const symbol *s = static_cast(rel.lhs().bp); const ex point = rel.rhs(); const symbol foo; - ex replarg = series(log(arg), *s==foo, order, false).subs(foo==point); + ex replarg = series(log(arg), *s==foo, order).subs(foo==point); epvector seq; seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0())); seq.push_back(expair(Order(_ex1()), order)); @@ -447,6 +453,7 @@ static ex tan_series(const ex &x, int order, unsigned options) { + GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol)); // method: // Taylor series where there is no pole falls back to tan_deriv. // On a pole simply expand sin(x)/cos(x). @@ -454,7 +461,7 @@ static ex tan_series(const ex &x, if (!(2*x_pt/Pi).info(info_flags::odd)) throw do_taylor(); // caught by function::series() // if we got here we have to care for a simple pole - return (sin(x)/cos(x)).series(rel, order+2); + return (sin(x)/cos(x)).series(rel, order+2, options); } REGISTER_FUNCTION(tan, eval_func(tan_eval). @@ -607,9 +614,53 @@ static ex atan_deriv(const ex & x, unsigned deriv_param) return power(_ex1()+power(x,_ex2()), _ex_1()); } +static ex atan_series(const ex &arg, + const relational &rel, + int order, + unsigned options) +{ + GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol)); + // method: + // Taylor series where there is no pole or cut falls back to atan_deriv. + // There are two branch cuts, one runnig from I up the imaginary axis and + // one running from -I down the imaginary axis. The points I and -I are + // poles. + // On the branch cuts and the poles series expand + // (log(1+I*x)-log(1-I*x))/(2*I) + // instead. + const ex arg_pt = arg.subs(rel); + if (!(I*arg_pt).info(info_flags::real)) + throw do_taylor(); // Re(x) != 0 + if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1()) + throw do_taylor(); // Re(x) == 0, but abs(x)<1 + // care for the poles, using the defining formula for atan()... + if (arg_pt.is_equal(I) || arg_pt.is_equal(-I)) + return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options); + if (!(options & series_options::suppress_branchcut)) { + // method: + // This is the branch cut: assemble the primitive series manually and + // then add the corresponding complex step function. + const symbol *s = static_cast(rel.lhs().bp); + const ex point = rel.rhs(); + const symbol foo; + ex replarg = series(atan(arg), *s==foo, order).subs(foo==point); + ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2(); + if ((I*arg_pt)<_ex0()) + Order0correction += log((I*arg_pt+_ex_1())/(I*arg_pt+_ex1()))*I*_ex_1_2(); + else + Order0correction += log((I*arg_pt+_ex1())/(I*arg_pt+_ex_1()))*I*_ex1_2(); + epvector seq; + seq.push_back(expair(Order0correction, _ex0())); + seq.push_back(expair(Order(_ex1()), order)); + return series(replarg - pseries(rel, seq), rel, order); + } + throw do_taylor(); +} + REGISTER_FUNCTION(atan, eval_func(atan_eval). evalf_func(atan_evalf). - derivative_func(atan_deriv)); + derivative_func(atan_deriv). + series_func(atan_series)); ////////// // inverse tangent (atan2(y,x)) @@ -815,6 +866,7 @@ static ex tanh_series(const ex &x, int order, unsigned options) { + GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol)); // method: // Taylor series where there is no pole falls back to tanh_deriv. // On a pole simply expand sinh(x)/cosh(x). @@ -822,7 +874,7 @@ static ex tanh_series(const ex &x, if (!(2*I*x_pt/Pi).info(info_flags::odd)) throw do_taylor(); // caught by function::series() // if we got here we have to care for a simple pole - return (sinh(x)/cosh(x)).series(rel, order+2); + return (sinh(x)/cosh(x)).series(rel, order+2, options); } REGISTER_FUNCTION(tanh, eval_func(tanh_eval). @@ -952,10 +1004,53 @@ static ex atanh_deriv(const ex & x, unsigned deriv_param) return power(_ex1()-power(x,_ex2()),_ex_1()); } +static ex atanh_series(const ex &arg, + const relational &rel, + int order, + unsigned options) +{ + GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol)); + // method: + // Taylor series where there is no pole or cut falls back to atanh_deriv. + // There are two branch cuts, one runnig from 1 up the real axis and one + // one running from -1 down the real axis. The points 1 and -1 are poles + // On the branch cuts and the poles series expand + // (log(1+x)-log(1-x))/2 + // instead. + const ex arg_pt = arg.subs(rel); + if (!(arg_pt).info(info_flags::real)) + throw do_taylor(); // Im(x) != 0 + if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1()) + throw do_taylor(); // Im(x) == 0, but abs(x)<1 + // care for the poles, using the defining formula for atanh()... + if (arg_pt.is_equal(_ex1()) || arg_pt.is_equal(_ex_1())) + return ((log(_ex1()+arg)-log(_ex1()-arg))*_ex1_2()).series(rel, order, options); + // ...and the branch cuts (the discontinuity at the cut being just I*Pi) + if (!(options & series_options::suppress_branchcut)) { + // method: + // This is the branch cut: assemble the primitive series manually and + // then add the corresponding complex step function. + const symbol *s = static_cast(rel.lhs().bp); + const ex point = rel.rhs(); + const symbol foo; + ex replarg = series(atanh(arg), *s==foo, order).subs(foo==point); + ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2(); + if (arg_pt<_ex0()) + Order0correction += log((arg_pt+_ex_1())/(arg_pt+_ex1()))*_ex1_2(); + else + Order0correction += log((arg_pt+_ex1())/(arg_pt+_ex_1()))*_ex_1_2(); + epvector seq; + seq.push_back(expair(Order0correction, _ex0())); + seq.push_back(expair(Order(_ex1()), order)); + return series(replarg - pseries(rel, seq), rel, order); + } + throw do_taylor(); +} + REGISTER_FUNCTION(atanh, eval_func(atanh_eval). evalf_func(atanh_evalf). - derivative_func(atanh_deriv)); + derivative_func(atanh_deriv). + series_func(atanh_series)); + -#ifndef NO_NAMESPACE_GINAC } // namespace GiNaC -#endif // ndef NO_NAMESPACE_GINAC