- GiNaC-devel - ginac.de

New bug in symmetrization?
by Jens Vollinga 12 Dec '06

12 Dec '06

Hi, 'make check' gives me a failure for the paranoia exam number 17. That exam is concerned with the symmetrization of indices. Can anybody reproduce that finding? It didn't fail before I did an cvs update (my local version was rather old), so either the bug must be really new or it is one of the bugs that appear only when the ordering of expressions is disadvantageous due to the memory layout. Regards, Jens

3 4

Numerical evaluation for tgamma, lgamma and beta.
by Chris Dams 12 Dec '06

12 Dec '06

Dear all, I added numerical evaluation for the functions tgamma, lgamma and beta to the CVS. The method employed is the Lanczos approximation as described on Wikepedia: http://en.wikipedia.org/wiki/Lanczos_approximation . (I also corrected some small errors on that page.) This method uses a fixed precision and tables are used to store coefficients. Currently the tables available are for maximally 20, 50, 100 and 200 digits. If Digits exceeds 200, no numerical evaluation is available. Best wishes, Chris

2 1

Bugfix for numerical integration.
by Chris Dams 11 Dec '06

11 Dec '06

Dear all, I discovered that the numerical integration routine that I wrote does not handle complex functions. The inequalities in the error checking are the problem. Therefore, I added some abs-function to the code. BTW, maybe it would be a good idea to include the numerical evalution for the gamma function that I comitted last friday in 1.3 as well? Maybe it could also go into release today? Bye, Chris

2 2

Upcoming release
by Jens Vollinga 11 Dec '06

11 Dec '06

Dear all, due to a newly discovered bug in the numerical evaluation of multiple polylogs I feel the strong urge to do a new GiNaC release very soon. Are there any unresolved issues or pending submissions that call for a delay of this release beyond the oncoming Monday? Regards, Jens

3 2

[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