]> www.ginac.de Git - ginac.git/blobdiff - ginac/mul.cpp
Changed naming convention for the library. Now, it is only libginac.so
[ginac.git] / ginac / mul.cpp
index c4cb2e3b90f72026c02a4b19cfd25c07efe66946..27add930cf56863ae7d635d3031e4dc2b1d7ea9f 100644 (file)
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's products of expressions. */
 
 /*
- *  GiNaC Copyright (C) 1999-2007 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2011 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
  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
-#include <iostream>
-#include <vector>
-#include <stdexcept>
-#include <limits>
-
 #include "mul.h"
 #include "add.h"
 #include "power.h"
 #include "symbol.h"
 #include "compiler.h"
 
+#include <iostream>
+#include <limits>
+#include <stdexcept>
+#include <vector>
+
 namespace GiNaC {
 
 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mul, expairseq,
@@ -53,7 +53,6 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mul, expairseq,
 
 mul::mul()
 {
-       tinfo_key = &mul::tinfo_static;
 }
 
 //////////
@@ -64,7 +63,6 @@ mul::mul()
 
 mul::mul(const ex & lh, const ex & rh)
 {
-       tinfo_key = &mul::tinfo_static;
        overall_coeff = _ex1;
        construct_from_2_ex(lh,rh);
        GINAC_ASSERT(is_canonical());
@@ -72,7 +70,6 @@ mul::mul(const ex & lh, const ex & rh)
 
 mul::mul(const exvector & v)
 {
-       tinfo_key = &mul::tinfo_static;
        overall_coeff = _ex1;
        construct_from_exvector(v);
        GINAC_ASSERT(is_canonical());
@@ -80,7 +77,6 @@ mul::mul(const exvector & v)
 
 mul::mul(const epvector & v)
 {
-       tinfo_key = &mul::tinfo_static;
        overall_coeff = _ex1;
        construct_from_epvector(v);
        GINAC_ASSERT(is_canonical());
@@ -88,7 +84,6 @@ mul::mul(const epvector & v)
 
 mul::mul(const epvector & v, const ex & oc, bool do_index_renaming)
 {
-       tinfo_key = &mul::tinfo_static;
        overall_coeff = oc;
        construct_from_epvector(v, do_index_renaming);
        GINAC_ASSERT(is_canonical());
@@ -96,7 +91,6 @@ mul::mul(const epvector & v, const ex & oc, bool do_index_renaming)
 
 mul::mul(std::auto_ptr<epvector> vp, const ex & oc, bool do_index_renaming)
 {
-       tinfo_key = &mul::tinfo_static;
        GINAC_ASSERT(vp.get()!=0);
        overall_coeff = oc;
        construct_from_epvector(*vp, do_index_renaming);
@@ -105,7 +99,6 @@ mul::mul(std::auto_ptr<epvector> vp, const ex & oc, bool do_index_renaming)
 
 mul::mul(const ex & lh, const ex & mh, const ex & rh)
 {
-       tinfo_key = &mul::tinfo_static;
        exvector factors;
        factors.reserve(3);
        factors.push_back(lh);
@@ -120,8 +113,6 @@ mul::mul(const ex & lh, const ex & mh, const ex & rh)
 // archiving
 //////////
 
-DEFAULT_ARCHIVING(mul)
-
 //////////
 // functions overriding virtual functions from base classes
 //////////
@@ -287,6 +278,16 @@ bool mul::info(unsigned inf) const
                case info_flags::integer_polynomial:
                case info_flags::cinteger_polynomial:
                case info_flags::rational_polynomial:
+               case info_flags::real:
+               case info_flags::rational:
+               case info_flags::integer:
+               case info_flags::crational:
+               case info_flags::cinteger:
+               case info_flags::positive:
+               case info_flags::nonnegative:
+               case info_flags::posint:
+               case info_flags::nonnegint:
+               case info_flags::even:
                case info_flags::crational_polynomial:
                case info_flags::rational_function: {
                        epvector::const_iterator i = seq.begin(), end = seq.end();
@@ -295,6 +296,8 @@ bool mul::info(unsigned inf) const
                                        return false;
                                ++i;
                        }
+                       if (overall_coeff.is_equal(*_num1_p) && inf == info_flags::even)
+                               return true;
                        return overall_coeff.info(inf);
                }
                case info_flags::algebraic: {
@@ -306,10 +309,55 @@ bool mul::info(unsigned inf) const
                        }
                        return false;
                }
+               case info_flags::negative: {
+                       bool neg = false;
+                       epvector::const_iterator i = seq.begin(), end = seq.end();
+                       while (i != end) {
+                               const ex& factor = recombine_pair_to_ex(*i++);
+                               if (factor.info(info_flags::positive))
+                                       continue;
+                               else if (factor.info(info_flags::negative))
+                                       neg = !neg;
+                               else
+                                       return false;
+                       }
+                       if (overall_coeff.info(info_flags::negative))
+                               neg = !neg;
+                       return neg;
+               }
+               case info_flags::negint: {
+                       bool neg = false;
+                       epvector::const_iterator i = seq.begin(), end = seq.end();
+                       while (i != end) {
+                               const ex& factor = recombine_pair_to_ex(*i++);
+                               if (factor.info(info_flags::posint))
+                                       continue;
+                               else if (factor.info(info_flags::negint))
+                                       neg = !neg;
+                               else
+                                       return false;
+                       }
+                       if (overall_coeff.info(info_flags::negint))
+                               neg = !neg;
+                       else if (!overall_coeff.info(info_flags::posint))
+                               return false;
+                       return neg;
+               }
        }
        return inherited::info(inf);
 }
 
+bool mul::is_polynomial(const ex & var) const
+{
+       for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+               if (!i->rest.is_polynomial(var) ||
+                   (i->rest.has(var) && !i->coeff.info(info_flags::integer))) {
+                       return false;
+               }
+       }
+       return true;
+}
+
 int mul::degree(const ex & s) const
 {
        // Sum up degrees of factors
@@ -469,8 +517,8 @@ ex mul::eval(int level) const
 
                        // XXX: What is the best way to check if the polynomial is a primitive? 
                        numeric c = i->rest.integer_content();
-                       const numeric& lead_coeff =
-                               ex_to<numeric>(ex_to<add>(i->rest).seq.begin()->coeff).div_dyn(c);
+                       const numeric lead_coeff =
+                               ex_to<numeric>(ex_to<add>(i->rest).seq.begin()->coeff).div(c);
                        const bool canonicalizable = lead_coeff.is_integer();
 
                        // XXX: The main variable is chosen in a random way, so this code 
@@ -637,7 +685,7 @@ ex mul::eval_ncmul(const exvector & v) const
        return inherited::eval_ncmul(v);
 }
 
-bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatches, lst & repls)
+bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatches, exmap& repls)
 {      
        ex origbase;
        int origexponent;
@@ -669,7 +717,7 @@ bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatch
                patternexpsign = 1;
        }
 
-       lst saverepls = repls;
+       exmap saverepls = repls;
        if (origexponent < patternexponent || origexpsign != patternexpsign || !origbase.match(patternbase,saverepls))
                return false;
        repls = saverepls;
@@ -688,17 +736,20 @@ bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatch
   * that already have been replaced by previous substitutions and matched[i]
   * is true for factors that have been matched by the current match.
   */
