X-Git-Url: https://ginac.de/ginac.git//ginac.git?a=blobdiff_plain;ds=sidebyside;f=ginac%2Fpseries.cpp;h=58f8d629c48c02ad535408cf2640fdc021e8817b;hb=cc65d5d5db813b87867255f32934c823e8ab1809;hp=459b46920b59793071185925a9efd61aa502c107;hpb=956eeb82c513a723e864edefa857133282cf692b;p=ginac.git diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index 459b4692..58f8d629 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -97,15 +97,17 @@ void pseries::destroy(bool call_parent) * the last coefficient can be Order(_ex1()) to represent a truncated, * non-terminating series. * - * @param var_ series variable (must hold a symbol) - * @param point_ expansion point + * @param rel__ expansion variable and point (must hold a relational) * @param ops_ vector of {coefficient, power} pairs (coefficient must not be zero) * @return newly constructed pseries */ -pseries::pseries(const ex &var_, const ex &point_, const epvector &ops_) - : basic(TINFO_pseries), seq(ops_), var(var_), point(point_) +pseries::pseries(const ex &rel_, const epvector &ops_) + : basic(TINFO_pseries), seq(ops_) { - debugmsg("pseries constructor from ex,ex,epvector", LOGLEVEL_CONSTRUCT); - GINAC_ASSERT(is_ex_exactly_of_type(var_, symbol)); + debugmsg("pseries constructor from rel,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(); + var = *static_cast(rel_.lhs().bp); } @@ -163,6 +165,9 @@ basic *pseries::duplicate() const void pseries::print(ostream &os, unsigned upper_precedence) const { debugmsg("pseries print", LOGLEVEL_PRINT); + // This could be made better, since series expansion at x==1 might print + // -1+2*x+Order((-1+x)^2) instead of 1+2*(-1+x)+Order((-1+x)^2), which is + // correct but can be rather confusing. convert_to_poly().print(os, upper_precedence); } @@ -176,6 +181,23 @@ void pseries::printraw(ostream &os) const os << ")"; } +void pseries::printtree(ostream & os, unsigned indent) const +{ + debugmsg("pseries printtree",LOGLEVEL_PRINT); + os << string(indent,' ') << "pseries " + << ", hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")" + << ", flags=" << flags << endl; + for (unsigned i=0; icoeff).to_int(); - if (pow == n) - return it->rest; - if (pow > n) - return _ex0(); - it++; - } - return _ex0(); + if (seq.size() == 0) + return _ex0(); + + // Binary search in sequence for given power + numeric looking_for = numeric(n); + int lo = 0, hi = seq.size() - 1; + while (lo <= hi) { + int mid = (lo + hi) / 2; + GINAC_ASSERT(is_ex_exactly_of_type(seq[mid].coeff, numeric)); + int cmp = ex_to_numeric(seq[mid].coeff).compare(looking_for); + switch (cmp) { + case -1: + lo = mid + 1; + break; + case 0: + return seq[mid].rest; + case 1: + hi = mid - 1; + break; + default: + throw(std::logic_error("pseries::coeff: compare() didn't return -1, 0 or 1")); + } + } + return _ex0(); } else return convert_to_poly().coeff(s, n); } +ex pseries::collect(const symbol &s) const +{ + return *this; +} + +/** Evaluate coefficients. */ ex pseries::eval(int level) const { if (level == 1) return this->hold(); + + if (level == -max_recursion_level) + throw (std::runtime_error("pseries::eval(): recursion limit exceeded")); // Construct a new series with evaluated coefficients epvector new_seq; @@ -269,13 +314,27 @@ ex pseries::eval(int level) const new_seq.push_back(expair(it->rest.eval(level-1), it->coeff)); it++; } - return (new pseries(var, point, new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated); + return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated); } -/** Evaluate numerically. The order term is dropped. */ +/** Evaluate coefficients numerically. */ ex pseries::evalf(int level) const { - return convert_to_poly().evalf(level); + if (level == 1) + return *this; + + if (level == -max_recursion_level) + throw (std::runtime_error("pseries::evalf(): recursion limit exceeded")); + + // Construct a new series with evaluated coefficients + epvector new_seq; + new_seq.reserve(seq.size()); + epvector::const_iterator it = seq.begin(), itend = seq.end(); + while (it != itend) { + new_seq.push_back(expair(it->rest.evalf(level-1), it->coeff)); + it++; + } + return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated); } ex pseries::subs(const lst & ls, const lst & lr) const @@ -295,7 +354,33 @@ ex pseries::subs(const lst & ls, const lst & lr) const new_seq.push_back(expair(it->rest.subs(ls, lr), it->coeff)); it++; } - return (new pseries(var, point.subs(ls, lr), new_seq))->setflag(status_flags::dynallocated); + return (new pseries(relational(var,point.subs(ls, lr)), new_seq))->setflag(status_flags::dynallocated); +} + +/** Implementation of ex::diff() for a power series. It treats the series as a + * polynomial. + * @see ex::diff */ +ex pseries::derivative(const symbol & s) const +{ + if (s == var) { + epvector new_seq; + epvector::const_iterator it = seq.begin(), itend = seq.end(); + + // FIXME: coeff might depend on var + while (it != itend) { + if (is_order_function(it->rest)) { + new_seq.push_back(expair(it->rest, it->coeff - 1)); + } else { + ex c = it->rest * it->coeff; + if (!c.is_zero()) + new_seq.push_back(expair(c, it->coeff - 1)); + } + it++; + } + return pseries(relational(var,point), new_seq); + } else { + return *this; + } } @@ -329,42 +414,48 @@ ex pseries::convert_to_poly(bool no_order) const /** Default implementation of ex::series(). This performs Taylor expansion. * @see ex::series */ -ex basic::series(const symbol & s, const ex & point, int order) const +ex basic::series(const relational & r, int order) const { epvector seq; numeric fac(1); ex deriv = *this; - ex coeff = deriv.subs(s == point); + ex coeff = deriv.subs(r); + const symbol *s = static_cast(r.lhs().bp); + if (!coeff.is_zero()) seq.push_back(expair(coeff, numeric(0))); int n; for (n=1; n(r.lhs().bp); + + if (this->is_equal(*s)) { if (order > 0 && !point.is_zero()) seq.push_back(expair(point, _ex0())); if (order > 1) @@ -373,7 +464,7 @@ ex symbol::series(const symbol & s, const ex & point, int order) const seq.push_back(expair(Order(_ex1()), numeric(order))); } else seq.push_back(expair(*this, _ex0())); - return pseries(s, point, seq); + return pseries(r, seq); } @@ -389,7 +480,7 @@ ex pseries::add_series(const pseries &other) const if (!is_compatible_to(other)) { epvector nul; nul.push_back(expair(Order(_ex1()), _ex0())); - return pseries(var, point, nul); + return pseries(relational(var,point), nul); } // Series addition @@ -447,19 +538,19 @@ ex pseries::add_series(const pseries &other) const } } } - return pseries(var, point, new_seq); + return pseries(relational(var,point), new_seq); } /** Implementation of ex::series() for sums. This performs series addition when * adding pseries objects. * @see ex::series */ -ex add::series(const symbol & s, const ex & point, int order) const +ex add::series(const relational & r, int order) const { ex acc; // Series accumulator // Get first term from overall_coeff - acc = overall_coeff.series(s, point, order); + acc = overall_coeff.series(r, order); // Add remaining terms epvector::const_iterator it = seq.begin(); @@ -469,7 +560,7 @@ ex add::series(const symbol & s, const ex & point, int order) const if (is_ex_exactly_of_type(it->rest, pseries)) op = it->rest; else - op = it->rest.series(s, point, order); + op = it->rest.series(r, order); if (!it->coeff.is_equal(_ex1())) op = ex_to_pseries(op).mul_const(ex_to_numeric(it->coeff)); @@ -498,7 +589,7 @@ ex pseries::mul_const(const numeric &other) const new_seq.push_back(*it); it++; } - return pseries(var, point, new_seq); + return pseries(relational(var,point), new_seq); } @@ -514,7 +605,7 @@ ex pseries::mul_series(const pseries &other) const if (!is_compatible_to(other)) { epvector nul; nul.push_back(expair(Order(_ex1()), _ex0())); - return pseries(var, point, nul); + return pseries(relational(var,point), nul); } // Series multiplication @@ -545,26 +636,26 @@ ex pseries::mul_series(const pseries &other) const ex a_coeff = coeff(*s, i); ex b_coeff = other.coeff(*s, cdeg-i); if (!is_order_function(a_coeff) && !is_order_function(b_coeff)) - co += coeff(*s, i) * other.coeff(*s, cdeg-i); + co += a_coeff * b_coeff; } if (!co.is_zero()) new_seq.push_back(expair(co, numeric(cdeg))); } if (higher_order_c < INT_MAX) new_seq.push_back(expair(Order(_ex1()), numeric(higher_order_c))); - return pseries(var, point, new_seq); + return pseries(relational(var,point), new_seq); } /** Implementation of ex::series() for product. This performs series * multiplication when multiplying series. * @see ex::series */ -ex mul::series(const symbol & s, const ex & point, int order) const +ex mul::series(const relational & r, int order) const { ex acc; // Series accumulator // Get first term from overall_coeff - acc = overall_coeff.series(s, point, order); + acc = overall_coeff.series(r, order); // Multiply with remaining terms epvector::const_iterator it = seq.begin(); @@ -577,7 +668,7 @@ ex mul::series(const symbol & s, const ex & point, int order) const acc = ex_to_pseries(acc).mul_const(ex_to_numeric(f)); continue; } else if (!is_ex_exactly_of_type(op, pseries)) - op = op.series(s, point, order); + op = op.series(r, order); if (!it->coeff.is_equal(_ex1())) op = ex_to_pseries(op).power_const(ex_to_numeric(it->coeff), order); @@ -632,27 +723,27 @@ ex pseries::power_const(const numeric &p, int deg) const } if (!higher_order && !all_sums_zero) new_seq.push_back(expair(Order(_ex1()), numeric(deg) + p * ldeg)); - return pseries(var, point, new_seq); + return pseries(relational(var,point), new_seq); } /** Implementation of ex::series() for powers. This performs Laurent expansion * of reciprocals of series at singularities. * @see ex::series */ -ex power::series(const symbol & s, const ex & point, int order) const +ex power::series(const relational & r, int order) const { ex e; if (!is_ex_exactly_of_type(basis, pseries)) { // Basis is not a series, may there be a singulary? if (!exponent.info(info_flags::negint)) - return basic::series(s, point, order); + return basic::series(r, order); // Expression is of type something^(-int), check for singularity - if (!basis.subs(s == point).is_zero()) - return basic::series(s, point, order); + if (!basis.subs(r).is_zero()) + return basic::series(r, order); // Singularity encountered, expand basis into series - e = basis.series(s, point, order); + e = basis.series(r, order); } else { // Basis is a series e = basis; @@ -663,19 +754,62 @@ ex power::series(const symbol & s, const ex & point, int order) const } +/** Re-expansion of a pseries object. */ +ex pseries::series(const relational & r, int order) const +{ + const ex p = r.rhs(); + GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol)); + const symbol *s = static_cast(r.lhs().bp); + + if (var.is_equal(*s) && point.is_equal(p)) { + if (order > degree(*s)) + return *this; + else { + epvector new_seq; + epvector::const_iterator it = seq.begin(), itend = seq.end(); + while (it != itend) { + int o = ex_to_numeric(it->coeff).to_int(); + if (o >= order) { + new_seq.push_back(expair(Order(_ex1()), o)); + break; + } + new_seq.push_back(*it); + it++; + } + return pseries(r, new_seq); + } + } else + return convert_to_poly().series(r, order); +} + + /** Compute the truncated series expansion of an expression. - * This function returns an expression containing an object of class pseries to - * represent the series. If the series does not terminate within the given + * This function returns an expression containing an object of class pseries + * to represent the series. If the series does not terminate within the given * truncation order, the last term of the series will be an order term. * - * @param s expansion variable - * @param point expansion point + * @param r expansion relation, lhs holds variable and rhs holds point * @param order truncation order of series calculations * @return an expression holding a pseries object */ -ex ex::series(const symbol &s, const ex &point, int order) const +ex ex::series(const ex & r, int order) const { GINAC_ASSERT(bp!=0); - return bp->series(s, point, order); + ex e; + relational rel_; + + if (is_ex_exactly_of_type(r,relational)) + rel_ = ex_to_relational(r); + else if (is_ex_exactly_of_type(r,symbol)) + rel_ = relational(r,_ex0()); + else + throw (std::logic_error("ex::series(): expansion point has unknown type")); + + try { + e = bp->series(rel_, order); + } catch (exception &x) { + throw (std::logic_error(string("unable to compute series (") + x.what() + ")")); + } + return e; }