From: Richard Kreckel Date: Wed, 7 Feb 2001 19:16:21 +0000 (+0000) Subject: trancendental function branch-point sanity X-Git-Tag: release_0-7-1~2 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=3d9039b1a02aedda5eff2392e9c76df986e94f52;ds=sidebyside trancendental function branch-point sanity --- diff --git a/check/exam_pseries.cpp b/check/exam_pseries.cpp index 54f4a806..cb65a533 100644 --- a/check/exam_pseries.cpp +++ b/check/exam_pseries.cpp @@ -220,6 +220,73 @@ static unsigned exam_series9(void) return check_series(e,2,d,5); } +// Series expansion of logarithms around branch points +static unsigned exam_series10(void) +{ + unsigned result = 0; + ex e, d; + symbol a("a"); + + e = log(x); + d = log(x); + result += check_series(e,0,d,5); + + e = log(3/x); + d = log(3)-log(x); + result += check_series(e,0,d,5); + + e = log(3*pow(x,2)); + d = log(3)+2*log(x); + result += check_series(e,0,d,5); + + // These ones must not be expanded because it would result in a branch cut + // running in the wrong direction. (Other systems tend to get this wrong.) + e = log(-x); + d = e; + result += check_series(e,0,d,5); + + e = log(I*(x-123)); + d = e; + result += check_series(e,123,d,5); + + e = log(a*x); + d = e; // we don't know anything about a! + result += check_series(e,0,d,5); + + e = log((1-x)/x); + d = log(1-x) - (x-1) + pow(x-1,2)/2 - pow(x-1,3)/3 + Order(pow(x-1,4)); + result += check_series(e,1,d,4); + + return result; +} + +// Series expansion of other functions around branch points +static unsigned exam_series11(void) +{ + unsigned result = 0; + ex e, d; + + // NB: Mma and Maple give different results, but they agree if one + // takes into account that by assumption |x|<1. + e = atan(x); + d = (I*log(2)/2-I*log(1+I*x)/2) + (x-I)/4 + I*pow(x-I,2)/16 + Order(pow(x-I,3)); + result += check_series(e,I,d,3); + + // NB: here, at -I, Mathematica disagrees, but it is wrong -- they + // pick up a complex phase by incorrectly expanding logarithms. + e = atan(x); + d = (-I*log(2)/2+I*log(1-I*x)/2) + (x+I)/4 - I*pow(x+I,2)/16 + Order(pow(x+I,3)); + result += check_series(e,-I,d,3); + + // This is basically the same as above, the branch point is at +/-1: + e = atanh(x); + d = (-log(2)/2+log(x+1)/2) + (x+1)/4 + pow(x+1,2)/16 + Order(pow(x+1,3)); + result += check_series(e,-1,d,3); + + return result; +} + + unsigned exam_pseries(void) { unsigned result = 0; @@ -236,6 +303,8 @@ unsigned exam_pseries(void) result += exam_series7(); cout << '.' << flush; result += exam_series8(); cout << '.' << flush; result += exam_series9(); cout << '.' << flush; + result += exam_series10(); cout << '.' << flush; + result += exam_series11(); cout << '.' << flush; if (!result) { cout << " passed " << endl; diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index f410e7aa..d5c90e1d 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -114,7 +114,8 @@ static ex csgn_series(const ex & arg, { const ex arg_pt = arg.subs(rel); if (arg_pt.info(info_flags::numeric) - && ex_to_numeric(arg_pt).real().is_zero()) + && ex_to_numeric(arg_pt).real().is_zero() + && !(options & series_options::suppress_branchcut)) throw (std::domain_error("csgn_series(): on imaginary axis")); epvector seq; diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp index c7dbfbe6..cceaad99 100644 --- a/ginac/inifcns_trans.cpp +++ b/ginac/inifcns_trans.cpp @@ -158,7 +158,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 +176,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)); @@ -608,7 +616,7 @@ static ex atan_deriv(const ex & x, unsigned deriv_param) return power(_ex1()+power(x,_ex2()), _ex_1()); } -static ex atan_series(const ex &x, +static ex atan_series(const ex &arg, const relational &rel, int order, unsigned options) @@ -620,16 +628,35 @@ static ex atan_series(const ex &x, // 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)/(1-I*x))/(2*I) + // (log(1+I*x)-log(1-I*x))/(2*I) // instead. - // (The constant term on the cut itself could be made simpler.) - const ex x_pt = x.subs(rel); - if (!(I*x_pt).info(info_flags::real)) + const ex arg_pt = arg.subs(rel); + if (!(I*arg_pt).info(info_flags::real)) throw do_taylor(); // Re(x) != 0 - if ((I*x_pt).info(info_flags::real) && abs(I*x_pt)<_ex1()) + if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1()) throw do_taylor(); // Re(x) == 0, but abs(x)<1 - // if we got here we have to care for cuts and poles - return (log((1+I*x)/(1-I*x))/(2*I)).series(rel, order, options); + // 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). @@ -979,27 +1006,47 @@ static ex atanh_deriv(const ex & x, unsigned deriv_param) return power(_ex1()-power(x,_ex2()),_ex_1()); } -static ex atanh_series(const ex &x, +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 atan_deriv. + // 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)/(1-x))/(2*I) + // (log(1+x)-log(1-x))/2 // instead. - // (The constant term on the cut itself could be made simpler.) - const ex x_pt = x.subs(rel); - if (!(x_pt).info(info_flags::real)) + const ex arg_pt = arg.subs(rel); + if (!(arg_pt).info(info_flags::real)) throw do_taylor(); // Im(x) != 0 - if ((x_pt).info(info_flags::real) && abs(x_pt)<_ex1()) + if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1()) throw do_taylor(); // Im(x) == 0, but abs(x)<1 - // if we got here we have to care for cuts and poles - return (log((1+x)/(1-x))/2).series(rel, order, options); + // 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).