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

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index c189c3d7..ef278417 100644
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -6,10 +6,31 @@
  *  provide implementation of arithmetic operators and numerical evaluation
  *  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
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
 #include <vector>
 #include <stdexcept>
 
-#include "ginac.h"
+#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:
@@ -19,6 +40,8 @@
 #include <cln.h>
 #endif
 
+namespace GiNaC {
+
 // linker has no problems finding text symbols for numerator or denominator
 //#define SANE_LINKER
 
@@ -30,7 +53,7 @@
 // public
 
 /** default ctor. Numerically it initializes to an integer zero. */
-numeric::numeric() : basic(TINFO_NUMERIC)
+numeric::numeric() : basic(TINFO_numeric)
 {
     debugmsg("numeric default constructor", LOGLEVEL_CONSTRUCT);
     value = new cl_N;
@@ -82,7 +105,7 @@ void numeric::destroy(bool call_parent)
 
 // public
 
-numeric::numeric(int i) : basic(TINFO_NUMERIC)
+numeric::numeric(int i) : basic(TINFO_numeric)
 {
     debugmsg("numeric constructor from int",LOGLEVEL_CONSTRUCT);
     // Not the whole int-range is available if we don't cast to long
@@ -94,7 +117,7 @@ numeric::numeric(int i) : basic(TINFO_NUMERIC)
             status_flags::hash_calculated);
 }
 
-numeric::numeric(unsigned int i) : basic(TINFO_NUMERIC)
+numeric::numeric(unsigned int i) : basic(TINFO_numeric)
 {
     debugmsg("numeric constructor from uint",LOGLEVEL_CONSTRUCT);
     // Not the whole uint-range is available if we don't cast to ulong
@@ -106,7 +129,7 @@ numeric::numeric(unsigned int i) : basic(TINFO_NUMERIC)
             status_flags::hash_calculated);
 }
 
-numeric::numeric(long i) : basic(TINFO_NUMERIC)
+numeric::numeric(long i) : basic(TINFO_numeric)
 {
     debugmsg("numeric constructor from long",LOGLEVEL_CONSTRUCT);
     value = new cl_I(i);
@@ -115,7 +138,7 @@ numeric::numeric(long i) : basic(TINFO_NUMERIC)
             status_flags::hash_calculated);
 }
 
-numeric::numeric(unsigned long i) : basic(TINFO_NUMERIC)
+numeric::numeric(unsigned long i) : basic(TINFO_numeric)
 {
     debugmsg("numeric constructor from ulong",LOGLEVEL_CONSTRUCT);
     value = new cl_I(i);
@@ -127,7 +150,7 @@ numeric::numeric(unsigned long i) : basic(TINFO_NUMERIC)
 /** Ctor for rational numerics a/b.
  *
  *  @exception overflow_error (division by zero) */
-numeric::numeric(long numer, long denom) : basic(TINFO_NUMERIC)
+numeric::numeric(long numer, long denom) : basic(TINFO_numeric)
 {
     debugmsg("numeric constructor from long/long",LOGLEVEL_CONSTRUCT);
     if (!denom)
@@ -139,7 +162,7 @@ numeric::numeric(long numer, long denom) : basic(TINFO_NUMERIC)
             status_flags::hash_calculated);
 }
 
-numeric::numeric(double d) : basic(TINFO_NUMERIC)
+numeric::numeric(double d) : basic(TINFO_numeric)
 {
     debugmsg("numeric constructor from double",LOGLEVEL_CONSTRUCT);
     // We really want to explicitly use the type cl_LF instead of the
@@ -152,7 +175,7 @@ numeric::numeric(double d) : basic(TINFO_NUMERIC)
             status_flags::hash_calculated);
 }
 
