X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmul.cpp;h=3385bcae2b072a65bf56db5aca883205bf294875;hp=6f99c515062ff86c64a253ae9ce80adf512680bb;hb=217bee1d939abf7f9ed3fe75928a62f7fdc3abcf;hpb=6d7bf9ee5a7ce05cb3a23dae664e781d7325d7b8 diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 6f99c515..3385bcae 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's products of expressions. */ /* - * GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2006 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 @@ -17,7 +17,7 @@ * * 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 @@ -27,9 +27,12 @@ #include "mul.h" #include "add.h" +#include "color.h" +#include "clifford.h" #include "power.h" #include "operators.h" #include "matrix.h" +#include "indexed.h" #include "lst.h" #include "archive.h" #include "utils.h" @@ -40,7 +43,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mul, expairseq, print_func(&mul::do_print). print_func(&mul::do_print_latex). print_func(&mul::do_print_csrc). - print_func(&inherited::do_print_tree). + print_func(&mul::do_print_tree). print_func(&mul::do_print_python_repr)) @@ -50,7 +53,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mul, expairseq, mul::mul() { - tinfo_key = TINFO_mul; + tinfo_key = &mul::tinfo_static; } ////////// @@ -61,7 +64,7 @@ mul::mul() mul::mul(const ex & lh, const ex & rh) { - tinfo_key = TINFO_mul; + tinfo_key = &mul::tinfo_static; overall_coeff = _ex1; construct_from_2_ex(lh,rh); GINAC_ASSERT(is_canonical()); @@ -69,7 +72,7 @@ mul::mul(const ex & lh, const ex & rh) mul::mul(const exvector & v) { - tinfo_key = TINFO_mul; + tinfo_key = &mul::tinfo_static; overall_coeff = _ex1; construct_from_exvector(v); GINAC_ASSERT(is_canonical()); @@ -77,7 +80,7 @@ mul::mul(const exvector & v) mul::mul(const epvector & v) { - tinfo_key = TINFO_mul; + tinfo_key = &mul::tinfo_static; overall_coeff = _ex1; construct_from_epvector(v); GINAC_ASSERT(is_canonical()); @@ -85,7 +88,7 @@ mul::mul(const epvector & v) mul::mul(const epvector & v, const ex & oc) { - tinfo_key = TINFO_mul; + tinfo_key = &mul::tinfo_static; overall_coeff = oc; construct_from_epvector(v); GINAC_ASSERT(is_canonical()); @@ -93,8 +96,8 @@ mul::mul(const epvector & v, const ex & oc) mul::mul(std::auto_ptr vp, const ex & oc) { - tinfo_key = TINFO_mul; - GINAC_ASSERT(vp!=0); + tinfo_key = &mul::tinfo_static; + GINAC_ASSERT(vp.get()!=0); overall_coeff = oc; construct_from_epvector(*vp); GINAC_ASSERT(is_canonical()); @@ -102,7 +105,7 @@ mul::mul(std::auto_ptr vp, const ex & oc) mul::mul(const ex & lh, const ex & mh, const ex & rh) { - tinfo_key = TINFO_mul; + tinfo_key = &mul::tinfo_static; exvector factors; factors.reserve(3); factors.push_back(lh); @@ -128,8 +131,8 @@ void mul::print_overall_coeff(const print_context & c, const char *mul_sym) cons const numeric &coeff = ex_to(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); @@ -421,7 +424,7 @@ ex mul::eval(int level) const return recombine_pair_to_ex(*(seq.begin())); } else if ((seq_size==1) && is_exactly_a((*seq.begin()).rest) && - ex_to((*seq.begin()).coeff).is_equal(_num1)) { + ex_to((*seq.begin()).coeff).is_equal(*_num1_p)) { // *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +()) const add & addref = ex_to((*seq.begin()).rest); std::auto_ptr distrseq(new epvector); @@ -558,6 +561,58 @@ 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 + * 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, lst &repls, + int factor, int &nummatches, const std::vector &subsed, + std::vector &matched) +{ + if (factor == pat.nops()) + return true; + + for (size_t i=0; i(pattern)) { + lst repls; + int nummatches = std::numeric_limits::max(); + std::vector subsed(seq.size(), false); + std::vector 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 subsed(seq.size(), false); @@ -566,29 +621,13 @@ ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) { if (is_exactly_a(it->first)) { - +retry1: int nummatches = std::numeric_limits::max(); std::vector currsubsed(seq.size(), false); bool succeed = true; lst repls; - - for (size_t j=0; jfirst.nops(); j++) { - bool found=false; - for (size_t k=0; kfirst.op(j), nummatches, repls)) { - currsubsed[k] = true; - found = true; - break; - } - } - if (!found) { - succeed = false; - break; - } - } - if (!succeed) + + if(!algebraic_match_mul_with_mul(*this, it->first, repls, 0, nummatches, subsed, currsubsed)) continue; bool foundfirstsubsedfactor = false; @@ -603,9 +642,10 @@ ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const subsed[j] = true; } } + goto retry1; } else { - +retry2: int nummatches = std::numeric_limits::max(); lst repls; @@ -613,6 +653,7 @@ ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const 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; } } } @@ -672,7 +713,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; } @@ -691,9 +732,22 @@ 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; + 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()) { + return return_types::noncommutative_composite; } } ++i; @@ -702,10 +756,10 @@ unsigned mul::return_type() const return all_commutative ? return_types::commutative : return_types::noncommutative; } -unsigned mul::return_type_tinfo() const +const basic* mul::return_type_tinfo() const { if (seq.empty()) - return tinfo_key; // mul without factors: should not happen + return this; // mul without factors: should not happen // return type_info of first noncommutative element epvector::const_iterator i = seq.begin(), end = seq.end(); @@ -715,7 +769,7 @@ unsigned mul::return_type_tinfo() const ++i; } // no noncommutative element found, should not happen - return tinfo_key; + return this; } ex mul::thisexpairseq(const epvector & v, const ex & oc) const @@ -766,7 +820,7 @@ expair mul::combine_pair_with_coeff_to_pair(const expair & p, ex mul::recombine_pair_to_ex(const expair & p) const { - if (ex_to(p.coeff).is_equal(_num1)) + if (ex_to(p.coeff).is_equal(*_num1_p)) return p.rest; else return (new power(p.rest,p.coeff))->setflag(status_flags::dynallocated); @@ -821,7 +875,21 @@ 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(p.coeff).is_equal(_num1); + return ex_to(p.coeff).is_equal(*_num1_p); +} + +bool mul::can_be_further_expanded(const ex & e) +{ + if (is_exactly_a(e)) { + for (epvector::const_iterator cit = ex_to(e).seq.begin(); cit != ex_to(e).seq.end(); ++cit) { + if (is_exactly_a(cit->rest) && cit->coeff.info(info_flags::posint)) + return true; + } + } else if (is_exactly_a(e)) { + if (is_exactly_a(e.op(0)) && e.op(1).info(info_flags::posint)) + return true; + } + return false; } ex mul::expand(unsigned options) const @@ -833,15 +901,14 @@ ex mul::expand(unsigned options) const // Now, look for all the factors that are sums and multiply each one out // with the next one that is found while collecting the factors which are // not sums - int number_of_adds = 0; ex last_expanded = _ex1; + epvector non_adds; non_adds.reserve(expanded_seq.size()); - epvector::const_iterator cit = expanded_seq.begin(), last = expanded_seq.end(); - while (cit != last) { + + for (epvector::const_iterator cit = expanded_seq.begin(); cit != expanded_seq.end(); ++cit) { if (is_exactly_a(cit->rest) && (cit->coeff.is_equal(_ex1))) { - ++number_of_adds; if (is_exactly_a(last_expanded)) { // Expand a product of two sums, aggressive version. @@ -859,6 +926,7 @@ ex mul::expand(unsigned options) const const epvector::const_iterator add2end = add2.seq.end(); epvector distrseq; distrseq.reserve(add1.seq.size()+add2.seq.size()); + // Multiply add2 with the overall coefficient of add1 and append it to distrseq: if (!add1.overall_coeff.is_zero()) { if (add1.overall_coeff.is_equal(_ex1)) @@ -867,6 +935,7 @@ ex mul::expand(unsigned options) const for (epvector::const_iterator i=add2begin; i!=add2end; ++i) distrseq.push_back(expair(i->rest, ex_to(i->coeff).mul_dyn(ex_to(add1.overall_coeff)))); } + // Multiply add1 with the overall coefficient of add2 and append it to distrseq: if (!add2.overall_coeff.is_zero()) { if (add2.overall_coeff.is_equal(_ex1)) @@ -875,56 +944,93 @@ ex mul::expand(unsigned options) const for (epvector::const_iterator i=add1begin; i!=add1end; ++i) distrseq.push_back(expair(i->rest, ex_to(i->coeff).mul_dyn(ex_to(add2.overall_coeff)))); } + // 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; + + for (epvector::const_iterator i=add1begin; i!=add1end; ++i) { + add_indices = get_all_dummy_indices(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); + 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()); + lst 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) { + ex i2_new = (dummy_subs.op(0).nops()>0? + i2->rest.subs((lst)dummy_subs.op(0), (lst)dummy_subs.op(1), subs_options::no_pattern) : i2->rest); + 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(rest)) + const ex rest = (new mul(i1->rest, i2_new))->setflag(status_flags::dynallocated); + if (is_exactly_a(rest)) { oc += ex_to(rest).mul(ex_to(i1->coeff).mul(ex_to(i2->coeff))); - else + } else { distrseq.push_back(expair(rest, ex_to(i1->coeff).mul_dyn(ex_to(i2->coeff)))); + } } tmp_accu += (new add(distrseq, oc))->setflag(status_flags::dynallocated); } last_expanded = tmp_accu; } else { - non_adds.push_back(split_ex_to_pair(last_expanded)); + if (!last_expanded.is_equal(_ex1)) + non_adds.push_back(split_ex_to_pair(last_expanded)); last_expanded = cit->rest; } + } else { non_adds.push_back(*cit); } - ++cit; } - + // Now the only remaining thing to do is to multiply the factors which // were not sums into the "last_expanded" sum if (is_exactly_a(last_expanded)) { - const add & finaladd = ex_to(last_expanded); + size_t n = last_expanded.nops(); exvector distrseq; - size_t n = finaladd.nops(); distrseq.reserve(n); + exvector va = get_all_dummy_indices(mul(non_adds)); + sort(va.begin(), va.end(), ex_is_less()); + for (size_t i=0; i - setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0))); + 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)) { + distrseq.push_back(term.expand()); + } else { + if (options == 0) + ex_to(term).setflag(status_flags::expanded); + distrseq.push_back(term); + } } + return ((new add(distrseq))-> setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0))); } + non_adds.push_back(split_ex_to_pair(last_expanded)); - return (new mul(non_adds, overall_coeff))-> - setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)); + ex result = (new mul(non_adds, overall_coeff))->setflag(status_flags::dynallocated); + if (can_be_further_expanded(result)) { + return result.expand(); + } else { + if (options == 0) + ex_to(result).setflag(status_flags::expanded); + return result; + } }