Fix some apparent typos and misindentations 2.
[ginac.git] / ginac / mul.cpp
index 19dd2bdf406a8daa4e2b6689f6e0fcc7479a2df6..2b404877c9840aa53a071e712aa222e51680533a 100644 (file)
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's products of expressions. */
 
 /*
- *  GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2015 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
  *
  *  You should have received a copy of the GNU General Public License
  *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *  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 "operators.h"
 #include "matrix.h"
+#include "indexed.h"
 #include "lst.h"
 #include "archive.h"
 #include "utils.h"
+#include "symbol.h"
+#include "compiler.h"
+
+#include <iostream>
+#include <limits>
+#include <stdexcept>
+#include <vector>
 
 namespace GiNaC {
 
@@ -40,7 +43,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mul, expairseq,
   print_func<print_context>(&mul::do_print).
   print_func<print_latex>(&mul::do_print_latex).
   print_func<print_csrc>(&mul::do_print_csrc).
-  print_func<print_tree>(&inherited::do_print_tree).
+  print_func<print_tree>(&mul::do_print_tree).
   print_func<print_python_repr>(&mul::do_print_python_repr))
 
 
@@ -50,7 +53,6 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mul, expairseq,
 
 mul::mul()
 {
-       tinfo_key = TINFO_mul;
 }
 
 //////////
@@ -61,7 +63,6 @@ mul::mul()
 
 mul::mul(const ex & lh, const ex & rh)
 {
-       tinfo_key = TINFO_mul;
        overall_coeff = _ex1;
        construct_from_2_ex(lh,rh);
        GINAC_ASSERT(is_canonical());
@@ -69,7 +70,6 @@ mul::mul(const ex & lh, const ex & rh)
 
 mul::mul(const exvector & v)
 {
-       tinfo_key = TINFO_mul;
        overall_coeff = _ex1;
        construct_from_exvector(v);
        GINAC_ASSERT(is_canonical());
@@ -77,32 +77,28 @@ mul::mul(const exvector & v)
 
 mul::mul(const epvector & v)
 {
-       tinfo_key = TINFO_mul;
        overall_coeff = _ex1;
        construct_from_epvector(v);
        GINAC_ASSERT(is_canonical());
 }
 
-mul::mul(const epvector & v, const ex & oc)
+mul::mul(const epvector & v, const ex & oc, bool do_index_renaming)
 {
-       tinfo_key = TINFO_mul;
        overall_coeff = oc;
-       construct_from_epvector(v);
+       construct_from_epvector(v, do_index_renaming);
        GINAC_ASSERT(is_canonical());
 }
 
-mul::mul(std::auto_ptr<epvector> vp, const ex & oc)
+mul::mul(std::auto_ptr<epvector> vp, const ex & oc, bool do_index_renaming)
 {
-       tinfo_key = TINFO_mul;
-       GINAC_ASSERT(vp!=0);
+       GINAC_ASSERT(vp.get()!=0);
        overall_coeff = oc;
-       construct_from_epvector(*vp);
+       construct_from_epvector(*vp, do_index_renaming);
        GINAC_ASSERT(is_canonical());
 }
 
 mul::mul(const ex & lh, const ex & mh, const ex & rh)
 {
-       tinfo_key = TINFO_mul;
        exvector factors;
        factors.reserve(3);
        factors.push_back(lh);
@@ -117,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
 //////////
@@ -128,8 +122,8 @@ void mul::print_overall_coeff(const print_context & c, const char *mul_sym) cons
        const numeric &coeff = ex_to<numeric>(overall_coeff);
        if (coeff.csgn() == -1)
                c.s << '-';
-       if (!coeff.is_equal(_num1) &&
-               !coeff.is_equal(_num_1)) {
+       if (!coeff.is_equal(*_num1_p) &&
+               !coeff.is_equal(*_num_1_p)) {
                if (coeff.is_rational()) {
                        if (coeff.is_negative())
                                (-coeff).print(c);
@@ -217,8 +211,12 @@ void mul::do_print_csrc(const print_csrc & c, unsigned level) const
                c.s << "(";
 
        if (!overall_coeff.is_equal(_ex1)) {
-               overall_coeff.print(c, precedence());
-               c.s << "*";
+               if (overall_coeff.is_equal(_ex_1))
+                       c.s << "-";
+               else {
+                       overall_coeff.print(c, precedence());
+                       c.s << "*";
+               }
        }
 
        // Print arguments, separated by "*" or "/"
@@ -280,6 +278,12 @@ 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::even:
                case info_flags::crational_polynomial:
                case info_flags::rational_function: {
                        epvector::const_iterator i = seq.begin(), end = seq.end();
@@ -288,6 +292,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: {
@@ -299,10 +305,114 @@ bool mul::info(unsigned inf) const
                        }
                        return false;
                }
+               case info_flags::positive:
+               case info_flags::negative: {
+                       if ((inf==info_flags::positive) && (flags & status_flags::is_positive))
+                               return true;
+                       else if ((inf==info_flags::negative) && (flags & status_flags::is_negative))
+                               return true;
+                       if (flags & status_flags::purely_indefinite)
+                               return false;
+
+                       bool pos = true;
+                       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))
+                                       pos = !pos;
+                               else
+                                       return false;
+                       }
+                       if (overall_coeff.info(info_flags::negative))
+                               pos = !pos;
+                       setflag(pos ? status_flags::is_positive : status_flags::is_negative);
+                       return (inf == info_flags::positive? pos : !pos);
+               }
+               case info_flags::nonnegative: {
+                       if  (flags & status_flags::is_positive)
+                               return true;
+                       bool pos = true;
+                       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::nonnegative) || factor.info(info_flags::positive))
+                                       continue;
+                               else if (factor.info(info_flags::negative))
+                                       pos = !pos;
+                               else
+                                       return false;
+                       }
+                       return (overall_coeff.info(info_flags::negative)? !pos : pos);
+               }
+               case info_flags::posint:
+               case info_flags::negint: {
+                       bool pos = true;
+                       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))
+                                       pos = !pos;
+                               else
+                                       return false;
+                       }
+                       if (overall_coeff.info(info_flags::negint))
+                               pos = !pos;
+                       else if (!overall_coeff.info(info_flags::posint))
+                               return false;
+                       return (inf ==info_flags::posint? pos : !pos); 
+               }
+               case info_flags::nonnegint: {
+                       bool pos = true;
+                       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::nonnegint) || factor.info(info_flags::posint))
+                                       continue;
+                               else if (factor.info(info_flags::negint))
+                                       pos = !pos;
+                               else
+                                       return false;
+                       }
+                       if (overall_coeff.info(info_flags::negint))
+                               pos = !pos;
+                       else if (!overall_coeff.info(info_flags::posint))
+                               return false;
+                       return pos; 
+               }
+               case info_flags::indefinite: {
+                       if (flags & status_flags::purely_indefinite)
+                               return true;
+                       if (flags & (status_flags::is_positive | status_flags::is_negative))
+                               return false;
+                       epvector::const_iterator i = seq.begin(), end = seq.end();
+                       while (i != end) {
+                               const ex& term = recombine_pair_to_ex(*i);
+                               if (term.info(info_flags::positive) || term.info(info_flags::negative))
+                                       return false;
+                               ++i;
+                       }
+                       setflag(status_flags::purely_indefinite);
+                       return true;
+               }
        }
        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::nonnegint))) {
+                       return false;
+               }
+       }
+       return true;
+}
+
 int mul::degree(const ex & s) const
 {
        // Sum up degrees of factors
@@ -310,7 +420,11 @@ int mul::degree(const ex & s) const
        epvector::const_iterator i = seq.begin(), end = seq.end();
        while (i != end) {
                if (ex_to<numeric>(i->coeff).is_integer())
-                       deg_sum += i->rest.degree(s) * ex_to<numeric>(i->coeff).to_int();
+                       deg_sum += recombine_pair_to_ex(*i).degree(s);
+               else {
+                       if (i->rest.has(s))
+                               throw std::runtime_error("mul::degree() undefined degree because of non-integer exponent");
+               }
                ++i;
        }
        return deg_sum;
@@ -323,7 +437,11 @@ int mul::ldegree(const ex & s) const
        epvector::const_iterator i = seq.begin(), end = seq.end();
        while (i != end) {
                if (ex_to<numeric>(i->coeff).is_integer())
-                       deg_sum += i->rest.ldegree(s) * ex_to<numeric>(i->coeff).to_int();
+                       deg_sum += recombine_pair_to_ex(*i).ldegree(s);
+               else {
+                       if (i->rest.has(s))
+                               throw std::runtime_error("mul::ldegree() undefined degree because of non-integer exponent");
+               }
                ++i;
        }
        return deg_sum;
@@ -385,31 +503,13 @@ ex mul::eval(int level) const
                           setflag(status_flags::dynallocated);
        }
        
-#ifdef DO_GINAC_ASSERT
-       epvector::const_iterator i = seq.begin(), end = seq.end();
-       while (i != end) {
-               GINAC_ASSERT((!is_exactly_a<mul>(i->rest)) ||
-                            (!(ex_to<numeric>(i->coeff).is_integer())));
-               GINAC_ASSERT(!(i->is_canonical_numeric()));
-               if (is_exactly_a<numeric>(recombine_pair_to_ex(*i)))
-                   print(print_tree(std::cerr));
-               GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(*i)));
-               /* for paranoia */
-               expair p = split_ex_to_pair(recombine_pair_to_ex(*i));
-               GINAC_ASSERT(p.rest.is_equal(i->rest));
-               GINAC_ASSERT(p.coeff.is_equal(i->coeff));
-               /* end paranoia */
-               ++i;
-       }
-#endif // def DO_GINAC_ASSERT
-       
        if (flags & status_flags::evaluated) {
                GINAC_ASSERT(seq.size()>0);
                GINAC_ASSERT(seq.size()>1 || !overall_coeff.is_equal(_ex1));
                return *this;
        }
        
