]> www.ginac.de Git - ginac.git/blobdiff - ginac/power.cpp
Remove info_flags::algebraic.
[ginac.git] / ginac / power.cpp
index e8f599d88abe08078f15faf2fb3db9a034347178..8e62af8b5e4f98c6b5e59f786f8ddb05efc0bd36 100644 (file)
@@ -42,6 +42,7 @@
 #include <limits>
 #include <stdexcept>
 #include <vector>
+#include <algorithm>
 
 namespace GiNaC {
 
@@ -233,9 +234,6 @@ bool power::info(unsigned inf) const
                case info_flags::rational_function:
                        return exponent.info(info_flags::integer) &&
                               basis.info(inf);
-               case info_flags::algebraic:
-                       return !exponent.info(info_flags::integer) ||
-                              basis.info(inf);
                case info_flags::expanded:
                        return (flags & status_flags::expanded);
                case info_flags::positive:
@@ -281,7 +279,7 @@ ex power::map(map_function & f) const
 
        if (!are_ex_trivially_equal(basis, mapped_basis)
         || !are_ex_trivially_equal(exponent, mapped_exponent))
-               return (new power(mapped_basis, mapped_exponent))->setflag(status_flags::dynallocated);
+               return dynallocate<power>(mapped_basis, mapped_exponent);
        else
                return *this;
 }
@@ -372,42 +370,36 @@ ex power::coeff(const ex & s, int n) const
  *  - ^(*(x,y,z),c) -> *(x^c,y^c,z^c)  (if c integer)
  *  - ^(*(x,c1),c2) -> ^(x,c2)*c1^c2  (c1>0)
  *  - ^(*(x,c1),c2) -> ^(-x,c2)*c1^c2  (c1<0)