-numeric::numeric(char const *s) : basic(TINFO_NUMERIC)
+numeric::numeric(char const *s) : basic(TINFO_numeric)
 {   // MISSING: treatment of complex and ints and rationals.
     debugmsg("numeric constructor from string",LOGLEVEL_CONSTRUCT);
     if (strchr(s, '.'))
@@ -166,7 +189,7 @@ numeric::numeric(char const *s) : basic(TINFO_NUMERIC)
 
 /** Ctor from CLN types.  This is for the initiated user or internal use
  *  only. */
-numeric::numeric(cl_N const & z) : basic(TINFO_NUMERIC)
+numeric::numeric(cl_N const & z) : basic(TINFO_numeric)
 {
     debugmsg("numeric constructor from cl_N", LOGLEVEL_CONSTRUCT);
     value = new cl_N(z);
@@ -197,13 +220,13 @@ void numeric::printraw(ostream & os) const
 }
 
 // The method print adds to the output so it blends more consistently together
-// with the other routines.
+// with the other routines and produces something compatible to Maple input.
 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 ((realpart(*value) < 0) && (precedence <= upper_precedence)) {
+        if ((precedence<=upper_precedence) && (!is_pos_integer())) {
             os << "(" << *value << ")";
         } else {
             os << *value;
@@ -211,7 +234,7 @@ void numeric::print(ostream & os, unsigned upper_precedence) const
     } else {
         // case 2, imaginary:  y*I  or  -y*I
         if (realpart(*value) == 0) {
-            if ((imagpart(*value) < 0) && (precedence <= upper_precedence)) {
+            if ((precedence<=upper_precedence) && (imagpart(*value) < 0)) {
                 if (imagpart(*value) == -1) {
                     os << "(-I)";
                 } else {
@@ -230,37 +253,22 @@ void numeric::print(ostream & os, unsigned upper_precedence) const
             }
         } else {
             // case 3, complex:  x+y*I  or  x-y*I  or  -x+y*I  or  -x-y*I
-            if ((realpart(*value) < 0) && (precedence <= upper_precedence)) {
-                os << "(" << realpart(*value);
-                if (imagpart(*value) < 0) {
-                    if (imagpart(*value) == -1) {
-                        os << "-I)";
-                    } else {
-                        os << imagpart(*value) << "*I)";
-                    }
+            if (precedence <= upper_precedence) os << "(";
+            os << realpart(*value);
+            if (imagpart(*value) < 0) {
+                if (imagpart(*value) == -1) {
+                    os << "-I";
                 } else {
-                    if (imagpart(*value) == 1) {
-                        os << "+I)";
-                    } else {
-                        os << "+" << imagpart(*value) << "*I)";
-                    }
+                    os << imagpart(*value) << "*I";
                 }
             } else {
-                os << realpart(*value);
-                if (imagpart(*value) < 0) {
-                    if (imagpart(*value) == -1) {
-                        os << "-I";
-                    } else {
-                        os << imagpart(*value) << "*I";
-                    }
+                if (imagpart(*value) == 1) {
+                    os << "+I";
                 } else {
-                    if (imagpart(*value) == 1) {
-                        os << "+I";
-                    } else {
-                        os << "+" << imagpart(*value) << "*I";
-                    }
+                    os << "+" << imagpart(*value) << "*I";
                 }
             }
+            if (precedence <= upper_precedence) os << ")";
         }
     }
 }
@@ -507,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?
@@ -846,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
@@ -899,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.
@@ -911,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).
@@ -919,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).
@@ -927,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).
@@ -935,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).
@@ -943,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).
@@ -951,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.
@@ -965,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.
@@ -976,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"));        
 }
@@ -986,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).
@@ -994,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).
@@ -1002,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).
@@ -1010,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).
@@ -1018,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).
@@ -1026,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);
 }
 
@@ -1046,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
@@ -1057,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
@@ -1118,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);
@@ -1129,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).
@@ -1142,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?
@@ -1157,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?
     }
@@ -1172,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?
@@ -1240,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. */
@@ -1248,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?
@@ -1261,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();
 }
@@ -1273,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;
 }
@@ -1339,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