November 2006 - GiNaC-devel

[PATCH] gcd: fix heuristics to chose main variable properly
by Alexei Sheplyakov 03 Dec '06

03 Dec '06

Hello! Heuristic code in gcd [very] often makes bad choice of the main variable. In particular, if one of polynomials (say, p2) does not contain the variable x which _is contained_ in another one (say, p1), than gcd(p1, p2) will use such a variable as a main, so subsequent calculation becomes unnecessary complicated. Fixing this substantially improves worst case performance (less than minute versus infinity). Best regards, Alexei -- All science is either physics or stamp collecting. --- ginac/normal.cpp | 46 +++++++++++++++++++++++++++++++++++++++++++--- 1 files changed, 43 insertions(+), 3 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index b2a4e8c..1cbf3c2 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1589,11 +1589,51 @@ factored_b: } } + if (is_exactly_a<numeric>(aex)) { + numeric bcont = bex.integer_content(); + numeric g = gcd(ex_to<numeric>(aex), bcont); + if (ca) + *ca = ex_to<numeric>(aex)/g; + if (cb) + *cb = bex/g; + return g; + } + + if (is_exactly_a<numeric>(bex)) { + numeric acont = aex.integer_content(); + numeric g = gcd(ex_to<numeric>(bex), acont); + if (ca) + *ca = aex/g; + if (cb) + *cb = ex_to<numeric>(bex)/g; + return g; + } + // Gather symbol statistics sym_desc_vec sym_stats; get_symbol_stats(a, b, sym_stats); - // The symbol with least degree is our main variable + // The symbol with least degree which is contained in both polynomials + // is our main variable + sym_desc_vec::iterator vari = sym_stats.begin(); + while ( (vari != sym_stats.end()) && + (((vari->ldeg_b == 0) && (vari->deg_b == 0)) || + ((vari->ldeg_a == 0) && (vari->deg_a == 0)))) + vari++; + + // No common symbols at all, just return 1: + if (vari == sym_stats.end()) { + // N.B: keep cofactors factored + if (ca) + *ca = a; + if (cb) + *cb = b; + return _ex1; + } + // move symbols which contained only in one of the polynomials + // to the end: + rotate(sym_stats.begin(), vari, sym_stats.end()); + sym_desc_vec::const_iterator var = sym_stats.begin(); const ex &x = var->sym; @@ -1607,14 +1647,14 @@ factored_b: } // Try to eliminate variables - if (var->deg_a == 0) { + if ((var->deg_a == 0) && (var->deg_b != 0)) { ex bex_u, bex_c, bex_p; bex.unitcontprim(x, bex_u, bex_c, bex_p); ex g = gcd(aex, bex_c, ca, cb, false); if (cb) *cb *= bex_u * bex_p; return g; - } else if (var->deg_b == 0) { + } else if ((var->deg_b == 0) && (var->deg_a != 0)) { ex aex_u, aex_c, aex_p; aex.unitcontprim(x, aex_u, aex_c, aex_p); ex g = gcd(aex_c, bex, ca, cb, false); -- 1.4.3.3

4 3

backported fix for reposition_dummy_indices to 1.3
by Chris Dams 22 Nov '06

22 Nov '06

Hello again, I backported the fixes to reposition_dummy_indices that were discussed on this list in October to 1.3. This was Alexeis patch, so thanks to Alexei! Best, Chris

1 0

Patches for collect_common_factors
by Chris Dams 22 Nov '06

22 Nov '06

Dear all, I applied the three patches that Alexei sent at the beginning of October for improving collect_common_factors to both 1.3 and head. Bye, Chris

1 0

Improvements for positive symbols.
by Chris Dams 22 Nov '06

22 Nov '06

Dear all, .info(info_flags::positive) did not yet return true for a postive symbol. I fixed this and also made atan2 handle pos. symbols better. Now atan2(0, x) with x a positive symbol returns 0 automatically. Best wishes, Chris

1 0

Re: [GiNaC-devel] [PATCH] gcd: fix heuristics to chose main variable properly
by Alexei Sheplyakov 21 Nov '06

21 Nov '06

[adding ginac-devel to CC] On Tue, Nov 21, 2006 at 09:28:13AM +0100, Bernard Parisse wrote: > BTW, I checked your example with the gcd code of my giac library, > it takes less than 3s. to normalize on a Intel(R) Pentium(R) M > processor 1.60GHz, how long does it take with your patched ginac version? [on a Pentium 4 3GHz box] $ g++ -O2 -march=pentium4 -finline-limit=1200 thepoly_another.cpp -o thepoly_another -lginac $ time ./thepoly_another >/dev/null 2>&1 real 0m7.846s user 0m7.836s sys 0m0.008s With unpatched GiNaC it takes way to long (after several days of calculation I get bored and killed the process). Best regards, Alexei -- All science is either physics or stamp collecting.

4 7

[PATCH] collect: improve distributed collect of sparse multivariate polynomials.
by Alexei Sheplyakov 21 Nov '06

21 Nov '06

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

1 0