X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=blobdiff_plain;ds=sidebyside;f=ginac%2Fpseries.cpp;h=fbb7583b3bc1b0b10e7c16eacff1b86f5049fc24;hb=4e3a4ac2bcb0837611ea31bc8fc05d84a20c33ac;hp=36fb4429294928b79dea975abb05f397e2931b3c;hpb=e58227e1112f989f3b5417e497a61d53fc2971fa;p=ginac.git diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index 36fb4429..fbb7583b 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -35,41 +35,17 @@ #include "utils.h" #include "debugmsg.h" -#ifndef NO_NAMESPACE_GINAC namespace GiNaC { -#endif // ndef NO_NAMESPACE_GINAC GINAC_IMPLEMENT_REGISTERED_CLASS(pseries, basic) /* - * Default constructor, destructor, copy constructor, assignment operator and helpers + * Default ctor, dtor, copy ctor, assignment operator and helpers */ pseries::pseries() : basic(TINFO_pseries) { - debugmsg("pseries default constructor", LOGLEVEL_CONSTRUCT); -} - -pseries::~pseries() -{ - debugmsg("pseries destructor", LOGLEVEL_DESTRUCT); - destroy(false); -} - -pseries::pseries(const pseries &other) -{ - debugmsg("pseries copy constructor", LOGLEVEL_CONSTRUCT); - copy(other); -} - -const pseries &pseries::operator=(const pseries & other) -{ - debugmsg("pseries operator=", LOGLEVEL_ASSIGNMENT); - if (this != &other) { - destroy(true); - copy(other); - } - return *this; + debugmsg("pseries default ctor", LOGLEVEL_CONSTRUCT); } void pseries::copy(const pseries &other) @@ -88,7 +64,7 @@ void pseries::destroy(bool call_parent) /* - * Other constructors + * Other ctors */ /** Construct pseries from a vector of coefficients and powers. @@ -102,7 +78,7 @@ void pseries::destroy(bool call_parent) * @return newly constructed pseries */ pseries::pseries(const ex &rel_, const epvector &ops_) : basic(TINFO_pseries), seq(ops_) { - debugmsg("pseries constructor from ex,epvector", LOGLEVEL_CONSTRUCT); + debugmsg("pseries ctor from ex,epvector", LOGLEVEL_CONSTRUCT); GINAC_ASSERT(is_ex_exactly_of_type(rel_, relational)); GINAC_ASSERT(is_ex_exactly_of_type(rel_.lhs(),symbol)); point = rel_.rhs(); @@ -117,7 +93,7 @@ pseries::pseries(const ex &rel_, const epvector &ops_) : basic(TINFO_pseries), s /** Construct object from archive_node. */ pseries::pseries(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) { - debugmsg("pseries constructor from archive_node", LOGLEVEL_CONSTRUCT); + debugmsg("pseries ctor from archive_node", LOGLEVEL_CONSTRUCT); for (unsigned int i=0; true; ++i) { ex rest; ex coeff; @@ -154,20 +130,15 @@ void pseries::archive(archive_node &n) const // functions overriding virtual functions from bases classes ////////// -basic *pseries::duplicate() const -{ - debugmsg("pseries duplicate", LOGLEVEL_DUPLICATE); - return new pseries(*this); -} - void pseries::print(std::ostream &os, unsigned upper_precedence) const { debugmsg("pseries print", LOGLEVEL_PRINT); if (precedence<=upper_precedence) os << "("; + // objects of type pseries must not have any zero entries, so the + // trivial (zero) pseries needs a special treatment here: + if (seq.size()==0) + os << '0'; for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) { - // omit zero terms - if (i->rest.is_zero()) - continue; // print a sign, if needed if (i!=seq.begin()) os << '+'; @@ -205,9 +176,8 @@ void pseries::printraw(std::ostream &os) const { debugmsg("pseries printraw", LOGLEVEL_PRINT); os << "pseries(" << var << ";" << point << ";"; - for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) { + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) os << "(" << (*i).rest << "," << (*i).coeff << "),"; - } os << ")"; } @@ -229,6 +199,31 @@ void pseries::printtree(std::ostream & os, unsigned indent) const point.printtree(os, indent+delta_indent); } +int pseries::compare_same_type(const basic & other) const +{ + GINAC_ASSERT(is_of_type(other, pseries)); + const pseries &o = static_cast(other); + + int cmpval = var.compare(o.var); + if (cmpval) + return cmpval; + cmpval = point.compare(o.point); + if (cmpval) + return cmpval; + + epvector::const_iterator it1 = seq.begin(), it2 = o.seq.begin(), it1end = seq.end(), it2end = o.seq.end(); + while ((it1 != it1end) && (it2 != it2end)) { + cmpval = it1->compare(*it2); + if (cmpval) + return cmpval; + it1++; it2++; + } + if (it1 == it1end) + return it2 == it2end ? 0 : -1; + + return 0; +} + /** Return the number of operands including a possible order term. */ unsigned pseries::nops(void) const { @@ -305,6 +300,13 @@ int pseries::ldegree(const symbol &s) const } } +/** Return coefficient of degree n in power series if s is the expansion + * variable. If the expansion point is nonzero, by definition the n=1 + * coefficient in s of a+b*(s-z)+c*(s-z)^2+Order((s-z)^3) is b (assuming + * the expansion took place in the s in the first place). + * If s is not the expansion variable, an attempt is made to convert the + * series to a polynomial and return the corresponding coefficient from + * there. */ ex pseries::coeff(const symbol &s, int n) const { if (var.is_equal(s)) { @@ -336,7 +338,7 @@ ex pseries::coeff(const symbol &s, int n) const return convert_to_poly().coeff(s, n); } - +/** Does nothing. */ ex pseries::collect(const symbol &s) const { return *this; @@ -407,14 +409,15 @@ ex pseries::subs(const lst & ls, const lst & lr) const /** Implementation of ex::expand() for a power series. It expands all the - * terms individually and returns the resulting series as a new pseries. - * @see ex::diff */ + * terms individually and returns the resulting series as a new pseries. */ ex pseries::expand(unsigned options) const { epvector newseq; - newseq.reserve(seq.size()); - for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) - newseq.push_back(expair(i->rest.expand(), i->coeff)); + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) { + ex restexp = i->rest.expand(); + if (!restexp.is_zero()) + newseq.push_back(expair(restexp, i->coeff)); + } return (new pseries(relational(var,point), newseq)) ->setflag(status_flags::dynallocated | status_flags::expanded); } @@ -447,10 +450,6 @@ ex pseries::derivative(const symbol & s) const } -/* - * Construct ordinary polynomial out of series - */ - /** Convert a pseries object to an ordinary polynomial. * * @param no_order flag: discard higher order terms */ @@ -470,6 +469,7 @@ ex pseries::convert_to_poly(bool no_order) const return e; } + /** Returns true if there is no order term, i.e. the series terminates and * false otherwise. */ bool pseries::is_terminating(void) const @@ -479,7 +479,7 @@ bool pseries::is_terminating(void) const /* - * Implementation of series expansion + * Implementations of series expansion */ /** Default implementation of ex::series(). This performs Taylor expansion. @@ -755,17 +755,46 @@ ex mul::series(const relational & r, int order, unsigned options) const * @param deg truncation order of series calculation */ ex pseries::power_const(const numeric &p, int deg) const { - int i; + // method: + // let A(x) be this series and for the time being let it start with a + // constant (later we'll generalize): + // A(x) = a_0 + a_1*x + a_2*x^2 + ... + // We want to compute + // C(x) = A(x)^p + // C(x) = c_0 + c_1*x + c_2*x^2 + ... + // Taking the derivative on both sides and multiplying with A(x) one + // immediately arrives at + // C'(x)*A(x) = p*C(x)*A'(x) + // Multiplying this out and comparing coefficients we get the recurrence + // formula + // c_i = (i*p*a_i*c_0 + ((i-1)*p-1)*a_{i-1}*c_1 + ... + // ... + (p-(i-1))*a_1*c_{i-1})/(a_0*i) + // which can easily be solved given the starting value c_0 = (a_0)^p. + // For the more general case where the leading coefficient of A(x) is not + // a constant, just consider A2(x) = A(x)*x^m, with some integer m and + // repeat the above derivation. The leading power of C2(x) = A2(x)^2 is + // then of course x^(p*m) but the recurrence formula still holds. + + if (seq.size()==0) { + // as a spacial case, handle the empty (zero) series honoring the + // usual power laws such as implemented in power::eval() + if (p.real().is_zero()) + throw (std::domain_error("pseries::power_const(): pow(0,I) is undefined")); + else if (p.real().is_negative()) + throw (pole_error("pseries::power_const(): division by zero",1)); + else + return *this; + } + const symbol *s = static_cast(var.bp); int ldeg = ldegree(*s); - // Calculate coefficients of powered series + // Compute coefficients of the powered series exvector co; co.reserve(deg); - ex co0; - co.push_back(co0 = power(coeff(*s, ldeg), p)); + co.push_back(power(coeff(*s, ldeg), p)); bool all_sums_zero = true; - for (i=1; i