Speed up power::real_part() and power::imag_part().

Add special case for real base and exponent (not just integer exponent).
Improve special case for integer exponent by explicitly constructing
the result using the Binomial expansion.
Add a test case for real/imaginary part expansion.

diff --git a/check/Makefile.am b/check/Makefile.am
index 2c1b5abd374d3bd6373d2db50ac22305179a9d9e..0f7e035e0d037dc8e10f7dc3256eafef01a615ee 100644 (file)
--- a/check/Makefile.am
+++ b/check/Makefile.am
@@ -34,7 +34,8 @@ EXAMS = exam_paranoia \
        factor_univariate_bug \
        pgcd_relatively_prime_bug \
        pgcd_infinite_loop \
        factor_univariate_bug \
        pgcd_relatively_prime_bug \
        pgcd_infinite_loop \

-       exam_cra
+       exam_cra \
+       exam_real_imag

 
 TIMES = time_dennyfliegner \
        time_gammaseries \
 
 TIMES = time_dennyfliegner \
        time_gammaseries \


@@ -157,6 +158,9 @@ exam_mod_gcd_LDADD = ../ginac/libginac.la
 exam_cra_SOURCES = exam_cra.cpp
 exam_cra_LDADD = ../ginac/libginac.la
 
 exam_cra_SOURCES = exam_cra.cpp
 exam_cra_LDADD = ../ginac/libginac.la
 

+exam_real_imag_SOURCES = exam_real_imag.cpp
+exam_real_imag_LDADD = ../ginac/libginac.la
+

 time_dennyfliegner_SOURCES = time_dennyfliegner.cpp \
                             randomize_serials.cpp timer.cpp timer.h
 time_dennyfliegner_LDADD = ../ginac/libginac.la
 time_dennyfliegner_SOURCES = time_dennyfliegner.cpp \
                             randomize_serials.cpp timer.cpp timer.h
 time_dennyfliegner_LDADD = ../ginac/libginac.la

diff --git a/check/exam_real_imag.cpp b/check/exam_real_imag.cpp
new file mode 100644 (file)
index 0000000..548321c
--- /dev/null
+++ b/check/exam_real_imag.cpp
@@ -0,0 +1,82 @@
+/** @file exam_misc.cpp
+ *
+ */
+
+/*
+ *  GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+#include "ginac.h"
+using namespace GiNaC;
+
+#include <iostream>
+using namespace std;
+
+/* Exam real/imaginary part of polynomials. */
+static unsigned exam_polynomials()
+{
+       realsymbol a("a"), b("b");
+       ex e = pow(a + I*b,3).expand() + I;
+
+       if (e.real_part() != pow(a,3)-3*a*pow(b,2) ||
+           e.imag_part() != 1+3*pow(a,2)*b-pow(b,3)) {
+               clog << "real / imaginary part miscomputed" << endl;
+               return 1;
+       }
+       return 0;
+}
+
+/* Exam symbolic expansion of nested expression. */
+static unsigned exam_monster()
+{
+       // This little monster is inspired by Sage's Symbench R1.
+       // It is much more aggressive that the original and covers more code.
+       struct {  // C++ doesn't have nested functions...
+               ex operator()(const ex & z) {
+                       return sqrt(ex(1)/3) * pow(z, 11) - I / pow(2*z, 3);
+               }
+       } f;
+       ex monster = f(f(f(f(I/2))));  // grows exponentially with number of nestings..
+       ex r = real_part(monster);
+       ex i = imag_part(monster);
+
+       // Check with precomputed result
+       double r_eps = ex_to<numeric>(evalf(r)).to_double() - 0.2000570104163233;
+       double i_eps = ex_to<numeric>(evalf(i)).to_double() - 0.5284320312415462;
+       if (abs(r_eps) > 1e-9  ||  abs(i_eps) > 1e-9) {
+               clog << "iterated function was miscomputed" << endl;
+               return 1;
+       }
+       return 0;
+}
+
+unsigned exam_real_imag()
+{
+       unsigned result = 0;
+
+       cout << "examining real/imaginary part separation" << flush;
+
+       result += exam_polynomials(); cout << '.' << flush;
+       result += exam_monster(); cout << '.' << flush;
+
+       return result;
+}
+
+int main(int argc, char** argv)
+{
+       return exam_real_imag();
+}

diff --git a/ginac/power.cpp b/ginac/power.cpp
index 99a5a3c4e635273d3eb1dd05da15e760eb4d70a8..bf72d97320b667d1bcee62cf6572adee24845678 100644 (file)
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -695,51 +695,71 @@ ex power::conjugate() const
 
 ex power::real_part() 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)) {
        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;
        }
                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
 {
        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)) {
        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;
        }
                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)));
 }
 
        return power(abs(basis),c)*exp(-d*atan2(b,a))*sin(c*atan2(b,a)+d*log(abs(basis)));
 }