- *
- *  @param level cut-off in recursive evaluation */
-ex power::eval(int level) const
+ */
+ex power::eval() const
 {
-       if ((level==1) && (flags & status_flags::evaluated))
+       if (flags & status_flags::evaluated)
                return *this;
-       else if (level == -max_recursion_level)
-               throw(std::runtime_error("max recursion level reached"));
-       
-       const ex & ebasis    = level==1 ? basis    : basis.eval(level-1);
-       const ex & eexponent = level==1 ? exponent : exponent.eval(level-1);
-       
+
        const numeric *num_basis = nullptr;
        const numeric *num_exponent = nullptr;
-       
-       if (is_exactly_a<numeric>(ebasis)) {
-               num_basis = &ex_to<numeric>(ebasis);
+
+       if (is_exactly_a<numeric>(basis)) {
+               num_basis = &ex_to<numeric>(basis);
        }
-       if (is_exactly_a<numeric>(eexponent)) {
-               num_exponent = &ex_to<numeric>(eexponent);
+       if (is_exactly_a<numeric>(exponent)) {
+               num_exponent = &ex_to<numeric>(exponent);
        }
        
        // ^(x,0) -> 1  (0^0 also handled here)
-       if (eexponent.is_zero()) {
-               if (ebasis.is_zero())
+       if (exponent.is_zero()) {
+               if (basis.is_zero())
                        throw (std::domain_error("power::eval(): pow(0,0) is undefined"));
                else
                        return _ex1;
        }
        
        // ^(x,1) -> x
-       if (eexponent.is_equal(_ex1))
-               return ebasis;
+       if (exponent.is_equal(_ex1))
+               return basis;
 
        // ^(0,c1) -> 0 or exception  (depending on real value of c1)
-       if ( ebasis.is_zero() && num_exponent ) {
+       if (basis.is_zero() && num_exponent) {
                if ((num_exponent->real()).is_zero())
                        throw (std::domain_error("power::eval(): pow(0,I) is undefined"));
                else if ((num_exponent->real()).is_negative())
@@ -417,16 +409,16 @@ ex power::eval(int level) const
        }
 
        // ^(1,x) -> 1
-       if (ebasis.is_equal(_ex1))
+       if (basis.is_equal(_ex1))
                return _ex1;
 
        // power of a function calculated by separate rules defined for this function
-       if (is_exactly_a<function>(ebasis))
-               return ex_to<function>(ebasis).power(eexponent);
+       if (is_exactly_a<function>(basis))
+               return ex_to<function>(basis).power(exponent);
 
        // Turn (x^c)^d into x^(c*d) in the case that x is positive and c is real.
-       if (is_exactly_a<power>(ebasis) && ebasis.op(0).info(info_flags::positive) && ebasis.op(1).info(info_flags::real))
-               return power(ebasis.op(0), ebasis.op(1) * eexponent);
+       if (is_exactly_a<power>(basis) && basis.op(0).info(info_flags::positive) && basis.op(1).info(info_flags::real))
+               return dynallocate<power>(basis.op(0), basis.op(1) * exponent);
 
        if ( num_exponent ) {
 
@@ -437,9 +429,7 @@ ex power::eval(int level) const
                        const bool exponent_is_crational = num_exponent->is_crational();
                        if (!basis_is_crational || !exponent_is_crational) {
                                // return a plain float
-                               return (new numeric(num_basis->power(*num_exponent)))->setflag(status_flags::dynallocated |
-                                                                                              status_flags::evaluated |
-                                                                                              status_flags::expanded);
+                               return dynallocate<numeric>(num_basis->power(*num_exponent));
                        }
 
                        const numeric res = num_basis->power(*num_exponent);
@@ -469,9 +459,9 @@ ex power::eval(int level) const
                                                const numeric res_bnum = bnum.power(*num_exponent);
                                                const numeric res_bden = bden.power(*num_exponent);
                                                if (res_bnum.is_integer())
-                                                       return (new mul(power(bden,-*num_exponent),res_bnum))->setflag(status_flags::dynallocated | status_flags::evaluated);
+                                                       return dynallocate<mul>(dynallocate<power>(bden,-*num_exponent),res_bnum).setflag(status_flags::evaluated);
                                                if (res_bden.is_integer())
-                                                       return (new mul(power(bnum,*num_exponent),res_bden.inverse()))->setflag(status_flags::dynallocated | status_flags::evaluated);
+                                                       return dynallocate<mul>(dynallocate<power>(bnum,*num_exponent),res_bden.inverse()).setflag(status_flags::evaluated);
                                        }
                                        return this->hold();
                                } else {
@@ -479,8 +469,7 @@ ex power::eval(int level) const
                                        // because otherwise we'll end up with something like
                                        //    (7/8)^(4/3)  ->  7/8*(1/2*7^(1/3))
                                        // instead of 7/16*7^(1/3).
-                                       ex prod = power(*num_basis,r.div(m));
-                                       return prod*power(*num_basis,q);
+                                       return pow(basis, r.div(m)) * pow(basis, q);
                                }
                        }
                }
@@ -488,8 +477,8 @@ ex power::eval(int level) const
                // ^(^(x,c1),c2) -> ^(x,c1*c2)
                // (c1, c2 numeric(), c2 integer or -1 < c1 <= 1 or (c1=-1 and c2>0),
                // case c1==1 should not happen, see below!)
-               if (is_exactly_a<power>(ebasis)) {
-                       const power & sub_power = ex_to<power>(ebasis);
+               if (is_exactly_a<power>(basis)) {
+                       const power & sub_power = ex_to<power>(basis);
                        const ex & sub_basis = sub_power.basis;
                        const ex & sub_exponent = sub_power.exponent;
                        if (is_exactly_a<numeric>(sub_exponent)) {
@@ -497,21 +486,21 @@ ex power::eval(int level) const
                                GINAC_ASSERT(num_sub_exponent!=numeric(1));
                                if (num_exponent->is_integer() || (abs(num_sub_exponent) - (*_num1_p)).is_negative() ||
                                    (num_sub_exponent == *_num_1_p && num_exponent->is_positive())) {
-                                       return power(sub_basis,num_sub_exponent.mul(*num_exponent));
+                                       return dynallocate<power>(sub_basis, num_sub_exponent.mul(*num_exponent));
                                }
                        }
                }
        
                // ^(*(x,y,z),c1) -> *(x^c1,y^c1,z^c1) (c1 integer)
-               if (num_exponent->is_integer() && is_exactly_a<mul>(ebasis)) {
-                       return expand_mul(ex_to<mul>(ebasis), *num_exponent, 0);
+               if (num_exponent->is_integer() && is_exactly_a<mul>(basis)) {
+                       return expand_mul(ex_to<mul>(basis), *num_exponent, false);
                }
 
                // (2*x + 6*y)^(-4) -> 1/16*(x + 3*y)^(-4)
-               if (num_exponent->is_integer() && is_exactly_a<add>(ebasis)) {
-                       numeric icont = ebasis.integer_content();
+               if (num_exponent->is_integer() && is_exactly_a<add>(basis)) {
+                       numeric icont = basis.integer_content();
                        const numeric lead_coeff = 
-                               ex_to<numeric>(ex_to<add>(ebasis).seq.begin()->coeff).div(icont);
+                               ex_to<numeric>(ex_to<add>(basis).seq.begin()->coeff).div(icont);
 
                        const bool canonicalizable = lead_coeff.is_integer();
                        const bool unit_normal = lead_coeff.is_pos_integer();
@@ -519,48 +508,43 @@ ex power::eval(int level) const
                                icont = icont.mul(*_num_1_p);
                        
                        if (canonicalizable && (icont != *_num1_p)) {
-                               const add& addref = ex_to<add>(ebasis);
-                               add* addp = new add(addref);
-                               addp->setflag(status_flags::dynallocated);
-                               addp->clearflag(status_flags::hash_calculated);
-                               addp->overall_coeff = ex_to<numeric>(addp->overall_coeff).div_dyn(icont);
-                               for (auto & i : addp->seq)
+                               const add& addref = ex_to<add>(basis);
+                               add & addp = dynallocate<add>(addref);
+                               addp.clearflag(status_flags::hash_calculated);
+                               addp.overall_coeff = ex_to<numeric>(addp.overall_coeff).div_dyn(icont);
+                               for (auto & i : addp.seq)
                                        i.coeff = ex_to<numeric>(i.coeff).div_dyn(icont);
 
                                const numeric c = icont.power(*num_exponent);
                                if (likely(c != *_num1_p))
-                                       return (new mul(power(*addp, *num_exponent), c))->setflag(status_flags::dynallocated);
+                                       return dynallocate<mul>(dynallocate<power>(addp, *num_exponent), c);
                                else
-                                       return power(*addp, *num_exponent);
+                                       return dynallocate<power>(addp, *num_exponent);
                        }
                }
 
                // ^(*(...,x;c1),c2) -> *(^(*(...,x;1),c2),c1^c2)  (c1, c2 numeric(), c1>0)
                // ^(*(...,x;c1),c2) -> *(^(*(...,x;-1),c2),(-c1)^c2)  (c1, c2 numeric(), c1<0)
-               if (is_exactly_a<mul>(ebasis)) {
+               if (is_exactly_a<mul>(basis)) {
                        GINAC_ASSERT(!num_exponent->is_integer()); // should have been handled above
-                       const mul & mulref = ex_to<mul>(ebasis);
+                       const mul & mulref = ex_to<mul>(basis);
                        if (!mulref.overall_coeff.is_equal(_ex1)) {
                                const numeric & num_coeff = ex_to<numeric>(mulref.overall_coeff);
                                if (num_coeff.is_real()) {
                                        if (num_coeff.is_positive()) {
-                                               mul *mulp = new mul(mulref);
-                                               mulp->overall_coeff = _ex1;
-                                               mulp->setflag(status_flags::dynallocated);
-                                               mulp->clearflag(status_flags::evaluated);
-                                               mulp->clearflag(status_flags::hash_calculated);
-                                               return (new mul(power(*mulp,exponent),
-                                                               power(num_coeff,*num_exponent)))->setflag(status_flags::dynallocated);
+                                               mul & mulp = dynallocate<mul>(mulref);
+                                               mulp.overall_coeff = _ex1;
+                                               mulp.clearflag(status_flags::evaluated | status_flags::hash_calculated);
+                                               return dynallocate<mul>(dynallocate<power>(mulp, exponent),
+                                                                       dynallocate<power>(num_coeff, *num_exponent));
                                        } else {
                                                GINAC_ASSERT(num_coeff.compare(*_num0_p)<0);
                                                if (!num_coeff.is_equal(*_num_1_p)) {
-                                                       mul *mulp = new mul(mulref);
-                                                       mulp->overall_coeff = _ex_1;
-                                                       mulp->setflag(status_flags::dynallocated);
-                                                       mulp->clearflag(status_flags::evaluated);
-                                                       mulp->clearflag(status_flags::hash_calculated);
-                                                       return (new mul(power(*mulp,exponent),
-                                                                       power(abs(num_coeff),*num_exponent)))->setflag(status_flags::dynallocated);
+                                                       mul & mulp = dynallocate<mul>(mulref);
+                                                       mulp.overall_coeff = _ex_1;
+                                                       mulp.clearflag(status_flags::evaluated | status_flags::hash_calculated);
+                                                       return dynallocate<mul>(dynallocate<power>(mulp, exponent),
+                                                                               dynallocate<power>(abs(num_coeff), *num_exponent));
                                                }
                                        }
                                }
@@ -569,18 +553,13 @@ ex power::eval(int level) const
 
                // ^(nc,c1) -> ncmul(nc,nc,...) (c1 positive integer, unless nc is a matrix)
                if (num_exponent->is_pos_integer() &&
-                   ebasis.return_type() != return_types::commutative &&
-                   !is_a<matrix>(ebasis)) {
-                       return ncmul(exvector(num_exponent->to_int(), ebasis));
+                   basis.return_type() != return_types::commutative &&
+                   !is_a<matrix>(basis)) {
+                       return ncmul(exvector(num_exponent->to_int(), basis));
                }
        }
-       
-       if (are_ex_trivially_equal(ebasis,basis) &&
-           are_ex_trivially_equal(eexponent,exponent)) {
-               return this->hold();
-       }
-       return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated |
-                                                      status_flags::evaluated);
+
+       return this->hold();
 }
 
 ex power::evalf(int level) const
@@ -601,7 +580,7 @@ ex power::evalf(int level) const
                        eexponent = exponent;
        }
 
-       return power(ebasis,eexponent);
+       return dynallocate<power>(ebasis, eexponent);
 }
 
 ex power::evalm() const
@@ -610,10 +589,10 @@ ex power::evalm() const
        const ex eexponent = exponent.evalm();
        if (is_a<matrix>(ebasis)) {
                if (is_exactly_a<numeric>(eexponent)) {
-                       return (new matrix(ex_to<matrix>(ebasis).pow(eexponent)))->setflag(status_flags::dynallocated);
+                       return dynallocate<matrix>(ex_to<matrix>(ebasis).pow(eexponent));
                }
        }
-       return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated);
+       return dynallocate<power>(ebasis, eexponent);
 }
 
 bool power::has(const ex & other, unsigned options) const
@@ -659,7 +638,7 @@ ex power::subs(const exmap & m, unsigned options) const
                if (tryfactsubs(*this, it.first, nummatches, repls)) {
                        ex anum = it.second.subs(repls, subs_options::no_pattern);
                        ex aden = it.first.subs(repls, subs_options::no_pattern);
-                       ex result = (*this)*power(anum/aden, nummatches);
+                       ex result = (*this) * pow(anum/aden, nummatches);
                        return (ex_to<basic>(result)).subs_one_level(m, options);
                }
        }
