In power::expand_add(), don't reserve excess monomial sizes.

[ginac.git] / ginac / power.cpp

diff --git a/ginac/power.cpp b/ginac/power.cpp
index e8f599d88abe08078f15faf2fb3db9a034347178..1a452e3fcc91f36a0e78ec767e6d3cf7efc410db 100644 (file)
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -42,6 +42,7 @@
 #include <limits>
 #include <stdexcept>
 #include <vector>
+#include <algorithm>
 
 namespace GiNaC {
 
@@ -1197,14 +1198,16 @@ ex power::expand_add(const add & a, long n, unsigned options)
                partition_generator partitions(k, a.seq.size());
                do {
                        const std::vector<int>& partition = partitions.current();
+                       // All monomials of this partition have the same number of terms and the same coefficient.
+                       const unsigned msize = std::count_if(partition.begin(), partition.end(), [](int i) { return i > 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<int>& 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;
@@ -1221,16 +1224,16 @@ ex power::expand_add(const add & a, long n, unsigned options)
                                                // 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());
        }