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()) {
}
// 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;
}
-/** 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() &&
}
-/** 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));
}