- 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
+ }
+ }