* methods for series expansion. */
/*
- * GiNaC Copyright (C) 1999-2019 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
* 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_)
{
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);
/** 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<pseries>(relational(var,point), std::move(new_seq)).setflag(status_flags::evaluated);
+ return hold();
}
/** Evaluate coefficients numerically. */
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<pseries>(relational(var,point), std::move(new_seq)).setflag(status_flags::evaluated);
}
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<pseries>(var==point, std::move(v));
}
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<pseries>(var==point, std::move(v));
}
std::unique_ptr<epvector> 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();
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));
}
}
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)) {
newseq.reserve(seq.size());
std::copy(seq.begin(), i, std::back_inserter<epvector>(newseq));
if (!newcoeff.is_zero())
- newseq.push_back(expair(newcoeff, i->coeff));
+ newseq.emplace_back(expair(newcoeff, i->coeff));
}
}
}
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<pseries>(relational(var,point.subs(m, options)), std::move(newseq));
}
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<pseries>(relational(var,point), std::move(newseq)).setflag(options == 0 ? status_flags::expanded : 0);
}
// 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));
}
}
} 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));
}
}
}
// 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));
}
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;
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));
}
if (this->is_equal_same_type(ex_to<symbol>(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));
}
} 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;
}
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);
}
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<int>::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));
}
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 ) {
// 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
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<pseries>(relational(var,point), std::move(epv));
+ return dynallocate<pseries>(relational(var, point),
+ epvector{{Order(_ex1), deg}});
}
// O(x^n)^(-m) is undefined
// 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<numcoeff; ++i) {
ex sum = _ex0;
for (int j=1; j<=i; ++j) {
- ex c = coeff(var, j + ldeg);
+ ex c = coeff(var, j + base_ldeg);
if (is_order_function(c)) {
co.push_back(Order(_ex1));
break;
} else
sum += (p * j - (i - j)) * co[i - j] * c;
}
- co.push_back(sum / coeff(var, ldeg) / i);
+ co.push_back(sum / coeff(var, base_ldeg) / i);
}
// Construct new series (of non-zero coefficients)
epvector new_seq;
bool higher_order = false;
for (int i=0; i<numcoeff; ++i) {
- if (!co[i].is_zero())
- new_seq.push_back(expair(co[i], p * ldeg + i));
+ if (!co[i].is_zero()) {
+ new_seq.emplace_back(expair(co[i], new_ldeg + i));
+ }
if (is_order_function(co[i])) {
higher_order = true;
break;
}
}
- if (!higher_order)
- new_seq.push_back(expair(Order(_ex1), p * ldeg + numcoeff));
+ if (!higher_order && new_deg == deg) {
+ new_seq.emplace_back(expair{Order(_ex1), new_deg});
+ }
return pseries(relational(var,point), std::move(new_seq));
}
bool must_expand_basis = false;
try {
basis.subs(r, subs_options::no_pattern);
- } catch (pole_error) {
+ } catch (pole_error &) {
must_expand_basis = true;
}
bool exponent_is_regular = true;
try {
exponent.subs(r, subs_options::no_pattern);
- } catch (pole_error) {
+ } catch (pole_error &) {
exponent_is_regular = false;
}
if (basis.is_equal(r.lhs() - r.rhs())) {
epvector new_seq;
if (ex_to<numeric>(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));
}
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<pseries>(e).power_const(numexp, order);
- } catch (pole_error) {
+ } catch (pole_error &) {
epvector ser { expair(Order(_ex1), order) };
result = pseries(r, std::move(ser));
}
for (auto & it : seq) {
int o = ex_to<numeric>(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);
? currcoeff
: integral(x, a.subs(r), b.subs(r), currcoeff);
if (currcoeff != 0)
- fexpansion.push_back(
+ fexpansion.emplace_back(
expair(currcoeff, ex_to<pseries>(fseries).exponop(i)));
}
git://www.ginac.de / ginac.git / blobdiff