Fix bug in power::expand() with the overall coefficient.

[ginac.git] / ginac / power.cpp

diff --git a/ginac/power.cpp b/ginac/power.cpp
index 3c718027bdaa530782d063620fe75b20b440f681..ba7a66f4bc77cb912947fa06975b2ee801e304d9 100644 (file)
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
 
 /*
- *  GiNaC Copyright (C) 1999-2017 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2020 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;
        }
@@ -807,9 +809,10 @@ ex power::expand(unsigned options) const
                ex coeff=(possign? _ex1 : _ex_1);
                if (m.overall_coeff.info(info_flags::positive) && m.overall_coeff != _ex1)
                        prodseq.push_back(pow(m.overall_coeff, exponent));
-               else if (m.overall_coeff.info(info_flags::negative) && m.overall_coeff != _ex_1)
+               else if (m.overall_coeff.info(info_flags::negative) && m.overall_coeff != _ex_1) {
                        prodseq.push_back(pow(-m.overall_coeff, exponent));
-               else
+                       coeff = -coeff;
+               } else
                        coeff *= m.overall_coeff;
 
                // If positive/negative factors are found, then extract them.
@@ -979,9 +982,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<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;
@@ -989,7 +992,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<int>& exponent = compositions.current();
+                               const std::vector<unsigned>& exponent = compositions.get();
                                epvector monomial;
                                monomial.reserve(msize);
                                numeric factor = coeff;
@@ -1008,16 +1011,16 @@ ex power::expand_add(const add & a, long n, unsigned options)
                                                // optimize away
                                        } else if (exponent[i] == 1) {
                                                // optimized
-                                               monomial.push_back(expair(r, _ex1));
+                                               monomial.emplace_back(expair(r, _ex1));
                                                if (c != *_num1_p)
                                                        factor = factor.mul(c);
                                        } else { // general case exponent[i] > 1
-                                               monomial.push_back(expair(r, exponent[i]));
+                                               monomial.emplace_back(expair(r, exponent[i]));
                                                if (c != *_num1_p)
                                                        factor = factor.mul(c.power(exponent[i]));
                                        }
                                }
-                               result.push_back(expair(mul(std::move(monomial)).expand(options), factor));
+                               result.emplace_back(expair(mul(std::move(monomial)).expand(options), factor));
                        } while (compositions.next());
                } while (partitions.next());
        }
@@ -1061,27 +1064,27 @@ ex power::expand_add_2(const add & a, unsigned options)
                
                if (c.is_equal(_ex1)) {
                        if (is_exactly_a<mul>(r)) {
-                               result.push_back(expair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
-                                                       _ex1));
+                               result.emplace_back(expair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
+                                                          _ex1));
                        } else {
-                               result.push_back(expair(dynallocate<power>(r, _ex2),
-                                                       _ex1));
+                               result.emplace_back(expair(dynallocate<power>(r, _ex2),
+                                                          _ex1));
                        }
                } else {
                        if (is_exactly_a<mul>(r)) {
-                               result.push_back(expair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
-                                                       ex_to<numeric>(c).power_dyn(*_num2_p)));
+                               result.emplace_back(expair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
+                                                          ex_to<numeric>(c).power_dyn(*_num2_p)));
                        } else {
-                               result.push_back(expair(dynallocate<power>(r, _ex2),
-                                                       ex_to<numeric>(c).power_dyn(*_num2_p)));
+                               result.emplace_back(expair(dynallocate<power>(r, _ex2),
+                                                          ex_to<numeric>(c).power_dyn(*_num2_p)));
                        }
                }
 
                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),
-                                               _num2_p->mul(ex_to<numeric>(c)).mul_dyn(ex_to<numeric>(c1))));
+                       result.emplace_back(expair(mul(r,r1).expand(options),
+                                                  _num2_p->mul(ex_to<numeric>(c)).mul_dyn(ex_to<numeric>(c1))));
                }
        }