Add some minor optimizations.

[ginac.git] / ginac / power.cpp

diff --git a/ginac/power.cpp b/ginac/power.cpp
index 99a5a3c4e635273d3eb1dd05da15e760eb4d70a8..28e5bd3a82acbef059c258882d22d44616319ede 100644 (file)
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -695,51 +695,71 @@ ex power::conjugate() const
 
 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)));
 }
 
@@ -1195,7 +1215,6 @@ ex power::expand_add(const add & a, int n, unsigned options) const
                                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) ||
@@ -1203,15 +1222,19 @@ ex power::expand_add(const add & a, int n, unsigned options) const
                                                     !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));
+                                               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]));
+                                               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));