From f16bf92e181b45c3caca15c6428bab0733accf1c Mon Sep 17 00:00:00 2001 From: Jens Vollinga Date: Wed, 10 Mar 2004 16:20:02 +0000 Subject: [PATCH] Improved series expansion synced to HEAD. --- ginac/inifcns_gamma.cpp | 4 +- ginac/inifcns_trans.cpp | 4 +- ginac/pseries.cpp | 89 +++++++++++++++++++++++++++++++++++------ 3 files changed, 80 insertions(+), 17 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 08758f56..c12df144 100644 --- a/ginac/inifcns_trans.cpp +++ b/ginac/inifcns_trans.cpp @@ -472,7 +472,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). @@ -882,7 +882,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 3528cbbd..607360e7 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" @@ -529,14 +529,23 @@ bool pseries::is_terminating() const ex basic::series(const relational & r, int order, unsigned options) const { epvector seq; + + // default for order-values that make no sense for Taylor expansion + if (order <= 0) { + seq.push_back(expair(Order(_ex1), order)); + return pseries(r, seq); + } + + // do Taylor expansion numeric fac = 1; ex deriv = *this; ex coeff = deriv.subs(r, subs_options::no_pattern); const symbol &s = ex_to(r.lhs()); - if (!coeff.is_zero()) + if (!coeff.is_zero()) { seq.push_back(expair(coeff, _ex0)); - + } + int n; for (n=1; n(r.lhs())); + const ex& sym = r.lhs(); + + // holds ldegrees of the series of individual factors + std::vector ldegrees; + + // find minimal degrees 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); + + ex buf = recombine_pair_to_ex(*it); + + int real_ldegree = buf.expand().ldegree(sym-r.rhs()); + if (real_ldegree == 0) { + int orderloop = 0; + do { + orderloop++; + real_ldegree = buf.series(r, orderloop, options).ldegree(sym); + } while (real_ldegree == orderloop); + } + + ldegrees.push_back(real_ldegree); + } + + int degsum = std::accumulate(ldegrees.begin(), ldegrees.end(), 0); + + // Multiply with remaining terms + std::vector::const_iterator itd = ldegrees.begin(); + for (epvector::const_iterator it=itbeg; it!=itend; ++it, ++itd) { + + // do series expansion with adjusted order + ex op = recombine_pair_to_ex(*it).series(r, order-degsum+(*itd), options); // Series multiplication if (it==itbeg) @@ -777,6 +814,7 @@ ex mul::series(const relational & r, int order, unsigned options) const else acc = ex_to(acc.mul_series(ex_to(op))); } + return acc.mul_const(ex_to(overall_coeff)); } @@ -823,6 +861,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); @@ -831,7 +872,6 @@ ex pseries::power_const(const numeric &p, int deg) const exvector co; co.reserve(deg); co.push_back(power(coeff(var, ldeg), p)); - bool all_sums_zero = true; for (int i=1; i(e).power_const(ex_to(exponent), order); + + int intexp = ex_to(exponent).to_int(); + const ex& sym = r.lhs(); + // find existing minimal degree + int real_ldegree = basis.expand().ldegree(sym-r.rhs()); + if (real_ldegree == 0) { + int orderloop = 0; + do { + orderloop++; + real_ldegree = basis.series(r, orderloop, options).ldegree(sym); + } while (real_ldegree == orderloop); + } + + ex e = basis.series(r, order + real_ldegree*(1-intexp), options); + + ex result; + try { + result = ex_to(e).power_const(intexp, order); + } + catch (pole_error) { + epvector ser; + ser.push_back(expair(Order(_ex1), order)); + result = pseries(r, ser); + } + + return result; } -- 2.44.0