* Methods of class ex which do absolutely nothing than type dispatch should

[ginac.git] / ginac / basic.cpp

diff --git a/ginac/basic.cpp b/ginac/basic.cpp
index c649f64b0e6dfd6e74b2c41f11393d8249447b7d..5b0b180378bc24df6bd95dbf671ece9f95bc6e56 100644 (file)
--- a/ginac/basic.cpp
+++ b/ginac/basic.cpp
@@ -22,6 +22,9 @@
 
 #include <iostream>
 #include <stdexcept>
+#ifdef DO_GINAC_ASSERT
+#  include <typeinfo>
+#endif
 
 #include "basic.h"
 #include "ex.h"
@@ -220,31 +223,103 @@ bool basic::has(const ex & other) const
        return false;
 }
 
-/** Return degree of highest power in symbol s. */
+/** Return degree of highest power in object s. */
 int basic::degree(const ex & s) const
 {
        return 0;
 }
 
-/** Return degree of lowest power in symbol s. */
+/** Return degree of lowest power in object s. */
 int basic::ldegree(const ex & s) const
 {
        return 0;
 }
 
-/** Return coefficient of degree n in symbol s. */
+/** Return coefficient of degree n in object s. */
 ex basic::coeff(const ex & s, int n) const
 {
        return n==0 ? *this : _ex0();
 }
 
-/** Sort expression in terms of powers of some symbol.
- *  @param s symbol to sort in. */
-ex basic::collect(const ex & s) const
+/** Sort expression in terms of powers of some object(s).
+ *  @param s object(s) to sort in
+ *  @param distributed recursive or distributed form (only used when s is a list) */
+ex basic::collect(const ex & s, bool distributed) const
 {
        ex x;
-       for (int n=this->ldegree(s); n<=this->degree(s); ++n)
-               x += this->coeff(s,n)*power(s,n);
+       if (is_ex_of_type(s, lst)) {
+
+               // List of objects specified
+               if (s.nops() == 1)
+                       return collect(s.op(0));
+
+               else if (distributed) {
+
+                       // Get lower/upper degree of all symbols in list
+                       int 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 (int 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 (int i=0; i<num; i++) {
+                                       int cnt = si[i].cnt;
+                                       y *= power(si[i].sym, cnt);
+                               }
+                               x += y * si[num - 1].coeff;
+
+                               // Increment counters
+                               int 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 (int 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--;
+                               }
+                       }
+
+done:          delete[] si;
+
+               } else {
+
+                       // Recursive form
+                       x = *this;
+                       for (int n=s.nops()-1; n>=0; n--)
+                               x = x.collect(s[n]);
+               }
+
+       } else {
+
+               // Only one object specified
+               for (int n=this->ldegree(s); n<=this->degree(s); ++n)
+                       x += this->coeff(s,n)*power(s,n);
+       }
        
        // correct for lost fractional arguments and return
        return x + (*this - x).expand();
@@ -467,7 +542,7 @@ ex basic::subs(const ex & e) const
        for (unsigned i=0; i<e.nops(); i++) {
                ex r = e.op(i);
                if (!r.info(info_flags::relation_equal)) {
-                       throw(std::invalid_argument("basic::subs(ex): argument must be a list or equations"));
+                       throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
                }
                ls.append(r.op(0));
                lr.append(r.op(1));