@@ -681,14 +660,14 @@ ex power::conjugate() const
                if (are_ex_trivially_equal(exponent, newexponent)) {
                        return *this;
                }
-               return (new power(basis, newexponent))->setflag(status_flags::dynallocated);
+               return dynallocate<power>(basis, newexponent);
        }
        if (exponent.info(info_flags::integer)) {
                ex newbasis = basis.conjugate();
                if (are_ex_trivially_equal(basis, newbasis)) {
                        return *this;
                }
-               return (new power(newbasis, exponent))->setflag(status_flags::dynallocated);
+               return dynallocate<power>(newbasis, exponent);
        }
        return conjugate_function(*this).hold();
 }
@@ -708,12 +687,12 @@ ex power::real_part() const
                // Re((a+I*b)^c)  w/  c ∈ ℤ
                long N = ex_to<numeric>(c).to_long();
                // Use real terms in Binomial expansion to construct
-               // Re(expand(power(a+I*b, N))).
+               // Re(expand(pow(a+I*b, N))).
                long NN = N > 0 ? N : -N;
-               ex numer = N > 0 ? _ex1 : power(power(a,2) + power(b,2), NN);
+               ex numer = N > 0 ? _ex1 : pow(pow(a,2) + pow(b,2), NN);
                ex result = 0;
                for (long n = 0; n <= NN; n += 2) {
-                       ex term = binomial(NN, n) * power(a, NN-n) * power(b, n) / numer;
+                       ex term = binomial(NN, n) * pow(a, NN-n) * pow(b, n) / numer;
                        if (n % 4 == 0) {
                                result += term;  // sign: I^n w/ n == 4*m
                        } else {
@@ -725,11 +704,12 @@ ex power::real_part() const
 
        // Re((a+I*b)^(c+I*d))
        const ex d = exponent.imag_part();
-       return power(abs(basis),c)*exp(-d*atan2(b,a))*cos(c*atan2(b,a)+d*log(abs(basis)));
+       return pow(abs(basis),c) * exp(-d*atan2(b,a)) * cos(c*atan2(b,a)+d*log(abs(basis)));
 }
 
 ex power::imag_part() const
 {
+       // basis == a+I*b, exponent == c+I*d
        const ex a = basis.real_part();
        const ex c = exponent.real_part();
        if (basis.is_equal(a) && exponent.is_equal(c)) {
@@ -742,13 +722,13 @@ ex power::imag_part() const
                // Im((a+I*b)^c)  w/  c ∈ ℤ
                long N = ex_to<numeric>(c).to_long();
                // Use imaginary terms in Binomial expansion to construct
-               // Im(expand(power(a+I*b, N))).
+               // Im(expand(pow(a+I*b, N))).
                long p = N > 0 ? 1 : 3;  // modulus for positive sign
                long NN = N > 0 ? N : -N;
-               ex numer = N > 0 ? _ex1 : power(power(a,2) + power(b,2), NN);
+               ex numer = N > 0 ? _ex1 : pow(pow(a,2) + pow(b,2), NN);
                ex result = 0;
                for (long n = 1; n <= NN; n += 2) {
-                       ex term = binomial(NN, n) * power(a, NN-n) * power(b, n) / numer;
+                       ex term = binomial(NN, n) * pow(a, NN-n) * pow(b, n) / numer;
                        if (n % 4 == p) {
                                result += term;  // sign: I^n w/ n == 4*m+p
                        } else {
@@ -760,7 +740,7 @@ ex power::imag_part() const
 
        // Im((a+I*b)^(c+I*d))
        const ex d = exponent.imag_part();
-       return power(abs(basis),c)*exp(-d*atan2(b,a))*sin(c*atan2(b,a)+d*log(abs(basis)));
+       return pow(abs(basis),c) * exp(-d*atan2(b,a)) * sin(c*atan2(b,a)+d*log(abs(basis)));
 }
 
 // protected
@@ -771,16 +751,11 @@ ex power::derivative(const symbol & s) const
 {
        if (is_a<numeric>(exponent)) {
                // D(b^r) = r * b^(r-1) * D(b) (faster than the formula below)
-               epvector newseq;
-               newseq.reserve(2);
-               newseq.push_back(expair(basis, exponent - _ex1));
-               newseq.push_back(expair(basis.diff(s), _ex1));
-               return mul(std::move(newseq), exponent);
+               const epvector newseq = {expair(basis, exponent - _ex1), expair(basis.diff(s), _ex1)};
+               return dynallocate<mul>(std::move(newseq), exponent);
        } else {
                // D(b^e) = b^e * (D(e)*ln(b) + e*D(b)/b)
-               return mul(*this,
-                          add(mul(exponent.diff(s), log(basis)),
-                          mul(mul(exponent, basis.diff(s)), power(basis, _ex_1))));
+               return *this * (exponent.diff(s)*log(basis) + exponent*basis.diff(s)*pow(basis, _ex_1));
        }
 }
 
@@ -839,9 +814,9 @@ ex power::expand(unsigned options) const
                // take care on the numeric coefficient
                ex coeff=(possign? _ex1 : _ex_1);
                if (m.overall_coeff.info(info_flags::positive) && m.overall_coeff != _ex1)
-                       prodseq.push_back(power(m.overall_coeff, exponent));
+                       prodseq.push_back(pow(m.overall_coeff, exponent));
                else if (m.overall_coeff.info(info_flags::negative) && m.overall_coeff != _ex_1)
-                       prodseq.push_back(power(-m.overall_coeff, exponent));
+                       prodseq.push_back(pow(-m.overall_coeff, exponent));
                else
                        coeff *= m.overall_coeff;
 
@@ -849,9 +824,9 @@ ex power::expand(unsigned options) const
                // In either case we set a flag to avoid the second run on a part
                // which does not have positive/negative terms.
                if (prodseq.size() > 0) {
-                       ex newbasis = coeff*mul(std::move(powseq));
+                       ex newbasis = dynallocate<mul>(std::move(powseq), coeff);
                        ex_to<basic>(newbasis).setflag(status_flags::purely_indefinite);
-                       return ((new mul(std::move(prodseq)))->setflag(status_flags::dynallocated)*(new power(newbasis, exponent))->setflag(status_flags::dynallocated).expand(options)).expand(options);
+                       return dynallocate<mul>(std::move(prodseq)) * pow(newbasis, exponent);
                } else
                        ex_to<basic>(basis).setflag(status_flags::purely_indefinite);
        }
@@ -865,7 +840,7 @@ ex power::expand(unsigned options) const
                exvector distrseq;
                distrseq.reserve(a.seq.size() + 1);
                for (auto & cit : a.seq) {
-                       distrseq.push_back(power(expanded_basis, a.recombine_pair_to_ex(cit)));
+                       distrseq.push_back(pow(expanded_basis, a.recombine_pair_to_ex(cit)));
                }
                
                // Make sure that e.g. (x+y)^(2+a) expands the (x+y)^2 factor
@@ -875,12 +850,12 @@ ex power::expand(unsigned options) const
                        if (int_exponent > 0 && is_exactly_a<add>(expanded_basis))
                                distrseq.push_back(expand_add(ex_to<add>(expanded_basis), int_exponent, options));
                        else
-                               distrseq.push_back(power(expanded_basis, a.overall_coeff));
+                               distrseq.push_back(pow(expanded_basis, a.overall_coeff));
                } else
-                       distrseq.push_back(power(expanded_basis, a.overall_coeff));
+                       distrseq.push_back(pow(expanded_basis, a.overall_coeff));
                
                // Make sure that e.g. (x+y)^(1+a) -> x*(x+y)^a + y*(x+y)^a
-               ex r = (new mul(distrseq))->setflag(status_flags::dynallocated);
+               ex r = dynallocate<mul>(distrseq);
                return r.expand(options);
        }
        
@@ -889,7 +864,7 @@ ex power::expand(unsigned options) const
                if (are_ex_trivially_equal(basis,expanded_basis) && are_ex_trivially_equal(exponent,expanded_exponent)) {
                        return this->hold();
                } else {
-                       return (new power(expanded_basis,expanded_exponent))->setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0));
+                       return dynallocate<power>(expanded_basis, expanded_exponent).setflag(options == 0 ? status_flags::expanded : 0);
                }
        }
        
