* Analytic continuation for atan2(y,x) for complex y and x.

diff --git a/NEWS b/NEWS
index d92cb439b2034d04390c007f58da48a2a9ed6a24..2345f5093b5cae5aa06e7a866c9a5be1a82bc95c 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -3,6 +3,7 @@ This file records noteworthy changes.
 1.3.8 (<date>)
 
 * Drop support of ginac-config and ginac.m4. Please use pkg-config instead.
+* atan2(y,x) branch cut correction.
 
 1.3.7 (05 February 2007)
 * Fixed bug in expansion of power.

diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp
index bf1d99ebc1f7e6ddaff04aeeac31e32afcb60e9c..7f3434621beaf1d218e64b465aa035378d7da410 100644 (file)
--- a/ginac/inifcns_trans.cpp
+++ b/ginac/inifcns_trans.cpp
@@ -770,9 +770,9 @@ static ex atan2_eval(const ex & y, const ex & x)
                        if (x.info(info_flags::positive))
                                return _ex0;
 
-                       // atan(0, x), x real and negative -> -Pi
+                       // atan(0, x), x real and negative -> Pi
                        if (x.info(info_flags::negative))
-                               return _ex_1*Pi;
+                               return Pi;
                }
 
                if (x.is_zero()) {
@@ -816,7 +816,7 @@ static ex atan2_eval(const ex & y, const ex & x)
                if (y.info(info_flags::real) && x.info(info_flags::real)) {
                        if (x.info(info_flags::positive))
                                return atan(y/x);
-                       else if(y.info(info_flags::positive))
+                       else if (y.info(info_flags::positive))
                                return atan(y/x)+Pi;
                        else
                                return atan(y/x)-Pi;

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index 9a567f929ead97201f5a7e018b827f6fd3d5c39c..81f285b5bbbce7f9cf6baf02ac0845744009450f 100644 (file)
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -1370,11 +1370,11 @@ const numeric acos(const numeric &x)
 }
        
 
-/** Arcustangent.
+/** Numeric arcustangent.
  *
  *  @param x complex number
  *  @return atan(x)
- *  @exception pole_error("atan(): logarithmic pole",0) */
+ *  @exception pole_error("atan(): logarithmic pole",0) if x==I or x==-I. */
 const numeric atan(const numeric &x)
 {
        if (!x.is_real() &&
@@ -1385,18 +1385,36 @@ const numeric atan(const numeric &x)
 }
 
 
-/** Arcustangent.
+/** Numeric arcustangent of two arguments, analytically continued in a suitable way.
  *
- *  @param x real number
- *  @param y real number
- *  @return atan(y/x) */
+ *  @param y complex number
+ *  @param x complex number
+ *  @return -I*log((x+I*y)/sqrt(x^2+y^2)), which is equal to atan(y/x) if y and
+ *    x are both real.
+ *  @exception pole_error("atan(): logarithmic pole",0) if y/x==+I or y/x==-I. */
 const numeric atan(const numeric &y, const numeric &x)
 {
+       if (x.is_zero() && y.is_zero())
+               return *_num0_p;
        if (x.is_real() && y.is_real())
                return cln::atan(cln::the<cln::cl_R>(x.to_cl_N()),
                                 cln::the<cln::cl_R>(y.to_cl_N()));
-       else
-               throw std::invalid_argument("atan(): complex argument");        
+
+       // Compute -I*log((x+I*y)/sqrt(x^2+y^2))
+       //      == -I*log((x+I*y)/sqrt((x+I*y)*(x-I*y)))
+       // Do not "simplify" this to -I/2*log((x+I*y)/(x-I*y))) or likewise.
+       // The branch cuts are easily messed up.
+       const cln::cl_N aux_p = x.to_cl_N()+cln::complex(0,1)*y.to_cl_N();
+       if (cln::zerop(aux_p)) {
+               // x+I*y==0 => y/x==I, so this is a pole (we have x!=0).
+               throw pole_error("atan(): logarithmic pole",0);
+       }
+       const cln::cl_N aux_m = x.to_cl_N()-cln::complex(0,1)*y.to_cl_N();
+       if (cln::zerop(aux_m)) {
+               // x-I*y==0 => y/x==-I, so this is a pole (we have x!=0).
+               throw pole_error("atan(): logarithmic pole",0);
+       }
+       return cln::complex(0,-1)*cln::log(aux_p/cln::sqrt(aux_p*aux_m));
 }