#include <limits>
#include <stdexcept>
#include <vector>
+#include <algorithm>
namespace GiNaC {
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)) {
+ // Re(a^c)
+ return *this;
+ }
+
+ const ex b = basis.imag_part();
if (exponent.info(info_flags::integer)) {
- ex basis_real = basis.real_part();
- if (basis_real == basis)
- return *this;
- realsymbol a("a"),b("b");
- ex result;
- if (exponent.info(info_flags::posint))
- result = power(a+I*b,exponent);
- else
- result = power(a/(a*a+b*b)-I*b/(a*a+b*b),-exponent);
- result = result.expand();
- result = result.real_part();
- result = result.subs(lst( a==basis_real, b==basis.imag_part() ));
+ // 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))).
+ long NN = N > 0 ? N : -N;
+ ex numer = N > 0 ? _ex1 : power(power(a,2) + power(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;
+ if (n % 4 == 0) {
+ result += term; // sign: I^n w/ n == 4*m
+ } else {
+ result -= term; // sign: I^n w/ n == 4*m+2
+ }
+ }
return result;
}
-
- ex a = basis.real_part();
- ex b = basis.imag_part();
- ex c = exponent.real_part();
- ex d = exponent.imag_part();
+
+ // 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)));
}
ex power::imag_part() const
{
+ const ex a = basis.real_part();
+ const ex c = exponent.real_part();
+ if (basis.is_equal(a) && exponent.is_equal(c)) {
+ // Im(a^c)
+ return 0;
+ }
+
+ const ex b = basis.imag_part();
if (exponent.info(info_flags::integer)) {
- ex basis_real = basis.real_part();
- if (basis_real == basis)
- return 0;
- realsymbol a("a"),b("b");
- ex result;
- if (exponent.info(info_flags::posint))
- result = power(a+I*b,exponent);
- else
- result = power(a/(a*a+b*b)-I*b/(a*a+b*b),-exponent);
- result = result.expand();
- result = result.imag_part();
- result = result.subs(lst( a==basis_real, b==basis.imag_part() ));
+ // 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))).
+ 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 result = 0;
+ for (long n = 1; n <= NN; n += 2) {
+ ex term = binomial(NN, n) * power(a, NN-n) * power(b, n) / numer;
+ if (n % 4 == p) {
+ result += term; // sign: I^n w/ n == 4*m+p
+ } else {
+ result -= term; // sign: I^n w/ n == 4*m+2+p
+ }
+ }
return result;
}
-
- ex a=basis.real_part();
- ex b=basis.imag_part();
- ex c=exponent.real_part();
- ex d=exponent.imag_part();
+
+ // 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)));
}
* where n = p1+p2+...+pk, i.e. p is a partition of n.
*/
const numeric
-multinomial_coefficient(const std::vector<int> p)
+multinomial_coefficient(const std::vector<int> & p)
{
numeric n = 0, d = 1;
std::vector<int>::const_iterator it = p.begin(), itend = p.end();
partition_generator partitions(k, a.seq.size());
do {
const std::vector<int>& partition = partitions.current();
+ // All monomials of this partition have the same number of terms and the same coefficient.
+ const unsigned msize = count_if(partition.begin(), partition.end(), bind2nd(std::greater<int>(), 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();
- exvector term;
- term.reserve(n);
+ exvector monomial;
+ monomial.reserve(msize);
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<add>(r));
GINAC_ASSERT(!is_exactly_a<power>(r) ||
!is_exactly_a<numeric>(ex_to<power>(r).exponent) ||
!is_exactly_a<add>(ex_to<power>(r).basis) ||
!is_exactly_a<mul>(ex_to<power>(r).basis) ||
!is_exactly_a<power>(ex_to<power>(r).basis));
+ GINAC_ASSERT(is_exactly_a<numeric>(a.seq[i].coeff));
+ const numeric & c = ex_to<numeric>(a.seq[i].coeff);
if (exponent[i] == 0) {
// optimize away
} else if (exponent[i] == 1) {
// optimized
- term.push_back(r);
- factor = factor.mul(ex_to<numeric>(c));
+ monomial.push_back(r);
+ if (c != *_num1_p)
+ factor = factor.mul(c);
} else { // general case exponent[i] > 1
- term.push_back((new power(r, exponent[i]))->setflag(status_flags::dynallocated));
- factor = factor.mul(ex_to<numeric>(c).power(exponent[i]));
+ monomial.push_back((new power(r, exponent[i]))->setflag(status_flags::dynallocated));
+ if (c != *_num1_p)
+ factor = factor.mul(c.power(exponent[i]));
}
}
- result.push_back(a.combine_ex_with_coeff_to_pair(mul(term).expand(options), factor));
+ result.push_back(a.combine_ex_with_coeff_to_pair(mul(monomial).expand(options), factor));
} while (compositions.next());
} while (partitions.next());
}
git://www.ginac.de / ginac.git / blobdiff