- Changed behaviour of .evalf() to not evalf any exponents.

[ginac.git] / ginac / power.cpp

diff --git a/ginac/power.cpp b/ginac/power.cpp
index 89ba78bcc9a79c9e6d184ced3eb0372e8152df8e..8cb6223860cdd2cfef628841cb3cb1fa86766d81 100644 (file)
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -29,6 +29,7 @@
 #include "add.h"
 #include "mul.h"
 #include "numeric.h"
+#include "inifcns.h"
 #include "relational.h"
 #include "symbol.h"
 #include "archive.h"
@@ -337,23 +338,26 @@ ex power::eval(int level) const
     const ex & ebasis    = level==1 ? basis    : basis.eval(level-1);
     const ex & eexponent = level==1 ? exponent : exponent.eval(level-1);
 
-    bool basis_is_numerical=0;
-    bool exponent_is_numerical=0;
+    bool basis_is_numerical = 0;
+    bool exponent_is_numerical = 0;
     numeric * num_basis;
     numeric * num_exponent;
 
     if (is_exactly_of_type(*ebasis.bp,numeric)) {
-        basis_is_numerical=1;
-        num_basis=static_cast<numeric *>(ebasis.bp);
+        basis_is_numerical = 1;
+        num_basis = static_cast<numeric *>(ebasis.bp);
     }
     if (is_exactly_of_type(*eexponent.bp,numeric)) {
-        exponent_is_numerical=1;
-        num_exponent=static_cast<numeric *>(eexponent.bp);
+        exponent_is_numerical = 1;
+        num_exponent = static_cast<numeric *>(eexponent.bp);
     }
 
     // ^(x,0) -> 1 (0^0 also handled here)
     if (eexponent.is_zero())
-        return _ex1();
+        if (ebasis.is_zero())
+            throw (std::domain_error("power::eval(): pow(0,0) is undefined"));
+        else
+            return _ex1();
 
     // ^(x,1) -> x
     if (eexponent.is_equal(_ex1()))
@@ -387,10 +391,10 @@ ex power::eval(int level) const
         if (basis_is_crational && exponent_is_crational
             && num_exponent->is_real()
             && !num_exponent->is_integer()) {
-            numeric r, q, n, m;
-            n = num_exponent->numer();
-            m = num_exponent->denom();
-            q = iquo(n, m, r);
+            numeric n = num_exponent->numer();
+            numeric m = num_exponent->denom();
+            numeric r;
+            numeric q = iquo(n, m, r);
             if (r.is_negative()) {
                 r = r.add(m);
                 q = q.sub(_num1());
@@ -398,29 +402,22 @@ ex power::eval(int level) const
             if (q.is_zero())  // the exponent was in the allowed range 0<(n/m)<1
                 return this->hold();
             else {
-                epvector res(2);
+                epvector res;
                 res.push_back(expair(ebasis,r.div(m)));
-                res.push_back(expair(ex(num_basis->power(q)),_ex1()));
-                return (new mul(res))->setflag(status_flags::dynallocated | status_flags::evaluated);
-                /*return mul(num_basis->power(q),
-                           power(ex(*num_basis),ex(r.div(m)))).hold();
-                */
-                /* return (new mul(num_basis->power(q),
-                   power(*num_basis,r.div(m)).hold()))->setflag(status_flags::dynallocated | status_flags::evaluated);
-                */
+                return (new mul(res,ex(num_basis->power(q))))->setflag(status_flags::dynallocated | status_flags::evaluated);
             }
         }
     }
 
     // ^(^(x,c1),c2) -> ^(x,c1*c2)
     // (c1, c2 numeric(), c2 integer or -1 < c1 <= 1,
-    // case c1=1 should not happen, see below!)
+    // case c1==1 should not happen, see below!)
     if (exponent_is_numerical && is_ex_exactly_of_type(ebasis,power)) {
-        const power & sub_power=ex_to_power(ebasis);
-        const ex & sub_basis=sub_power.basis;
-        const ex & sub_exponent=sub_power.exponent;
+        const power & sub_power = ex_to_power(ebasis);
+        const ex & sub_basis = sub_power.basis;
+        const ex & sub_exponent = sub_power.exponent;
         if (is_ex_exactly_of_type(sub_exponent,numeric)) {
-            const numeric & num_sub_exponent=ex_to_numeric(sub_exponent);
+            const numeric & num_sub_exponent = ex_to_numeric(sub_exponent);
             GINAC_ASSERT(num_sub_exponent!=numeric(1));
             if (num_exponent->is_integer() || abs(num_sub_exponent)<1) {
                 return power(sub_basis,num_sub_exponent.mul(*num_exponent));
@@ -482,13 +479,16 @@ ex power::evalf(int level) const
     ex eexponent;
     
     if (level==1) {
-        ebasis=basis;
-        eexponent=exponent;
+        ebasis = basis;
+        eexponent = exponent;
     } else if (level == -max_recursion_level) {
         throw(std::runtime_error("max recursion level reached"));
     } else {
-        ebasis=basis.evalf(level-1);
-        eexponent=exponent.evalf(level-1);
+        ebasis = basis.evalf(level-1);
+        if (!is_ex_exactly_of_type(eexponent,numeric))
+            eexponent = exponent.evalf(level-1);
+        else
+            eexponent = exponent;
     }
 
     return power(ebasis,eexponent);
@@ -514,6 +514,21 @@ ex power::simplify_ncmul(const exvector & v) const
 
 // protected
 
+/** Implementation of ex::diff() for a power.
+ *  @see ex::diff */
+ex power::derivative(const symbol & s) const
+{
+    if (exponent.info(info_flags::real)) {
+        // D(b^r) = r * b^(r-1) * D(b) (faster than the formula below)
+        return mul(mul(exponent, power(basis, exponent - _ex1())), basis.diff(s));
+    } else {
+        // D(b^e) = b^e * (D(e)*ln(b) + e*D(b)/b)
+        return mul(power(basis, exponent),
+                   add(mul(exponent.diff(s), log(basis)),
+                       mul(mul(exponent, basis.diff(s)), power(basis, -1))));
+    }
+}
+
 int power::compare_same_type(const basic & other) const
 {
     GINAC_ASSERT(is_exactly_of_type(other, power));