]> www.ginac.de Git - ginac.git/blobdiff - ginac/basic.cpp
Improved collecting for sparse multivariate polynomials.
[ginac.git] / ginac / basic.cpp
index 17f265debbe27f9b5c35d07b0599ed108c6469c5..c7a33ca88f870bb5e84e4c52dd266afa4ee36d2e 100644 (file)
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's ABC. */
 
 /*
- *  GiNaC Copyright (C) 1999-2004 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 <iostream>
@@ -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 {
 
@@ -57,7 +59,6 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(basic, void,
  *  tinfo_key and the hash value. */
 basic::basic(const basic & other) : tinfo_key(other.tinfo_key), flags(other.flags & ~status_flags::dynallocated), hashvalue(other.hashvalue)
 {
-       GINAC_ASSERT(typeid(*this) == typeid(other));
 }
 
 /** basic assignment operator: the other object might be of a derived class. */
@@ -208,7 +209,7 @@ void basic::do_print_python_repr(const print_python_repr & c, unsigned level) co
  *  @see basic::dbgprinttree */
 void basic::dbgprint() const
 {
-       this->print(std::cerr);
+       this->print(print_dflt(std::cerr));
        std::cerr << std::endl;
 }
 
@@ -287,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;
@@ -326,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
 {
@@ -363,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 {
 
@@ -483,6 +466,20 @@ ex basic::evalm() const
                return map(map_evalm);
 }
 
+/** Function object to be applied by basic::eval_integ(). */
+struct eval_integ_map_function : public map_function {
+       ex operator()(const ex & e) { return eval_integ(e); }
+} map_eval_integ;
+
+/** Evaluate integrals, if result is known. */
+ex basic::eval_integ() const
+{
+       if (nops() == 0)
+               return *this;
+       else
+               return map(map_eval_integ);
+}
+
 /** Perform automatic symbolic evaluations on indexed expression that
  *  contains this object as the base expression. */
 ex basic::eval_indexed(const basic & i) const
@@ -497,7 +494,7 @@ ex basic::eval_indexed(const basic & i) const
  *  (or a subclass) and their indices are compatible. This function is used
  *  internally by simplify_indexed().
  *
- *  @param self First indexed expression; it's base object is *this
+ *  @param self First indexed expression; its base object is *this
  *  @param other Second indexed expression
  *  @return sum of self and other 
  *  @see ex::simplify_indexed() */
@@ -509,7 +506,7 @@ ex basic::add_indexed(const ex & self, const ex & other) const
 /** Multiply an indexed expression with a scalar. This function is used
  *  internally by simplify_indexed().
  *
- *  @param self Indexed expression; it's base object is *this
+ *  @param self Indexed expression; its base object is *this
  *  @param other Numeric value
  *  @return product of self and other
  *  @see ex::simplify_indexed() */
@@ -525,7 +522,7 @@ ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
  *  and that at least one dummy index has been found. This functions is
  *  used internally by simplify_indexed().
  *
- *  @param self Pointer to first indexed expression; it's base object is *this
+ *  @param self Pointer to first indexed expression; its base object is *this
  *  @param other Pointer to second indexed expression
  *  @param v The complete vector of factors
  *  @return true if the contraction was successful, false otherwise
@@ -604,9 +601,11 @@ ex basic::subs_one_level(const exmap & m, unsigned options) const
        exmap::const_iterator it;
 
        if (options & subs_options::no_pattern) {
-               it = m.find(*this);
+               ex thisex = *this;
+               it = m.find(thisex);
                if (it != m.end())
                        return it->second;
+               return thisex;
        } else {
                for (it = m.begin(); it != m.end(); ++it) {
                        lst repl_lst;
@@ -690,6 +689,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);
@@ -760,9 +769,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
@@ -773,7 +782,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();
@@ -830,8 +839,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);