* 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
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 dynallocate<power>(ebasis, eexponent);
}
// 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),
git://www.ginac.de / ginac.git / blobdiff