-bool algebraic_match_mul_with_mul(const mul &e, const ex &pat, lst &repls,
+bool algebraic_match_mul_with_mul(const mul &e, const ex &pat, exmap& repls,
                int factor, int &nummatches, const std::vector<bool> &subsed,
                std::vector<bool> &matched)
 {
-       if (factor == pat.nops())
+       GINAC_ASSERT(subsed.size() == e.nops());
+       GINAC_ASSERT(matched.size() == e.nops());
+
+       if (factor == (int)pat.nops())
                return true;
 
        for (size_t i=0; i<e.nops(); ++i) {
                if(subsed[i] || matched[i])
                        continue;
-               lst newrepls = repls;
+               exmap newrepls = repls;
                int newnummatches = nummatches;
                if (tryfactsubs(e.op(i), pat.op(factor), newnummatches, newrepls)) {
                        matched[i] = true;
@@ -721,10 +772,10 @@ bool mul::has(const ex & pattern, unsigned options) const
        if(!(options&has_options::algebraic))
                return basic::has(pattern,options);
        if(is_a<mul>(pattern)) {
-               lst repls;
+               exmap repls;
                int nummatches = std::numeric_limits<int>::max();
-               std::vector<bool> subsed(seq.size(), false);
-               std::vector<bool> matched(seq.size(), false);
+               std::vector<bool> subsed(nops(), false);
+               std::vector<bool> matched(nops(), false);
                if(algebraic_match_mul_with_mul(*this, pattern, repls, 0, nummatches,
                                subsed, matched))
                        return true;
@@ -734,8 +785,7 @@ bool mul::has(const ex & pattern, unsigned options) const
 
 ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const
 {      
-       std::vector<bool> subsed(seq.size(), false);
-       exvector subsresult(seq.size());
+       std::vector<bool> subsed(nops(), false);
        ex divide_by = 1;
        ex multiply_by = 1;
 
@@ -744,8 +794,8 @@ ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const
                if (is_exactly_a<mul>(it->first)) {
 retry1:
                        int nummatches = std::numeric_limits<int>::max();
-                       std::vector<bool> currsubsed(seq.size(), false);
-                       lst repls;
+                       std::vector<bool> currsubsed(nops(), false);
+                       exmap repls;
                        
                        if(!algebraic_match_mul_with_mul(*this, it->first, repls, 0, nummatches, subsed, currsubsed))
                                continue;
@@ -754,10 +804,10 @@ retry1:
                                if (currsubsed[j])
                                        subsed[j] = true;
                        ex subsed_pattern
-                               = it->first.subs(ex(repls), subs_options::no_pattern);
+                               = it->first.subs(repls, subs_options::no_pattern);
                        divide_by *= power(subsed_pattern, nummatches);
                        ex subsed_result
-                               = it->second.subs(ex(repls), subs_options::no_pattern);
+                               = it->second.subs(repls, subs_options::no_pattern);
                        multiply_by *= power(subsed_result, nummatches);
                        goto retry1;
 
@@ -765,14 +815,14 @@ retry1:
 
                        for (size_t j=0; j<this->nops(); j++) {
                                int nummatches = std::numeric_limits<int>::max();
-                               lst repls;
+                               exmap repls;
                                if (!subsed[j] && tryfactsubs(op(j), it->first, nummatches, repls)){
                                        subsed[j] = true;
                                        ex subsed_pattern
-                                               = it->first.subs(ex(repls), subs_options::no_pattern);
+                                               = it->first.subs(repls, subs_options::no_pattern);
                                        divide_by *= power(subsed_pattern, nummatches);
                                        ex subsed_result
-                                               = it->second.subs(ex(repls), subs_options::no_pattern);
+                                               = it->second.subs(repls, subs_options::no_pattern);
                                        multiply_by *= power(subsed_result, nummatches);
                                }
                        }
@@ -792,6 +842,38 @@ retry1:
        return ((*this)/divide_by)*multiply_by;
 }
 
+ex mul::conjugate() const
+{
+       // The base class' method is wrong here because we have to be careful at
+       // branch cuts. power::conjugate takes care of that already, so use it.
+       epvector *newepv = 0;
+       for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+               if (newepv) {
+                       newepv->push_back(split_ex_to_pair(recombine_pair_to_ex(*i).conjugate()));
+                       continue;
+               }
+               ex x = recombine_pair_to_ex(*i);
+               ex c = x.conjugate();
+               if (c.is_equal(x)) {
+                       continue;
+               }
+               newepv = new epvector;
+               newepv->reserve(seq.size());
+               for (epvector::const_iterator j=seq.begin(); j!=i; ++j) {
+                       newepv->push_back(*j);
+               }
+               newepv->push_back(split_ex_to_pair(c));
+       }
+       ex x = overall_coeff.conjugate();
+       if (!newepv && are_ex_trivially_equal(x, overall_coeff)) {
+               return *this;
+       }
+       ex result = thisexpairseq(newepv ? *newepv : seq, x);
+       delete newepv;
+       return result;
+}
+
+
 // protected
 
 /** Implementation of ex::diff() for a product.  It applies the product rule.
@@ -855,10 +937,10 @@ unsigned mul::return_type() const
        return all_commutative ? return_types::commutative : return_types::noncommutative;
 }
    
-tinfo_t mul::return_type_tinfo() const
+return_type_t mul::return_type_tinfo() const
 {
        if (seq.empty())
-               return this;  // mul without factors: should not happen
+               return make_return_type_t<mul>(); // mul without factors: should not happen
        
        // return type_info of first noncommutative element
        epvector::const_iterator i = seq.begin(), end = seq.end();
@@ -868,7 +950,7 @@ tinfo_t mul::return_type_tinfo() const
                ++i;
        }
        // no noncommutative element found, should not happen
-       return this;
+       return make_return_type_t<mul>();
 }
 
 ex mul::thisexpairseq(const epvector & v, const ex & oc, bool do_index_renaming) const
@@ -1217,4 +1299,6 @@ std::auto_ptr<epvector> mul::expandchildren(unsigned options) const
        return std::auto_ptr<epvector>(0); // nothing has changed
 }
 
+GINAC_BIND_UNARCHIVER(mul);
+
 } // namespace GiNaC