]> www.ginac.de Git - ginac.git/blobdiff - ginac/basic.cpp
match: don't modify subexpression list if expression doesn't match the pattern.
[ginac.git] / ginac / basic.cpp
index 9cb480f9f56ba06c9e4486316593da6ad570ebe7..7a0633de19a18221f241292b21742416cdb29dce 100644 (file)
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's ABC. */
 
 /*
- *  GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2008 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
@@ -30,6 +30,7 @@
 #include "ex.h"
 #include "numeric.h"
 #include "power.h"
+#include "add.h"
 #include "symbol.h"
 #include "lst.h"
 #include "ncmul.h"
@@ -38,6 +39,7 @@
 #include "wildcard.h"
 #include "archive.h"
 #include "utils.h"
+#include "inifcns.h"
 
 namespace GiNaC {
 
@@ -286,13 +288,13 @@ ex & basic::operator[](size_t i)
  *  the pattern itself or one of the children 'has' it.  As a consequence
  *  (according to the definition of children) given e=x+y+z, e.has(x) is true
  *  but e.has(x+y) is false. */
-bool basic::has(const ex & pattern) const
+bool basic::has(const ex & pattern, unsigned options) const
 {
        lst repl_lst;
        if (match(pattern, repl_lst))
                return true;
        for (size_t i=0; i<nops(); i++)
-               if (op(i).has(pattern))
+               if (op(i).has(pattern, options))
                        return true;
        
        return false;
@@ -325,6 +327,12 @@ ex basic::map(map_function & f) const
                return *this;
 }
 
+/** Check whether this is a polynomial in the given variables. */
+bool basic::is_polynomial(const ex & var) const
+{
+       return !has(var) || is_equal(ex_to<basic>(var));
+}
+
 /** Return degree of highest power in object s. */
 int basic::degree(const ex & s) const
 {
@@ -362,56 +370,32 @@ ex basic::collect(const ex & s, bool distributed) const
 
                else if (distributed) {
 
-                       // Get lower/upper degree of all symbols in list
-                       size_t num = s.nops();
-                       struct sym_info {
-                               ex sym;
-                               int ldeg, deg;
-                               int cnt;  // current degree, 'counter'
-                               ex coeff; // coefficient for degree 'cnt'
-                       };
-                       sym_info *si = new sym_info[num];
-                       ex c = *this;
-                       for (size_t i=0; i<num; i++) {
-                               si[i].sym = s.op(i);
-                               si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
-                               si[i].deg = this->degree(si[i].sym);
-                               c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
-                       }
-
-                       while (true) {
-
-                               // Calculate coeff*x1^c1*...*xn^cn
-                               ex y = _ex1;
-                               for (size_t i=0; i<num; i++) {
-                                       int cnt = si[i].cnt;
-                                       y *= power(si[i].sym, cnt);
-                               }
-                               x += y * si[num - 1].coeff;
-
-                               // Increment counters
-                               size_t n = num - 1;
-                               while (true) {
-                                       ++si[n].cnt;
-                                       if (si[n].cnt <= si[n].deg) {
-                                               // Update coefficients
-                                               ex c;
-                                               if (n == 0)
-                                                       c = *this;
-                                               else
-                                                       c = si[n - 1].coeff;
-                                               for (size_t i=n; i<num; i++)
-                                                       c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
-                                               break;
-                                       }
-                                       if (n == 0)
-                                               goto done;
-                                       si[n].cnt = si[n].ldeg;
-                                       n--;
+                       x = this->expand();
+                       if (! is_a<add>(x))
+                               return x; 
+                       const lst& l(ex_to<lst>(s));
+
+                       exmap cmap;
+                       cmap[_ex1] = _ex0;
+                       for (const_iterator xi=x.begin(); xi!=x.end(); ++xi) {
+                               ex key = _ex1;
+                               ex pre_coeff = *xi;
+                               for (lst::const_iterator li=l.begin(); li!=l.end(); ++li) {
+                                       int cexp = pre_coeff.degree(*li);
+                                       pre_coeff = pre_coeff.coeff(*li, cexp);
+                                       key *= pow(*li, cexp);
                                }
+                               exmap::iterator ci = cmap.find(key);
+                               if (ci != cmap.end())
+                                       ci->second += pre_coeff;
+                               else
+                                       cmap.insert(exmap::value_type(key, pre_coeff));
                        }
 
-done:          delete[] si;
+                       exvector resv;
+                       for (exmap::const_iterator mi=cmap.begin(); mi != cmap.end(); ++mi)
+                               resv.push_back((mi->first)*(mi->second));
+                       return (new add(resv))->setflag(status_flags::dynallocated);
 
                } else {
 
@@ -601,12 +585,18 @@ bool basic::match(const ex & pattern, lst & repl_lst) const
                if (!match_same_type(ex_to<basic>(pattern)))
                        return false;
 
+               // Even if the expression does not match the pattern, some of
+               // its subexpressions could match it. For example, x^5*y^(-1)
+               // does not match the pattern $0^5, but its subexpression x^5
+               // does. So, save repl_lst in order to not add bogus entries.
+               lst tmp_repl = repl_lst;
                // Otherwise the subexpressions must match one-to-one
                for (size_t i=0; i<nops(); i++)
-                       if (!op(i).match(pattern.op(i), repl_lst))
+                       if (!op(i).match(pattern.op(i), tmp_repl))
                                return false;
 
                // Looks similar enough, match found
+               repl_lst = tmp_repl;
                return true;
        }
 }
@@ -705,6 +695,16 @@ ex basic::conjugate() const
        return *this;
 }
 
+ex basic::real_part() const
+{
+       return real_part_function(*this).hold();
+}
+
+ex basic::imag_part() const
+{
+       return imag_part_function(*this).hold();
+}
+
 ex basic::eval_ncmul(const exvector & v) const
 {
        return hold_ncmul(v);
@@ -775,9 +775,9 @@ unsigned basic::return_type() const
        return return_types::commutative;
 }
 
-unsigned basic::return_type_tinfo() const
+tinfo_t basic::return_type_tinfo() const
 {
-       return tinfo();
+       return tinfo_key;
 }
 
 /** Compute the hash value of an object and if it makes sense to store it in
@@ -788,7 +788,7 @@ unsigned basic::return_type_tinfo() const
  *  would all end up with the same hashvalue. */
 unsigned basic::calchash() const
 {
-       unsigned v = golden_ratio_hash(tinfo());
+       unsigned v = golden_ratio_hash((p_int)tinfo());
        for (size_t i=0; i<nops(); i++) {
                v = rotate_left(v);
                v ^= this->op(i).gethash();
@@ -845,8 +845,8 @@ int basic::compare(const basic & other) const
        compare_statistics.compare_same_hashvalue++;
 #endif
 
-       const unsigned typeid_this = tinfo();
-       const unsigned typeid_other = other.tinfo();
+       const tinfo_t typeid_this = tinfo();
+       const tinfo_t typeid_other = other.tinfo();
        if (typeid_this==typeid_other) {
                GINAC_ASSERT(typeid(*this)==typeid(other));
 //             int cmpval = compare_same_type(other);