From: Richard Kreckel Date: Mon, 3 Apr 2000 22:03:23 +0000 (+0000) Subject: - Some corrections in series expansion of several functions. X-Git-Tag: release_0-6-0~34 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=ca32e8cce5576051183fc6be17ae0c624a185a6d - Some corrections in series expansion of several functions. --- diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index 96f4c9a5..c03f4a15 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -65,6 +65,51 @@ static ex abs_eval(const ex & x) REGISTER_FUNCTION(abs, eval_func(abs_eval). evalf_func(abs_evalf)); + +////////// +// Complex sign +////////// + +static ex csgn_evalf(const ex & x) +{ + BEGIN_TYPECHECK + TYPECHECK(x,numeric) + END_TYPECHECK(csgn(x)) + + return csgn(ex_to_numeric(x)); +} + +static ex csgn_eval(const ex & x) +{ + if (is_ex_exactly_of_type(x, numeric)) + return csgn(ex_to_numeric(x)); + + if (is_ex_exactly_of_type(x, mul)) { + numeric oc = ex_to_numeric(x.op(x.nops()-1)); + if (oc.is_real()) { + if (oc > 0) + // csgn(42*x) -> csgn(x) + return csgn(x/oc).hold(); + else + // csgn(-42*x) -> -csgn(x) + return -csgn(x/oc).hold(); + } + if (oc.real().is_zero()) { + if (oc.imag() > 0) + // csgn(42*I*x) -> csgn(I*x) + return csgn(I*x/oc).hold(); + else + // csgn(-42*I*x) -> -csgn(I*x) + return -csgn(I*x/oc).hold(); + } + } + + return csgn(x).hold(); +} + +REGISTER_FUNCTION(csgn, eval_func(csgn_eval). + evalf_func(csgn_evalf)); + ////////// // dilogarithm ////////// diff --git a/ginac/inifcns.h b/ginac/inifcns.h index 75599ea9..9781d935 100644 --- a/ginac/inifcns.h +++ b/ginac/inifcns.h @@ -32,6 +32,9 @@ namespace GiNaC { /** Absolute value. */ DECLARE_FUNCTION_1P(abs) + +/** Complex sign. */ +DECLARE_FUNCTION_1P(csgn) /** Sine. */ DECLARE_FUNCTION_1P(sin) diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp index b28c2462..85275f62 100644 --- a/ginac/inifcns_gamma.cpp +++ b/ginac/inifcns_gamma.cpp @@ -84,11 +84,33 @@ static ex lgamma_deriv(const ex & x, unsigned deriv_param) return psi(x); } -// need to implement lgamma_series. + +static ex lgamma_series(const ex & x, const relational & rel, int order) +{ + // method: + // Taylor series where there is no pole falls back to psi function + // evaluation. + // On a pole at -m use the recurrence relation + // lgamma(x) == lgamma(x+1)-log(x) + // from which follows + // series(lgamma(x),x==-m,order) == + // series(lgamma(x+m+1)-log(x)...-log(x+m)),x==-m,order+1); + const ex x_pt = x.subs(rel); + if (!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive)) + throw do_taylor(); // caught by function::series() + // if we got here we have to care for a simple pole at -m: + numeric m = -ex_to_numeric(x_pt); + ex ser_sub = _ex0(); + for (numeric p; p<=m; ++p) + ser_sub += log(x+p); + return (lgamma(x+m+_ex1())-ser_sub).series(rel, order); +} + REGISTER_FUNCTION(lgamma, eval_func(lgamma_eval). evalf_func(lgamma_evalf). - derivative_func(lgamma_deriv)); + derivative_func(lgamma_deriv). + series_func(lgamma_series)); ////////// @@ -156,7 +178,7 @@ static ex tgamma_deriv(const ex & x, unsigned deriv_param) } -static ex tgamma_series(const ex & x, const relational & r, int order) +static ex tgamma_series(const ex & x, const relational & rel, int order) { // method: // Taylor series where there is no pole falls back to psi function @@ -164,9 +186,9 @@ static ex tgamma_series(const ex & x, const relational & r, int order) // On a pole at -m use the recurrence relation // tgamma(x) == tgamma(x+1) / x // from which follows - // series(tgamma(x),x,-m,order) == - // series(tgamma(x+m+1)/(x*(x+1)*...*(x+m)),x,-m,order+1); - const ex x_pt = x.subs(r); + // series(tgamma(x),x==-m,order) == + // series(tgamma(x+m+1)/(x*(x+1)*...*(x+m)),x==-m,order+1); + const ex x_pt = x.subs(rel); if (!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive)) throw do_taylor(); // caught by function::series() // if we got here we have to care for a simple pole at -m: @@ -174,7 +196,7 @@ static ex tgamma_series(const ex & x, const relational & r, int order) ex ser_denom = _ex1(); for (numeric p; p<=m; ++p) ser_denom *= x+p; - return (tgamma(x+m+_ex1())/ser_denom).series(r, order+1); + return (tgamma(x+m+_ex1())/ser_denom).series(rel, order+1); } @@ -251,38 +273,37 @@ static ex beta_deriv(const ex & x, const ex & y, unsigned deriv_param) } -static ex beta_series(const ex & x, const ex & y, const relational & r, int order) +static ex beta_series(const ex & x, const ex & y, const relational & rel, int order) { // method: // Taylor series where there is no pole of one of the tgamma functions // falls back to beta function evaluation. Otherwise, fall back to // tgamma series directly. - // FIXME: this could need some testing, maybe it's wrong in some cases? - const ex x_pt = x.subs(r); - const ex y_pt = y.subs(r); - GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol)); - const symbol *s = static_cast(r.lhs().bp); + const ex x_pt = x.subs(rel); + const ex y_pt = y.subs(rel); + GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol)); + const symbol *s = static_cast(rel.lhs().bp); ex x_ser, y_ser, xy_ser; if ((!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive)) && (!y_pt.info(info_flags::integer) || y_pt.info(info_flags::positive))) throw do_taylor(); // caught by function::series() // trap the case where x is on a pole directly: if (x.info(info_flags::integer) && !x.info(info_flags::positive)) - x_ser = tgamma(x+*s).series(r,order); + x_ser = tgamma(x+*s).series(rel,order); else - x_ser = tgamma(x).series(r,order); + x_ser = tgamma(x).series(rel,order); // trap the case where y is on a pole directly: if (y.info(info_flags::integer) && !y.info(info_flags::positive)) - y_ser = tgamma(y+*s).series(r,order); + y_ser = tgamma(y+*s).series(rel,order); else - y_ser = tgamma(y).series(r,order); + y_ser = tgamma(y).series(rel,order); // trap the case where y is on a pole directly: if ((x+y).info(info_flags::integer) && !(x+y).info(info_flags::positive)) - xy_ser = tgamma(y+x+*s).series(r,order); + xy_ser = tgamma(y+x+*s).series(rel,order); else - xy_ser = tgamma(y+x).series(r,order); + xy_ser = tgamma(y+x).series(rel,order); // compose the result: - return (x_ser*y_ser/xy_ser).series(r,order); + return (x_ser*y_ser/xy_ser).series(rel,order); } @@ -358,7 +379,7 @@ static ex psi1_deriv(const ex & x, unsigned deriv_param) return psi(_ex1(), x); } -static ex psi1_series(const ex & x, const relational & r, int order) +static ex psi1_series(const ex & x, const relational & rel, int order) { // method: // Taylor series where there is no pole falls back to polygamma function @@ -366,9 +387,9 @@ static ex psi1_series(const ex & x, const relational & r, int order) // On a pole at -m use the recurrence relation // psi(x) == psi(x+1) - 1/z // from which follows - // series(psi(x),x,-m,order) == - // series(psi(x+m+1) - 1/x - 1/(x+1) - 1/(x+m)),x,-m,order); - const ex x_pt = x.subs(r); + // series(psi(x),x==-m,order) == + // series(psi(x+m+1) - 1/x - 1/(x+1) - 1/(x+m)),x==-m,order); + const ex x_pt = x.subs(rel); if (!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive)) throw do_taylor(); // caught by function::series() // if we got here we have to care for a simple pole at -m: @@ -376,7 +397,7 @@ static ex psi1_series(const ex & x, const relational & r, int order) ex recur; for (numeric p; p<=m; ++p) recur += power(x+p,_ex_1()); - return (psi(x+m+_ex1())-recur).series(r, order); + return (psi(x+m+_ex1())-recur).series(rel, order); } const unsigned function_index_psi1 = @@ -478,7 +499,7 @@ static ex psi2_deriv(const ex & n, const ex & x, unsigned deriv_param) return psi(n+_ex1(), x); } -static ex psi2_series(const ex & n, const ex & x, const relational & r, int order) +static ex psi2_series(const ex & n, const ex & x, const relational & rel, int order) { // method: // Taylor series where there is no pole falls back to polygamma function @@ -486,10 +507,10 @@ static ex psi2_series(const ex & n, const ex & x, const relational & r, int orde // On a pole at -m use the recurrence relation // psi(n,x) == psi(n,x+1) - (-)^n * n! / x^(n+1) // from which follows - // series(psi(x),x,-m,order) == + // series(psi(x),x==-m,order) == // series(psi(x+m+1) - (-1)^n * n! * ((x)^(-n-1) + (x+1)^(-n-1) + ... - // ... + (x+m)^(-n-1))),x,-m,order); - const ex x_pt = x.subs(r); + // ... + (x+m)^(-n-1))),x==-m,order); + const ex x_pt = x.subs(rel); if (!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive)) throw do_taylor(); // caught by function::series() // if we got here we have to care for a pole of order n+1 at -m: @@ -498,7 +519,7 @@ static ex psi2_series(const ex & n, const ex & x, const relational & r, int orde for (numeric p; p<=m; ++p) recur += power(x+p,-n+_ex_1()); recur *= factorial(n)*power(_ex_1(),n); - return (psi(n, x+m+_ex1())-recur).series(r, order); + return (psi(n, x+m+_ex1())-recur).series(rel, order); } const unsigned function_index_psi2 = diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp index 12fabc0e..87059014 100644 --- a/ginac/inifcns_trans.cpp +++ b/ginac/inifcns_trans.cpp @@ -144,6 +144,29 @@ static ex log_deriv(const ex & x, unsigned deriv_param) return power(x, _ex_1()); } +/*static ex log_series(const ex &x, const relational &rel, int order) +{ + const ex x_pt = x.subs(rel); + if (!x_pt.info(info_flags::negative) && !x_pt.is_zero()) + throw do_taylor(); // caught by function::series() + // now we either have to care for the branch cut or the branch point: + if (x_pt.is_zero()) { // branch point: return a plain log(x). + epvector seq; + seq.push_back(expair(log(x), _ex0())); + return pseries(rel, seq); + } // the branch cut: + const ex point = rel.rhs(); + const symbol *s = static_cast(rel.lhs().bp); + const symbol foo; + // compute the formal series: + ex replx = series(log(x),*s==foo,order).subs(foo==point); + epvector seq; + // FIXME: this is probably off by 2 or so: + seq.push_back(expair(-I*csgn(x*I)*Pi,_ex0())); + seq.push_back(expair(Order(_ex1()),order)); + return series(replx + pseries(rel, seq),rel,order); +}*/ + static ex log_series(const ex &x, const relational &r, int order) { if (x.subs(r).is_zero()) { @@ -395,16 +418,16 @@ static ex tan_deriv(const ex & x, unsigned deriv_param) return (_ex1()+power(tan(x),_ex2())); } -static ex tan_series(const ex &x, const relational &r, int order) +static ex tan_series(const ex &x, const relational &rel, int order) { // method: // Taylor series where there is no pole falls back to tan_deriv. // On a pole simply expand sin(x)/cos(x). - const ex x_pt = x.subs(r); + const ex x_pt = x.subs(rel); 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(r, order+2); + return (sin(x)/cos(x)).series(rel, order+2); } REGISTER_FUNCTION(tan, eval_func(tan_eval). @@ -752,16 +775,16 @@ static ex tanh_deriv(const ex & x, unsigned deriv_param) return _ex1()-power(tanh(x),_ex2()); } -static ex tanh_series(const ex &x, const relational &r, int order) +static ex tanh_series(const ex &x, const relational &rel, int order) { // method: // Taylor series where there is no pole falls back to tanh_deriv. // On a pole simply expand sinh(x)/cosh(x). - const ex x_pt = x.subs(r); + const ex x_pt = x.subs(rel); 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(r, order+2); + return (sinh(x)/cosh(x)).series(rel, order+2); } REGISTER_FUNCTION(tanh, eval_func(tanh_eval).