From cc308418b3e239e05a7f1f0f6cfc5919672fccfd Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Thu, 5 Nov 2015 13:30:48 +0100 Subject: [PATCH] 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. --- check/Makefile.am | 6 ++- check/exam_real_imag.cpp | 82 +++++++++++++++++++++++++++++++++++++ ginac/power.cpp | 88 ++++++++++++++++++++++++---------------- 3 files changed, 141 insertions(+), 35 deletions(-) create mode 100644 check/exam_real_imag.cpp diff --git a/check/Makefile.am b/check/Makefile.am index 2c1b5abd..0f7e035e 100644 --- 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 \ - exam_cra + exam_cra \ + exam_real_imag 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_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 diff --git a/check/exam_real_imag.cpp b/check/exam_real_imag.cpp new file mode 100644 index 00000000..548321c8 --- /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 +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(evalf(r)).to_double() - 0.2000570104163233; + double i_eps = ex_to(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 f951e4fd..b2622d60 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(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(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))); } -- 2.44.0