-       int seq_size = seq.size();
+       size_t seq_size = seq.size();
        if (overall_coeff.is_zero()) {
                // *(...,x;0) -> 0
                return _ex0;
@@ -421,7 +521,7 @@ ex mul::eval(int level) const
                return recombine_pair_to_ex(*(seq.begin()));
        } else if ((seq_size==1) &&
                   is_exactly_a<add>((*seq.begin()).rest) &&
-                  ex_to<numeric>((*seq.begin()).coeff).is_equal(_num1)) {
+                  ex_to<numeric>((*seq.begin()).coeff).is_equal(*_num1_p)) {
                // *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
                const add & addref = ex_to<add>((*seq.begin()).rest);
                std::auto_ptr<epvector> distrseq(new epvector);
@@ -433,9 +533,80 @@ ex mul::eval(int level) const
                }
                return (new add(distrseq,
                                ex_to<numeric>(addref.overall_coeff).
-                               mul_dyn(ex_to<numeric>(overall_coeff))))
-                     ->setflag(status_flags::dynallocated | status_flags::evaluated);
+                               mul_dyn(ex_to<numeric>(overall_coeff)))
+                      )->setflag(status_flags::dynallocated | status_flags::evaluated);
+       } else if ((seq_size >= 2) && (! (flags & status_flags::expanded))) {
+               // Strip the content and the unit part from each term. Thus
+               // things like (-x+a)*(3*x-3*a) automagically turn into - 3*(x-a)^2
+
+               epvector::const_iterator last = seq.end();
+               epvector::const_iterator i = seq.begin();
+               epvector::const_iterator j = seq.begin();
+               std::auto_ptr<epvector> s(new epvector);
+               numeric oc = *_num1_p;
+               bool something_changed = false;
+               while (i!=last) {
+                       if (likely(! (is_a<add>(i->rest) && i->coeff.is_equal(_ex1)))) {
+                               // power::eval has such a rule, no need to handle powers here
+                               ++i;
+                               continue;
+                       }
+
+                       // 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(c);
+                       const bool canonicalizable = lead_coeff.is_integer();
+
+                       // XXX: The main variable is chosen in a random way, so this code 
+                       // does NOT transform the term into the canonical form (thus, in some
+                       // very unlucky event it can even loop forever). Hopefully the main
+                       // variable will be the same for all terms in *this
+                       const bool unit_normal = lead_coeff.is_pos_integer();
+                       if (likely((c == *_num1_p) && ((! canonicalizable) || unit_normal))) {
+                               ++i;
+                               continue;
+                       }
+
+                       if (! something_changed) {
+                               s->reserve(seq_size);
+                               something_changed = true;
+                       }
+
+                       while ((j!=i) && (j!=last)) {
+                               s->push_back(*j);
+                               ++j;
+                       }
+
+                       if (! unit_normal)
+                               c = c.mul(*_num_1_p);
+
+                       oc = oc.mul(c);
+
+                       // divide add by the number in place to save at least 2 .eval() calls
+                       const add& addref = ex_to<add>(i->rest);
+                       add* primitive = new add(addref);
+                       primitive->setflag(status_flags::dynallocated);
+                       primitive->clearflag(status_flags::hash_calculated);
+                       primitive->overall_coeff = ex_to<numeric>(primitive->overall_coeff).div_dyn(c);
+                       for (epvector::iterator ai = primitive->seq.begin(); ai != primitive->seq.end(); ++ai)
+                               ai->coeff = ex_to<numeric>(ai->coeff).div_dyn(c);
+                       
+                       s->push_back(expair(*primitive, _ex1));
+
+                       ++i;
+                       ++j;
+               }
+               if (something_changed) {
+                       while (j!=last) {
+                               s->push_back(*j);
+                               ++j;
+                       }
+                       return (new mul(s, ex_to<numeric>(overall_coeff).mul_dyn(oc))
+                              )->setflag(status_flags::dynallocated);
+               }
        }
