}
}
+template<typename T1, typename T2>
+static inline bool coerce(T1& dst, const T2& arg);
+
+/**
+ * @brief Check if CLN integer can be converted into int
+ *
+ * @sa http://www.ginac.de/pipermail/cln-list/2006-October/000248.html
+ */
+template<>
+static inline bool coerce<int, cln::cl_I>(int& dst, const cln::cl_I& arg)
+{
+ static const cln::cl_I cl_max_int =
+ (cln::cl_I)(long)(std::numeric_limits<int>::max());
+ static const cln::cl_I cl_min_int =
+ (cln::cl_I)(long)(std::numeric_limits<int>::min());
+ if ((arg >= cl_min_int) && (arg <= cl_max_int)) {
+ dst = cl_I_to_int(arg);
+ return true;
+ }
+ return false;
+}
+
+template<>
+static inline bool coerce<unsigned int, cln::cl_I>(unsigned int& dst, const cln::cl_I& arg)
+{
+ static const cln::cl_I cl_max_uint =
+ (cln::cl_I)(unsigned long)(std::numeric_limits<unsigned int>::max());
+ if ((! minusp(arg)) && (arg <= cl_max_uint)) {
+ dst = cl_I_to_uint(arg);
+ return true;
+ }
+ return false;
+}
+
/** Helper function to print real number in C++ source format using cl_N types.
*
* @see numeric::print() */
{
if (cln::instanceof(x, cln::cl_I_ring)) {
- // Integer number
- c.s << "cln::cl_I(\"";
- print_real_number(c, x);
- c.s << "\")";
-
+ int dst;
+ // fixnum
+ if (coerce(dst, cln::the<cln::cl_I>(x))) {
+ // can be converted to native int
+ if (dst < 0)
+ c.s << "(-" << dst << ")";
+ else
+ c.s << dst;
+ } else {
+ // bignum
+ c.s << "cln::cl_I(\"";
+ print_real_number(c, x);
+ c.s << "\")";
+ }
} else if (cln::instanceof(x, cln::cl_RA_ring)) {
// Rational number
case info_flags::numeric:
case info_flags::polynomial:
case info_flags::rational_function:
+ case info_flags::expanded:
return true;
case info_flags::real:
return is_real();
}
-/** 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));
}