]> www.ginac.de Git - ginac.git/blobdiff - ginac/mul.cpp
Improved dummy index renaming.
[ginac.git] / ginac / mul.cpp
index dd929f7e136ffc90546550ee57e87dfdb5069dab..e5b38fd6485c1022f1707f75c02d1d20e317fb88 100644 (file)
@@ -84,20 +84,20 @@ mul::mul(const 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 = &mul::tinfo_static;
        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 = &mul::tinfo_static;
        GINAC_ASSERT(vp.get()!=0);
        overall_coeff = oc;
-       construct_from_epvector(*vp);
+       construct_from_epvector(*vp, do_index_renaming);
        GINAC_ASSERT(is_canonical());
 }
 
@@ -461,6 +461,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
@@ -559,7 +594,7 @@ bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatch
        return true;
 }
 
-/** Checks wheter e matches to the pattern pat and the (possibly to be updated
+/** Checks wheter 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
@@ -615,6 +650,8 @@ ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const
 {      
        std::vector<bool> subsed(seq.size(), false);
        exvector subsresult(seq.size());
+       ex divide_by = 1;
+       ex multiply_by = 1;
 
        for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
 
@@ -628,30 +665,30 @@ retry1:
                        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(ex(repls), subs_options::no_pattern);
+                       divide_by *= power(subsed_pattern, nummatches);
+                       ex subsed_result
+                               = it->second.subs(ex(repls), subs_options::no_pattern);
+                       multiply_by *= power(subsed_result, nummatches);
                        goto retry1;
 
                } else {
-retry2:
-                       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();
+                               lst 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);
-                                       goto retry2;
+                                       ex subsed_pattern
+                                               = it->first.subs(ex(repls), subs_options::no_pattern);
+                                       divide_by *= power(subsed_pattern, nummatches);
+                                       ex subsed_result
+                                               = it->second.subs(ex(repls), subs_options::no_pattern);
+                                       multiply_by *= power(subsed_result, nummatches);
                                }
                        }
                }
@@ -667,15 +704,7 @@ retry2:
        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 (new mul(ev))->setflag(status_flags::dynallocated);
+       return ((*this)/divide_by)*multiply_by;
 }
 
 // protected
@@ -757,14 +786,14 @@ tinfo_t mul::return_type_tinfo() const
        return this;
 }
 
-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
@@ -936,11 +965,11 @@ ex mul::expand(unsigned options) const
                                exvector add1_dummy_indices, add2_dummy_indices, add_indices;
 
                                for (epvector::const_iterator i=add1begin; i!=add1end; ++i) {
-                                       add_indices = get_all_dummy_indices(i->rest);
+                                       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(i->rest);
+                                       add_indices = get_all_dummy_indices_safely(i->rest);
                                        add2_dummy_indices.insert(add2_dummy_indices.end(), add_indices.begin(), add_indices.end());
                                }
 
@@ -987,7 +1016,7 @@ ex mul::expand(unsigned options) const
                size_t n = last_expanded.nops();
                exvector distrseq;
                distrseq.reserve(n);
-               exvector va = get_all_dummy_indices(mul(non_adds));
+               exvector 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) {