From: Richard Kreckel Date: Sat, 22 May 2010 20:34:42 +0000 (+0200) Subject: Be more careful with conjugate(f(x)) -> f(conjugate(x)). X-Git-Tag: release_1-5-8~4 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=5fb83676210401cc08fae91831514bac44502209 Be more careful with conjugate(f(x)) -> f(conjugate(x)). That identity is correct for holomorphic functions, but we have to be careful with some functions not being holomorphic everywhere. On branch cuts, it is wrong. For a discussion, see: . This patch changes the default behavior of user-defined functions. From now on, a user-defined conjugate_func must be registered, in order to evaluate conjugate(f(x)) -> f(conjugate(x)). This patch also adds such functions for most initially known functions. --- diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index a628ff17..b3befd3f 100644 --- a/doc/tutorial/ginac.texi +++ b/doc/tutorial/ginac.texi @@ -6149,12 +6149,13 @@ For example, @} @end example -If you declare your own GiNaC functions, then they will conjugate themselves -by conjugating their arguments. This is the default strategy. If you want to -change this behavior, you have to supply a specialized conjugation method -for your function (see @ref{Symbolic functions} and the GiNaC source-code -for @code{abs} as an example). Also, specialized methods can be provided -to take real and imaginary parts of user-defined functions. +If you declare your own GiNaC functions and you want to conjugate them, you +will have to supply a specialized conjugation method for them (see +@ref{Symbolic functions} and the GiNaC source-code for @code{abs} as an +example). GiNaC does not automatically conjugate user-supplied functions +by conjugating their arguments because this would be incorrect on branch +cuts. Also, specialized methods can be provided to take real and imaginary +parts of user-defined functions. @node Solving linear systems of equations, Input/output, Complex expressions, Methods and functions @c node-name, next, previous, up diff --git a/ginac/function.pl b/ginac/function.pl index aa55c9bc..1ca923d4 100644 --- a/ginac/function.pl +++ b/ginac/function.pl @@ -1124,7 +1124,7 @@ ex function::conjugate() const const function_options & opt = registered_functions()[serial]; if (opt.conjugate_f==0) { - return exprseq::conjugate(); + return conjugate_function(*this).hold(); } if (opt.conjugate_use_exvector_args) { diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index 63c2d74f..981e7b25 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -319,8 +319,8 @@ REGISTER_FUNCTION(step, eval_func(step_eval). evalf_func(step_evalf). series_func(step_series). conjugate_func(step_conjugate). - real_part_func(step_real_part). - imag_part_func(step_imag_part)); + real_part_func(step_real_part). + imag_part_func(step_imag_part)); ////////// // Complex sign @@ -643,10 +643,25 @@ static ex Li2_series(const ex &x, const relational &rel, int order, unsigned opt throw do_taylor(); // caught by function::series() } +static ex Li2_conjugate(const ex & x) +{ + // conjugate(Li2(x))==Li2(conjugate(x)) unless on the branch cuts which + // run along the positive real axis beginning at 1. + if (x.info(info_flags::negative)) { + return Li2(x); + } + if (is_exactly_a(x) && + (!x.imag_part().is_zero() || x < *_num1_p)) { + return Li2(x.conjugate()); + } + return conjugate_function(Li2(x)).hold(); +} + REGISTER_FUNCTION(Li2, eval_func(Li2_eval). evalf_func(Li2_evalf). derivative_func(Li2_deriv). series_func(Li2_series). + conjugate_func(Li2_conjugate). latex_name("\\mathrm{Li}_2")); ////////// diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp index 0fabbca7..46c8ea9e 100644 --- a/ginac/inifcns_gamma.cpp +++ b/ginac/inifcns_gamma.cpp @@ -110,10 +110,26 @@ static ex lgamma_series(const ex & arg, } +static ex lgamma_conjugate(const ex & x) +{ + // conjugate(lgamma(x))==lgamma(conjugate(x)) unless on the branch cut + // which runs along the negative real axis. + if (x.info(info_flags::positive)) { + return lgamma(x); + } + if (is_exactly_a(x) && + !x.imag_part().is_zero()) { + return lgamma(x.conjugate()); + } + return conjugate_function(lgamma(x)).hold(); +} + + REGISTER_FUNCTION(lgamma, eval_func(lgamma_eval). evalf_func(lgamma_evalf). derivative_func(lgamma_deriv). series_func(lgamma_series). + conjugate_func(lgamma_conjugate). latex_name("\\log \\Gamma")); @@ -207,10 +223,18 @@ static ex tgamma_series(const ex & arg, } +static ex tgamma_conjugate(const ex & x) +{ + // conjugate(tgamma(x))==tgamma(conjugate(x)) + return tgamma(x.conjugate()); +} + + REGISTER_FUNCTION(tgamma, eval_func(tgamma_eval). evalf_func(tgamma_evalf). derivative_func(tgamma_deriv). series_func(tgamma_series). + conjugate_func(tgamma_conjugate). latex_name("\\Gamma")); @@ -329,7 +353,7 @@ REGISTER_FUNCTION(beta, eval_func(beta_eval). derivative_func(beta_deriv). series_func(beta_series). latex_name("\\mathrm{B}"). - set_symmetry(sy_symm(0, 1))); + set_symmetry(sy_symm(0, 1))); ////////// diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp index c8a0a821..3884d79c 100644 --- a/ginac/inifcns_trans.cpp +++ b/ginac/inifcns_trans.cpp @@ -99,11 +99,18 @@ static ex exp_imag_part(const ex & x) return exp(GiNaC::real_part(x))*sin(GiNaC::imag_part(x)); } +static ex exp_conjugate(const ex & x) +{ + // conjugate(exp(x))==exp(conjugate(x)) + return exp(x.conjugate()); +} + REGISTER_FUNCTION(exp, eval_func(exp_eval). evalf_func(exp_evalf). derivative_func(exp_deriv). real_part_func(exp_real_part). imag_part_func(exp_imag_part). + conjugate_func(exp_conjugate). latex_name("\\exp")); ////////// @@ -258,12 +265,27 @@ static ex log_imag_part(const ex & x) return atan2(GiNaC::imag_part(x), GiNaC::real_part(x)); } +static ex log_conjugate(const ex & x) +{ + // conjugate(log(x))==log(conjugate(x)) unless on the branch cut which + // runs along the negative real axis. + if (x.info(info_flags::positive)) { + return log(x); + } + if (is_exactly_a(x) && + !x.imag_part().is_zero()) { + return log(x.conjugate()); + } + return conjugate_function(log(x)).hold(); +} + REGISTER_FUNCTION(log, eval_func(log_eval). evalf_func(log_evalf). derivative_func(log_deriv). series_func(log_series). real_part_func(log_real_part). imag_part_func(log_imag_part). + conjugate_func(log_conjugate). latex_name("\\ln")); ////////// @@ -359,11 +381,18 @@ static ex sin_imag_part(const ex & x) return sinh(GiNaC::imag_part(x))*cos(GiNaC::real_part(x)); } +static ex sin_conjugate(const ex & x) +{ + // conjugate(sin(x))==sin(conjugate(x)) + return sin(x.conjugate()); +} + REGISTER_FUNCTION(sin, eval_func(sin_eval). evalf_func(sin_evalf). derivative_func(sin_deriv). real_part_func(sin_real_part). imag_part_func(sin_imag_part). + conjugate_func(sin_conjugate). latex_name("\\sin")); ////////// @@ -459,11 +488,18 @@ static ex cos_imag_part(const ex & x) return -sinh(GiNaC::imag_part(x))*sin(GiNaC::real_part(x)); } +static ex cos_conjugate(const ex & x) +{ + // conjugate(cos(x))==cos(conjugate(x)) + return cos(x.conjugate()); +} + REGISTER_FUNCTION(cos, eval_func(cos_eval). evalf_func(cos_evalf). derivative_func(cos_deriv). real_part_func(cos_real_part). imag_part_func(cos_imag_part). + conjugate_func(cos_conjugate). latex_name("\\cos")); ////////// @@ -576,12 +612,19 @@ static ex tan_series(const ex &x, return (sin(x)/cos(x)).series(rel, order, options); } +static ex tan_conjugate(const ex & x) +{ + // conjugate(tan(x))==tan(conjugate(x)) + return tan(x.conjugate()); +} + REGISTER_FUNCTION(tan, eval_func(tan_eval). evalf_func(tan_evalf). derivative_func(tan_deriv). series_func(tan_series). real_part_func(tan_real_part). imag_part_func(tan_imag_part). + conjugate_func(tan_conjugate). latex_name("\\tan")); ////////// @@ -640,9 +683,21 @@ static ex asin_deriv(const ex & x, unsigned deriv_param) return power(1-power(x,_ex2),_ex_1_2); } +static ex asin_conjugate(const ex & x) +{ + // conjugate(asin(x))==asin(conjugate(x)) unless on the branch cuts which + // run along the real axis outside the interval [-1, +1]. + if (is_exactly_a(x) && + (!x.imag_part().is_zero() || (x > *_num_1_p && x < *_num1_p))) { + return asin(x.conjugate()); + } + return conjugate_function(asin(x)).hold(); +} + REGISTER_FUNCTION(asin, eval_func(asin_eval). evalf_func(asin_evalf). derivative_func(asin_deriv). + conjugate_func(asin_conjugate). latex_name("\\arcsin")); ////////// @@ -701,9 +756,21 @@ static ex acos_deriv(const ex & x, unsigned deriv_param) return -power(1-power(x,_ex2),_ex_1_2); } +static ex acos_conjugate(const ex & x) +{ + // conjugate(acos(x))==acos(conjugate(x)) unless on the branch cuts which + // run along the real axis outside the interval [-1, +1]. + if (is_exactly_a(x) && + (!x.imag_part().is_zero() || (x > *_num_1_p && x < *_num1_p))) { + return acos(x.conjugate()); + } + return conjugate_function(acos(x)).hold(); +} + REGISTER_FUNCTION(acos, eval_func(acos_eval). evalf_func(acos_evalf). derivative_func(acos_deriv). + conjugate_func(acos_conjugate). latex_name("\\arccos")); ////////// @@ -800,10 +867,27 @@ static ex atan_series(const ex &arg, throw do_taylor(); } +static ex atan_conjugate(const ex & x) +{ + // conjugate(atan(x))==atan(conjugate(x)) unless on the branch cuts which + // run along the imaginary axis outside the interval [-I, +I]. + if (x.info(info_flags::real)) + return atan(x); + if (is_exactly_a(x)) { + const numeric x_re = ex_to(x.real_part()); + const numeric x_im = ex_to(x.imag_part()); + if (!x_re.is_zero() || + (x_im > *_num_1_p && x_im < *_num1_p)) + return atan(x.conjugate()); + } + return conjugate_function(atan(x)).hold(); +} + REGISTER_FUNCTION(atan, eval_func(atan_eval). evalf_func(atan_evalf). derivative_func(atan_deriv). series_func(atan_series). + conjugate_func(atan_conjugate). latex_name("\\arctan")); ////////// @@ -1197,9 +1281,26 @@ static ex asinh_deriv(const ex & x, unsigned deriv_param) return power(_ex1+power(x,_ex2),_ex_1_2); } +static ex asinh_conjugate(const ex & x) +{ + // conjugate(asinh(x))==asinh(conjugate(x)) unless on the branch cuts which + // run along the imaginary axis outside the interval [-I, +I]. + if (x.info(info_flags::real)) + return asinh(x); + if (is_exactly_a(x)) { + const numeric x_re = ex_to(x.real_part()); + const numeric x_im = ex_to(x.imag_part()); + if (!x_re.is_zero() || + (x_im > *_num_1_p && x_im < *_num1_p)) + return asinh(x.conjugate()); + } + return conjugate_function(asinh(x)).hold(); +} + REGISTER_FUNCTION(asinh, eval_func(asinh_eval). evalf_func(asinh_evalf). - derivative_func(asinh_deriv)); + derivative_func(asinh_deriv). + conjugate_func(asinh_conjugate)); ////////// // inverse hyperbolic cosine (trigonometric function) @@ -1249,9 +1350,21 @@ static ex acosh_deriv(const ex & x, unsigned deriv_param) return power(x+_ex_1,_ex_1_2)*power(x+_ex1,_ex_1_2); } +static ex acosh_conjugate(const ex & x) +{ + // conjugate(acosh(x))==acosh(conjugate(x)) unless on the branch cut + // which runs along the real axis from +1 to -inf. + if (is_exactly_a(x) && + (!x.imag_part().is_zero() || x > *_num1_p)) { + return acosh(x.conjugate()); + } + return conjugate_function(acosh(x)).hold(); +} + REGISTER_FUNCTION(acosh, eval_func(acosh_eval). evalf_func(acosh_evalf). - derivative_func(acosh_deriv)); + derivative_func(acosh_deriv). + conjugate_func(acosh_conjugate)); ////////// // inverse hyperbolic tangent (trigonometric function) @@ -1340,10 +1453,22 @@ static ex atanh_series(const ex &arg, throw do_taylor(); } +static ex atanh_conjugate(const ex & x) +{ + // conjugate(atanh(x))==atanh(conjugate(x)) unless on the branch cuts which + // run along the real axis outside the interval [-1, +1]. + if (is_exactly_a(x) && + (!x.imag_part().is_zero() || (x > *_num_1_p && x < *_num1_p))) { + return atanh(x.conjugate()); + } + return conjugate_function(atanh(x)).hold(); +} + REGISTER_FUNCTION(atanh, eval_func(atanh_eval). evalf_func(atanh_evalf). derivative_func(atanh_deriv). - series_func(atanh_series)); + series_func(atanh_series). + conjugate_func(atanh_conjugate)); } // namespace GiNaC diff --git a/ginac/power.cpp b/ginac/power.cpp index 4349ab7a..d24c969a 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -667,12 +667,23 @@ ex power::eval_ncmul(const exvector & v) const ex power::conjugate() const { - ex newbasis = basis.conjugate(); - ex newexponent = exponent.conjugate(); - if (are_ex_trivially_equal(basis, newbasis) && are_ex_trivially_equal(exponent, newexponent)) { - return *this; + // conjugate(pow(x,y))==pow(conjugate(x),conjugate(y)) unless on the + // branch cut which runs along the negative real axis. + if (basis.info(info_flags::positive)) { + ex newexponent = exponent.conjugate(); + if (are_ex_trivially_equal(exponent, newexponent)) { + return *this; + } + return (new power(basis, newexponent))->setflag(status_flags::dynallocated); + } + if (exponent.info(info_flags::integer)) { + ex newbasis = basis.conjugate(); + if (are_ex_trivially_equal(basis, newbasis)) { + return *this; + } + return (new power(newbasis, exponent))->setflag(status_flags::dynallocated); } - return (new power(newbasis, newexponent))->setflag(status_flags::dynallocated); + return conjugate_function(*this).hold(); } ex power::real_part() const