@@ -909,7 +884,7 @@ ex power::expand(unsigned options) const
        if (are_ex_trivially_equal(basis,expanded_basis) && are_ex_trivially_equal(exponent,expanded_exponent))
                return this->hold();
        else
-               return (new power(expanded_basis,expanded_exponent))->setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0));
+               return dynallocate<power>(expanded_basis, expanded_exponent).setflag(options == 0 ? status_flags::expanded : 0);
 }
 
 //////////
@@ -1197,14 +1172,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);
+                               epvector monomial;
+                               monomial.reserve(msize);
                                numeric factor = coeff;
                                for (unsigned i = 0; i < exponent.size(); ++i) {
                                        const ex & r = a.seq[i].rest;
@@ -1221,28 +1198,25 @@ 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(expair(r, _ex1));
                                                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(expair(r, exponent[i]));
                                                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(expair(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(std::move(result)))->setflag(status_flags::dynallocated |
-                                                            status_flags::expanded);
+               return dynallocate<add>(std::move(result)).setflag(status_flags::expanded);
        } else {
-               return (new add(std::move(result), ex_to<numeric>(a.overall_coeff).power(n)))->setflag(status_flags::dynallocated |
-                                                                                                      status_flags::expanded);
+               return dynallocate<add>(std::move(result), ex_to<numeric>(a.overall_coeff).power(n)).setflag(status_flags::expanded);
        }
 }
 
