X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fpseries.cpp;h=2ea4365f57324d868b6e956e9e6d940b92c52fa6;hp=31a92002958399818659b65e83094b498a66dbc4;hb=c12c8ec3c5cf0c75f061f6c52d04206277bbdcca;hpb=6c946d4c762f5a0d6a3b742f03556dd018d63886 diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index 31a92002..2ea4365f 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -4,7 +4,7 @@ * methods for series expansion. */ /* - * GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2016 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 @@ -71,6 +71,19 @@ pseries::pseries() { } pseries::pseries(const ex &rel_, const epvector &ops_) : seq(ops_) { +#ifdef DO_GINAC_ASSERT + auto i = seq.begin(); + while (i != seq.end()) { + auto ip1 = i+1; + if (ip1 != seq.end()) + GINAC_ASSERT(!is_order_function(i->rest)); + else + break; + GINAC_ASSERT(is_a(i->coeff)); + GINAC_ASSERT(ex_to(i->coeff) < ex_to(ip1->coeff)); + ++i; + } +#endif // def DO_GINAC_ASSERT GINAC_ASSERT(is_a(rel_)); GINAC_ASSERT(is_a(rel_.lhs())); point = rel_.rhs(); @@ -79,6 +92,19 @@ pseries::pseries(const ex &rel_, const epvector &ops_) pseries::pseries(const ex &rel_, epvector &&ops_) : seq(std::move(ops_)) { +#ifdef DO_GINAC_ASSERT + auto i = seq.begin(); + while (i != seq.end()) { + auto ip1 = i+1; + if (ip1 != seq.end()) + GINAC_ASSERT(!is_order_function(i->rest)); + else + break; + GINAC_ASSERT(is_a(i->coeff)); + GINAC_ASSERT(ex_to(i->coeff) < ex_to(ip1->coeff)); + ++i; + } +#endif // def DO_GINAC_ASSERT GINAC_ASSERT(is_a(rel_)); GINAC_ASSERT(is_a(rel_.lhs())); point = rel_.rhs(); @@ -177,7 +203,7 @@ void pseries::print_series(const print_context & c, const char *openbrace, const } } } else - Order(power(var-point,i->coeff)).print(c); + Order(pow(var - point, i->coeff)).print(c); ++i; } @@ -281,8 +307,8 @@ ex pseries::op(size_t i) const 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); + return Order(pow(var-point, seq[i].coeff)); + return seq[i].rest * pow(var - point, seq[i].coeff); } /** Return degree of highest power of the series. This is usually the exponent @@ -290,25 +316,17 @@ ex pseries::op(size_t i) const * series is examined termwise. */ int pseries::degree(const ex &s) const { - if (var.is_equal(s)) { - // Return last exponent - if (seq.size()) - return ex_to((seq.end()-1)->coeff).to_int(); - else - return 0; - } else { - epvector::const_iterator it = seq.begin(), itend = seq.end(); - if (it == itend) - return 0; - int max_pow = std::numeric_limits::min(); - while (it != itend) { - int pow = it->rest.degree(s); - if (pow > max_pow) - max_pow = pow; - ++it; - } - return max_pow; - } + if (seq.empty()) + return 0; + + if (var.is_equal(s)) + // Return last/greatest exponent + return ex_to((seq.end()-1)->coeff).to_int(); + + int max_pow = std::numeric_limits::min(); + for (auto & it : seq) + max_pow = std::max(max_pow, it.rest.degree(s)); + return max_pow; } /** Return degree of lowest power of the series. This is usually the exponent @@ -318,25 +336,17 @@ int pseries::degree(const ex &s) const * I.e.: (1-x) + (1-x)^2 + Order((1-x)^3) has ldegree(x) 1, not 0. */ int pseries::ldegree(const ex &s) const { - if (var.is_equal(s)) { - // Return first exponent - if (seq.size()) - return ex_to((seq.begin())->coeff).to_int(); - else - return 0; - } else { - epvector::const_iterator it = seq.begin(), itend = seq.end(); - if (it == itend) - return 0; - int min_pow = std::numeric_limits::max(); - while (it != itend) { - int pow = it->rest.ldegree(s); - if (pow < min_pow) - min_pow = pow; - ++it; - } - return min_pow; - } + if (seq.empty()) + return 0; + + if (var.is_equal(s)) + // Return first/smallest exponent + return ex_to((seq.begin())->coeff).to_int(); + + int min_pow = std::numeric_limits::max(); + for (auto & it : seq) + min_pow = std::min(min_pow, it.rest.degree(s)); + return min_pow; } /** Return coefficient of degree n in power series if s is the expansion @@ -389,35 +399,25 @@ 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()); - epvector::const_iterator it = seq.begin(), itend = seq.end(); - while (it != itend) { - new_seq.push_back(expair(it->rest, it->coeff)); - ++it; - } + 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); } /** Evaluate coefficients numerically. */ -ex pseries::evalf(int level) const +ex pseries::evalf() const { - 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; - } + 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); } @@ -500,8 +500,7 @@ ex pseries::evalm() const ex newcoeff = i->rest.evalm(); if (!newcoeff.is_zero()) newseq.push_back(expair(newcoeff, i->coeff)); - } - else { + } else { ex newcoeff = i->rest.evalm(); if (!are_ex_trivially_equal(newcoeff, i->rest)) { something_changed = true; @@ -589,9 +588,9 @@ ex pseries::convert_to_poly(bool no_order) const for (auto & it : seq) { if (is_order_function(it.rest)) { if (!no_order) - e += Order(power(var - point, it.coeff)); + e += Order(pow(var - point, it.coeff)); } else - e += it.rest * power(var - point, it.coeff); + e += it.rest * pow(var - point, it.coeff); } return e; } @@ -644,7 +643,7 @@ ex basic::series(const relational & r, int order, unsigned options) const int n; for (n=1; n::max(); @@ -837,18 +836,30 @@ ex pseries::mul_series(const pseries &other) const higher_order_a = a_max + b_min; if (is_order_function(other.coeff(var, b_max))) higher_order_b = b_max + a_min; - int higher_order_c = std::min(higher_order_a, higher_order_b); + const int higher_order_c = std::min(higher_order_a, higher_order_b); if (cdeg_max >= higher_order_c) cdeg_max = higher_order_c - 1; - + + std::map rest_map_a, rest_map_b; + for (const auto& it : seq) + rest_map_a[ex_to(it.coeff).to_int()] = it.rest; + + if (other.var.is_equal(var)) + for (const auto& it : other.seq) + rest_map_b[ex_to(it.coeff).to_int()] = it.rest; + for (int cdeg=cdeg_min; cdeg<=cdeg_max; ++cdeg) { ex co = _ex0; // c(i)=a(0)b(i)+...+a(i)b(0) for (int i=a_min; cdeg-i>=b_min; ++i) { - ex a_coeff = coeff(var, i); - ex b_coeff = other.coeff(var, cdeg-i); - if (!is_order_function(a_coeff) && !is_order_function(b_coeff)) - co += a_coeff * b_coeff; + const auto& ita = rest_map_a.find(i); + if (ita == rest_map_a.end()) + continue; + const auto& itb = rest_map_b.find(cdeg-i); + if (itb == rest_map_b.end()) + continue; + if (!is_order_function(ita->second) && !is_order_function(itb->second)) + co += ita->second * itb->second; } if (!co.is_zero()) new_seq.push_back(expair(co, numeric(cdeg))); @@ -1027,7 +1038,7 @@ ex pseries::power_const(const numeric &p, int deg) const // Compute coefficients of the powered series exvector co; co.reserve(numcoeff); - co.push_back(power(coeff(var, ldeg), p)); + co.push_back(pow(coeff(var, ldeg), p)); for (int i=1; i