X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fpseries.cpp;h=dcc25a0493a767d45980c2d1bc0be27b3f6d1ce4;hp=cc756ac43e9c70b0cb3bfec01af080aa8a63ed90;hb=d0ff428fb5b7bb565a0aea69e05e5705d840c16d;hpb=82df718524319471d3a92fb051329aa8cd529c22 diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index cc756ac4..dcc25a04 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-2018 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(); @@ -255,7 +281,7 @@ int pseries::compare_same_type(const basic & other) const return cmpval; // ...and if that failed the individual elements - epvector::const_iterator it = seq.begin(), o_it = o.seq.begin(); + auto it = seq.begin(), o_it = o.seq.begin(); while (it!=seq.end() && o_it!=o.seq.end()) { cmpval = it->compare(*o_it); if (cmpval) @@ -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); } @@ -823,11 +823,11 @@ ex pseries::mul_series(const pseries &other) const // Series multiplication epvector new_seq; - int a_max = degree(var); - int b_max = other.degree(var); - int a_min = ldegree(var); - int b_min = other.ldegree(var); - int cdeg_min = a_min + b_min; + const int a_max = degree(var); + const int b_max = other.degree(var); + const int a_min = ldegree(var); + const int b_min = other.ldegree(var); + const int cdeg_min = a_min + b_min; int cdeg_max = a_max + b_max; int higher_order_a = std::numeric_limits::max(); @@ -836,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)));