@@ -1277,27 +1251,27 @@ ex power::expand_add_2(const add & a, unsigned options)
                
                if (c.is_equal(_ex1)) {
                        if (is_exactly_a<mul>(r)) {
-                               result.push_back(a.combine_ex_with_coeff_to_pair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
-                                                                                _ex1));
+                               result.push_back(expair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
+                                                       _ex1));
                        } else {
-                               result.push_back(a.combine_ex_with_coeff_to_pair((new power(r,_ex2))->setflag(status_flags::dynallocated),
-                                                                                _ex1));
+                               result.push_back(expair(dynallocate<power>(r, _ex2),
+                                                       _ex1));
                        }
                } else {
                        if (is_exactly_a<mul>(r)) {
-                               result.push_back(a.combine_ex_with_coeff_to_pair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
-                                                                                ex_to<numeric>(c).power_dyn(*_num2_p)));
+                               result.push_back(expair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
+                                                       ex_to<numeric>(c).power_dyn(*_num2_p)));
                        } else {
-                               result.push_back(a.combine_ex_with_coeff_to_pair((new power(r,_ex2))->setflag(status_flags::dynallocated),
-                                                                                ex_to<numeric>(c).power_dyn(*_num2_p)));
+                               result.push_back(expair(dynallocate<power>(r, _ex2),
+                                                       ex_to<numeric>(c).power_dyn(*_num2_p)));
                        }
                }
 
                for (epvector::const_iterator cit1=cit0+1; cit1!=last; ++cit1) {
                        const ex & r1 = cit1->rest;
                        const ex & c1 = cit1->coeff;
-                       result.push_back(a.combine_ex_with_coeff_to_pair(mul(r,r1).expand(options),
-                                                                        _num2_p->mul(ex_to<numeric>(c)).mul_dyn(ex_to<numeric>(c1))));
+                       result.push_back(expair(mul(r,r1).expand(options),
+                                               _num2_p->mul(ex_to<numeric>(c)).mul_dyn(ex_to<numeric>(c1))));
                }
        }
        
@@ -1310,11 +1284,9 @@ ex power::expand_add_2(const add & a, unsigned options)
        GINAC_ASSERT(result.size() == result_size);
 
        if (a.overall_coeff.is_zero()) {
-               return (new add(std::move(result)))->setflag(status_flags::dynallocated |
-                                                            status_flags::expanded);
+               return dynallocate<add>(std::move(result)).setflag(status_flags::expanded);
        } else {
-               return (new add(std::move(result), ex_to<numeric>(a.overall_coeff).power(2)))->setflag(status_flags::dynallocated |
-                                                                                                      status_flags::expanded);
+               return dynallocate<add>(std::move(result), ex_to<numeric>(a.overall_coeff).power(2)).setflag(status_flags::expanded);
        }
 }
 
@@ -1358,7 +1330,7 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr
                distrseq.push_back(p);
        }
 
-       const mul & result = static_cast<const mul &>((new mul(distrseq, ex_to<numeric>(m.overall_coeff).power_dyn(n)))->setflag(status_flags::dynallocated));
+       const mul & result = dynallocate<mul>(std::move(distrseq), ex_to<numeric>(m.overall_coeff).power_dyn(n));
        if (need_reexpand)
                return ex(result).expand(options);
        if (from_expand)