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:
<http://www.ginac.de/pipermail/ginac-list/2010-April/001601.html>.

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 a628ff17a4c2cc9c6e420d039417863af5f8b614..b3befd3fb7d170510ca7f20f501e9f3be01e546a 100644 (file)
--- 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 a866b7e35f557a8deb06da4b39cf47dd9fc3274b..07e82e19d2578c73e6c10e3e9cf76caa6a60e797 100644 (file)
--- a/ginac/function.pl
+++ b/ginac/function.pl
@@ -1122,7 +1122,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 63c2d74f556d74c2c9cce124729107695a3bdf43..981e7b25f0bfda83d1d9ad542966ce4512b3ea64 100644 (file)
--- 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<numeric>(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 0fabbca7157a7bfc5363a203ebbdd4a3fc7fd362..46c8ea9e00e717652b188a0979f18224bc639fa0 100644 (file)
--- 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<numeric>(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 c8a0a8214ca7f276da625096c36899cdd1be9d78..3884d79c5275abda52f8917901769739c8ae9e64 100644 (file)
--- 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<numeric>(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<numeric>(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<numeric>(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<numeric>(x)) {
+               const numeric x_re = ex_to<numeric>(x.real_part());
+               const numeric x_im = ex_to<numeric>(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<numeric>(x)) {
+               const numeric x_re = ex_to<numeric>(x.real_part());
+               const numeric x_im = ex_to<numeric>(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<numeric>(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<numeric>(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 4349ab7a59a47959d0d600607e7148dc13084a74..d24c969a2145d5a755d574c475a4ea57fd05d139 100644 (file)
--- 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