#include <limits>
#include <stdexcept>
#include <vector>
+#include <algorithm>
namespace GiNaC {
static void print_sym_pow(const print_context & c, const symbol &x, int exp)
{
// Optimal output of integer powers of symbols to aid compiler CSE.
- // C.f. ISO/IEC 14882:1998, section 1.9 [intro execution], paragraph 15
+ // C.f. ISO/IEC 14882:2011, section 1.9 [intro execution], paragraph 15
// to learn why such a parenthesation is really necessary.
if (exp == 1) {
x.print(c);
void power::do_print_csrc(const print_csrc & c, unsigned level) const
{
// Integer powers of symbols are printed in a special, optimized way
- if (exponent.info(info_flags::integer)
- && (is_a<symbol>(basis) || is_a<constant>(basis))) {
+ if (exponent.info(info_flags::integer) &&
+ (is_a<symbol>(basis) || is_a<constant>(basis))) {
int exp = ex_to<numeric>(exponent).to_int();
if (exp > 0)
c.s << '(';
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;
}
const ex & ebasis = level==1 ? basis : basis.eval(level-1);
const ex & eexponent = level==1 ? exponent : exponent.eval(level-1);
- const numeric *num_basis = NULL;
- const numeric *num_exponent = NULL;
+ const numeric *num_basis = nullptr;
+ const numeric *num_exponent = nullptr;
if (is_exactly_a<numeric>(ebasis)) {
num_basis = &ex_to<numeric>(ebasis);
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);
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 {
if (is_exactly_a<numeric>(sub_exponent)) {
const numeric & num_sub_exponent = ex_to<numeric>(sub_exponent);
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())) {
+ 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));
}
}
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 (epvector::iterator i = addp->seq.begin(); i != addp->seq.end(); ++i)
- i->coeff = ex_to<numeric>(i->coeff).div_dyn(icont);
+ 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);
}
}
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));
}
}
}
if (num_exponent->is_pos_integer() &&
ebasis.return_type() != return_types::commutative &&
!is_a<matrix>(ebasis)) {
- return ncmul(exvector(num_exponent->to_int(), ebasis), true);
+ return ncmul(exvector(num_exponent->to_int(), ebasis));
}
}
are_ex_trivially_equal(eexponent,exponent)) {
return this->hold();
}
- return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated |
- status_flags::evaluated);
+ return dynallocate<power>(ebasis, eexponent).setflag(status_flags::evaluated);
}
ex power::evalf(int level) 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
return basic::has(other, options);
if (!is_a<power>(other))
return basic::has(other, options);
- if (!exponent.info(info_flags::integer)
- || !other.op(1).info(info_flags::integer))
+ if (!exponent.info(info_flags::integer) ||
+ !other.op(1).info(info_flags::integer))
return basic::has(other, options);
- if (exponent.info(info_flags::posint)
- && other.op(1).info(info_flags::posint)
- && ex_to<numeric>(exponent).to_int()
- > ex_to<numeric>(other.op(1)).to_int()
- && basis.match(other.op(0)))
+ if (exponent.info(info_flags::posint) &&
+ other.op(1).info(info_flags::posint) &&
+ ex_to<numeric>(exponent) > ex_to<numeric>(other.op(1)) &&
+ basis.match(other.op(0)))
return true;
- if (exponent.info(info_flags::negint)
- && other.op(1).info(info_flags::negint)
- && ex_to<numeric>(exponent).to_int()
- < ex_to<numeric>(other.op(1)).to_int()
- && basis.match(other.op(0)))
+ if (exponent.info(info_flags::negint) &&
+ other.op(1).info(info_flags::negint) &&
+ ex_to<numeric>(exponent) < ex_to<numeric>(other.op(1)) &&
+ basis.match(other.op(0)))
return true;
return basic::has(other, options);
}
if (!(options & subs_options::algebraic))
return subs_one_level(m, options);
- for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
+ for (auto & it : m) {
int nummatches = std::numeric_limits<int>::max();
exmap repls;
- 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);
+ 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);
return (ex_to<basic>(result)).subs_one_level(m, options);
}
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();
}
ex power::real_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)) {
+ // Re(a^c)
+ return *this;
+ }
+
+ const ex b = basis.imag_part();
if (exponent.info(info_flags::integer)) {
- ex basis_real = basis.real_part();
- if (basis_real == basis)
- return *this;
- realsymbol a("a"),b("b");
- ex result;
- if (exponent.info(info_flags::posint))
- result = power(a+I*b,exponent);
- else
- result = power(a/(a*a+b*b)-I*b/(a*a+b*b),-exponent);
- result = result.expand();
- result = result.real_part();
- result = result.subs(lst( a==basis_real, b==basis.imag_part() ));
+ // 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))).
+ long NN = N > 0 ? N : -N;
+ ex numer = N > 0 ? _ex1 : power(power(a,2) + power(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;
+ if (n % 4 == 0) {
+ result += term; // sign: I^n w/ n == 4*m
+ } else {
+ result -= term; // sign: I^n w/ n == 4*m+2
+ }
+ }
return result;
}
-
- ex a = basis.real_part();
- ex b = basis.imag_part();
- ex c = exponent.real_part();
- ex d = exponent.imag_part();
+
+ // 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)));
}
ex power::imag_part() const
{
+ const ex a = basis.real_part();
+ const ex c = exponent.real_part();
+ if (basis.is_equal(a) && exponent.is_equal(c)) {
+ // Im(a^c)
+ return 0;
+ }
+
+ const ex b = basis.imag_part();
if (exponent.info(info_flags::integer)) {
- ex basis_real = basis.real_part();
- if (basis_real == basis)
- return 0;
- realsymbol a("a"),b("b");
- ex result;
- if (exponent.info(info_flags::posint))
- result = power(a+I*b,exponent);
- else
- result = power(a/(a*a+b*b)-I*b/(a*a+b*b),-exponent);
- result = result.expand();
- result = result.imag_part();
- result = result.subs(lst( a==basis_real, b==basis.imag_part() ));
+ // 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))).
+ 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 result = 0;
+ for (long n = 1; n <= NN; n += 2) {
+ ex term = binomial(NN, n) * power(a, NN-n) * power(b, n) / numer;
+ if (n % 4 == p) {
+ result += term; // sign: I^n w/ n == 4*m+p
+ } else {
+ result -= term; // sign: I^n w/ n == 4*m+2+p
+ }
+ }
return result;
}
-
- ex a=basis.real_part();
- ex b=basis.imag_part();
- ex c=exponent.real_part();
- ex d=exponent.imag_part();
+
+ // 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)));
}
newseq.reserve(2);
newseq.push_back(expair(basis, exponent - _ex1));
newseq.push_back(expair(basis.diff(s), _ex1));
- return mul(newseq, exponent);
+ return mul(std::move(newseq), exponent);
} else {
// D(b^e) = b^e * (D(e)*ln(b) + e*D(b)/b)
return mul(*this,
epvector powseq;
prodseq.reserve(m.seq.size() + 1);
powseq.reserve(m.seq.size() + 1);
- epvector::const_iterator last = m.seq.end();
- epvector::const_iterator cit = m.seq.begin();
bool possign = true;
// search for positive/negative factors
- while (cit!=last) {
- ex e=m.recombine_pair_to_ex(*cit);
+ for (auto & cit : m.seq) {
+ ex e=m.recombine_pair_to_ex(cit);
if (e.info(info_flags::positive))
prodseq.push_back(pow(e, exponent).expand(options));
else if (e.info(info_flags::negative)) {
prodseq.push_back(pow(-e, exponent).expand(options));
possign = !possign;
} else
- powseq.push_back(*cit);
- ++cit;
+ powseq.push_back(cit);
}
// take care on the numeric coefficient
// 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(powseq);
+ ex newbasis = coeff*mul(std::move(powseq));
ex_to<basic>(newbasis).setflag(status_flags::purely_indefinite);
- return ((new mul(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);
}
const add &a = ex_to<add>(expanded_exponent);
exvector distrseq;
distrseq.reserve(a.seq.size() + 1);
- epvector::const_iterator last = a.seq.end();
- epvector::const_iterator cit = a.seq.begin();
- while (cit!=last) {
- distrseq.push_back(power(expanded_basis, a.recombine_pair_to_ex(*cit)));
- ++cit;
+ for (auto & cit : a.seq) {
+ distrseq.push_back(power(expanded_basis, a.recombine_pair_to_ex(cit)));
}
// Make sure that e.g. (x+y)^(2+a) expands the (x+y)^2 factor
if (ex_to<numeric>(a.overall_coeff).is_integer()) {
const numeric &num_exponent = ex_to<numeric>(a.overall_coeff);
- int int_exponent = num_exponent.to_int();
+ long int_exponent = num_exponent.to_int();
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));
// 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);
}
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);
}
}
// integer numeric exponent
const numeric & num_exponent = ex_to<numeric>(expanded_exponent);
- int int_exponent = num_exponent.to_int();
+ long int_exponent = num_exponent.to_long();
// (x+y)^n, n>0
if (int_exponent > 0 && is_exactly_a<add>(expanded_basis))
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);
}
//////////
element *head, *i, *after_i;
// NB: Partition must be sorted in non-decreasing order.
explicit coolmulti(const std::vector<int>& partition)
+ : head(nullptr), i(nullptr), after_i(nullptr)
{
- head = NULL;
for (unsigned n = 0; n < partition.size(); ++n) {
head = new element(partition[n], head);
if (n <= 1)
void next_permutation()
{
element *before_k;
- if (after_i->next != NULL && i->value >= after_i->next->value)
+ if (after_i->next != nullptr && i->value >= after_i->next->value)
before_k = after_i;
else
before_k = i;
}
bool finished() const
{
- return after_i->next == NULL && after_i->value >= head->value;
+ return after_i->next == nullptr && after_i->value >= head->value;
}
} cmgen;
bool atend; // needed for simplifying iteration over permutations
{
coolmulti::element* it = cmgen.head;
size_t i = 0;
- while (it != NULL) {
+ while (it != nullptr) {
composition[i] = it->value;
it = it->next;
++i;
* where n = p1+p2+...+pk, i.e. p is a partition of n.
*/
const numeric
-multinomial_coefficient(const std::vector<int> p)
+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;
+ for (auto & it : p) {
+ n += numeric(it);
+ d *= factorial(numeric(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
+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)
// 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();
+ 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;
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;
- const ex & c = a.seq[i].coeff;
GINAC_ASSERT(!is_exactly_a<add>(r));
GINAC_ASSERT(!is_exactly_a<power>(r) ||
!is_exactly_a<numeric>(ex_to<power>(r).exponent) ||
!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
- term.push_back(r);
- factor = factor.mul(ex_to<numeric>(c));
+ 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));
- factor = factor.mul(ex_to<numeric>(c).power(exponent[i]));
+ 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(result))->setflag(status_flags::dynallocated |
- status_flags::expanded);
+ return dynallocate<add>(std::move(result)).setflag(status_flags::expanded);
} else {
- return (new add(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);
}
}
/** Special case of power::expand_add. Expands a^2 where a is an add.
* @see power::expand_add */
-ex power::expand_add_2(const add & a, unsigned options) const
+ex power::expand_add_2(const add & a, unsigned options)
{
- epvector sum;
- size_t a_nops = a.nops();
- sum.reserve((a_nops*(a_nops+1))/2);
+ epvector result;
+ size_t result_size = (a.nops() * (a.nops()+1)) / 2;
+ if (!a.overall_coeff.is_zero()) {
+ // the result's overall_coeff is one of the terms
+ --result_size;
+ }
+ result.reserve(result_size);
+
epvector::const_iterator last = a.seq.end();
// power(+(x,...,z;c),2)=power(+(x,...,z;0),2)+2*c*+(x,...,z;0)+c*c
if (c.is_equal(_ex1)) {
if (is_exactly_a<mul>(r)) {
- sum.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 {
- sum.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)) {
- sum.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 {
- sum.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;
- sum.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))));
}
}
- GINAC_ASSERT(sum.size()==(a.seq.size()*(a.seq.size()+1))/2);
-
// second part: add terms coming from overall_coeff (if != 0)
if (!a.overall_coeff.is_zero()) {
- epvector::const_iterator i = a.seq.begin(), end = a.seq.end();
- while (i != end) {
- sum.push_back(a.combine_pair_with_coeff_to_pair(*i, ex_to<numeric>(a.overall_coeff).mul_dyn(*_num2_p)));
- ++i;
- }
- sum.push_back(expair(ex_to<numeric>(a.overall_coeff).power_dyn(*_num2_p),_ex1));
+ for (auto & i : a.seq)
+ result.push_back(a.combine_pair_with_coeff_to_pair(i, ex_to<numeric>(a.overall_coeff).mul_dyn(*_num2_p)));
+ }
+
+ GINAC_ASSERT(result.size() == result_size);
+
+ if (a.overall_coeff.is_zero()) {
+ return dynallocate<add>(std::move(result)).setflag(status_flags::expanded);
+ } else {
+ return dynallocate<add>(std::move(result), ex_to<numeric>(a.overall_coeff).power(2)).setflag(status_flags::expanded);
}
-
- GINAC_ASSERT(sum.size()==(a_nops*(a_nops+1))/2);
-
- return (new add(sum))->setflag(status_flags::dynallocated | status_flags::expanded);
}
/** Expand factors of m in m^n where m is a mul and n is an integer.
* @see power::expand */
-ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool from_expand) const
+ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool from_expand)
{
GINAC_ASSERT(n.is_integer());
}
// do not bother to rename indices if there are no any.
- if ((!(options & expand_options::expand_rename_idx))
- && m.info(info_flags::has_indices))
+ if (!(options & expand_options::expand_rename_idx) &&
+ m.info(info_flags::has_indices))
options |= expand_options::expand_rename_idx;
// Leave it to multiplication since dummy indices have to be renamed
if ((options & expand_options::expand_rename_idx) &&
- (get_all_dummy_indices(m).size() > 0) && n.is_positive()) {
+ (get_all_dummy_indices(m).size() > 0) && n.is_positive()) {
ex result = m;
exvector va = get_all_dummy_indices(m);
sort(va.begin(), va.end(), ex_is_less());
distrseq.reserve(m.seq.size());
bool need_reexpand = false;
- epvector::const_iterator last = m.seq.end();
- epvector::const_iterator cit = m.seq.begin();
- while (cit!=last) {
- expair p = m.combine_pair_with_coeff_to_pair(*cit, n);
- if (from_expand && is_exactly_a<add>(cit->rest) && ex_to<numeric>(p.coeff).is_pos_integer()) {
+ for (auto & cit : m.seq) {
+ expair p = m.combine_pair_with_coeff_to_pair(cit, n);
+ if (from_expand && is_exactly_a<add>(cit.rest) && ex_to<numeric>(p.coeff).is_pos_integer()) {
// this happens when e.g. (a+b)^(1/2) gets squared and
// the resulting product needs to be reexpanded
need_reexpand = true;
}
distrseq.push_back(p);
- ++cit;
}
- 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)
git://www.ginac.de / ginac.git / blobdiff