* Avoid bad choice of main variable in heuristic gcd [Alexei Sheplyakov].

[ginac.git] / ginac / power.cpp

diff --git a/ginac/power.cpp b/ginac/power.cpp
index 339b48c1a3cfd067ef692ea968082e0ff3102b7b..75fd2ff51a9d12bd09d867d30fdf50b860ec6240 100644 (file)
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
 
 /*
- *  GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2006 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
@@ -40,6 +40,7 @@
 #include "lst.h"
 #include "archive.h"
 #include "utils.h"
+#include "relational.h"
 
 namespace GiNaC {
 
@@ -56,7 +57,7 @@ typedef std::vector<int> intvector;
 // default constructor
 //////////
 
-power::power() : inherited(TINFO_power) { }
+power::power() : inherited(&power::tinfo_static) { }
 
 //////////
 // other constructors
@@ -258,6 +259,15 @@ ex power::map(map_function & f) const
                return *this;
 }
 
+bool power::is_polynomial(const ex & var) const
+{
+       if (exponent.has(var))
+               return false;
+       if (!exponent.info(info_flags::nonnegint))
+               return false;
+       return basis.is_polynomial(var);
+}
+
 int power::degree(const ex & s) const
 {
        if (is_equal(ex_to<basic>(s)))
@@ -381,6 +391,10 @@ ex power::eval(int level) const
        if (ebasis.is_equal(_ex1))
                return _ex1;
 
+       // power of a function calculated by separate rules defined for this function
+       if (is_exactly_a<function>(ebasis))
+               return ex_to<function>(ebasis).power(eexponent);
+
        if (exponent_is_numerical) {
 
                // ^(c1,c2) -> c1^c2  (c1, c2 numeric(),
@@ -537,6 +551,30 @@ ex power::evalm() const
        return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated);
 }
 
+bool power::has(const ex & other, unsigned options) const
+{
+       if (!(options & has_options::algebraic))
+               return basic::has(other, options);
+       if (!is_a<power>(other))
+               return basic::has(other, options);
+       if (!exponent.info(info_flags::integer)
+                       || !other.op(1).info(info_flags::integer))
+               return basic::has(other, options);
+       if (exponent.info(info_flags::posint)
+                       && other.op(1).info(info_flags::posint)
+                       && ex_to<numeric>(exponent).to_int()
+                                       > ex_to<numeric>(other.op(1)).to_int()
+                       && basis.match(other.op(0)))
+               return true;
+       if (exponent.info(info_flags::negint)
+                       && other.op(1).info(info_flags::negint)
+                       && ex_to<numeric>(exponent).to_int()
+                                       < ex_to<numeric>(other.op(1)).to_int()
+                       && basis.match(other.op(0)))
+               return true;
+       return basic::has(other, options);
+}
+
 // from mul.cpp
 extern bool tryfactsubs(const ex &, const ex &, int &, lst &);
 
@@ -577,13 +615,66 @@ ex power::conjugate() const
        return (new power(newbasis, newexponent))->setflag(status_flags::dynallocated);
 }
 
+ex power::real_part() const
+{
+       if (exponent.info(info_flags::integer)) {
+               ex basis_real = basis.real_part();
+               if (basis_real == basis)
+                       return *this;
+               realsymbol a("a"),b("b");
+               ex result;
+               if (exponent.info(info_flags::posint))
+                       result = power(a+I*b,exponent);
+               else
+                       result = power(a/(a*a+b*b)-I*b/(a*a+b*b),-exponent);
+               result = result.expand();
+               result = result.real_part();
+               result = result.subs(lst( a==basis_real, b==basis.imag_part() ));
+               return result;
+       }
+       
+       ex a = basis.real_part();
+       ex b = basis.imag_part();
+       ex c = exponent.real_part();
+       ex d = exponent.imag_part();
+       return power(abs(basis),c)*exp(-d*atan2(b,a))*cos(c*atan2(b,a)+d*log(abs(basis)));
+}
+
+ex power::imag_part() const
+{
+       if (exponent.info(info_flags::integer)) {
+               ex basis_real = basis.real_part();
+               if (basis_real == basis)
+                       return 0;
+               realsymbol a("a"),b("b");
+               ex result;
+               if (exponent.info(info_flags::posint))
+                       result = power(a+I*b,exponent);
+               else
+                       result = power(a/(a*a+b*b)-I*b/(a*a+b*b),-exponent);
+               result = result.expand();
+               result = result.imag_part();
+               result = result.subs(lst( a==basis_real, b==basis.imag_part() ));
+               return result;
+       }
+       
+       ex a=basis.real_part();
+       ex b=basis.imag_part();
+       ex c=exponent.real_part();
+       ex d=exponent.imag_part();
+       return
+               power(abs(basis),c)*exp(-d*atan2(b,a))*sin(c*atan2(b,a)+d*log(abs(basis)));
+}
+
+// protected
+
 // protected
 
 /** Implementation of ex::diff() for a power.
  *  @see ex::diff */
 ex power::derivative(const symbol & s) const
 {
-       if (exponent.info(info_flags::real)) {
+       if (is_a<numeric>(exponent)) {
                // D(b^r) = r * b^(r-1) * D(b) (faster than the formula below)
                epvector newseq;
                newseq.reserve(2);
@@ -615,7 +706,7 @@ unsigned power::return_type() const
        return basis.return_type();
 }
 
-unsigned power::return_type_tinfo() const
+tinfo_t power::return_type_tinfo() const
 {
        return basis.return_type_tinfo();
 }
@@ -860,8 +951,11 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr
        // Leave it to multiplication since dummy indices have to be renamed
        if (get_all_dummy_indices(m).size() > 0 && n.is_positive()) {
                ex result = m;
+               exvector va = get_all_dummy_indices(m);
+               sort(va.begin(), va.end(), ex_is_less());
+
                for (int i=1; i < n.to_int(); i++)
-                       result *= rename_dummy_indices_uniquely(m,m);
+                       result *= rename_dummy_indices_uniquely(va, m);
                return result;
        }