From ae6c004bd31e02dda37357d74b641c101b116c73 Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Wed, 16 Dec 2015 13:00:30 +0100 Subject: [PATCH] Make add::eval(), mul::eval() work without compromise. If a user asks to evaluate an object, the expectation is that the library returns an evaluated, canonical expression. Everything else is a bug. (It doesn't matter whether the user asks explicitly or by assigning to an ex.) It turns out that it was possible to construct objects which didn't fully evaluate. This patch fixes this by making eval() a little bit more careful. This obsoletes the need to go through combine_ex_with_coeff_to_pair() when constructing an epvector that is then used to construct an add or mul. (Alas, this was omitted frequently enough...) --- check/exam_paranoia.cpp | 28 ++++++++++++++++ ginac/add.cpp | 28 ++++++++-------- ginac/expairseq.cpp | 73 +++++++++++++++++------------------------ ginac/function.cppy | 4 +++ ginac/mul.cpp | 30 ++++++++++------- ginac/normal.cpp | 4 +-- ginac/power.cpp | 29 ++++++++-------- 7 files changed, 112 insertions(+), 84 deletions(-) diff --git a/check/exam_paranoia.cpp b/check/exam_paranoia.cpp index 711c64a6..6ab452ba 100644 --- a/check/exam_paranoia.cpp +++ b/check/exam_paranoia.cpp @@ -568,6 +568,33 @@ static unsigned exam_paranoia22() return 0; } +// Bug in expairseq::evalchildren(). +static unsigned exam_paranoia23() +{ + unsigned result = 0; + symbol x("x"); + + epvector v1; + v1.push_back(expair(1, 1)); + v1.push_back(expair(2*x, -1)); + ex e1 = add(v1); // Should be e==1-2*x, + if (!e1.is_equal(1-2*x)) { + clog << "Failure constructing " << e1 << " from add.\n"; + ++result; + } + + epvector v2; + v2.push_back(expair(x, 1)); + v2.push_back(expair(1,-1)); + ex e2 = mul(v2); // Should be e==x; + if (!e2.is_equal(x)) { + clog << "Failure constructing " << e2 << " from mul.\n"; + ++result; + } + + return result; +} + unsigned exam_paranoia() { unsigned result = 0; @@ -597,6 +624,7 @@ unsigned exam_paranoia() result += is_polynomial_false_positive(); cout << '.' << flush; result += exam_paranoia21(); cout << '.' << flush; result += exam_paranoia22(); cout << '.' << flush; + result += exam_paranoia23(); cout << '.' << flush; return result; } diff --git a/ginac/add.cpp b/ginac/add.cpp index 978757fb..ed587309 100644 --- a/ginac/add.cpp +++ b/ginac/add.cpp @@ -307,13 +307,13 @@ ex add::coeff(const ex & s, int n) const if (!restcoeff.is_zero()) { if (do_clifford) { if (clifford_max_label(restcoeff) == -1) { - coeffseq_cliff.push_back(combine_ex_with_coeff_to_pair(ncmul(restcoeff, dirac_ONE(rl)), i.coeff)); + coeffseq_cliff.push_back(expair(ncmul(restcoeff, dirac_ONE(rl)), i.coeff)); } else { - coeffseq_cliff.push_back(combine_ex_with_coeff_to_pair(restcoeff, i.coeff)); + coeffseq_cliff.push_back(expair(restcoeff, i.coeff)); nonscalar = true; } } - coeffseq.push_back(combine_ex_with_coeff_to_pair(restcoeff, i.coeff)); + coeffseq.push_back(expair(restcoeff, i.coeff)); } } @@ -330,9 +330,15 @@ ex add::coeff(const ex & s, int n) const * @param level cut-off in recursive evaluation */ ex add::eval(int level) const { - epvector evaled = evalchildren(level); + if ((level == 1) && (flags & status_flags::evaluated)) { + GINAC_ASSERT(seq.size()>0); + GINAC_ASSERT(seq.size()>1 || !overall_coeff.is_zero()); + return *this; + } + + const epvector evaled = evalchildren(level); if (unlikely(!evaled.empty())) { - // do more evaluation later + // start over evaluating a new object return dynallocate(std::move(evaled), overall_coeff); } @@ -341,13 +347,7 @@ ex add::eval(int level) const GINAC_ASSERT(!is_exactly_a(i.rest)); } #endif // def DO_GINAC_ASSERT - - if (flags & status_flags::evaluated) { - GINAC_ASSERT(seq.size()>0); - GINAC_ASSERT(seq.size()>1 || !overall_coeff.is_zero()); - return *this; - } - + int seq_size = seq.size(); if (seq_size == 0) { // +(;c) -> c @@ -491,7 +491,7 @@ ex add::derivative(const symbol & y) const // than the default implementation in basic::derivative() although // if performs the same function (differentiate each term). for (auto & it : seq) - s.push_back(combine_ex_with_coeff_to_pair(it.rest.diff(y), it.coeff)); + s.push_back(expair(it.rest.diff(y), it.coeff)); return dynallocate(std::move(s), _ex0); } @@ -551,7 +551,7 @@ expair add::combine_ex_with_coeff_to_pair(const ex & e, if (is_exactly_a(e)) { const mul &mulref(ex_to(e)); const ex &numfactor = mulref.overall_coeff; - if (numfactor.is_equal(_ex1)) + if (likely(numfactor.is_equal(_ex1))) return expair(e, c); mul & mulcopy = dynallocate(mulref); mulcopy.overall_coeff = _ex1; diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp index 84b756f9..9ca40709 100644 --- a/ginac/expairseq.cpp +++ b/ginac/expairseq.cpp @@ -244,7 +244,7 @@ ex expairseq::eval(int level) const if ((level==1) && (flags &status_flags::evaluated)) return *this; - epvector evaled = evalchildren(level); + const epvector evaled = evalchildren(level); if (!evaled.empty()) return dynallocate(std::move(evaled), overall_coeff).setflag(status_flags::evaluated); else @@ -1014,7 +1014,7 @@ epvector expairseq::expandchildren(unsigned options) const { auto cit = seq.begin(), last = seq.end(); while (cit!=last) { - const ex &expanded_ex = cit->rest.expand(options); + const ex expanded_ex = cit->rest.expand(options); if (!are_ex_trivially_equal(cit->rest,expanded_ex)) { // something changed, copy seq, eval and return it @@ -1022,25 +1022,20 @@ epvector expairseq::expandchildren(unsigned options) const s.reserve(seq.size()); // copy parts of seq which are known not to have changed - auto cit2 = seq.begin(); - while (cit2!=cit) { - s.push_back(*cit2); - ++cit2; - } + s.insert(s.begin(), seq.begin(), cit); // copy first changed element - s.push_back(combine_ex_with_coeff_to_pair(expanded_ex, - cit2->coeff)); - ++cit2; + s.push_back(expair(expanded_ex, cit->coeff)); + ++cit; // copy rest - while (cit2!=last) { - s.push_back(combine_ex_with_coeff_to_pair(cit2->rest.expand(options), - cit2->coeff)); - ++cit2; + while (cit != last) { + s.push_back(expair(cit->rest.expand(options), cit->coeff)); + ++cit; } return s; } + ++cit; } @@ -1055,42 +1050,35 @@ epvector expairseq::expandchildren(unsigned options) const * had to be changed. */ epvector expairseq::evalchildren(int level) const { - if (likely(level==1)) - return epvector(); // nothing had to be evaluated - if (level == -max_recursion_level) throw(std::runtime_error("max recursion level reached")); - - --level; + auto cit = seq.begin(), last = seq.end(); while (cit!=last) { - const ex evaled_ex = cit->rest.eval(level); - if (!are_ex_trivially_equal(cit->rest,evaled_ex)) { - - // something changed, copy seq, eval and return it + const ex evaled_rest = level==1 ? cit->rest : cit->rest.eval(level-1); + const expair evaled_pair = combine_ex_with_coeff_to_pair(evaled_rest, cit->coeff); + if (unlikely(!evaled_pair.is_equal(*cit))) { + + // something changed: copy seq, eval and return it epvector s; s.reserve(seq.size()); - + // copy parts of seq which are known not to have changed - auto cit2 = seq.begin(); - while (cit2!=cit) { - s.push_back(*cit2); - ++cit2; - } + s.insert(s.begin(), seq.begin(), cit); // copy first changed element - s.push_back(combine_ex_with_coeff_to_pair(evaled_ex, - cit2->coeff)); - ++cit2; + s.push_back(evaled_pair); + ++cit; // copy rest - while (cit2!=last) { - s.push_back(combine_ex_with_coeff_to_pair(cit2->rest.eval(level), - cit2->coeff)); - ++cit2; + while (cit != last) { + const ex evaled_rest = level==1 ? cit->rest : cit->rest.eval(level-1); + s.push_back(combine_ex_with_coeff_to_pair(evaled_rest, cit->coeff)); + ++cit; } - return std::move(s); + return s; } + ++cit; } @@ -1130,7 +1118,7 @@ epvector expairseq::subschildren(const exmap & m, unsigned options) const const ex &subsed_ex = orig_ex.subs(m, options); if (!are_ex_trivially_equal(orig_ex, subsed_ex)) { - // Something changed, copy seq, subs and return it + // Something changed: copy seq, subs and return it epvector s; s.reserve(seq.size()); @@ -1158,10 +1146,11 @@ epvector expairseq::subschildren(const exmap & m, unsigned options) const auto cit = seq.begin(), last = seq.end(); while (cit != last) { - const ex &subsed_ex = cit->rest.subs(m, options); - if (!are_ex_trivially_equal(cit->rest, subsed_ex)) { + const ex subsed_rest = cit->rest.subs(m, options); + const expair subsed_pair = combine_ex_with_coeff_to_pair(subsed_rest, cit->coeff); + if (!subsed_pair.is_equal(*cit)) { - // Something changed, copy seq, subs and return it + // Something changed: copy seq, subs and return it epvector s; s.reserve(seq.size()); @@ -1169,7 +1158,7 @@ epvector expairseq::subschildren(const exmap & m, unsigned options) const s.insert(s.begin(), seq.begin(), cit); // Copy first changed element - s.push_back(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff)); + s.push_back(subsed_pair); ++cit; // Copy rest diff --git a/ginac/function.cppy b/ginac/function.cppy index 00e4b2a6..789ca056 100644 --- a/ginac/function.cppy +++ b/ginac/function.cppy @@ -364,6 +364,10 @@ next_context: ex function::eval(int level) const { + if ((level == 1) && (flags & status_flags::evaluated)) { + return *this; + } + if (level>1) { // first evaluate children, then we will end up here again return function(serial,evalchildren(level)); diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 48d64787..a83b771d 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -468,18 +468,18 @@ ex mul::coeff(const ex & s, int n) const * @param level cut-off in recursive evaluation */ ex mul::eval(int level) const { - epvector evaled = evalchildren(level); - if (unlikely(!evaled.empty())) { - // do more evaluation later - return dynallocate(std::move(evaled), overall_coeff); - } - - if (flags & status_flags::evaluated) { + if ((level == 1) && (flags & status_flags::evaluated)) { GINAC_ASSERT(seq.size()>0); GINAC_ASSERT(seq.size()>1 || !overall_coeff.is_equal(_ex1)); return *this; } - + + const epvector evaled = evalchildren(level); + if (unlikely(!evaled.empty())) { + // start over evaluating a new object + return dynallocate(std::move(evaled), overall_coeff); + } + size_t seq_size = seq.size(); if (overall_coeff.is_zero()) { // *(...,x;0) -> 0 @@ -588,10 +588,9 @@ ex mul::evalf(int level) const --level; for (auto & it : seq) { - s.push_back(combine_ex_with_coeff_to_pair(it.rest.evalf(level), - it.coeff)); + s.push_back(expair(it.rest.evalf(level), it.coeff)); } - return mul(std::move(s), overall_coeff.evalf(level)); + return dynallocate(std::move(s), overall_coeff.evalf(level)); } void mul::find_real_imag(ex & rp, ex & ip) const @@ -966,6 +965,12 @@ expair mul::split_ex_to_pair(const ex & e) const expair mul::combine_ex_with_coeff_to_pair(const ex & e, const ex & c) const { + GINAC_ASSERT(is_exactly_a(c)); + + // First, try a common shortcut: + if (is_exactly_a(e)) + return expair(e, c); + // to avoid duplication of power simplification rules, // we create a temporary power object // otherwise it would be hard to correctly evaluate @@ -979,6 +984,9 @@ expair mul::combine_ex_with_coeff_to_pair(const ex & e, expair mul::combine_pair_with_coeff_to_pair(const expair & p, const ex & c) const { + GINAC_ASSERT(is_exactly_a(p.coeff)); + GINAC_ASSERT(is_exactly_a(c)); + // to avoid duplication of power simplification rules, // we create a temporary power object // otherwise it would be hard to correctly evaluate diff --git a/ginac/normal.cpp b/ginac/normal.cpp index ccf11e8d..b2e8d109 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -2555,7 +2555,7 @@ ex expairseq::to_rational(exmap & repl) const if (oc.info(info_flags::numeric)) return thisexpairseq(std::move(s), overall_coeff); else - s.push_back(combine_ex_with_coeff_to_pair(oc, _ex1)); + s.push_back(expair(oc, _ex1)); return thisexpairseq(std::move(s), default_overall_coeff()); } @@ -2571,7 +2571,7 @@ ex expairseq::to_polynomial(exmap & repl) const if (oc.info(info_flags::numeric)) return thisexpairseq(std::move(s), overall_coeff); else - s.push_back(combine_ex_with_coeff_to_pair(oc, _ex1)); + s.push_back(expair(oc, _ex1)); return thisexpairseq(std::move(s), default_overall_coeff()); } diff --git a/ginac/power.cpp b/ginac/power.cpp index 2aedc9bf..afe2f1ef 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -1198,7 +1198,7 @@ ex power::expand_add(const add & a, long n, unsigned options) composition_generator compositions(partition); do { const std::vector& exponent = compositions.current(); - exvector monomial; + epvector monomial; monomial.reserve(msize); numeric factor = coeff; for (unsigned i = 0; i < exponent.size(); ++i) { @@ -1216,22 +1216,21 @@ ex power::expand_add(const add & a, long n, unsigned options) // optimize away } else if (exponent[i] == 1) { // optimized - monomial.push_back(r); + monomial.push_back(expair(r, _ex1)); if (c != *_num1_p) factor = factor.mul(c); } else { // general case exponent[i] > 1 - monomial.push_back(dynallocate(r, exponent[i])); + monomial.push_back(expair(r, exponent[i])); if (c != *_num1_p) factor = factor.mul(c.power(exponent[i])); } } - result.push_back(a.combine_ex_with_coeff_to_pair(mul(monomial).expand(options), factor)); + result.push_back(expair(mul(monomial).expand(options), factor)); } while (compositions.next()); } while (partitions.next()); } GINAC_ASSERT(result.size() == result_size); - if (a.overall_coeff.is_zero()) { return dynallocate(std::move(result)).setflag(status_flags::expanded); } else { @@ -1270,27 +1269,27 @@ ex power::expand_add_2(const add & a, unsigned options) if (c.is_equal(_ex1)) { if (is_exactly_a(r)) { - result.push_back(a.combine_ex_with_coeff_to_pair(expand_mul(ex_to(r), *_num2_p, options, true), - _ex1)); + result.push_back(expair(expand_mul(ex_to(r), *_num2_p, options, true), + _ex1)); } else { - result.push_back(a.combine_ex_with_coeff_to_pair(dynallocate(r, _ex2), - _ex1)); + result.push_back(expair(dynallocate(r, _ex2), + _ex1)); } } else { if (is_exactly_a(r)) { - result.push_back(a.combine_ex_with_coeff_to_pair(expand_mul(ex_to(r), *_num2_p, options, true), - ex_to(c).power_dyn(*_num2_p))); + result.push_back(expair(expand_mul(ex_to(r), *_num2_p, options, true), + ex_to(c).power_dyn(*_num2_p))); } else { - result.push_back(a.combine_ex_with_coeff_to_pair(dynallocate(r, _ex2), - ex_to(c).power_dyn(*_num2_p))); + result.push_back(expair(dynallocate(r, _ex2), + ex_to(c).power_dyn(*_num2_p))); } } for (epvector::const_iterator cit1=cit0+1; cit1!=last; ++cit1) { const ex & r1 = cit1->rest; const ex & c1 = cit1->coeff; - result.push_back(a.combine_ex_with_coeff_to_pair(mul(r,r1).expand(options), - _num2_p->mul(ex_to(c)).mul_dyn(ex_to(c1)))); + result.push_back(expair(mul(r,r1).expand(options), + _num2_p->mul(ex_to(c)).mul_dyn(ex_to(c1)))); } } -- 2.49.0