From 2c9ea31bba5cd8a4079a69715445debe84305549 Mon Sep 17 00:00:00 2001 From: Jens Vollinga Date: Fri, 13 Feb 2004 21:23:04 +0000 Subject: [PATCH] Modification of series expansion --- ginac/inifcns_gamma.cpp | 4 +-- ginac/inifcns_trans.cpp | 4 +-- ginac/pseries.cpp | 55 ++++++++++++++++++++++++++++++++++++++--- 3 files changed, 56 insertions(+), 7 deletions(-) diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp index 545fcb84..6c80ed6d 100644 --- a/ginac/inifcns_gamma.cpp +++ b/ginac/inifcns_gamma.cpp @@ -193,7 +193,7 @@ static ex tgamma_series(const ex & arg, // 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); + // series(tgamma(x+m+1)/(x*(x+1)*...*(x+m)),x==-m,order); const ex arg_pt = arg.subs(rel, subs_options::no_pattern); if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive)) throw do_taylor(); // caught by function::series() @@ -202,7 +202,7 @@ static ex tgamma_series(const ex & arg, ex ser_denom = _ex1; for (numeric p; p<=m; ++p) ser_denom *= arg+p; - return (tgamma(arg+m+_ex1)/ser_denom).series(rel, order+1, options); + return (tgamma(arg+m+_ex1)/ser_denom).series(rel, order, options); } diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp index 73fd4f48..806ed3aa 100644 --- a/ginac/inifcns_trans.cpp +++ b/ginac/inifcns_trans.cpp @@ -469,7 +469,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, options); + return (sin(x)/cos(x)).series(rel, order, options); } REGISTER_FUNCTION(tan, eval_func(tan_eval). @@ -879,7 +879,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, options); + return (sinh(x)/cosh(x)).series(rel, order, options); } REGISTER_FUNCTION(tanh, eval_func(tanh_eval). diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index a0250478..dd4f51ab 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -21,7 +21,7 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ -#include +#include #include #include "pseries.h" @@ -747,19 +747,59 @@ ex mul::series(const relational & r, int order, unsigned options) const { pseries acc; // Series accumulator + GINAC_ASSERT(is_a(r.lhs())); + const ex& sym = r.lhs(); + + // holds ldegrees of the series of individual factors + std::vector ldegrees; + // flag if series have to be re-calculated + bool negldegree = false; + // Multiply with remaining terms const epvector::const_iterator itbeg = seq.begin(); const epvector::const_iterator itend = seq.end(); for (epvector::const_iterator it=itbeg; it!=itend; ++it) { ex op = recombine_pair_to_ex(*it).series(r, order, options); + int ldeg = op.ldegree(sym); + ldegrees.push_back(ldeg); + if (ldeg < 0) { + negldegree = true; + } + // Series multiplication if (it==itbeg) acc = ex_to(op); else acc = ex_to(acc.mul_series(ex_to(op))); } - return acc.mul_const(ex_to(overall_coeff)); + + if (negldegree && (seq.size() > 1)) { + + // re-calculation of series + + pseries newacc; + + const int degsum = std::accumulate(ldegrees.begin(), ldegrees.end(), 0); + + // Multiply with remaining terms + const epvector::const_reverse_iterator itbeg = seq.rbegin(); + const epvector::const_reverse_iterator itend = seq.rend(); + std::vector::const_reverse_iterator itd = ldegrees.rbegin(); + for (epvector::const_reverse_iterator it=itbeg; it!=itend; ++it, ++itd) { + + ex op = recombine_pair_to_ex(*it).series(r, order-degsum+(*itd), options); + + // Series multiplication + if (it==itbeg) + newacc = ex_to(op); + else + newacc = ex_to(newacc.mul_series(ex_to(op))); + } + return newacc.mul_const(ex_to(overall_coeff)); + } else { + return acc.mul_const(ex_to(overall_coeff)); + } } @@ -805,6 +845,9 @@ ex pseries::power_const(const numeric &p, int deg) const if (!(p*ldeg).is_integer()) throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series"); + // adjust number of coefficients + deg = deg - p.to_int()*ldeg; + // O(x^n)^(-m) is undefined if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative()) throw pole_error("pseries::power_const(): division by zero",1); @@ -842,6 +885,7 @@ ex pseries::power_const(const numeric &p, int deg) const } if (!higher_order && !all_sums_zero) new_seq.push_back(expair(Order(_ex1), p * ldeg + deg)); + return pseries(relational(var,point), new_seq); } @@ -895,8 +939,13 @@ ex power::series(const relational & r, int order, unsigned options) const } // No, expand basis into series + int intexp = ex_to(exponent).to_int(); ex e = basis.series(r, order, options); - return ex_to(e).power_const(ex_to(exponent), order); + int ldeg = ex_to(e).ldegree(r.lhs()); + if (intexp * ldeg < 0) { + e = basis.series(r, order + ldeg*(1-intexp), options); + } + return ex_to(e).power_const(intexp, order); } -- 2.50.0