fixed a problem with expand() not fully expanding

diff --git a/ginac/mul.cpp b/ginac/mul.cpp
index 2352729c1adbf74159a29d88b400ef4231b61ade..4af3bb3ea81eccdf4d462baaa16ad23520a474be 100644 (file)
--- a/ginac/mul.cpp
+++ b/ginac/mul.cpp
@@ -839,8 +839,11 @@ ex mul::expand(unsigned options) const
        // not sums
        int number_of_adds = 0;
        ex last_expanded = _ex1;
+
        epvector non_adds;
        non_adds.reserve(expanded_seq.size());
+       bool non_adds_has_sums = false; // Look for sums or powers of sums in the non_adds (we need this later)
+
        epvector::const_iterator cit = expanded_seq.begin(), last = expanded_seq.end();
        while (cit != last) {
                if (is_exactly_a<add>(cit->rest) &&
@@ -890,7 +893,7 @@ ex mul::expand(unsigned options) const
                                        for (epvector::const_iterator i2=add2begin; i2!=add2end; ++i2) {
                                                // Don't push_back expairs which might have a rest that evaluates to a numeric,
                                                // since that would violate an invariant of expairseq:
-                                               const ex rest = (new mul(i1->rest, i2->rest))->setflag(status_flags::dynallocated);
+                                               const ex rest = ex((new mul(i1->rest, i2->rest))->setflag(status_flags::dynallocated)).expand();
                                                if (is_exactly_a<numeric>(rest))
                                                        oc += ex_to<numeric>(rest).mul(ex_to<numeric>(i1->coeff).mul(ex_to<numeric>(i2->coeff)));
                                                else
@@ -905,25 +908,38 @@ ex mul::expand(unsigned options) const
                                last_expanded = cit->rest;
                        }
                } else {
+                       if (is_exactly_a<add>(cit->rest))
+                               non_adds_has_sums = true;
                        non_adds.push_back(*cit);
                }
                ++cit;
        }
        if (expanded_seqp)
                delete expanded_seqp;
-       
+
        // Now the only remaining thing to do is to multiply the factors which
        // were not sums into the "last_expanded" sum
        if (is_exactly_a<add>(last_expanded)) {
                const add & finaladd = ex_to<add>(last_expanded);
-               exvector distrseq;
+
                size_t n = finaladd.nops();
+               exvector distrseq;
                distrseq.reserve(n);
+
                for (size_t i=0; i<n; ++i) {
                        epvector factors = non_adds;
-                       factors.push_back(split_ex_to_pair(finaladd.op(i)));
-                       distrseq.push_back((new mul(factors, overall_coeff))->
-                                           setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)));
+                       expair new_factor = split_ex_to_pair(finaladd.op(i).expand());
+                       factors.push_back(new_factor);
+
+                       const mul & term = static_cast<const mul &>((new mul(factors, overall_coeff))->setflag(status_flags::dynallocated));
+
+                       // The new term may have sums in it if e.g. a sqrt() of a sum in
+                       // the non_adds meets a sqrt() of a sum in the factor from
+                       // last_expanded. In this case we should re-expand the term.
+                       if (non_adds_has_sums || is_exactly_a<add>(new_factor.rest))
+                               distrseq.push_back(ex(term).expand());
+                       else
+                               distrseq.push_back(term.setflag(options == 0 ? status_flags::expanded : 0));
                }
                return ((new add(distrseq))->
                        setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)));

diff --git a/ginac/power.cpp b/ginac/power.cpp
index fb67a504de19f3d4ac162b3935cf1a4246f45105..b286db365ad8bac7bb6e54e2f57da929031737f3 100644 (file)
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -837,6 +837,8 @@ ex power::expand_mul(const mul & m, const numeric & n) const
 
        epvector distrseq;
        distrseq.reserve(m.seq.size());
+       bool need_reexpand = false;
+
        epvector::const_iterator last = m.seq.end();
        epvector::const_iterator cit = m.seq.begin();
        while (cit!=last) {
@@ -845,11 +847,22 @@ ex power::expand_mul(const mul & m, const numeric & n) const
                } else {
                        // it is safe not to call mul::combine_pair_with_coeff_to_pair()
                        // since n is an integer
-                       distrseq.push_back(expair(cit->rest, ex_to<numeric>(cit->coeff).mul(n)));
+                       numeric new_coeff = ex_to<numeric>(cit->coeff).mul(n);
+                       if (is_exactly_a<add>(cit->rest) && new_coeff.is_pos_integer()) {
+                               // this happens when e.g. (a+b)^(1/2) gets squared and
+                               // the resulting product needs to be reexpanded
+                               need_reexpand = true;
+                       }
+                       distrseq.push_back(expair(cit->rest, new_coeff));
                }
                ++cit;
        }
-       return (new mul(distrseq, ex_to<numeric>(m.overall_coeff).power_dyn(n)))->setflag(status_flags::dynallocated);
+
+       const mul & result = static_cast<const mul &>((new mul(distrseq, ex_to<numeric>(m.overall_coeff).power_dyn(n)))->setflag(status_flags::dynallocated));
+       if (need_reexpand)
+               return ex(result).expand();
+       else
+               return result.setflag(status_flags::expanded);
 }
 
 } // namespace GiNaC