X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=ef2784172107f22965381b271e2fafcd92d7bbca;hp=886fac72ccba04eca9ab78582750cd085c7dd984;hb=8d16780dbd9f4da9397e638aca213745589818c0;hpb=66c0f31c678e6c1938d637636b230ea376c157c1

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index 886fac72..ef278417 100644
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -4,8 +4,9 @@
  *  Its most important design principle is to completely hide the inner
  *  working of that other package from the user of GiNaC.  It must either 
  *  provide implementation of arithmetic operators and numerical evaluation
- *  of special functions or implement the interface to the bignum package.
- *
+ *  of special functions or implement the interface to the bignum package. */
+
+/*
  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
@@ -29,6 +30,7 @@
 #include "numeric.h"
 #include "ex.h"
 #include "config.h"
+#include "debugmsg.h"
 
 // CLN should not pollute the global namespace, hence we include it here
 // instead of in some header file where it would propagate to other parts:
@@ -38,6 +40,8 @@
 #include <cln.h>
 #endif
 
+namespace GiNaC {
+
 // linker has no problems finding text symbols for numerator or denominator
 //#define SANE_LINKER
 
@@ -220,7 +224,7 @@ void numeric::printraw(ostream & os) const
 void numeric::print(ostream & os, unsigned upper_precedence) const
 {
     debugmsg("numeric print", LOGLEVEL_PRINT);
-    if (is_real()) {  
+    if (is_real()) {
         // case 1, real:  x  or  -x
         if ((precedence<=upper_precedence) && (!is_pos_integer())) {
             os << "(" << *value << ")";
@@ -511,10 +515,35 @@ numeric const & numeric::operator=(char const * s)
     return operator=(numeric(s));
 }
 
+/** Return the complex half-plane (left or right) in which the number lies.
+ *  csgn(x)==0 for x==0, csgn(x)==1 for Re(x)>0 or Re(x)=0 and Im(x)>0,
+ *  csgn(x)==-1 for Re(x)<0 or Re(x)=0 and Im(x)<0.
+ *
+ *  @see numeric::compare(numeric const & other) */
+int numeric::csgn(void) const
+{
+    if (is_zero())
+        return 0;
+    if (!zerop(realpart(*value))) {
+        if (plusp(realpart(*value)))
+            return 1;
+        else
+            return -1;
+    } else {
+        if (plusp(imagpart(*value)))
+            return 1;
+        else
+            return -1;
+    }
+}
+
 /** This method establishes a canonical order on all numbers.  For complex
  *  numbers this is not possible in a mathematically consistent way but we need
  *  to establish some order and it ought to be fast.  So we simply define it
- *  similar to Maple's csgn. */
+ *  to be compatible with our method csgn.
+ *
+ *  @return csgn(*this-other)
+ *  @see numeric::csgn(void) */
 int numeric::compare(numeric const & other) const
 {
     // Comparing two real numbers?
@@ -850,7 +879,7 @@ const numeric some_numeric;
 type_info const & typeid_numeric=typeid(some_numeric);
 /** Imaginary unit.  This is not a constant but a numeric since we are
  *  natively handing complex numbers anyways. */
-const numeric I = (complex(cl_I(0),cl_I(1)));
+const numeric I = numeric(complex(cl_I(0),cl_I(1)));
 
 //////////
 // global functions
@@ -903,7 +932,7 @@ numeric const & numHALF(void)
  *  @return  arbitrary precision numerical exp(x). */
 numeric exp(numeric const & x)
 {
-    return exp(*x.value);  // -> CLN
+    return ::exp(*x.value);  // -> CLN
 }
 
 /** Natural logarithm.
@@ -915,7 +944,7 @@ numeric log(numeric const & z)
 {
     if (z.is_zero())
         throw (std::overflow_error("log(): logarithmic singularity"));
-    return log(*z.value);  // -> CLN
+    return ::log(*z.value);  // -> CLN
 }
 
 /** Numeric sine (trigonometric function).
@@ -923,7 +952,7 @@ numeric log(numeric const & z)
  *  @return  arbitrary precision numerical sin(x). */
 numeric sin(numeric const & x)
 {
-    return sin(*x.value);  // -> CLN
+    return ::sin(*x.value);  // -> CLN
 }
 
 /** Numeric cosine (trigonometric function).
@@ -931,7 +960,7 @@ numeric sin(numeric const & x)
  *  @return  arbitrary precision numerical cos(x). */
 numeric cos(numeric const & x)
 {
-    return cos(*x.value);  // -> CLN
+    return ::cos(*x.value);  // -> CLN
 }
     
 /** Numeric tangent (trigonometric function).
@@ -939,7 +968,7 @@ numeric cos(numeric const & x)
  *  @return  arbitrary precision numerical tan(x). */
 numeric tan(numeric const & x)
 {
-    return tan(*x.value);  // -> CLN
+    return ::tan(*x.value);  // -> CLN
 }
     
 /** Numeric inverse sine (trigonometric function).
@@ -947,7 +976,7 @@ numeric tan(numeric const & x)
  *  @return  arbitrary precision numerical asin(x). */
 numeric asin(numeric const & x)
 {
-    return asin(*x.value);  // -> CLN
+    return ::asin(*x.value);  // -> CLN
 }
     
 /** Numeric inverse cosine (trigonometric function).
@@ -955,7 +984,7 @@ numeric asin(numeric const & x)
  *  @return  arbitrary precision numerical acos(x). */
 numeric acos(numeric const & x)
 {
-    return acos(*x.value);  // -> CLN
+    return ::acos(*x.value);  // -> CLN
 }
     
 /** Arcustangents.
@@ -969,7 +998,7 @@ numeric atan(numeric const & x)
         x.real().is_zero() &&
         !abs(x.imag()).is_equal(numONE()))
         throw (std::overflow_error("atan(): logarithmic singularity"));
-    return atan(*x.value);  // -> CLN
+    return ::atan(*x.value);  // -> CLN
 }
 
 /** Arcustangents.
@@ -980,7 +1009,7 @@ numeric atan(numeric const & x)
 numeric atan(numeric const & y, numeric const & x)
 {
     if (x.is_real() && y.is_real())
-        return atan(realpart(*x.value), realpart(*y.value));  // -> CLN
+        return ::atan(realpart(*x.value), realpart(*y.value));  // -> CLN
     else
         throw (std::invalid_argument("numeric::atan(): complex argument"));        
 }
@@ -990,7 +1019,7 @@ numeric atan(numeric const & y, numeric const & x)
  *  @return  arbitrary precision numerical sinh(x). */
 numeric sinh(numeric const & x)
 {
-    return sinh(*x.value);  // -> CLN
+    return ::sinh(*x.value);  // -> CLN
 }
 
 /** Numeric hyperbolic cosine (trigonometric function).
@@ -998,7 +1027,7 @@ numeric sinh(numeric const & x)
  *  @return  arbitrary precision numerical cosh(x). */
 numeric cosh(numeric const & x)
 {
-    return cosh(*x.value);  // -> CLN
+    return ::cosh(*x.value);  // -> CLN
 }
     
 /** Numeric hyperbolic tangent (trigonometric function).
@@ -1006,7 +1035,7 @@ numeric cosh(numeric const & x)
  *  @return  arbitrary precision numerical tanh(x). */
 numeric tanh(numeric const & x)
 {
-    return tanh(*x.value);  // -> CLN
+    return ::tanh(*x.value);  // -> CLN
 }
     
 /** Numeric inverse hyperbolic sine (trigonometric function).
@@ -1014,7 +1043,7 @@ numeric tanh(numeric const & x)
  *  @return  arbitrary precision numerical asinh(x). */
 numeric asinh(numeric const & x)
 {
-    return asinh(*x.value);  // -> CLN
+    return ::asinh(*x.value);  // -> CLN
 }
 
 /** Numeric inverse hyperbolic cosine (trigonometric function).
@@ -1022,7 +1051,7 @@ numeric asinh(numeric const & x)
  *  @return  arbitrary precision numerical acosh(x). */
 numeric acosh(numeric const & x)
 {
-    return acosh(*x.value);  // -> CLN
+    return ::acosh(*x.value);  // -> CLN
 }
 
 /** Numeric inverse hyperbolic tangent (trigonometric function).
@@ -1030,14 +1059,22 @@ numeric acosh(numeric const & x)
  *  @return  arbitrary precision numerical atanh(x). */
 numeric atanh(numeric const & x)
 {
-    return atanh(*x.value);  // -> CLN
+    return ::atanh(*x.value);  // -> CLN
 }
 
 /** The gamma function.
- *  stub stub stub stub stub stub! */
+ *  This is only a stub! */
 numeric gamma(numeric const & x)
 {
-    clog << "gamma(): Nobody expects the Spanish inquisition" << endl;
+    clog << "gamma(): Does anybody know good way to calculate this numerically?" << endl;
+    return numeric(0);
+}
+
+/** The psi function (aka polygamma function).
+ *  This is only a stub! */
+numeric psi(numeric const & n, numeric const & x)
+{
+    clog << "psi(): Does anybody know good way to calculate this numerically?" << endl;
     return numeric(0);
 }
 
@@ -1050,7 +1087,7 @@ numeric factorial(numeric const & nn)
         throw (std::range_error("numeric::factorial(): argument must be integer >= 0"));
     }
     
