return result;
}
+// Expanding products containing powers of sums could return results that
+// were not fully expanded. Fixed on Dec 10, 2003.
+static unsigned exam_paranoia16()
+{
+ unsigned result = 0;
+ symbol a("a"), b("b"), c("c"), d("d");
+ ex e1, e2, e3;
+
+ e1 = pow(1+a*sqrt(b+c), 2);
+ e2 = e1.expand();
+
+ if (e2.has(pow(a, 2)*(b+c))) {
+ clog << "expand(" << e1 << ") didn't fully expand\n";
+ ++result;
+ }
+
+ e1 = (d*sqrt(a+b)+a*sqrt(c+d))*(b*sqrt(a+b)+a*sqrt(c+d));
+ e2 = e1.expand();
+
+ if (e2.has(pow(a, 2)*(c+d))) {
+ clog << "expand(" << e1 << ") didn't fully expand\n";
+ ++result;
+ }
+
+ e1 = (a+sqrt(b+c))*sqrt(b+c)*(d+sqrt(b+c));
+ e2 = e1.expand();
+
+ if (e2.has(a*(b+c))) {
+ clog << "expand(" << e1 << ") didn't fully expand\n";
+ ++result;
+ }
+
+ return result;
+}
+
unsigned exam_paranoia()
{
unsigned result = 0;
result += exam_paranoia13(); cout << '.' << flush;
result += exam_paranoia14(); cout << '.' << flush;
result += exam_paranoia15(); cout << '.' << flush;
+ result += exam_paranoia16(); cout << '.' << flush;
if (!result) {
cout << " passed " << endl;
// 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) &&
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
last_expanded = cit->rest;
}
} else {
+ if (is_exactly_a<add>(cit->rest))
+ non_adds_has_sums = true;
non_adds.push_back(*cit);
}
++cit;
}
-
+
// 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)));
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) {
} 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