]> www.ginac.de Git - ginac.git/blobdiff - ginac/power.cpp
- Complete revamp of methods in class matrix. Some redundant (and poor)
[ginac.git] / ginac / power.cpp
index 32038209f51e69e53c0d93afab773c3ca03b7750..53cafc52e6b64f15e4c73eb730e2482a00f0bd61 100644 (file)
@@ -42,7 +42,7 @@ namespace GiNaC {
 
 GINAC_IMPLEMENT_REGISTERED_CLASS(power, basic)
 
-typedef vector<int> intvector;
+typedef std::vector<int> intvector;
 
 //////////
 // default constructor, destructor, copy constructor assignment operator and helpers
@@ -147,7 +147,7 @@ basic * power::duplicate() const
     return new power(*this);
 }
 
-void power::print(ostream & os, unsigned upper_precedence) const
+void power::print(std::ostream & os, unsigned upper_precedence) const
 {
     debugmsg("power print",LOGLEVEL_PRINT);
     if (exponent.is_equal(_ex1_2())) {
@@ -161,7 +161,7 @@ void power::print(ostream & os, unsigned upper_precedence) const
     }
 }
 
-void power::printraw(ostream & os) const
+void power::printraw(std::ostream & os) const
 {
     debugmsg("power printraw",LOGLEVEL_PRINT);
 
@@ -172,18 +172,19 @@ void power::printraw(ostream & os) const
     os << ",hash=" << hashvalue << ",flags=" << flags << ")";
 }
 
-void power::printtree(ostream & os, unsigned indent) const
+void power::printtree(std::ostream & os, unsigned indent) const
 {
     debugmsg("power printtree",LOGLEVEL_PRINT);
 
-    os << string(indent,' ') << "power: "
-       << "hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
-       << ", flags=" << flags << endl;
-    basis.printtree(os,indent+delta_indent);
-    exponent.printtree(os,indent+delta_indent);
+    os << std::string(indent,' ') << "power: "
+       << "hash=" << hashvalue
+       << " (0x" << std::hex << hashvalue << std::dec << ")"
+       << ", flags=" << flags << std::endl;
+    basis.printtree(os, indent+delta_indent);
+    exponent.printtree(os, indent+delta_indent);
 }
 
-static void print_sym_pow(ostream & os, unsigned type, const symbol &x, int exp)
+static void print_sym_pow(std::ostream & os, unsigned type, const symbol &x, int exp)
 {
     // Optimal output of integer powers of symbols to aid compiler CSE
     if (exp == 1) {
@@ -205,7 +206,7 @@ static void print_sym_pow(ostream & os, unsigned type, const symbol &x, int exp)
     }
 }
 
-void power::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
+void power::printcsrc(std::ostream & os, unsigned type, unsigned upper_precedence) const
 {
     debugmsg("power print csrc", LOGLEVEL_PRINT);
     
@@ -250,17 +251,20 @@ void power::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) co
 
 bool power::info(unsigned inf) const
 {
-    if (inf==info_flags::polynomial ||
-        inf==info_flags::integer_polynomial ||
-        inf==info_flags::cinteger_polynomial ||
-        inf==info_flags::rational_polynomial ||
-        inf==info_flags::crational_polynomial) {
-        return exponent.info(info_flags::nonnegint);
-    } else if (inf==info_flags::rational_function) {
-        return exponent.info(info_flags::integer);
-    } else {
-        return inherited::info(inf);
+    switch (inf) {
+        case info_flags::polynomial:
+        case info_flags::integer_polynomial:
+        case info_flags::cinteger_polynomial:
+        case info_flags::rational_polynomial:
+        case info_flags::crational_polynomial:
+            return exponent.info(info_flags::nonnegint);
+        case info_flags::rational_function:
+            return exponent.info(info_flags::integer);
+        case info_flags::algebraic:
+            return (!exponent.info(info_flags::integer) ||
+                    basis.info(inf));
     }
+    return inherited::info(inf);
 }
 
 unsigned power::nops() const
@@ -319,7 +323,7 @@ ex power::eval(int level) const
 {
     // simplifications: ^(x,0) -> 1 (0^0 handled here)
     //                  ^(x,1) -> x
-    //                  ^(0,x) -> 0 (except if the realpart of x is non-positive, in which case an exception is thrown)
+    //                  ^(0,c1) -> 0 or exception (depending on real value of c1)
     //                  ^(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!)
@@ -328,7 +332,7 @@ ex power::eval(int level) const
     //                  ^(*(x,c1),c2) -> ^(-x,c2)*c1^c2 (c1, c2 numeric(), c1<0)
     
     debugmsg("power eval",LOGLEVEL_MEMBER_FUNCTION);
-
+    
     if ((level==1) && (flags & status_flags::evaluated))
         return *this;
     else if (level == -max_recursion_level)
@@ -336,12 +340,12 @@ 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;
     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);
@@ -350,35 +354,32 @@ ex power::eval(int level) const
         exponent_is_numerical = 1;
         num_exponent = static_cast<numeric *>(eexponent.bp);
     }
