+ return composition;
+ }
+ bool next()
+ {
+ // This ugly contortion is needed because the original coolmulti
+ // algorithm requires code duplication of the payload procedure,
+ // one before the loop and one inside it.
+ if (trivial || atend)
+ return false;
+ cmgen.next_permutation();
+ atend = cmgen.finished();
+ return true;
+ }
+};
+
+/** Helper function to compute the multinomial coefficient n!/(p1!*p2!*...*pk!)
+ * where n = p1+p2+...+pk, i.e. p is a partition of n.
+ */
+const numeric
+multinomial_coefficient(const std::vector<int> & p)
+{
+ numeric n = 0, d = 1;
+ std::vector<int>::const_iterator it = p.begin(), itend = p.end();
+ while (it != itend) {
+ n += numeric(*it);
+ d *= factorial(numeric(*it));
+ ++it;
+ }
+ return factorial(numeric(n)) / d;
+}
+
+} // anonymous namespace
+
+/** expand a^n where a is an add and n is a positive integer.
+ * @see power::expand */
+ex power::expand_add(const add & a, int n, unsigned options) const
+{
+ // The special case power(+(x,...y;x),2) can be optimized better.
+ if (n==2)
+ return expand_add_2(a, options);
+
+ // method:
+ //
+ // Consider base as the sum of all symbolic terms and the overall numeric
+ // coefficient and apply the binomial theorem:
+ // S = power(+(x,...,z;c),n)
+ // = power(+(+(x,...,z;0);c),n)
+ // = sum(binomial(n,k)*power(+(x,...,z;0),k)*c^(n-k), k=1..n) + c^n
+ // Then, apply the multinomial theorem to expand all power(+(x,...,z;0),k):
+ // The multinomial theorem is computed by an outer loop over all
+ // partitions of the exponent and an inner loop over all compositions of
+ // that partition. This method makes the expansion a combinatorial
+ // problem and allows us to directly construct the expanded sum and also
+ // to re-use the multinomial coefficients (since they depend only on the
+ // partition, not on the composition).
+ //
+ // multinomial power(+(x,y,z;0),3) example:
+ // partition : compositions : multinomial coefficient
+ // [0,0,3] : [3,0,0],[0,3,0],[0,0,3] : 3!/(3!*0!*0!) = 1
+ // [0,1,2] : [2,1,0],[1,2,0],[2,0,1],... : 3!/(2!*1!*0!) = 3
+ // [1,1,1] : [1,1,1] : 3!/(1!*1!*1!) = 6
+ // => (x + y + z)^3 =
+ // x^3 + y^3 + z^3
+ // + 3*x^2*y + 3*x*y^2 + 3*y^2*z + 3*y*z^2 + 3*x*z^2 + 3*x^2*z
+ // + 6*x*y*z
+ //
+ // multinomial power(+(x,y,z;0),4) example:
+ // partition : compositions : multinomial coefficient
+ // [0,0,4] : [4,0,0],[0,4,0],[0,0,4] : 4!/(4!*0!*0!) = 1
+ // [0,1,3] : [3,1,0],[1,3,0],[3,0,1],... : 4!/(3!*1!*0!) = 4
+ // [0,2,2] : [2,2,0],[2,0,2],[0,2,2] : 4!/(2!*2!*0!) = 6
+ // [1,1,2] : [2,1,1],[1,2,1],[1,1,2] : 4!/(2!*1!*1!) = 12
+ // (no [1,1,1,1] partition since it has too many parts)
+ // => (x + y + z)^4 =
+ // x^4 + y^4 + z^4
+ // + 4*x^3*y + 4*x*y^3 + 4*y^3*z + 4*y*z^3 + 4*x*z^3 + 4*x^3*z
+ // + 6*x^2*y^2 + 6*y^2*z^2 + 6*x^2*z^2
+ // + 12*x^2*y*z + 12*x*y^2*z + 12*x*y*z^2
+ //
+ // Summary:
+ // r = 0
+ // for k from 0 to n:
+ // f = c^(n-k)*binomial(n,k)
+ // for p in all partitions of n with m parts (including zero parts):
+ // h = f * multinomial coefficient of p
+ // for c in all compositions of p:
+ // t = 1
+ // for e in all elements of c:
+ // t = t * a[e]^e
+ // r = r + h*t
+ // return r
+
+ epvector result;
+ // The number of terms will be the number of combinatorial compositions,
+ // i.e. the number of unordered arrangements of m nonnegative integers
+ // which sum up to n. It is frequently written as C_n(m) and directly
+ // related with binomial coefficients: binomial(n+m-1,m-1).
+ size_t result_size = binomial(numeric(n+a.nops()-1), numeric(a.nops()-1)).to_int();
+ if (!a.overall_coeff.is_zero()) {
+ // the result's overall_coeff is one of the terms
+ --result_size;
+ }
+ result.reserve(result_size);
+
+ // Iterate over all terms in binomial expansion of
+ // S = power(+(x,...,z;c),n)
+ // = sum(binomial(n,k)*power(+(x,...,z;0),k)*c^(n-k), k=1..n) + c^n
+ for (int k = 1; k <= n; ++k) {
+ numeric binomial_coefficient; // binomial(n,k)*c^(n-k)
+ if (a.overall_coeff.is_zero()) {
+ // degenerate case with zero overall_coeff:
+ // apply multinomial theorem directly to power(+(x,...z;0),n)
+ binomial_coefficient = 1;
+ if (k < n) {
+ continue;
+ }
+ } else {
+ binomial_coefficient = binomial(numeric(n), numeric(k)) * pow(ex_to<numeric>(a.overall_coeff), numeric(n-k));
+ }
+
+ // Multinomial expansion of power(+(x,...,z;0),k)*c^(n-k):
+ // Iterate over all partitions of k with exactly as many parts as
+ // there are symbolic terms in the basis (including zero parts).
+ 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 = count_if(partition.begin(), partition.end(), bind2nd(std::greater<int>(), 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 monomial;
+ monomial.reserve(msize);
+ numeric factor = coeff;
+ for (unsigned i = 0; i < exponent.size(); ++i) {
+ const ex & r = a.seq[i].rest;
+ GINAC_ASSERT(!is_exactly_a<add>(r));
+ GINAC_ASSERT(!is_exactly_a<power>(r) ||
+ !is_exactly_a<numeric>(ex_to<power>(r).exponent) ||
+ !ex_to<numeric>(ex_to<power>(r).exponent).is_pos_integer() ||
+ !is_exactly_a<add>(ex_to<power>(r).basis) ||
+ !is_exactly_a<mul>(ex_to<power>(r).basis) ||
+ !is_exactly_a<power>(ex_to<power>(r).basis));
+ GINAC_ASSERT(is_exactly_a<numeric>(a.seq[i].coeff));
+ const numeric & c = ex_to<numeric>(a.seq[i].coeff);
+ if (exponent[i] == 0) {
+ // optimize away
+ } else if (exponent[i] == 1) {
+ // optimized
+ monomial.push_back(r);
+ if (c != *_num1_p)
+ factor = factor.mul(c);
+ } else { // general case exponent[i] > 1
+ 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(monomial).expand(options), factor));
+ } while (compositions.next());
+ } while (partitions.next());
+ }
+
+ GINAC_ASSERT(result.size() == result_size);
+
+ if (a.overall_coeff.is_zero()) {
+ return (new add(result))->setflag(status_flags::dynallocated |
+ status_flags::expanded);
+ } else {
+ return (new add(result, ex_to<numeric>(a.overall_coeff).power(n)))->setflag(status_flags::dynallocated |
+ status_flags::expanded);