X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fpseries.cpp;h=6fe4e54c1871c015cf0b4284c90481e85814f97c;hp=dcc25a0493a767d45980c2d1bc0be27b3f6d1ce4;hb=HEAD;hpb=63f3e977f92d51ea173382a9b7c4c3b18bda7b8e diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index dcc25a04..b3ec633e 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -4,7 +4,7 @@ * methods for series expansion. */ /* - * GiNaC Copyright (C) 1999-2018 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2024 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -66,8 +66,7 @@ pseries::pseries() { } * non-terminating series. * * @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 */ + * @param ops_ vector of {coefficient, power} pairs (coefficient must not be zero) */ pseries::pseries(const ex &rel_, const epvector &ops_) : seq(ops_) { @@ -119,17 +118,15 @@ pseries::pseries(const ex &rel_, epvector &&ops_) void pseries::read_archive(const archive_node &n, lst &sym_lst) { inherited::read_archive(n, sym_lst); - auto first = n.find_first("coeff"); - auto last = n.find_last("power"); - ++last; - seq.reserve((last-first)/2); + auto range = n.find_property_range("coeff", "power"); + seq.reserve((range.end-range.begin)/2); - for (auto loc = first; loc < last;) { + for (auto loc = range.begin; loc < range.end;) { ex rest; ex coeff; n.find_ex_by_loc(loc++, rest, sym_lst); n.find_ex_by_loc(loc++, coeff, sym_lst); - seq.push_back(expair(rest, coeff)); + seq.emplace_back(expair(rest, coeff)); } n.find_ex("var", var, sym_lst); @@ -396,17 +393,7 @@ ex pseries::collect(const ex &s, bool distributed) const /** Perform coefficient-wise automatic term rewriting rules in this class. */ ex pseries::eval() const { - if (flags & status_flags::evaluated) { - return *this; - } - - // Construct a new series with evaluated coefficients - epvector new_seq; - new_seq.reserve(seq.size()); - for (auto & it : seq) - new_seq.push_back(expair(it.rest, it.coeff)); - - return dynallocate(relational(var,point), std::move(new_seq)).setflag(status_flags::evaluated); + return hold(); } /** Evaluate coefficients numerically. */ @@ -416,7 +403,7 @@ ex pseries::evalf() const epvector new_seq; new_seq.reserve(seq.size()); for (auto & it : seq) - new_seq.push_back(expair(it.rest, it.coeff)); + new_seq.emplace_back(expair(it.rest.evalf(), it.coeff)); return dynallocate(relational(var,point), std::move(new_seq)).setflag(status_flags::evaluated); } @@ -447,7 +434,7 @@ ex pseries::real_part() const epvector v; v.reserve(seq.size()); for (auto & it : seq) - v.push_back(expair((it.rest).real_part(), it.coeff)); + v.emplace_back(expair(it.rest.real_part(), it.coeff)); return dynallocate(var==point, std::move(v)); } @@ -462,7 +449,7 @@ ex pseries::imag_part() const epvector v; v.reserve(seq.size()); for (auto & it : seq) - v.push_back(expair((it.rest).imag_part(), it.coeff)); + v.emplace_back(expair(it.rest.imag_part(), it.coeff)); return dynallocate(var==point, std::move(v)); } @@ -471,7 +458,7 @@ ex pseries::eval_integ() const std::unique_ptr newseq(nullptr); for (auto i=seq.begin(); i!=seq.end(); ++i) { if (newseq) { - newseq->push_back(expair(i->rest.eval_integ(), i->coeff)); + newseq->emplace_back(expair(i->rest.eval_integ(), i->coeff)); continue; } ex newterm = i->rest.eval_integ(); @@ -480,7 +467,7 @@ ex pseries::eval_integ() const newseq->reserve(seq.size()); for (auto j=seq.begin(); j!=i; ++j) newseq->push_back(*j); - newseq->push_back(expair(newterm, i->coeff)); + newseq->emplace_back(expair(newterm, i->coeff)); } } @@ -499,7 +486,7 @@ ex pseries::evalm() const if (something_changed) { ex newcoeff = i->rest.evalm(); if (!newcoeff.is_zero()) - newseq.push_back(expair(newcoeff, i->coeff)); + newseq.emplace_back(expair(newcoeff, i->coeff)); } else { ex newcoeff = i->rest.evalm(); if (!are_ex_trivially_equal(newcoeff, i->rest)) { @@ -507,7 +494,7 @@ ex pseries::evalm() const newseq.reserve(seq.size()); std::copy(seq.begin(), i, std::back_inserter(newseq)); if (!newcoeff.is_zero()) - newseq.push_back(expair(newcoeff, i->coeff)); + newseq.emplace_back(expair(newcoeff, i->coeff)); } } } @@ -530,7 +517,7 @@ ex pseries::subs(const exmap & m, unsigned options) const epvector newseq; newseq.reserve(seq.size()); for (auto & it : seq) - newseq.push_back(expair(it.rest.subs(m, options), it.coeff)); + newseq.emplace_back(expair(it.rest.subs(m, options), it.coeff)); return dynallocate(relational(var,point.subs(m, options)), std::move(newseq)); } @@ -542,7 +529,7 @@ ex pseries::expand(unsigned options) const for (auto & it : seq) { ex restexp = it.rest.expand(); if (!restexp.is_zero()) - newseq.push_back(expair(restexp, it.coeff)); + newseq.emplace_back(expair(restexp, it.coeff)); } return dynallocate(relational(var,point), std::move(newseq)).setflag(options == 0 ? status_flags::expanded : 0); } @@ -558,11 +545,11 @@ ex pseries::derivative(const symbol & s) const // FIXME: coeff might depend on var for (auto & it : seq) { if (is_order_function(it.rest)) { - new_seq.push_back(expair(it.rest, it.coeff - 1)); + new_seq.emplace_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)); + new_seq.emplace_back(expair(c, it.coeff - 1)); } } @@ -574,7 +561,7 @@ ex pseries::derivative(const symbol & s) const } else { ex c = it.rest.diff(s); if (!c.is_zero()) - new_seq.push_back(expair(c, it.coeff)); + new_seq.emplace_back(expair(c, it.coeff)); } } } @@ -628,7 +615,7 @@ ex basic::series(const relational & r, int order, unsigned options) const // default for order-values that make no sense for Taylor expansion if ((order <= 0) && this->has(s)) { - seq.push_back(expair(Order(_ex1), order)); + seq.emplace_back(expair(Order(_ex1), order)); return pseries(r, std::move(seq)); } @@ -638,7 +625,7 @@ ex basic::series(const relational & r, int order, unsigned options) const ex coeff = deriv.subs(r, subs_options::no_pattern); if (!coeff.is_zero()) { - seq.push_back(expair(coeff, _ex0)); + seq.emplace_back(expair(coeff, _ex0)); } int n; @@ -653,13 +640,13 @@ ex basic::series(const relational & r, int order, unsigned options) const coeff = deriv.subs(r, subs_options::no_pattern); if (!coeff.is_zero()) - seq.push_back(expair(fac * coeff, n)); + seq.emplace_back(expair(fac * coeff, n)); } // Higher-order terms, if present deriv = deriv.diff(s); if (!deriv.expand().is_zero()) - seq.push_back(expair(Order(_ex1), n)); + seq.emplace_back(expair(Order(_ex1), n)); return pseries(r, std::move(seq)); } @@ -674,13 +661,13 @@ ex symbol::series(const relational & r, int order, unsigned options) const if (this->is_equal_same_type(ex_to(r.lhs()))) { if (order > 0 && !point.is_zero()) - seq.push_back(expair(point, _ex0)); + seq.emplace_back(expair(point, _ex0)); if (order > 1) - seq.push_back(expair(_ex1, _ex1)); + seq.emplace_back(expair(_ex1, _ex1)); else - seq.push_back(expair(Order(_ex1), numeric(order))); + seq.emplace_back(expair(Order(_ex1), numeric(order))); } else - seq.push_back(expair(*this, _ex0)); + seq.emplace_back(expair(*this, _ex0)); return pseries(r, std::move(seq)); } @@ -741,12 +728,12 @@ ex pseries::add_series(const pseries &other) const } else { // Add coefficient of a and b if (is_order_function((*a).rest) || is_order_function((*b).rest)) { - new_seq.push_back(expair(Order(_ex1), (*a).coeff)); + new_seq.emplace_back(expair(Order(_ex1), (*a).coeff)); break; // Order term ends the sequence } else { ex sum = (*a).rest + (*b).rest; if (!(sum.is_zero())) - new_seq.push_back(expair(sum, numeric(pow_a))); + new_seq.emplace_back(expair(sum, numeric(pow_a))); ++a; ++b; } @@ -795,7 +782,7 @@ ex pseries::mul_const(const numeric &other) const for (auto & it : seq) { if (!is_order_function(it.rest)) - new_seq.push_back(expair(it.rest * other, it.coeff)); + new_seq.emplace_back(expair(it.rest * other, it.coeff)); else new_seq.push_back(it); } @@ -862,10 +849,10 @@ ex pseries::mul_series(const pseries &other) const co += ita->second * itb->second; } if (!co.is_zero()) - new_seq.push_back(expair(co, numeric(cdeg))); + new_seq.emplace_back(expair(co, numeric(cdeg))); } if (higher_order_c < std::numeric_limits::max()) - new_seq.push_back(expair(Order(_ex1), numeric(higher_order_c))); + new_seq.emplace_back(expair(Order(_ex1), numeric(higher_order_c))); return pseries(relational(var, point), std::move(new_seq)); } @@ -903,7 +890,7 @@ ex mul::series(const relational & r, int order, unsigned options) const bool flag_redo = false; try { real_ldegree = buf.expand().ldegree(sym-r.rhs()); - } catch (std::runtime_error) {} + } catch (std::runtime_error &) {} if (real_ldegree == 0) { if ( factor < 0 ) { @@ -959,9 +946,8 @@ ex mul::series(const relational & r, int order, unsigned options) const int degsum = std::accumulate(ldegrees.begin(), ldegrees.end(), 0); - if (degsum >= order) { - epvector epv { expair(Order(_ex1), order) }; - return dynallocate(r, std::move(epv)); + if (degsum > order) { + return dynallocate(r, epvector{{Order(_ex1), order}}); } // Multiply with remaining terms @@ -1006,8 +992,8 @@ ex pseries::power_const(const numeric &p, int deg) const // 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. + // repeat the above derivation. The leading power of C2(x) = A2(x)^p is + // then of course a_0^p*x^(p*m) but the recurrence formula still holds. if (seq.empty()) { // as a special case, handle the empty (zero) series honoring the @@ -1019,16 +1005,27 @@ ex pseries::power_const(const numeric &p, int deg) const else return *this; } - - const int ldeg = ldegree(var); - if (!(p*ldeg).is_integer()) + + const int base_ldeg = ldegree(var); + if (!(p*base_ldeg).is_integer()) throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series"); + int new_ldeg = (p*base_ldeg).to_int(); + + const int base_deg = degree(var); + int new_deg = deg; + if (p.is_pos_integer()) { + // No need to compute beyond p*base_deg. + new_deg = std::min((p*base_deg).to_int(), deg); + } // adjust number of coefficients - int numcoeff = deg - (p*ldeg).to_int(); + int numcoeff = new_deg - new_ldeg; + if (new_deg < deg) + numcoeff += 1; + if (numcoeff <= 0) { - epvector epv { expair(Order(_ex1), deg) }; - return dynallocate(relational(var,point), std::move(epv)); + return dynallocate(relational(var, point), + epvector{{Order(_ex1), deg}}); } // O(x^n)^(-m) is undefined @@ -1038,33 +1035,35 @@ ex pseries::power_const(const numeric &p, int deg) const // Compute coefficients of the powered series exvector co; co.reserve(numcoeff); - co.push_back(pow(coeff(var, ldeg), p)); + co.push_back(pow(coeff(var, base_ldeg), p)); for (int i=1; i(exponent).to_int() < order) - new_seq.push_back(expair(_ex1, exponent)); + new_seq.emplace_back(expair(_ex1, exponent)); else - new_seq.push_back(expair(Order(_ex1), exponent)); + new_seq.emplace_back(expair(Order(_ex1), exponent)); return pseries(r, std::move(new_seq)); } @@ -1164,12 +1163,13 @@ ex power::series(const relational & r, int order, unsigned options) const if (!(real_ldegree*numexp).is_integer()) throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series"); - ex e = basis.series(r, (order + real_ldegree*(1-numexp)).to_int(), options); + int extra_terms = (real_ldegree*(1-numexp)).to_int(); + ex e = basis.series(r, order + std::max(0, extra_terms), options); ex result; try { result = ex_to(e).power_const(numexp, order); - } catch (pole_error) { + } catch (pole_error &) { epvector ser { expair(Order(_ex1), order) }; result = pseries(r, std::move(ser)); } @@ -1193,7 +1193,7 @@ ex pseries::series(const relational & r, int order, unsigned options) const for (auto & it : seq) { int o = ex_to(it.coeff).to_int(); if (o >= order) { - new_seq.push_back(expair(Order(_ex1), o)); + new_seq.emplace_back(expair(Order(_ex1), o)); break; } new_seq.push_back(it); @@ -1219,7 +1219,7 @@ ex integral::series(const relational & r, int order, unsigned options) const ? currcoeff : integral(x, a.subs(r), b.subs(r), currcoeff); if (currcoeff != 0) - fexpansion.push_back( + fexpansion.emplace_back( expair(currcoeff, ex_to(fseries).exponop(i))); }