+/** expand a^n where a is an add and n is a positive integer.
+ * @see power::expand */
+ex power::expand_add(const add & a, long n, unsigned options)
+{
+ // 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_long();
+ 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));