-
+    
     // ^(x,0) -> 1 (0^0 also handled here)
     if (eexponent.is_zero())
         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()))
         return ebasis;
-
-    // ^(0,x) -> 0 (except if the realpart of x is non-positive)
-    if (ebasis.is_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
+    
+    // ^(0,c1) -> 0 or exception (depending on real value of c1)
+    if (ebasis.is_zero() && 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 (pole_error("power::eval(): division by zero",1));
+        else
             return _ex0();
     }
-
+    
     // ^(1,x) -> 1
     if (ebasis.is_equal(_ex1()))
         return _ex1();
-
+    
     if (basis_is_numerical && exponent_is_numerical) {
         // ^(c1,c2) -> c1^c2 (c1, c2 numeric(),
         // except if c1,c2 are rational, but c1^c2 is not)
@@ -412,7 +413,7 @@ ex power::eval(int level) const
             }
         }
     }
-
+    
     // ^(^(x,c1),c2) -> ^(x,c1*c2)
     // (c1, c2 numeric(), c2 integer or -1 < c1 <= 1,
     // case c1==1 should not happen, see below!)
@@ -434,7 +435,7 @@ ex power::eval(int level) const
         is_ex_exactly_of_type(ebasis,mul)) {
         return expand_mul(ex_to_mul(ebasis), *num_exponent);
     }
-
+    
     // ^(*(...,x;c1),c2) -> ^(*(...,x;1),c2)*c1^c2 (c1, c2 numeric(), c1>0)
     // ^(*(...,x,c1),c2) -> ^(*(...,x;-1),c2)*(-c1)^c2 (c1, c2 numeric(), c1<0)
     if (exponent_is_numerical && is_ex_exactly_of_type(ebasis,mul)) {
@@ -524,7 +525,11 @@ 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));
+        epvector newseq;
+        newseq.reserve(2);
+        newseq.push_back(expair(basis, exponent - _ex1()));
+        newseq.push_back(expair(basis.diff(s), _ex1()));
+        return mul(newseq, exponent);
     } else {
         // D(b^e) = b^e * (D(e)*ln(b) + e*D(b)/b)
         return mul(power(basis, exponent),
@@ -635,7 +640,10 @@ ex power::expand_add(const add & a, int n) const
             GINAC_ASSERT(!is_ex_exactly_of_type(b,add));
             GINAC_ASSERT(!is_ex_exactly_of_type(b,power)||
                          !is_ex_exactly_of_type(ex_to_power(b).exponent,numeric)||
-                         !ex_to_numeric(ex_to_power(b).exponent).is_pos_integer());
+                         !ex_to_numeric(ex_to_power(b).exponent).is_pos_integer()||
+                         !is_ex_exactly_of_type(ex_to_power(b).basis,add)||
+                         !is_ex_exactly_of_type(ex_to_power(b).basis,mul)||
+                         !is_ex_exactly_of_type(ex_to_power(b).basis,power));
             if (is_ex_exactly_of_type(b,mul)) {
                 term.push_back(expand_mul(ex_to_mul(b),numeric(k[l])));
             } else {
@@ -647,7 +655,10 @@ ex power::expand_add(const add & a, int n) const
         GINAC_ASSERT(!is_ex_exactly_of_type(b,add));
         GINAC_ASSERT(!is_ex_exactly_of_type(b,power)||
                      !is_ex_exactly_of_type(ex_to_power(b).exponent,numeric)||
-                     !ex_to_numeric(ex_to_power(b).exponent).is_pos_integer());
+                     !ex_to_numeric(ex_to_power(b).exponent).is_pos_integer()||
+                     !is_ex_exactly_of_type(ex_to_power(b).basis,add)||
+                     !is_ex_exactly_of_type(ex_to_power(b).basis,mul)||
+                     !is_ex_exactly_of_type(ex_to_power(b).basis,power));
         if (is_ex_exactly_of_type(b,mul)) {
             term.push_back(expand_mul(ex_to_mul(b),numeric(n-k_cum[m-2])));
         } else {