From: Richard Kreckel <kreckel@ginac.de>
Date: Wed, 25 Nov 2015 10:28:10 +0000 (+0100)
Subject: In power::expand_add(), don't reserve excess monomial sizes.
X-Git-Tag: release_1-7-0~7^2~48
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=b236efe23093bf2c4b5e7702b60d45505724d915

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.
---

diff --git a/ginac/power.cpp b/ginac/power.cpp
index e8f599d8..1a452e3f 100644
--- 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());
 	}