]> www.ginac.de Git - ginac.git/blobdiff - ginac/mul.cpp
Print -x as -x instead of -1.0*x when printing C-source.
[ginac.git] / ginac / mul.cpp
index 962947898baf696c245ab816866fb286a39c5fdf..63cc7e72a39e62c29e3e2f7444e458c80e3a5b82 100644 (file)
@@ -27,8 +27,6 @@
 
 #include "mul.h"
 #include "add.h"
-#include "color.h"
-#include "clifford.h"
 #include "power.h"
 #include "operators.h"
 #include "matrix.h"
@@ -86,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());
 }
 
@@ -220,8 +218,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 "/"
@@ -463,6 +465,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
@@ -561,7 +598,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
@@ -597,10 +634,28 @@ bool algebraic_match_mul_with_mul(const mul &e, const ex &pat, lst &repls,
        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)) {
+               lst repls;
+               int nummatches = std::numeric_limits<int>::max();
+               std::vector<bool> subsed(seq.size(), false);
+               std::vector<bool> matched(seq.size(), 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());
+       ex divide_by = 1;
+       ex multiply_by = 1;
 
        for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
 
@@ -614,30 +669,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);
                                }
                        }
                }
@@ -653,15 +708,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
@@ -716,21 +763,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()->tinfo() == &clifford::tinfo_static) {
-                               if (i->rest.return_type_tinfo()->tinfo() != &clifford::tinfo_static ||
-                                   ((clifford*)(noncommutative_element->rest.return_type_tinfo()))->get_representation_label() !=
-                                   ((clifford*)(i->rest.return_type_tinfo()))->get_representation_label()) {
-                                       // diffent types -> mul is ncc
-                                       return return_types::noncommutative_composite;
-                               }
-                       } else if (noncommutative_element->rest.return_type_tinfo()->tinfo() == &color::tinfo_static) {
-                               if (i->rest.return_type_tinfo()->tinfo() != &color::tinfo_static ||
-                                   ((color*)(noncommutative_element->rest.return_type_tinfo()))->get_representation_label() !=
-                                   ((color*)(i->rest.return_type_tinfo()))->get_representation_label()) {
-                                       // diffent types -> mul is ncc
-                                       return return_types::noncommutative_composite;
-                               }
-                       } else if (noncommutative_element->rest.return_type_tinfo()->tinfo() != i->rest.return_type_tinfo()->tinfo()) {
+                       if (noncommutative_element->rest.return_type_tinfo() != i->rest.return_type_tinfo()) {
+                                       // different types -> mul is ncc
                                        return return_types::noncommutative_composite;
                        }
                }
@@ -740,7 +774,7 @@ unsigned mul::return_type() const
        return all_commutative ? return_types::commutative : return_types::noncommutative;
 }
    
-const basic* mul::return_type_tinfo() const
+tinfo_t mul::return_type_tinfo() const
 {
        if (seq.empty())
                return this;  // mul without factors: should not happen
@@ -756,14 +790,14 @@ const basic* 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
@@ -935,11 +969,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());
                                }
 
@@ -986,7 +1020,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) {