]> www.ginac.de Git - ginac.git/blobdiff - ginac/power.cpp
expand() always expands the exponent and transforms x^(a+b) -> x^a*x^b
[ginac.git] / ginac / power.cpp
index 2b6049ba4f3eb90dbf787ddcb6818c6d889de7c1..efddc39338b3f4d88276e6ba7a78d979829bab3b 100644 (file)
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
 
 /*
- *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2001 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
@@ -58,7 +58,7 @@ power::power() : basic(TINFO_power)
 power::~power()
 {
        debugmsg("power destructor",LOGLEVEL_DESTRUCT);
-       destroy(0);
+       destroy(false);
 }
 
 power::power(const power & other)
@@ -71,7 +71,7 @@ const power & power::operator=(const power & other)
 {
        debugmsg("power operator=",LOGLEVEL_ASSIGNMENT);
        if (this != &other) {
-               destroy(1);
+               destroy(true);
                copy(other);
        }
        return *this;
@@ -211,9 +211,8 @@ void power::printcsrc(std::ostream & os, unsigned type, unsigned upper_precedenc
        debugmsg("power print csrc", LOGLEVEL_PRINT);
        
        // Integer powers of symbols are printed in a special, optimized way
-       if (exponent.info(info_flags::integer) &&
-               (is_ex_exactly_of_type(basis, symbol) ||
-                is_ex_exactly_of_type(basis, constant))) {
+       if (exponent.info(info_flags::integer)
+        && (is_ex_exactly_of_type(basis, symbol) || is_ex_exactly_of_type(basis, constant))) {
                int exp = ex_to_numeric(exponent).to_int();
                if (exp > 0)
                        os << "(";
@@ -283,9 +282,12 @@ ex & power::let_op(int i)
 int power::degree(const symbol & s) const
 {
        if (is_exactly_of_type(*exponent.bp,numeric)) {
-               if ((*basis.bp).compare(s)==0)
-                       return ex_to_numeric(exponent).to_int();
-               else
+               if ((*basis.bp).compare(s)==0) {
+                       if (ex_to_numeric(exponent).is_integer())
+                               return ex_to_numeric(exponent).to_int();
+                       else
+                               return 0;
+               } else
                        return basis.degree(s) * ex_to_numeric(exponent).to_int();
        }
        return 0;
@@ -294,9 +296,12 @@ int power::degree(const symbol & s) const
 int power::ldegree(const symbol & s) const 
 {
        if (is_exactly_of_type(*exponent.bp,numeric)) {
-               if ((*basis.bp).compare(s)==0)
-                       return ex_to_numeric(exponent).to_int();
-               else
+               if ((*basis.bp).compare(s)==0) {
+                       if (ex_to_numeric(exponent).is_integer())
+                               return ex_to_numeric(exponent).to_int();
+                       else
+                               return 0;
+               } else
                        return basis.ldegree(s) * ex_to_numeric(exponent).to_int();
        }
        return 0;
@@ -306,17 +311,27 @@ ex power::coeff(const symbol & s, int n) const
 {
        if ((*basis.bp).compare(s)!=0) {
                // basis not equal to s
-               if (n==0) {
+               if (n == 0)
                        return *this;
-               } else {
+               else
                        return _ex0();
+       } else {
+               // basis equal to s
+               if (is_exactly_of_type(*exponent.bp, numeric) && ex_to_numeric(exponent).is_integer()) {
+                       // integer exponent
+                       int int_exp = ex_to_numeric(exponent).to_int();
+                       if (n == int_exp)
+                               return _ex1();
+                       else
+                               return _ex0();
+               } else {
+                       // non-integer exponents are treated as zero
+                       if (n == 0)
+                               return *this;
+                       else
+                               return _ex0();
                }
-       } else if (is_exactly_of_type(*exponent.bp,numeric)&&
-                          (static_cast<const numeric &>(*exponent.bp).compare(numeric(n))==0)) {
-               return _ex1();
        }
-
-       return _ex0();
 }
 
 ex power::eval(int level) const
@@ -450,8 +465,7 @@ ex power::eval(int level) const
                                        mulp->clearflag(status_flags::evaluated);
                                        mulp->clearflag(status_flags::hash_calculated);
                                        return (new mul(power(*mulp,exponent),
-                                                                       power(num_coeff,*num_exponent)))->
-                                               setflag(status_flags::dynallocated);
+                                                       power(num_coeff,*num_exponent)))->setflag(status_flags::dynallocated);
                                } else {
                                        GINAC_ASSERT(num_coeff.compare(_num0())<0);
                                        if (num_coeff.compare(_num_1())!=0) {
@@ -460,8 +474,7 @@ ex power::eval(int level) const
                                                mulp->clearflag(status_flags::evaluated);
                                                mulp->clearflag(status_flags::hash_calculated);
                                                return (new mul(power(*mulp,exponent),
-                                                                               power(abs(num_coeff),*num_exponent)))->
-                                                       setflag(status_flags::dynallocated);
+                                                               power(abs(num_coeff),*num_exponent)))->setflag(status_flags::dynallocated);
                                        }
                                }
                        }
@@ -532,9 +545,9 @@ ex power::derivative(const symbol & s) const
                return mul(newseq, exponent);
        } 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))));
+               return mul(*this,
+                          add(mul(exponent.diff(s), log(basis)),
+                          mul(mul(exponent, basis.diff(s)), power(basis, _ex_1()))));
        }
 }
 
@@ -567,37 +580,64 @@ ex power::expand(unsigned options) const
                return *this;
        
        ex expanded_basis = basis.expand(options);
-       
-       if (!is_ex_exactly_of_type(exponent,numeric) ||
-               !ex_to_numeric(exponent).is_integer()) {
-               if (are_ex_trivially_equal(basis,expanded_basis)) {
+       ex expanded_exponent = exponent.expand(options);
+
+       // x^(a+b) -> x^a * x^b
+       if (is_ex_exactly_of_type(expanded_exponent, add)) {
+               const add &a = ex_to_add(expanded_exponent);
+               exvector distrseq;
+               distrseq.reserve(a.seq.size() + 1);
+               epvector::const_iterator last = a.seq.end();
+               epvector::const_iterator cit = a.seq.begin();
+               while (cit!=last) {
+                       distrseq.push_back(power(expanded_basis, a.recombine_pair_to_ex(*cit)));
+                       cit++;
+               }
+
+               // Make sure that e.g. (x+y)^(2+a) expands the (x+y)^2 factor
+               if (ex_to_numeric(a.overall_coeff).is_integer()) {
+                       const numeric &num_exponent = ex_to_numeric(a.overall_coeff);
+                       int int_exponent = num_exponent.to_int();
+                       if (int_exponent > 0 && is_ex_exactly_of_type(expanded_basis, add))
+                               distrseq.push_back(expand_add(ex_to_add(expanded_basis), int_exponent));
+                       else
+                               distrseq.push_back(power(expanded_basis, a.overall_coeff));
+               } else
+                       distrseq.push_back(power(expanded_basis, a.overall_coeff));
+
+               // Make sure that e.g. (x+y)^(1+a) -> x*(x+y)^a + y*(x+y)^a
+               ex r = (new mul(distrseq))->setflag(status_flags::dynallocated);
+               return r.expand();
+       }
+
+       if (!is_ex_exactly_of_type(expanded_exponent, numeric) ||
+               !ex_to_numeric(expanded_exponent).is_integer()) {
+               if (are_ex_trivially_equal(basis,expanded_basis) && are_ex_trivially_equal(exponent,expanded_exponent)) {
                        return this->hold();
                } else {
-                       return (new power(expanded_basis,exponent))->
-                               setflag(status_flags::dynallocated |
-                                               status_flags::expanded);
+                       return (new power(expanded_basis,expanded_exponent))->setflag(status_flags::dynallocated | status_flags::expanded);
                }
        }
        
        // integer numeric exponent
-       const numeric & num_exponent = ex_to_numeric(exponent);
+       const numeric & num_exponent = ex_to_numeric(expanded_exponent);
        int int_exponent = num_exponent.to_int();
        
+       // (x+y)^n, n>0
        if (int_exponent > 0 && is_ex_exactly_of_type(expanded_basis,add)) {
                return expand_add(ex_to_add(expanded_basis), int_exponent);
        }
        
+       // (x*y)^n -> x^n * y^n
        if (is_ex_exactly_of_type(expanded_basis,mul)) {
                return expand_mul(ex_to_mul(expanded_basis), num_exponent);
        }
        
        // cannot expand further
-       if (are_ex_trivially_equal(basis,expanded_basis)) {
+       if (are_ex_trivially_equal(basis,expanded_basis) && are_ex_trivially_equal(exponent,expanded_exponent)) {
                return this->hold();
        } else {
-               return (new power(expanded_basis,exponent))->
-                          setflag(status_flags::dynallocated |
-                                          status_flags::expanded);
+               return (new power(expanded_basis,expanded_exponent))->setflag(status_flags::dynallocated | status_flags::expanded);
        }
 }
 
@@ -638,12 +678,12 @@ ex power::expand_add(const add & a, int n) const
                for (l=0; l<m-1; l++) {
                        const ex & b = a.op(l);
                        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()||
-                                                !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));
+                       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() ||
+                                    !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 {
@@ -653,12 +693,12 @@ ex power::expand_add(const add & a, int n) const
                
                const ex & b = a.op(l);
                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()||
-                                        !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));
+               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() ||
+                            !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 {
@@ -730,27 +770,28 @@ ex power::expand_add_2(const add & a) const
                const ex & c=(*cit0).coeff;
                
                GINAC_ASSERT(!is_ex_exactly_of_type(r,add));
-               GINAC_ASSERT(!is_ex_exactly_of_type(r,power)||
-                          !is_ex_exactly_of_type(ex_to_power(r).exponent,numeric)||
-                          !ex_to_numeric(ex_to_power(r).exponent).is_pos_integer()||
-                          !is_ex_exactly_of_type(ex_to_power(r).basis,add)||
-                          !is_ex_exactly_of_type(ex_to_power(r).basis,mul)||
-                          !is_ex_exactly_of_type(ex_to_power(r).basis,power));
+               GINAC_ASSERT(!is_ex_exactly_of_type(r,power) ||
+                            !is_ex_exactly_of_type(ex_to_power(r).exponent,numeric) ||
+                            !ex_to_numeric(ex_to_power(r).exponent).is_pos_integer() ||
+                            !is_ex_exactly_of_type(ex_to_power(r).basis,add) ||
+                            !is_ex_exactly_of_type(ex_to_power(r).basis,mul) ||
+                            !is_ex_exactly_of_type(ex_to_power(r).basis,power));
 
                if (are_ex_trivially_equal(c,_ex1())) {
                        if (is_ex_exactly_of_type(r,mul)) {
-                               sum.push_back(expair(expand_mul(ex_to_mul(r),_num2()),_ex1()));
+                               sum.push_back(expair(expand_mul(ex_to_mul(r),_num2()),
+                                                    _ex1()));
                        } else {
                                sum.push_back(expair((new power(r,_ex2()))->setflag(status_flags::dynallocated),
-                                                                        _ex1()));
+                                                    _ex1()));
                        }
                } else {
                        if (is_ex_exactly_of_type(r,mul)) {
                                sum.push_back(expair(expand_mul(ex_to_mul(r),_num2()),
-                                                                        ex_to_numeric(c).power_dyn(_num2())));
+                                                    ex_to_numeric(c).power_dyn(_num2())));
                        } else {
                                sum.push_back(expair((new power(r,_ex2()))->setflag(status_flags::dynallocated),
-                                                                        ex_to_numeric(c).power_dyn(_num2())));
+                                                    ex_to_numeric(c).power_dyn(_num2())));
                        }
                }
                        
@@ -758,7 +799,7 @@ ex power::expand_add_2(const add & a) const
                        const ex & r1=(*cit1).rest;
                        const ex & c1=(*cit1).coeff;
                        sum.push_back(a.combine_ex_with_coeff_to_pair((new mul(r,r1))->setflag(status_flags::dynallocated),
-                                                                                                                 _num2().mul(ex_to_numeric(c)).mul_dyn(ex_to_numeric(c1))));
+                                                                     _num2().mul(ex_to_numeric(c)).mul_dyn(ex_to_numeric(c1))));
                }
        }
 
@@ -774,8 +815,7 @@ ex power::expand_add_2(const add & a) const
                
        GINAC_ASSERT(sum.size()==(a_nops*(a_nops+1))/2);
        
-       return (new add(sum))->setflag(status_flags::dynallocated |
-                                                                  status_flags::expanded );
+       return (new add(sum))->setflag(status_flags::dynallocated | status_flags::expanded);
 }
 
 /** Expand factors of m in m^n where m is a mul and n is and integer
@@ -795,18 +835,16 @@ ex power::expand_mul(const mul & m, const numeric & n) const
                } else {
                        // it is safe not to call mul::combine_pair_with_coeff_to_pair()
                        // since n is an integer
-                       distrseq.push_back(expair((*cit).rest,
-                                                                         ex_to_numeric((*cit).coeff).mul(n)));
+                       distrseq.push_back(expair((*cit).rest, ex_to_numeric((*cit).coeff).mul(n)));
                }
                ++cit;
        }
-       return (new mul(distrseq,ex_to_numeric(m.overall_coeff).power_dyn(n)))
-               ->setflag(status_flags::dynallocated);
+       return (new mul(distrseq,ex_to_numeric(m.overall_coeff).power_dyn(n)))->setflag(status_flags::dynallocated);
 }
 
 /*
 ex power::expand_commutative_3(const ex & basis, const numeric & exponent,
-                                                        unsigned options) const
+                               unsigned options) const
 {
        // obsolete
 
@@ -822,12 +860,10 @@ ex power::expand_commutative_3(const ex & basis, const numeric & exponent,
        
        int n=exponent.to_int();
        for (int k=0; k<=n; k++) {
-               distrseq.push_back(binomial(n,k)*power(first_operands,numeric(k))*
-                                                  power(last_operand,numeric(n-k)));
+               distrseq.push_back(binomial(n,k) * power(first_operands,numeric(k))
+                                                * power(last_operand,numeric(n-k)));
        }
-       return ex((new add(distrseq))->setflag(status_flags::expanded |
-                                                                                  status_flags::dynallocated )).
-                  expand(options);
+       return ex((new add(distrseq))->setflag(status_flags::expanded | status_flags::dynallocated)).expand(options);
 }
 */
 
@@ -835,11 +871,11 @@ ex power::expand_commutative_3(const ex & basis, const numeric & exponent,
 ex power::expand_noncommutative(const ex & basis, const numeric & exponent,
                                                                unsigned options) const
 {
-       ex rest_power=ex(power(basis,exponent.add(_num_1()))).
-                                 expand(options | expand_options::internal_do_not_expand_power_operands);
+       ex rest_power = ex(power(basis,exponent.add(_num_1()))).
+                       expand(options | expand_options::internal_do_not_expand_power_operands);
 
        return ex(mul(rest_power,basis),0).
-                  expand(options | expand_options::internal_do_not_expand_mul_operands);
+              expand(options | expand_options::internal_do_not_expand_mul_operands);
 }
 */
 
@@ -856,7 +892,7 @@ unsigned power::precedence = 60;
 //////////
 
 const power some_power;
-const type_info & typeid_power=typeid(some_power);
+const std::type_info & typeid_power=typeid(some_power);
 
 // helper function