From f9ce2b39c8429897d9e29087828df68260659290 Mon Sep 17 00:00:00 2001
From: Richard Kreckel <Richard.Kreckel@uni-mainz.de>
Date: Wed, 25 Jul 2007 20:59:02 +0000
Subject: [PATCH] * Analytic continuation for atan2(y,x) for complex y and x.

---
 ginac/inifcns_trans.cpp | 10 +++++-----
 ginac/numeric.cpp       | 34 ++++++++++++++++++++++++++--------
 2 files changed, 31 insertions(+), 13 deletions(-)

diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp
index f0d78590..90380efe 100644
--- a/ginac/inifcns_trans.cpp
+++ b/ginac/inifcns_trans.cpp
@@ -830,9 +830,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()) {
@@ -869,15 +869,15 @@ static ex atan2_eval(const ex & y, const ex & x)
 	}
 
 	// atan(float, float) -> float
-	if (is_a<numeric>(y) && is_a<numeric>(x) && !y.info(info_flags::crational)
-			&& !x.info(info_flags::crational))
+	if (is_a<numeric>(y) && !y.info(info_flags::crational) &&
+	    is_a<numeric>(x) && !x.info(info_flags::crational))
 		return atan(ex_to<numeric>(y), ex_to<numeric>(x));
 
 	// atan(real, real) -> atan(y/x) +/- Pi
 	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 e4b1f229..1858cbd9 100644
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -1397,11 +1397,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() &&
@@ -1412,18 +1412,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));
 }
 
 
-- 
2.49.0