X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Fpower.cpp;h=62fc3a8ddf070bd5c02a3fe0e2b75ad0006330e6;hb=e432b3743d3d32f6060600af580e49d7033dcae9;hp=4ba3877844038d07958869f7936159cc84eb624d;hpb=9413cd14faaf2980de3884915e22a5beda068ecc;p=ginac.git diff --git a/ginac/power.cpp b/ginac/power.cpp index 4ba38778..62fc3a8d 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-2016 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2018 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 @@ -669,7 +669,8 @@ 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)) { + if (basis.is_equal(a) && exponent.is_equal(c) && + (a.info(info_flags::nonnegative) || c.info(info_flags::integer))) { // Re(a^c) return *this; } @@ -704,7 +705,8 @@ 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; } @@ -889,188 +891,6 @@ 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-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 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(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 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 != 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 & p) -{ - numeric n = 0, d = 1; - for (auto & it : p) { - n += numeric(it); - d *= factorial(numeric(it)); - } - return factorial(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) @@ -1161,9 +981,9 @@ 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& partition = partitions.current(); + const std::vector& 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; @@ -1171,7 +991,7 @@ ex power::expand_add(const add & a, long n, unsigned options) // Iterate over all compositions of the current partition. composition_generator compositions(partition); do { - const std::vector& exponent = compositions.current(); + const std::vector& exponent = compositions.get(); epvector monomial; monomial.reserve(msize); numeric factor = coeff;