-    return numeric(factorial(nn.to_int()));  // -> CLN
+    return numeric(::factorial(nn.to_int()));  // -> CLN
 }
 
 /** The double factorial combinatorial function.  (Scarcely used, but still
@@ -1061,6 +1098,10 @@ numeric factorial(numeric const & nn)
  *  @exception range_error (argument must be integer >= -1) */
 numeric doublefactorial(numeric const & nn)
 {
+    // META-NOTE:  The whole shit here will become obsolete and may be moved
+    // out once CLN learns about double factorial, which should be as soon as
+    // 1.0.3 rolls out.
+    
     // We store the results separately for even and odd arguments.  This has
     // the advantage that we don't have to compute any even result at all if
     // the function is always called with odd arguments and vice versa.  There
@@ -1122,7 +1163,7 @@ numeric doublefactorial(numeric const & nn)
 numeric binomial(numeric const & n, numeric const & k)
 {
     if (n.is_nonneg_integer() && k.is_nonneg_integer()) {
-        return numeric(binomial(n.to_int(),k.to_int()));  // -> CLN
+        return numeric(::binomial(n.to_int(),k.to_int()));  // -> CLN
     } else {
         // should really be gamma(n+1)/(gamma(r+1)/gamma(n-r+1)
         return numeric(0);
@@ -1133,7 +1174,7 @@ numeric binomial(numeric const & n, numeric const & k)
 /** Absolute value. */
 numeric abs(numeric const & x)
 {
-    return abs(*x.value);  // -> CLN
+    return ::abs(*x.value);  // -> CLN
 }
 
 /** Modulus (in positive representation).
@@ -1146,7 +1187,7 @@ numeric abs(numeric const & x)
 numeric mod(numeric const & a, numeric const & b)
 {
     if (a.is_integer() && b.is_integer()) {
-        return mod(The(cl_I)(*a.value), The(cl_I)(*b.value));  // -> CLN
+        return ::mod(The(cl_I)(*a.value), The(cl_I)(*b.value));  // -> CLN
     }
     else {
         return numZERO();  // Throw?
@@ -1161,7 +1202,7 @@ numeric smod(numeric const & a, numeric const & b)
 {
     if (a.is_integer() && b.is_integer()) {
         cl_I b2 = The(cl_I)(ceiling1(The(cl_I)(*b.value) / 2)) - 1;
-        return mod(The(cl_I)(*a.value) + b2, The(cl_I)(*b.value)) - b2;
+        return ::mod(The(cl_I)(*a.value) + b2, The(cl_I)(*b.value)) - b2;
     } else {
         return numZERO();  // Throw?
     }
@@ -1176,7 +1217,7 @@ numeric smod(numeric const & a, numeric const & b)
 numeric irem(numeric const & a, numeric const & b)
 {
     if (a.is_integer() && b.is_integer()) {
-        return rem(The(cl_I)(*a.value), The(cl_I)(*b.value));  // -> CLN
+        return ::rem(The(cl_I)(*a.value), The(cl_I)(*b.value));  // -> CLN
     }
     else {
         return numZERO();  // Throw?
@@ -1244,7 +1285,7 @@ numeric iquo(numeric const & a, numeric const & b, numeric & r)
  *  where imag(z)>0. */
 numeric sqrt(numeric const & z)
 {
-    return sqrt(*z.value);  // -> CLN
+    return ::sqrt(*z.value);  // -> CLN
 }
 
 /** Integer numeric square root. */
@@ -1252,7 +1293,7 @@ numeric isqrt(numeric const & x)
 {
 	if (x.is_integer()) {
 		cl_I root;
-		isqrt(The(cl_I)(*x.value), &root);	// -> CLN
+		::isqrt(The(cl_I)(*x.value), &root);	// -> CLN
 		return root;
 	} else
 		return numZERO();  // Throw?
@@ -1265,7 +1306,7 @@ numeric isqrt(numeric const & x)
 numeric gcd(numeric const & a, numeric const & b)
 {
     if (a.is_integer() && b.is_integer())
-        return gcd(The(cl_I)(*a.value), The(cl_I)(*b.value));	// -> CLN
+        return ::gcd(The(cl_I)(*a.value), The(cl_I)(*b.value));	// -> CLN
     else
         return numONE();
 }
@@ -1277,7 +1318,7 @@ numeric gcd(numeric const & a, numeric const & b)
 numeric lcm(numeric const & a, numeric const & b)
 {
     if (a.is_integer() && b.is_integer())
-        return lcm(The(cl_I)(*a.value), The(cl_I)(*b.value));	// -> CLN
+        return ::lcm(The(cl_I)(*a.value), The(cl_I)(*b.value));	// -> CLN
     else
         return *a.value * *b.value;
 }
@@ -1343,3 +1384,5 @@ bool _numeric_digits::too_late = false;
 /** Accuracy in decimal digits.  Only object of this type!  Can be set using
  *  assignment from C++ unsigned ints and evaluated like any built-in type. */
 _numeric_digits Digits;
+
+} // namespace GiNaC