Make add::eval(), mul::eval() work without compromise.
authorRichard Kreckel <kreckel@ginac.de>
Wed, 16 Dec 2015 12:00:30 +0000 (13:00 +0100)
committerRichard Kreckel <kreckel@ginac.de>
Wed, 16 Dec 2015 12:00:30 +0000 (13:00 +0100)
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
ginac/add.cpp
ginac/expairseq.cpp
ginac/function.cppy
ginac/mul.cpp
ginac/normal.cpp
ginac/power.cpp

index 711c64a..6ab452b 100644 (file)
@@ -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;
 }
index 978757f..ed58730 100644 (file)
@@ -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<add>(std::move(evaled), overall_coeff);
        }
 
@@ -341,13 +347,7 @@ ex add::eval(int level) const
                GINAC_ASSERT(!is_exactly_a<add>(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<add>(std::move(s), _ex0);
 }
@@ -551,7 +551,7 @@ expair add::combine_ex_with_coeff_to_pair(const ex & e,
        if (is_exactly_a<mul>(e)) {
                const mul &mulref(ex_to<mul>(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<mul>(mulref);
                mulcopy.overall_coeff = _ex1;
index 84b756f..9ca4070 100644 (file)
@@ -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<expairseq>(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
index 00e4b2a..789ca05 100644 (file)
@@ -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));
index 48d6478..a83b771 100644 (file)
@@ -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<mul>(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<mul>(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<mul>(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<numeric>(c));
+
+       // First, try a common shortcut:
+       if (is_exactly_a<symbol>(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<numeric>(p.coeff));
+       GINAC_ASSERT(is_exactly_a<numeric>(c));
+
        // to avoid duplication of power simplification rules,
        // we create a temporary power object
        // otherwise it would be hard to correctly evaluate
index ccf11e8..b2e8d10 100644 (file)
@@ -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());
 }
 
index 2aedc9b..afe2f1e 100644 (file)
@@ -1198,7 +1198,7 @@ ex power::expand_add(const add & a, long n, unsigned options)
                        composition_generator compositions(partition);
                        do {
                                const std::vector<int>& 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<power>(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<add>(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<mul>(r)) {
-                               result.push_back(a.combine_ex_with_coeff_to_pair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
-                                                                                _ex1));
+                               result.push_back(expair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
+                                                       _ex1));
                        } else {
-                               result.push_back(a.combine_ex_with_coeff_to_pair(dynallocate<power>(r, _ex2),
-                                                                                _ex1));
+                               result.push_back(expair(dynallocate<power>(r, _ex2),
+                                                       _ex1));
                        }
                } else {
                        if (is_exactly_a<mul>(r)) {
-                               result.push_back(a.combine_ex_with_coeff_to_pair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
-                                                                                ex_to<numeric>(c).power_dyn(*_num2_p)));
+                               result.push_back(expair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
+                                                       ex_to<numeric>(c).power_dyn(*_num2_p)));
                        } else {
-                               result.push_back(a.combine_ex_with_coeff_to_pair(dynallocate<power>(r, _ex2),
-                                                                                ex_to<numeric>(c).power_dyn(*_num2_p)));
+                               result.push_back(expair(dynallocate<power>(r, _ex2),
+                                                       ex_to<numeric>(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<numeric>(c)).mul_dyn(ex_to<numeric>(c1))));
+                       result.push_back(expair(mul(r,r1).expand(options),
+                                               _num2_p->mul(ex_to<numeric>(c)).mul_dyn(ex_to<numeric>(c1))));
                }
        }