X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fpower.cpp;h=c3fe14a132e14f98b79a24f19ec9ad5d0ef95d71;hp=53f7bb8700e9d27ce386ef24b6a5f13f38db4ff0;hb=d56a0f74afa5380a1730599c3b1b21f34be2f061;hpb=4c47ecd9caa39ba1a31f5294e395fcbdf2006431 diff --git a/ginac/power.cpp b/ginac/power.cpp index 53f7bb87..c3fe14a1 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's symbolic exponentiation (basis^exponent). */ /* - * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2015 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 @@ -20,11 +20,6 @@ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ -#include -#include -#include -#include - #include "power.h" #include "expairseq.h" #include "add.h" @@ -43,6 +38,11 @@ #include "relational.h" #include "compiler.h" +#include +#include +#include +#include + namespace GiNaC { GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(power, basic, @@ -53,8 +53,6 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(power, basic, print_func(&power::do_print_python_repr). print_func(&power::do_print_csrc_cl_N)) -typedef std::vector intvector; - ////////// // default constructor ////////// @@ -240,6 +238,11 @@ bool power::info(unsigned inf) const basis.info(inf); 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)); case info_flags::has_indices: { if (flags & status_flags::has_indices) return true; @@ -285,11 +288,16 @@ ex power::map(map_function & f) const bool power::is_polynomial(const ex & var) const { - if (exponent.has(var)) - return false; - if (!exponent.info(info_flags::nonnegint)) - return false; - return basis.is_polynomial(var); + if (basis.is_polynomial(var)) { + if (basis.has(var)) + // basis is non-constant polynomial in var + return exponent.info(info_flags::nonnegint); + else + // basis is constant in var + return !exponent.has(var); + } + // basis is a non-polynomial function of var + return false; } int power::degree(const ex & s) const @@ -360,7 +368,7 @@ ex power::coeff(const ex & s, int n) const * - ^(1,x) -> 1 * - ^(c1,c2) -> *(c1^n,c1^(c2-n)) (so that 0<(c2-n)<1, try to evaluate roots, possibly in numerator and denominator of c1) * - ^(^(x,c1),c2) -> ^(x,c1*c2) if x is positive and c1 is real. - * - ^(^(x,c1),c2) -> ^(x,c1*c2) (c2 integer or -1 < c1 <= 1, case c1=1 should not happen, see below!) + * - ^(^(x,c1),c2) -> ^(x,c1*c2) (c2 integer or -1 < c1 <= 1 or (c1=-1 and c2>0), case c1=1 should not happen, see below!) * - ^(*(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) @@ -478,7 +486,7 @@ ex power::eval(int level) const } // ^(^(x,c1),c2) -> ^(x,c1*c2) - // (c1, c2 numeric(), c2 integer or -1 < c1 <= 1, + // (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(ebasis)) { const power & sub_power = ex_to(ebasis); @@ -487,7 +495,8 @@ ex power::eval(int level) const if (is_exactly_a(sub_exponent)) { const numeric & num_sub_exponent = ex_to(sub_exponent); GINAC_ASSERT(num_sub_exponent!=numeric(1)); - if (num_exponent->is_integer() || (abs(num_sub_exponent) - (*_num1_p)).is_negative()) { + 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)); } } @@ -537,6 +546,7 @@ ex power::eval(int level) const 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), @@ -546,6 +556,7 @@ ex power::eval(int level) const 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), @@ -665,12 +676,23 @@ ex power::eval_ncmul(const exvector & v) const ex power::conjugate() const { - ex newbasis = basis.conjugate(); - ex newexponent = exponent.conjugate(); - if (are_ex_trivially_equal(basis, newbasis) && are_ex_trivially_equal(exponent, newexponent)) { - return *this; + // conjugate(pow(x,y))==pow(conjugate(x),conjugate(y)) unless on the + // branch cut which runs along the negative real axis. + if (basis.info(info_flags::positive)) { + ex newexponent = exponent.conjugate(); + if (are_ex_trivially_equal(exponent, newexponent)) { + return *this; + } + return (new power(basis, newexponent))->setflag(status_flags::dynallocated); + } + 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 (new power(newbasis, newexponent))->setflag(status_flags::dynallocated); + return conjugate_function(*this).hold(); } ex power::real_part() const @@ -720,14 +742,11 @@ ex power::imag_part() const ex b=basis.imag_part(); ex c=exponent.real_part(); 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 power(abs(basis),c)*exp(-d*atan2(b,a))*sin(c*atan2(b,a)+d*log(abs(basis))); } // protected -// protected - /** Implementation of ex::diff() for a power. * @see ex::diff */ ex power::derivative(const symbol & s) const @@ -777,6 +796,51 @@ ex power::expand(unsigned options) const return *this; } + // (x*p)^c -> x^c * p^c, if p>0 + // makes sense before expanding the basis + if (is_exactly_a(basis) && !basis.info(info_flags::indefinite)) { + const mul &m = ex_to(basis); + exvector prodseq; + 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); + 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; + } + + // 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)); + else if (m.overall_coeff.info(info_flags::negative) && m.overall_coeff != _ex_1) + prodseq.push_back(power(-m.overall_coeff, exponent)); + else + coeff *= m.overall_coeff; + + // If positive/negative factors are found, then extract them. + // 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_to(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); + } else + ex_to(basis).setflag(status_flags::purely_indefinite); + } + const ex expanded_basis = basis.expand(options); const ex expanded_exponent = exponent.expand(options); @@ -846,96 +910,326 @@ ex power::expand(unsigned options) const // 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-decreaing 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 x; + int n; // n>0 + int m; // 0 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& 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." ) + 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& partition) + { + head = NULL; + 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 != NULL && 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 == NULL && 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 composition; // current compositions +public: + explicit composition_generator(const std::vector& partition) + : cmgen(partition), atend(false), trivial(true), composition(partition.size()) + { + for (unsigned i=1; i& current() const + { + coolmulti::element* it = cmgen.head; + size_t i = 0; + while (it != NULL) { + 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 p) +{ + numeric n = 0, d = 1; + std::vector::const_iterator it = p.begin(), itend = p.end(); + while (it != itend) { + n += numeric(*it); + d *= factorial(numeric(*it)); + ++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 { + // The special case power(+(x,...y;x),2) can be optimized better. if (n==2) return expand_add_2(a, options); - const size_t m = a.nops(); - exvector result; + // method: + // + // Consider base as the sum of all symbolic terms and the overall numeric + // coefficient and apply the binomial theorem: + // S = power(+(x,...,z;c),n) + // = power(+(+(x,...,z;0);c),n) + // = sum(binomial(n,k)*power(+(x,...,z;0),k)*c^(n-k), k=1..n) + c^n + // Then, apply the multinomial theorem to expand all power(+(x,...,z;0),k): + // The multinomial theorem is computed by an outer loop over all + // partitions of the exponent and an inner loop over all compositions of + // that partition. This method makes the expansion a combinatorial + // problem and allows us to directly construct the expanded sum and also + // to re-use the multinomial coefficients (since they depend only on the + // partition, not on the composition). + // + // multinomial power(+(x,y,z;0),3) example: + // partition : compositions : multinomial coefficient + // [0,0,3] : [3,0,0],[0,3,0],[0,0,3] : 3!/(3!*0!*0!) = 1 + // [0,1,2] : [2,1,0],[1,2,0],[2,0,1],... : 3!/(2!*1!*0!) = 3 + // [1,1,1] : [1,1,1] : 3!/(1!*1!*1!) = 6 + // => (x + y + z)^3 = + // x^3 + y^3 + z^3 + // + 3*x^2*y + 3*x*y^2 + 3*y^2*z + 3*y*z^2 + 3*x*z^2 + 3*x^2*z + // + 6*x*y*z + // + // multinomial power(+(x,y,z;0),4) example: + // partition : compositions : multinomial coefficient + // [0,0,4] : [4,0,0],[0,4,0],[0,0,4] : 4!/(4!*0!*0!) = 1 + // [0,1,3] : [3,1,0],[1,3,0],[3,0,1],... : 4!/(3!*1!*0!) = 4 + // [0,2,2] : [2,2,0],[2,0,2],[0,2,2] : 4!/(2!*2!*0!) = 6 + // [1,1,2] : [2,1,1],[1,2,1],[1,1,2] : 4!/(2!*1!*1!) = 12 + // (no [1,1,1,1] partition since it has too many parts) + // => (x + y + z)^4 = + // x^4 + y^4 + z^4 + // + 4*x^3*y + 4*x*y^3 + 4*y^3*z + 4*y*z^3 + 4*x*z^3 + 4*x^3*z + // + 6*x^2*y^2 + 6*y^2*z^2 + 6*x^2*z^2 + // + 12*x^2*y*z + 12*x*y^2*z + 12*x*y*z^2 + // + // Summary: + // r = 0 + // for k from 0 to n: + // f = c^(n-k)*binomial(n,k) + // for p in all partitions of n with m parts (including zero parts): + // h = f * multinomial coefficient of p + // for c in all compositions of p: + // t = 1 + // for e in all elements of c: + // t = t * a[e]^e + // r = r + h*t + // return r + + epvector result; // The number of terms will be the number of combinatorial compositions, // i.e. the number of unordered arrangements of m nonnegative integers // which sum up to n. It is frequently written as C_n(m) and directly - // related with binomial coefficients: - result.reserve(binomial(numeric(n+m-1), numeric(m-1)).to_int()); - intvector k(m-1); - intvector k_cum(m-1); // k_cum[l]:=sum(i=0,l,k[l]); - intvector upper_limit(m-1); - - for (size_t l=0; l(b)); - GINAC_ASSERT(!is_exactly_a(b) || - !is_exactly_a(ex_to(b).exponent) || - !ex_to(ex_to(b).exponent).is_pos_integer() || - !is_exactly_a(ex_to(b).basis) || - !is_exactly_a(ex_to(b).basis) || - !is_exactly_a(ex_to(b).basis)); - if (is_exactly_a(b)) - term.push_back(expand_mul(ex_to(b), numeric(k[l]), options, true)); - else - term.push_back(power(b,k[l])); - } - - const ex & b = a.op(m - 1); - GINAC_ASSERT(!is_exactly_a(b)); - GINAC_ASSERT(!is_exactly_a(b) || - !is_exactly_a(ex_to(b).exponent) || - !ex_to(ex_to(b).exponent).is_pos_integer() || - !is_exactly_a(ex_to(b).basis) || - !is_exactly_a(ex_to(b).basis) || - !is_exactly_a(ex_to(b).basis)); - if (is_exactly_a(b)) - term.push_back(expand_mul(ex_to(b), numeric(n-k_cum[m-2]), options, true)); - else - term.push_back(power(b,n-k_cum[m-2])); - - numeric f = binomial(numeric(n),numeric(k[0])); - for (std::size_t l = 1; l < m - 1; ++l) - f *= binomial(numeric(n-k_cum[l-1]),numeric(k[l])); - - term.push_back(f); - - result.push_back(ex((new mul(term))->setflag(status_flags::dynallocated)).expand(options)); - - // increment k[] - bool done = false; - std::size_t l = m - 2; - while ((++k[l]) > upper_limit[l]) { - k[l] = 0; - if (l != 0) - --l; - else { - done = true; - break; + result.reserve(result_size); + + // Iterate over all terms in binomial expansion of + // S = power(+(x,...,z;c),n) + // = sum(binomial(n,k)*power(+(x,...,z;0),k)*c^(n-k), k=1..n) + c^n + for (int k = 1; k <= n; ++k) { + numeric binomial_coefficient; // binomial(n,k)*c^(n-k) + if (a.overall_coeff.is_zero()) { + // degenerate case with zero overall_coeff: + // apply multinomial theorem directly to power(+(x,...z;0),n) + binomial_coefficient = 1; + if (k < n) { + continue; } + } else { + binomial_coefficient = binomial(numeric(n), numeric(k)) * pow(ex_to(a.overall_coeff), numeric(n-k)); } - if (done) - break; - // recalc k_cum[] and upper_limit[] - k_cum[l] = (l==0 ? k[0] : k_cum[l-1]+k[l]); + // 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()); + do { + const std::vector& partition = partitions.current(); + const numeric coeff = multinomial_coefficient(partition) * binomial_coefficient; + + // Iterate over all compositions of the current partition. + composition_generator compositions(partition); + do { + const std::vector& exponent = compositions.current(); + exvector term; + term.reserve(n); + 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(r)); + GINAC_ASSERT(!is_exactly_a(r) || + !is_exactly_a(ex_to(r).exponent) || + !ex_to(ex_to(r).exponent).is_pos_integer() || + !is_exactly_a(ex_to(r).basis) || + !is_exactly_a(ex_to(r).basis) || + !is_exactly_a(ex_to(r).basis)); + if (exponent[i] == 0) { + // optimize away + } else if (exponent[i] == 1) { + // optimized + term.push_back(r); + factor = factor.mul(ex_to(c)); + } else { // general case exponent[i] > 1 + term.push_back((new power(r, exponent[i]))->setflag(status_flags::dynallocated)); + factor = factor.mul(ex_to(c).power(exponent[i])); + } + } + result.push_back(a.combine_ex_with_coeff_to_pair(mul(term).expand(options), factor)); + } while (compositions.next()); + } while (partitions.next()); + } - for (size_t i=l+1; isetflag(status_flags::dynallocated | + status_flags::expanded); + } else { + return (new add(result, ex_to(a.overall_coeff).power(n)))->setflag(status_flags::dynallocated | + status_flags::expanded); } - - return (new add(result))->setflag(status_flags::dynallocated | - status_flags::expanded); } @@ -964,11 +1258,11 @@ ex power::expand_add_2(const add & a, unsigned options) const if (c.is_equal(_ex1)) { if (is_exactly_a(r)) { - sum.push_back(expair(expand_mul(ex_to(r), *_num2_p, options, true), - _ex1)); + sum.push_back(a.combine_ex_with_coeff_to_pair(expand_mul(ex_to(r), *_num2_p, options, true), + _ex1)); } else { - sum.push_back(expair((new power(r,_ex2))->setflag(status_flags::dynallocated), - _ex1)); + sum.push_back(a.combine_ex_with_coeff_to_pair((new power(r,_ex2))->setflag(status_flags::dynallocated), + _ex1)); } } else { if (is_exactly_a(r)) { @@ -983,14 +1277,14 @@ ex power::expand_add_2(const add & a, unsigned options) const 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((new mul(r,r1))->setflag(status_flags::dynallocated), + sum.push_back(a.combine_ex_with_coeff_to_pair(mul(r,r1).expand(options), _num2_p->mul(ex_to(c)).mul_dyn(ex_to(c1)))); } } GINAC_ASSERT(sum.size()==(a.seq.size()*(a.seq.size()+1))/2); - // second part: add terms coming from overall_factor (if != 0) + // 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) {