From: Richard Kreckel <Richard.Kreckel@uni-mainz.de>
Date: Wed, 15 Mar 2000 22:28:41 +0000 (+0000)
Subject: - According to CLTL 0^I is undefined, 0^(I+epsilon) is 0 and 0^(I-epsilon)
X-Git-Tag: release_0-5-4~4
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=f2ec274c38bc881d4c52a0e5eb215fd78c730b4d;hp=f6f99c4d47762da9e3d73a2f5ec6f062e82505b8

- According to CLTL 0^I is undefined, 0^(I+epsilon) is 0 and 0^(I-epsilon)
  is an overflow.  Now, this should be honored.
---

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index 94012384..b11761f5 100644
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -685,13 +685,15 @@ numeric numeric::div(const numeric & other) const
 
 numeric numeric::power(const numeric & other) const
 {
-    static const numeric * _num1p=&_num1();
+    static const numeric * _num1p = &_num1();
     if (&other==_num1p)
         return *this;
     if (::zerop(*value)) {
         if (::zerop(*other.value))
             throw (std::domain_error("numeric::eval(): pow(0,0) is undefined"));
-        else if (other.is_real() && !::plusp(::realpart(*other.value)))
+        else if (::zerop(::realpart(*other.value)))
+            throw (std::domain_error("numeric::eval(): pow(0,I) is undefined"));
+        else if (::minusp(::realpart(*other.value)))
             throw (std::overflow_error("numeric::eval(): division by zero"));
         else
             return _num0();
@@ -745,7 +747,9 @@ const numeric & numeric::power_dyn(const numeric & other) const
     if (::zerop(*value)) {
         if (::zerop(*other.value))
             throw (std::domain_error("numeric::eval(): pow(0,0) is undefined"));
-        else if (other.is_real() && !::plusp(::realpart(*other.value)))
+        else if (::zerop(::realpart(*other.value)))
+            throw (std::domain_error("numeric::eval(): pow(0,I) is undefined"));
+        else if (::minusp(::realpart(*other.value)))
             throw (std::overflow_error("numeric::eval(): division by zero"));
         else
             return _num0();
@@ -1016,13 +1020,13 @@ double numeric::to_double(void) const
 }
 
 /** Real part of a number. */
-numeric numeric::real(void) const
+const numeric numeric::real(void) const
 {
     return numeric(::realpart(*value));  // -> CLN
 }
 
 /** Imaginary part of a number. */
-numeric numeric::imag(void) const
+const numeric numeric::imag(void) const
 {
     return numeric(::imagpart(*value));  // -> CLN
 }
@@ -1044,7 +1048,7 @@ inline cl_heap_ratio* TheRatio (const cl_N& obj)
  *  numerator of complex if real and imaginary part are both rational numbers
  *  (i.e numer(4/3+5/6*I) == 8+5*I), the number carrying the sign in all other
  *  cases. */
-numeric numeric::numer(void) const
+const numeric numeric::numer(void) const
 {
     if (this->is_integer()) {
         return numeric(*this);
@@ -1095,7 +1099,7 @@ numeric numeric::numer(void) const
 /** Denominator.  Computes the denominator of rational numbers, common integer
  *  denominator of complex if real and imaginary part are both rational numbers
  *  (i.e denom(4/3+5/6*I) == 6), one in all other cases. */
-numeric numeric::denom(void) const
+const numeric numeric::denom(void) const
 {
     if (this->is_integer()) {
         return _num1();
diff --git a/ginac/numeric.h b/ginac/numeric.h
index 30091760..95ba4419 100644
--- a/ginac/numeric.h
+++ b/ginac/numeric.h
@@ -188,10 +188,10 @@ public:
     int to_int(void) const;
     long to_long(void) const;
     double to_double(void) const;
-    numeric real(void) const;
-    numeric imag(void) const;
-    numeric numer(void) const;
-    numeric denom(void) const;
+    const numeric real(void) const;
+    const numeric imag(void) const;
+    const numeric numer(void) const;
+    const numeric denom(void) const;
     int int_length(void) const;
 
 // member variables
diff --git a/ginac/power.cpp b/ginac/power.cpp
index 2baa08d9..32038209 100644
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -319,7 +319,7 @@ ex power::eval(int level) const
 {
     // simplifications: ^(x,0) -> 1 (0^0 handled here)
     //                  ^(x,1) -> x
-    //                  ^(0,x) -> 0 (except if x is real and negative, in which case an exception is thrown)
+    //                  ^(0,x) -> 0 (except if the realpart of x is non-positive, in which case an exception is thrown)
     //                  ^(1,x) -> 1
     //                  ^(c1,c2) -> *(c1^n,c1^(c2-n)) (c1, c2 numeric(), 0<(c2-n)<1 except if c1,c2 are rational, but c1^c2 is not)
     //                  ^(^(x,c1),c2) -> ^(x,c1*c2) (c1, c2 numeric(), c2 integer or -1 < c1 <= 1, case c1=1 should not happen, see below!)
@@ -362,10 +362,15 @@ ex power::eval(int level) const
     if (eexponent.is_equal(_ex1()))
         return ebasis;
 
-    // ^(0,x) -> 0 (except if x is real and negative)
+    // ^(0,x) -> 0 (except if the realpart of x is non-positive)
     if (ebasis.is_zero()) {
-        if (exponent_is_numerical && num_exponent->is_negative()) {
-            throw(std::overflow_error("power::eval(): division by zero"));
+        if (exponent_is_numerical) {
+            if ((num_exponent->real()).is_zero())
+                throw (std::domain_error("power::eval(): pow(0,I) is undefined"));
+            else if ((num_exponent->real()).is_negative())
+                throw (std::overflow_error("power::eval(): division by zero"));
+            else
+                return _ex0();
         } else
             return _ex0();
     }