From 2639812c0ba4e1f9620660bbba1f12bf5b865e29 Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Wed, 25 Nov 2015 11:28:10 +0100 Subject: [PATCH 1/1] In power::expand_add(), don't reserve excess monomial sizes. There is no need to reserve n terms in each of the monomials of the result of power(+(x,y,z...;0),n): We can compute it exactly as the number of nonzero exponents in the multinomial expansion. The good thing is that this counting is the same for each composition of a partition, so it can be hoisted out of the loop over compositions. --- ginac/power.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/ginac/power.cpp b/ginac/power.cpp index 7298d569..6fdc3fcd 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -42,6 +42,7 @@ #include #include #include +#include namespace GiNaC { @@ -1204,14 +1205,16 @@ ex power::expand_add(const add & a, int n, unsigned options) const partition_generator partitions(k, a.seq.size()); do { const std::vector& partition = partitions.current(); + // All monomials of this partition have the same number of terms and the same coefficient. + const unsigned msize = count_if(partition.begin(), partition.end(), bind2nd(std::greater(), 0)); const numeric coeff = multinomial_coefficient(partition) * binomial_coefficient; // Iterate over all compositions of the current partition. composition_generator compositions(partition); do { const std::vector& exponent = compositions.current(); - exvector term; - term.reserve(n); + exvector monomial; + monomial.reserve(msize); numeric factor = coeff; for (unsigned i = 0; i < exponent.size(); ++i) { const ex & r = a.seq[i].rest; @@ -1228,16 +1231,16 @@ ex power::expand_add(const add & a, int n, unsigned options) const // optimize away } else if (exponent[i] == 1) { // optimized - term.push_back(r); + monomial.push_back(r); if (c != *_num1_p) factor = factor.mul(c); } else { // general case exponent[i] > 1 - term.push_back((new power(r, exponent[i]))->setflag(status_flags::dynallocated)); + monomial.push_back((new power(r, exponent[i]))->setflag(status_flags::dynallocated)); if (c != *_num1_p) factor = factor.mul(c.power(exponent[i])); } } - result.push_back(a.combine_ex_with_coeff_to_pair(mul(term).expand(options), factor)); + result.push_back(a.combine_ex_with_coeff_to_pair(mul(monomial).expand(options), factor)); } while (compositions.next()); } while (partitions.next()); } -- 2.44.0