[GiNaC-devel] [PATCH] collect: improve distributed collect of sparse multivariate polynomials.

Alexei Sheplyakov varg at theor.jinr.ru
Tue Nov 21 20:27:46 CET 2006


Hi,
as I promised here is a patch for basic::collect()

The old algorithm seems to be particularly inefficient when *this
is sparse multivariate polynomial. New algorithm handles such polynomials
much better (this is yet another "zero or infinity evaluation time"
issue). It should handle dense polynomials fairly well too. Another bonus
is that new code is simpler.

Best regards,
 Alexei

-- 
All science is either physics or stamp collecting.

---
 ginac/basic.cpp |   71 ++++++++++++++++++------------------------------------
 1 files changed, 24 insertions(+), 47 deletions(-)

diff --git a/ginac/basic.cpp b/ginac/basic.cpp
index f86483e..c7a33ca 100644
--- a/ginac/basic.cpp
+++ b/ginac/basic.cpp
@@ -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"
@@ -369,56 +370,32 @@ ex basic::collect(const ex & s, bool dis
 
 		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 {
 
-- 
1.4.3.3



-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 827 bytes
Desc: Digital signature
Url : http://www.cebix.net/pipermail/ginac-devel/attachments/20061121/b1b1b777/attachment.pgp


More information about the GiNaC-devel mailing list