+
        return this->hold();
 }
 
@@ -460,6 +631,41 @@ ex mul::evalf(int level) const
        return mul(s, overall_coeff.evalf(level));
 }
 
+void mul::find_real_imag(ex & rp, ex & ip) const
+{
+       rp = overall_coeff.real_part();
+       ip = overall_coeff.imag_part();
+       for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+               ex factor = recombine_pair_to_ex(*i);
+               ex new_rp = factor.real_part();
+               ex new_ip = factor.imag_part();
+               if(new_ip.is_zero()) {
+                       rp *= new_rp;
+                       ip *= new_rp;
+               } else {
+                       ex temp = rp*new_rp - ip*new_ip;
+                       ip = ip*new_rp + rp*new_ip;
+                       rp = temp;
+               }
+       }
+       rp = rp.expand();
+       ip = ip.expand();
+}
+
+ex mul::real_part() const
+{
+       ex rp, ip;
+       find_real_imag(rp, ip);
+       return rp;
+}
+
+ex mul::imag_part() const
+{
+       ex rp, ip;
+       find_real_imag(rp, ip);
+       return ip;
+}
+
 ex mul::evalm() const
 {
        // numeric*matrix
@@ -515,7 +721,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;
@@ -547,7 +753,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;
@@ -558,61 +764,102 @@ bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatch
        return true;
 }
 
