Improved CLN output [Sheplyakov].

[ginac.git] / ginac / numeric.cpp

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index e4b1f22942d7e4e86ec242bcd2fcd2a2b1d719ca..3ca4b58c5ef31487c1e7a654a47885e67df7c4aa 100644 (file)
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -399,6 +399,40 @@ static void print_real_csrc(const print_context & c, const cln::cl_R & x)
        }
 }
 
        }
 }
 

+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() */
 /** Helper function to print real number in C++ source format using cl_N types.
  *
  *  @see numeric::print() */


@@ -406,11 +440,20 @@ static void print_real_cl_N(const print_context & c, const cln::cl_R & x)
 {
        if (cln::instanceof(x, cln::cl_I_ring)) {
 
 {
        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
        } else if (cln::instanceof(x, cln::cl_RA_ring)) {
 
                // Rational number


@@ -576,6 +619,7 @@ bool numeric::info(unsigned inf) const
                case info_flags::numeric:
                case info_flags::polynomial:
                case info_flags::rational_function:
                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();
                        return true;
                case info_flags::real:
                        return is_real();


@@ -1397,11 +1441,11 @@ const numeric acos(const numeric &x)
 }
        
 
 }
        
 

-/** Arcustangent.
+/** Numeric arcustangent.

  *
  *  @param x complex number
  *  @return atan(x)
  *
  *  @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() &&
 const numeric atan(const numeric &x)
 {
        if (!x.is_real() &&


@@ -1412,18 +1456,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)
 {
 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()));
        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));

 }
 
 
 }