* Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
/*
- * GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2019 Johannes Gutenberg University Mainz, Germany
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
case info_flags::cinteger_polynomial:
case info_flags::rational_polynomial:
case info_flags::crational_polynomial:
- return exponent.info(info_flags::nonnegint) &&
- basis.info(inf);
+ return basis.info(inf) && exponent.info(info_flags::nonnegint);
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);
+ return basis.info(inf) && exponent.info(info_flags::integer);
+ case info_flags::real:
+ return basis.info(inf) && exponent.info(info_flags::integer);
case info_flags::expanded:
return (flags & status_flags::expanded);
case info_flags::positive:
return basis.info(info_flags::positive) && exponent.info(info_flags::real);
case info_flags::nonnegative:
return (basis.info(info_flags::positive) && exponent.info(info_flags::real)) ||
- (basis.info(info_flags::real) && exponent.info(info_flags::integer) && exponent.info(info_flags::even));
+ (basis.info(info_flags::real) && exponent.info(info_flags::even));
case info_flags::has_indices: {
if (flags & status_flags::has_indices)
return true;
// Turn (x^c)^d into x^(c*d) in the case that x is positive and c is real.
if (is_exactly_a<power>(basis) && basis.op(0).info(info_flags::positive) && basis.op(1).info(info_flags::real))
- return power(basis.op(0), basis.op(1) * exponent);
+ return dynallocate<power>(basis.op(0), basis.op(1) * exponent);
if ( num_exponent ) {
// 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);
}
}
}
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>(basis)) {
- return expand_mul(ex_to<mul>(basis), *num_exponent, 0);
+ return expand_mul(ex_to<mul>(basis), *num_exponent, false);
}
// (2*x + 6*y)^(-4) -> 1/16*(x + 3*y)^(-4)
return this->hold();
}
-ex power::evalf(int level) const
+ex power::evalf() const
{
- ex ebasis;
+ ex ebasis = basis.evalf();
ex eexponent;
- if (level==1) {
- ebasis = basis;
+ if (!is_exactly_a<numeric>(exponent))
+ eexponent = exponent.evalf();
+ else
eexponent = exponent;
- } else if (level == -max_recursion_level) {
- throw(std::runtime_error("max recursion level reached"));
- } else {
- ebasis = basis.evalf(level-1);
- if (!is_exactly_a<numeric>(exponent))
- eexponent = exponent.evalf(level-1);
- else
- eexponent = exponent;
- }
- return power(ebasis,eexponent);
+ return dynallocate<power>(ebasis, eexponent);
}
ex power::evalm() 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);
}
}
// 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)) {
+ if (basis.is_equal(a) && exponent.is_equal(c) &&
+ (a.info(info_flags::nonnegative) || c.info(info_flags::integer))) {
// Re(a^c)
return *this;
}
// 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 {
// 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)) {
+ if (basis.is_equal(a) && exponent.is_equal(c) &&
+ (a.info(info_flags::nonnegative) || c.info(info_flags::integer))) {
// Im(a^c)
return 0;
}
// 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 {
// 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
{
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));
}
}
// 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;
// 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 dynallocate<mul>(std::move(prodseq)) * pow(newbasis, exponent);
} else
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
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 = dynallocate<mul>(distrseq);
// non-virtual functions in this class
//////////
-namespace { // anonymous namespace for power::expand_add() helpers
-
-/** Helper class to generate all bounded combinatorial partitions of an integer
- * n with exactly m parts (including zero parts) in non-decreasing order.
- */
-class partition_generator {
-private:
- // Partitions n into m parts, not including zero parts.
- // (Cf. OEIS sequence A008284; implementation adapted from Jörg Arndt's
- // FXT library)
- struct mpartition2
- {
- // partition: x[1] + x[2] + ... + x[m] = n and sentinel x[0] == 0
- std::vector<int> x;
- int n; // n>0
- int m; // 0<m<=n
- mpartition2(unsigned n_, unsigned m_)
- : x(m_+1), n(n_), m(m_)
- {
- for (int k=1; k<m; ++k)
- x[k] = 1;
- x[m] = n - m + 1;
- }
- bool next_partition()
- {
- int u = x[m]; // last element
- int k = m;
- int s = u;
- while (--k) {
- s += x[k];
- if (x[k] + 2 <= u)
- break;
- }
- if (k==0)
- return false; // current is last
- int f = x[k] + 1;
- while (k < m) {
- x[k] = f;
- s -= f;
- ++k;
- }
- x[m] = s;
- return true;
- }
- } mpgen;
- int m; // number of parts 0<m<=n
- mutable std::vector<int> partition; // current partition
-public:
- partition_generator(unsigned n_, unsigned m_)
- : mpgen(n_, 1), m(m_), partition(m_)
- { }
- // returns current partition in non-decreasing order, padded with zeros
- const std::vector<int>& current() const
- {
- for (int i = 0; i < m - mpgen.m; ++i)
- partition[i] = 0; // pad with zeros
-
- for (int i = m - mpgen.m; i < m; ++i)
- partition[i] = mpgen.x[i - m + mpgen.m + 1];
-
- return partition;
- }
- bool next()
- {
- if (!mpgen.next_partition()) {
- if (mpgen.m == m || mpgen.m == mpgen.n)
- return false; // current is last
- // increment number of parts
- mpgen = mpartition2(mpgen.n, mpgen.m + 1);
- }
- return true;
- }
-};
-
-/** Helper class to generate all compositions of a partition of an integer n,
- * starting with the compositions which has non-decreasing order.
- */
-class composition_generator {
-private:
- // Generates all distinct permutations of a multiset.
- // (Based on Aaron Williams' algorithm 1 from "Loopless Generation of
- // Multiset Permutations using a Constant Number of Variables by Prefix
- // Shifts." <http://webhome.csc.uvic.ca/~haron/CoolMulti.pdf>)
- struct coolmulti {
- // element of singly linked list
- struct element {
- int value;
- element* next;
- element(int val, element* n)
- : value(val), next(n) {}
- ~element()
- { // recurses down to the end of the singly linked list
- delete next;
- }
- };
- 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)
- {
- for (unsigned n = 0; n < partition.size(); ++n) {
- head = new element(partition[n], head);
- if (n <= 1)
- i = head;
- }
- after_i = i->next;
- }
- ~coolmulti()
- { // deletes singly linked list
- delete head;
- }
- void next_permutation()
- {
- element *before_k;
- if (after_i->next != nullptr && i->value >= after_i->next->value)
- before_k = after_i;
- else
- before_k = i;
- element *k = before_k->next;
- before_k->next = k->next;
- k->next = head;
- if (k->value < head->value)
- i = k;
- after_i = i->next;
- head = k;
- }
- bool finished() const
- {
- return after_i->next == nullptr && after_i->value >= head->value;
- }
- } cmgen;
- bool atend; // needed for simplifying iteration over permutations
- bool trivial; // likewise, true if all elements are equal
- mutable std::vector<int> composition; // current compositions
-public:
- explicit composition_generator(const std::vector<int>& partition)
- : cmgen(partition), atend(false), trivial(true), composition(partition.size())
- {
- for (unsigned i=1; i<partition.size(); ++i)
- trivial = trivial && (partition[0] == partition[i]);
- }
- const std::vector<int>& current() const
- {
- coolmulti::element* it = cmgen.head;
- size_t i = 0;
- while (it != nullptr) {
- composition[i] = it->value;
- it = it->next;
- ++i;
- }
- 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;
- 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, long n, unsigned options)
// 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());
+ partition_with_zero_parts_generator partitions(k, a.seq.size());
do {
- const std::vector<int>& partition = partitions.current();
+ const std::vector<unsigned>& partition = partitions.get();
// 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();
+ const std::vector<unsigned>& exponent = compositions.get();
epvector monomial;
monomial.reserve(msize);
numeric factor = coeff;
factor = factor.mul(c.power(exponent[i]));
}
}
- result.push_back(expair(mul(monomial).expand(options), factor));
+ result.push_back(expair(mul(std::move(monomial)).expand(options), factor));
} while (compositions.next());
} while (partitions.next());
}
}
result.reserve(result_size);
- epvector::const_iterator last = a.seq.end();
+ auto last = a.seq.end();
// power(+(x,...,z;c),2)=power(+(x,...,z;0),2)+2*c*+(x,...,z;0)+c*c
// first part: ignore overall_coeff and expand other terms
- for (epvector::const_iterator cit0=a.seq.begin(); cit0!=last; ++cit0) {
+ for (auto cit0=a.seq.begin(); cit0!=last; ++cit0) {
const ex & r = cit0->rest;
const ex & c = cit0->coeff;
}
}
- for (epvector::const_iterator cit1=cit0+1; cit1!=last; ++cit1) {
+ for (auto cit1=cit0+1; cit1!=last; ++cit1) {
const ex & r1 = cit1->rest;
const ex & c1 = cit1->coeff;
result.push_back(expair(mul(r,r1).expand(options),