Added methods for taking real and imaginary parts.

[ginac.git] / ginac / power.cpp

diff --git a/ginac/power.cpp b/ginac/power.cpp
index 95a01acec2b5d31c4312ced86020990d99dc4fd7..75fd2ff51a9d12bd09d867d30fdf50b860ec6240 100644 (file)
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -40,6 +40,7 @@
 #include "lst.h"
 #include "archive.h"
 #include "utils.h"
+#include "relational.h"
 
 namespace GiNaC {
 
@@ -614,13 +615,66 @@ ex power::conjugate() const
        return (new power(newbasis, newexponent))->setflag(status_flags::dynallocated);
 }
 
+ex power::real_part() const
+{
+       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() ));
+               return result;
+       }
+       
+       ex a = basis.real_part();
+       ex b = basis.imag_part();
+       ex c = exponent.real_part();
+       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
+{
+       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() ));
+               return result;
+       }
+       
+       ex a=basis.real_part();
+       ex b=basis.imag_part();
+       ex c=exponent.real_part();
+       ex d=exponent.imag_part();
+       return
+               power(abs(basis),c)*exp(-d*atan2(b,a))*sin(c*atan2(b,a)+d*log(abs(basis)));
+}
+
+// protected
+
 // protected
 
 /** Implementation of ex::diff() for a power.
  *  @see ex::diff */
 ex power::derivative(const symbol & s) const
 {
-       if (exponent.info(info_flags::real)) {
+       if (is_a<numeric>(exponent)) {
                // D(b^r) = r * b^(r-1) * D(b) (faster than the formula below)
                epvector newseq;
                newseq.reserve(2);