* Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
/*
- * GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2017 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;
* - ^(*(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())
}
// ^(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 ) {
// 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);
}
}
}
// ^(^(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)) {
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();
icont = icont.mul(*_num_1_p);
if (canonicalizable && (icont != *_num1_p)) {
- const add& addref = ex_to<add>(ebasis);
+ 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);
// ^(*(...,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()) {
// ^(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 dynallocate<power>(ebasis, eexponent).setflag(status_flags::evaluated);
+
+ 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);
}
}
// 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)) {
// 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)
composition_generator compositions(partition);
do {
const std::vector<int>& exponent = compositions.current();
- exvector monomial;
+ epvector monomial;
monomial.reserve(msize);
numeric factor = coeff;
for (unsigned i = 0; i < exponent.size(); ++i) {
// optimize away
} else if (exponent[i] == 1) {
// optimized
- monomial.push_back(r);
+ monomial.push_back(expair(r, _ex1));
if (c != *_num1_p)
factor = factor.mul(c);
} else { // general case exponent[i] > 1
- monomial.push_back(dynallocate<power>(r, 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(monomial).expand(options), factor));
+ result.push_back(expair(mul(std::move(monomial)).expand(options), factor));
} while (compositions.next());
} while (partitions.next());
}
GINAC_ASSERT(result.size() == result_size);
-
if (a.overall_coeff.is_zero()) {
return dynallocate<add>(std::move(result)).setflag(status_flags::expanded);
} else {
}
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;
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(dynallocate<power>(r, _ex2),
- _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(dynallocate<power>(r, _ex2),
- 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) {
+ for (auto 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))));
}
}
git://www.ginac.de / ginac.git / blobdiff