X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fpower.cpp;h=bf72d97320b667d1bcee62cf6572adee24845678;hp=99a5a3c4e635273d3eb1dd05da15e760eb4d70a8;hb=6cac49558b75dce07f607e26ba74aa9148f92720;hpb=1fa7d19af08892800dfdb0b700bc0a9362a3fae8;ds=sidebyside

diff --git a/ginac/power.cpp b/ginac/power.cpp
index 99a5a3c4..bf72d973 100644
--- 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)));
 }