- 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 94012384f6377aea088ae4253bbb61ad99260aa5..b11761f58dc4e2dc1c48d4dda011995721925af3 100644 (file)
--- 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 30091760507018a5943819d9a1e35253f9c3c32c..95ba4419073d4b75b2d30e2b3bc3bc4176231883 100644 (file)
--- 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 2baa08d98bfc573abbe9b089eaab444a2dc78473..32038209f51e69e53c0d93afab773c3ca03b7750 100644 (file)
--- 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();
     }