From: Christian Bauer Date: Fri, 29 Oct 2004 15:11:39 +0000 (+0000) Subject: improvements to pseries, esp. wrt series expansion of integrals [Chris Dams] X-Git-Tag: release_1-4-0~206 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=dc0683db2b9fcb2b3533942cbfd497f17bda0122 improvements to pseries, esp. wrt series expansion of integrals [Chris Dams] --- diff --git a/ginac/integral.h b/ginac/integral.h index ef2e6626..22f3ad1c 100644 --- a/ginac/integral.h +++ b/ginac/integral.h @@ -56,6 +56,7 @@ public: ex eval_integ() const; protected: ex derivative(const symbol & s) const; + ex series(const relational & r, int order, unsigned options = 0) const; // new virtual functions which can be overridden by derived classes // none diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index 4ec5bcef..bc04b186 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -33,6 +33,7 @@ #include "relational.h" #include "operators.h" #include "symbol.h" +#include "integral.h" #include "archive.h" #include "utils.h" @@ -266,6 +267,8 @@ ex pseries::op(size_t i) const if (i >= seq.size()) throw (std::out_of_range("op() out of range")); + if (is_order_function(seq[i].rest)) + return Order(power(var-point, seq[i].coeff)); return seq[i].rest * power(var - point, seq[i].coeff); } @@ -424,6 +427,31 @@ ex pseries::conjugate() const return result; } +ex pseries::eval_integ() const +{ + epvector *newseq = NULL; + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) { + if (newseq) { + newseq->push_back(expair(i->rest.eval_integ(), i->coeff)); + continue; + } + ex newterm = i->rest.eval_integ(); + if (!are_ex_trivially_equal(newterm, i->rest)) { + newseq = new epvector; + newseq->reserve(seq.size()); + for (epvector::const_iterator j=seq.begin(); j!=i; ++j) + newseq->push_back(*j); + newseq->push_back(expair(newterm, i->coeff)); + } + } + + ex newpoint = point.eval_integ(); + if (newseq || !are_ex_trivially_equal(newpoint, point)) + return (new pseries(var==newpoint, *newseq)) + ->setflag(status_flags::dynallocated); + return *this; +} + ex pseries::subs(const exmap & m, unsigned options) const { // If expansion variable is being substituted, convert the series to a @@ -519,6 +547,20 @@ bool pseries::is_terminating() const return seq.empty() || !is_order_function((seq.end()-1)->rest); } +ex pseries::coeffop(size_t i) const +{ + if (i >=nops()) + throw (std::out_of_range("coeffop() out of range")); + return seq[i].rest; +} + +ex pseries::exponop(size_t i) const +{ + if (i >= nops()) + throw (std::out_of_range("exponop() out of range")); + return seq[i].coeff; +} + /* * Implementations of series expansion @@ -729,6 +771,11 @@ ex pseries::mul_series(const pseries &other) const nul.push_back(expair(Order(_ex1), _ex0)); return pseries(relational(var,point), nul); } + + if (seq.empty() || other.seq.empty()) { + return (new pseries(var==point, epvector())) + ->setflag(status_flags::dynallocated); + } // Series multiplication epvector new_seq; @@ -880,7 +927,14 @@ ex pseries::power_const(const numeric &p, int deg) const throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series"); // adjust number of coefficients - deg = deg - (p*ldeg).to_int(); + int numcoeff = deg - (p*ldeg).to_int(); + if (numcoeff <= 0) { + epvector epv; + epv.reserve(1); + epv.push_back(expair(Order(_ex1), deg)); + return (new pseries(relational(var,point), epv)) + ->setflag(status_flags::dynallocated); + } // O(x^n)^(-m) is undefined if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative()) @@ -888,9 +942,9 @@ ex pseries::power_const(const numeric &p, int deg) const // Compute coefficients of the powered series exvector co; - co.reserve(deg); + co.reserve(numcoeff); co.push_back(power(coeff(var, ldeg), p)); - for (int i=1; i(fseries).coeffop(i); + fexpansion.push_back(expair( + currcoeff == Order(_ex1) + ? currcoeff + : integral(x, a.subs(r), b.subs(r), currcoeff), + ex_to(fseries).exponop(i) + )); + } + + // Expanding lower boundary + ex result = (new pseries(r, fexpansion))->setflag(status_flags::dynallocated); + ex aseries = (a-a.subs(r)).series(r, order, options); + fseries = f.series(x == (a.subs(r)), order, options); + for (size_t i=0; i(fseries).coeffop(i); + if (is_order_function(currcoeff)) + break; + ex currexpon = ex_to(fseries).exponop(i); + int orderforf = order-ex_to(currexpon).to_int()-1; + currcoeff=currcoeff.series(r, orderforf); + ex term = ex_to(aseries).power_const(ex_to(currexpon+1),order); + term = ex_to(term).mul_const(ex_to(-1/(currexpon+1))); + term = ex_to(term).mul_series(ex_to(currcoeff)); + result = ex_to(result).add_series(ex_to(term)); + } + + // Expanding upper boundary + ex bseries = (b-b.subs(r)).series(r, order, options); + fseries = f.series(x == (b.subs(r)), order, options); + for (size_t i=0; i(fseries).coeffop(i); + if (is_order_function(currcoeff)) + break; + ex currexpon = ex_to(fseries).exponop(i); + int orderforf = order-ex_to(currexpon).to_int()-1; + currcoeff=currcoeff.series(r, orderforf); + ex term = ex_to(bseries).power_const(ex_to(currexpon+1),order); + term = ex_to(term).mul_const(ex_to(1/(currexpon+1))); + term = ex_to(term).mul_series(ex_to(currcoeff)); + result = ex_to(result).add_series(ex_to(term)); + } + + return result; +} + /** Compute the truncated series expansion of an expression. * This function returns an expression containing an object of class pseries diff --git a/ginac/pseries.h b/ginac/pseries.h index fd2a0428..7222535c 100644 --- a/ginac/pseries.h +++ b/ginac/pseries.h @@ -56,6 +56,7 @@ public: ex normal(exmap & repl, exmap & rev_lookup, int level = 0) const; ex expand(unsigned options = 0) const; ex conjugate() const; + ex eval_integ() const; protected: ex derivative(const symbol & s) const; @@ -83,6 +84,10 @@ public: * false otherwise. */ bool is_terminating() const; + /** Get coefficients and exponents. */ + ex coeffop(size_t i) const; + ex exponop(size_t i) const; + ex add_series(const pseries &other) const; ex mul_const(const numeric &other) const; ex mul_series(const pseries &other) const;