+/** Checks whether e matches to the pattern pat and the (possibly to be updated)
+  * list of replacements repls. This matching is in the sense of algebraic
+  * substitutions. Matching starts with pat.op(factor) of the pattern because
+  * the factors before this one have already been matched. The (possibly
+  * updated) number of matches is in nummatches. subsed[i] is true for factors
+  * 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, exmap& repls,
+                                  int factor, int &nummatches, const std::vector<bool> &subsed,
+                                  std::vector<bool> &matched)
+{
+       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;
+               exmap newrepls = repls;
+               int newnummatches = nummatches;
+               if (tryfactsubs(e.op(i), pat.op(factor), newnummatches, newrepls)) {
+                       matched[i] = true;
+                       if (algebraic_match_mul_with_mul(e, pat, newrepls, factor+1,
+                                       newnummatches, subsed, matched)) {
+                               repls = newrepls;
+                               nummatches = newnummatches;
+                               return true;
+                       }
+                       else
+                               matched[i] = false;
+               }
+       }
+
+       return false;
+}
+
+bool mul::has(const ex & pattern, unsigned options) const
+{
+       if(!(options & has_options::algebraic))
+               return basic::has(pattern,options);
+       if(is_a<mul>(pattern)) {
+               exmap repls;
+               int nummatches = std::numeric_limits<int>::max();
+               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;
+       }
+       return basic::has(pattern, options);
+}
+
 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;
 
        for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
 
                if (is_exactly_a<mul>(it->first)) {
-
+retry1:
                        int nummatches = std::numeric_limits<int>::max();
-                       std::vector<bool> currsubsed(seq.size(), false);
-                       bool succeed = true;
-                       lst repls;
-
-                       for (size_t j=0; j<it->first.nops(); j++) {
-                               bool found=false;
-                               for (size_t k=0; k<nops(); k++) {
-                                       if (currsubsed[k] || subsed[k])
-                                               continue;
-                                       if (tryfactsubs(op(k), it->first.op(j), nummatches, repls)) {
-                                               currsubsed[k] = true;
-                                               found = true;
-                                               break;
-                                       }
-                               }
-                               if (!found) {
-                                       succeed = false;
-                                       break;
-                               }
-                       }
-                       if (!succeed)
+                       std::vector<bool> currsubsed(nops(), false);
+                       exmap repls;
+                       
+                       if(!algebraic_match_mul_with_mul(*this, it->first, repls, 0, nummatches, subsed, currsubsed))
                                continue;
 
-                       bool foundfirstsubsedfactor = false;
-                       for (size_t j=0; j<subsed.size(); j++) {
-                               if (currsubsed[j]) {
-                                       if (foundfirstsubsedfactor)
-                                               subsresult[j] = op(j);
-                                       else {
-                                               foundfirstsubsedfactor = true;
-                                               subsresult[j] = op(j) * power(it->second.subs(ex(repls), subs_options::no_pattern) / it->first.subs(ex(repls), subs_options::no_pattern), nummatches);
-                                       }
+                       for (size_t j=0; j<subsed.size(); j++)
+                               if (currsubsed[j])
                                        subsed[j] = true;
-                               }
-                       }
+                       ex subsed_pattern
+                               = it->first.subs(repls, subs_options::no_pattern);
+                       divide_by *= power(subsed_pattern, nummatches);
+                       ex subsed_result
+                               = it->second.subs(repls, subs_options::no_pattern);
+                       multiply_by *= power(subsed_result, nummatches);
+                       goto retry1;
 
                } else {
 
-                       int nummatches = std::numeric_limits<int>::max();
-                       lst repls;
-
                        for (size_t j=0; j<this->nops(); j++) {
-                               if (!subsed[j] && tryfactsubs(op(j), it->first, nummatches, repls)) {
+                               int nummatches = std::numeric_limits<int>::max();
+                               exmap repls;
+                               if (!subsed[j] && tryfactsubs(op(j), it->first, nummatches, repls)){
                                        subsed[j] = true;
-                                       subsresult[j] = op(j) * power(it->second.subs(ex(repls), subs_options::no_pattern) / it->first.subs(ex(repls), subs_options::no_pattern), nummatches);
+                                       ex subsed_pattern
+                                               = it->first.subs(repls, subs_options::no_pattern);
+                                       divide_by *= power(subsed_pattern, nummatches);
+                                       ex subsed_result
+                                               = it->second.subs(repls, subs_options::no_pattern);
+                                       multiply_by *= power(subsed_result, nummatches);
                                }
                        }
                }
@@ -628,17 +875,41 @@ ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const
        if (!subsfound)
                return subs_one_level(m, options | subs_options::algebraic);
 
-       exvector ev; ev.reserve(nops());
-       for (size_t i=0; i<nops(); i++) {
-               if (subsed[i])
-                       ev.push_back(subsresult[i]);
-               else
-                       ev.push_back(op(i));
-       }
+       return ((*this)/divide_by)*multiply_by;
+}
 
-       return (new mul(ev))->setflag(status_flags::dynallocated);
+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.
@@ -672,7 +943,7 @@ int mul::compare_same_type(const basic & other) const
 unsigned mul::return_type() const
 {
        if (seq.empty()) {
-               // mul without factors: should not happen, but commutes
+               // mul without factors: should not happen, but commutates
                return return_types::commutative;
        }
        
@@ -692,8 +963,8 @@ unsigned mul::return_type() const
                if ((rt == return_types::noncommutative) && (!all_commutative)) {
                        // another nc element found, compare type_infos
                        if (noncommutative_element->rest.return_type_tinfo() != i->rest.return_type_tinfo()) {
-                               // diffent types -> mul is ncc
-                               return return_types::noncommutative_composite;
+                                       // different types -> mul is ncc
+                                       return return_types::noncommutative_composite;
                        }
                }
                ++i;
@@ -701,11 +972,11 @@ unsigned mul::return_type() const
        // all factors checked
        return all_commutative ? return_types::commutative : return_types::noncommutative;
 }
-   
-unsigned mul::return_type_tinfo() const
+
+return_type_t mul::return_type_tinfo() const
 {
        if (seq.empty())
-               return tinfo_key;  // 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();
@@ -715,17 +986,17 @@ unsigned mul::return_type_tinfo() const
                ++i;
        }
        // no noncommutative element found, should not happen
-       return tinfo_key;
+       return make_return_type_t<mul>();
 }
 
-ex mul::thisexpairseq(const epvector & v, const ex & oc) const
+ex mul::thisexpairseq(const epvector & v, const ex & oc, bool do_index_renaming) const
 {
-       return (new mul(v, oc))->setflag(status_flags::dynallocated);
+       return (new mul(v, oc, do_index_renaming))->setflag(status_flags::dynallocated);
 }
 
-ex mul::thisexpairseq(std::auto_ptr<epvector> vp, const ex & oc) const
+ex mul::thisexpairseq(std::auto_ptr<epvector> vp, const ex & oc, bool do_index_renaming) const
 {
-       return (new mul(vp, oc))->setflag(status_flags::dynallocated);
+       return (new mul(vp, oc, do_index_renaming))->setflag(status_flags::dynallocated);
 }
 
 expair mul::split_ex_to_pair(const ex & e) const
@@ -737,7 +1008,7 @@ expair mul::split_ex_to_pair(const ex & e) const
        }
        return expair(e,_ex1);
 }
-       
+
 expair mul::combine_ex_with_coeff_to_pair(const ex & e,
                                           const ex & c) const
 {
@@ -750,7 +1021,7 @@ expair mul::combine_ex_with_coeff_to_pair(const ex & e,
 
        return split_ex_to_pair(power(e,c));
 }
-       
+
 expair mul::combine_pair_with_coeff_to_pair(const expair & p,
                                             const ex & c) const
 {
@@ -763,10 +1034,10 @@ expair mul::combine_pair_with_coeff_to_pair(const expair & p,
 
        return split_ex_to_pair(power(recombine_pair_to_ex(p),c));
 }
-       
+
 ex mul::recombine_pair_to_ex(const expair & p) const
 {
-       if (ex_to<numeric>(p.coeff).is_equal(_num1)) 
+       if (ex_to<numeric>(p.coeff).is_equal(*_num1_p)) 
                return p.rest;
        else
                return (new power(p.rest,p.coeff))->setflag(status_flags::dynallocated);
@@ -775,22 +1046,22 @@ ex mul::recombine_pair_to_ex(const expair & p) const
 bool mul::expair_needs_further_processing(epp it)
 {
        if (is_exactly_a<mul>(it->rest) &&
-               ex_to<numeric>(it->coeff).is_integer()) {
+           ex_to<numeric>(it->coeff).is_integer()) {
                // combined pair is product with integer power -> expand it
                *it = split_ex_to_pair(recombine_pair_to_ex(*it));
                return true;
        }
        if (is_exactly_a<numeric>(it->rest)) {
+               if (it->coeff.is_equal(_ex1)) {
+                       // pair has coeff 1 and must be moved to the end
+                       return true;
+               }
                expair ep = split_ex_to_pair(recombine_pair_to_ex(*it));
                if (!ep.is_equal(*it)) {
                        // combined pair is a numeric power which can be simplified
                        *it = ep;
                        return true;
                }
-               if (it->coeff.is_equal(_ex1)) {
-                       // combined pair has coeff 1 and must be moved to the end
-                       return true;
-               }
        }
        return false;
 }       
@@ -821,7 +1092,7 @@ bool mul::can_make_flat(const expair & p) const
        // this assertion will probably fail somewhere
        // it would require a more careful make_flat, obeying the power laws
        // probably should return true only if p.coeff is integer
-       return ex_to<numeric>(p.coeff).is_equal(_num1);
+       return ex_to<numeric>(p.coeff).is_equal(*_num1_p);
 }
 
 bool mul::can_be_further_expanded(const ex & e)
@@ -840,6 +1111,24 @@ bool mul::can_be_further_expanded(const ex & e)
 
 ex mul::expand(unsigned options) const
 {
+       {
+       // trivial case: expanding the monomial (~ 30% of all calls)
+               epvector::const_iterator i = seq.begin(), seq_end = seq.end();
+               while ((i != seq.end()) &&  is_a<symbol>(i->rest) && i->coeff.info(info_flags::integer))
+                       ++i;
+               if (i == seq_end) {
+                       setflag(status_flags::expanded);
+                       return *this;
+               }
+       }
+
+       // do not rename indices if the object has no indices at all
+       if ((!(options & expand_options::expand_rename_idx)) && 
+                       this->info(info_flags::has_indices))
+               options |= expand_options::expand_rename_idx;
+
+       const bool skip_idx_rename = !(options & expand_options::expand_rename_idx);
+
        // First, expand the children
        std::auto_ptr<epvector> expanded_seqp = expandchildren(options);
        const epvector & expanded_seq = (expanded_seqp.get() ? *expanded_seqp : seq);
@@ -848,7 +1137,6 @@ ex mul::expand(unsigned options) const
        // with the next one that is found while collecting the factors which are
        // not sums
        ex last_expanded = _ex1;
-       bool need_reexpand = false;
 
        epvector non_adds;
        non_adds.reserve(expanded_seq.size());
@@ -895,25 +1183,48 @@ ex mul::expand(unsigned options) const
                                // Compute the new overall coefficient and put it together:
                                ex tmp_accu = (new add(distrseq, add1.overall_coeff*add2.overall_coeff))->setflag(status_flags::dynallocated);
 
+                               exvector add1_dummy_indices, add2_dummy_indices, add_indices;
+                               lst dummy_subs;
+
+                               if (!skip_idx_rename) {
+                                       for (epvector::const_iterator i=add1begin; i!=add1end; ++i) {
+                                               add_indices = get_all_dummy_indices_safely(i->rest);
+                                               add1_dummy_indices.insert(add1_dummy_indices.end(), add_indices.begin(), add_indices.end());
+                                       }
+                                       for (epvector::const_iterator i=add2begin; i!=add2end; ++i) {
+                                               add_indices = get_all_dummy_indices_safely(i->rest);
+                                               add2_dummy_indices.insert(add2_dummy_indices.end(), add_indices.begin(), add_indices.end());
+                                       }
+
+                                       sort(add1_dummy_indices.begin(), add1_dummy_indices.end(), ex_is_less());
+                                       sort(add2_dummy_indices.begin(), add2_dummy_indices.end(), ex_is_less());
+                                       dummy_subs = rename_dummy_indices_uniquely(add1_dummy_indices, add2_dummy_indices);
+                               }
+
                                // Multiply explicitly all non-numeric terms of add1 and add2:
-                               for (epvector::const_iterator i1=add1begin; i1!=add1end; ++i1) {
+                               for (epvector::const_iterator i2=add2begin; i2!=add2end; ++i2) {
                                        // We really have to combine terms here in order to compactify
                                        // the result.  Otherwise it would become waayy tooo bigg.
-                                       numeric oc;
-                                       distrseq.clear();
-                                       for (epvector::const_iterator i2=add2begin; i2!=add2end; ++i2) {
+                                       numeric oc(*_num0_p);
+                                       epvector distrseq2;
+                                       distrseq2.reserve(add1.seq.size());
+                                       const ex i2_new = (skip_idx_rename || (dummy_subs.op(0).nops() == 0) ?
+                                                       i2->rest :
+                                                       i2->rest.subs(ex_to<lst>(dummy_subs.op(0)), 
+                                                               ex_to<lst>(dummy_subs.op(1)), subs_options::no_pattern));
+                                       for (epvector::const_iterator i1=add1begin; i1!=add1end; ++i1) {
                                                // Don't push_back expairs which might have a rest that evaluates to a numeric,
                                                // since that would violate an invariant of expairseq:
-                                               const ex rest = (new mul(i1->rest, i2->rest))->setflag(status_flags::dynallocated);
-                                               if (is_exactly_a<numeric>(rest))
+                                               const ex rest = (new mul(i1->rest, i2_new))->setflag(status_flags::dynallocated);
+                                               if (is_exactly_a<numeric>(rest)) {
                                                        oc += ex_to<numeric>(rest).mul(ex_to<numeric>(i1->coeff).mul(ex_to<numeric>(i2->coeff)));
-                                               else
-                                                       distrseq.push_back(expair(rest, ex_to<numeric>(i1->coeff).mul_dyn(ex_to<numeric>(i2->coeff))));
+                                               } else {
+                                                       distrseq2.push_back(expair(rest, ex_to<numeric>(i1->coeff).mul_dyn(ex_to<numeric>(i2->coeff))));
+                                               }
                                        }
-                                       tmp_accu += (new add(distrseq, oc))->setflag(status_flags::dynallocated);
-                               }
+                                       tmp_accu += (new add(distrseq2, oc))->setflag(status_flags::dynallocated);
+                               } 
                                last_expanded = tmp_accu;
-
                        } else {
                                if (!last_expanded.is_equal(_ex1))
                                        non_adds.push_back(split_ex_to_pair(last_expanded));
@@ -931,14 +1242,22 @@ ex mul::expand(unsigned options) const
                size_t n = last_expanded.nops();
                exvector distrseq;
                distrseq.reserve(n);
+               exvector va;
+               if (! skip_idx_rename) {
+                       va = get_all_dummy_indices_safely(mul(non_adds));
+                       sort(va.begin(), va.end(), ex_is_less());
+               }
 
                for (size_t i=0; i<n; ++i) {
                        epvector factors = non_adds;
-                       factors.push_back(split_ex_to_pair(last_expanded.op(i)));
+                       if (skip_idx_rename)
+                               factors.push_back(split_ex_to_pair(last_expanded.op(i)));
+                       else
+                               factors.push_back(split_ex_to_pair(rename_dummy_indices_uniquely(va, last_expanded.op(i))));
                        ex term = (new mul(factors, overall_coeff))->setflag(status_flags::dynallocated);
-                       if (can_be_further_expanded(term))
+                       if (can_be_further_expanded(term)) {
                                distrseq.push_back(term.expand());
-                       else {
+                       else {
                                if (options == 0)
                                        ex_to<basic>(term).setflag(status_flags::expanded);
                                distrseq.push_back(term);
@@ -974,7 +1293,7 @@ ex mul::expand(unsigned options) const
 
 /** Member-wise expand the expairs representing this sequence.  This must be
  *  overridden from expairseq::expandchildren() and done iteratively in order
- *  to allow for early cancallations and thus safe memory.
+ *  to allow for early cancellations and thus safe memory.
  *
  *  @see mul::expand()
  *  @return pointer to epvector containing expanded representation or zero
@@ -1016,4 +1335,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