In power::expand_add(), don't reserve excess monomial sizes.
authorRichard Kreckel <kreckel@ginac.de>
Wed, 25 Nov 2015 10:28:10 +0000 (11:28 +0100)
committerRichard Kreckel <kreckel@ginac.de>
Wed, 25 Nov 2015 11:40:00 +0000 (12:40 +0100)
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

index e8f599d..1a452e3 100644 (file)
@@ -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());
        }