August 2008 - GiNaC-devel

factorization and stuff
by Jens Vollinga 10 Sep '08

10 Sep '08

Hi, I am re-sending this email because it looks like the server problems this morning have swallowed the original posting. Hi, I recently checked in code for polynomial factorization. The multivariate case and maybe algebraic extensions are in the making. Some notes on the code: - The code is absolutely not optimized. Don't look at the code, please! - Alternative methods are missing yet. It is doing the factorization with only one combination of methods and there are a lot of (well, infinitely many) use cases where the code performs amazingly poor. - The code is debugged and it should give the right answer for a proper expressions. But handling of non-proper expressions is still missing. Proper ~ Polynomial with symbols. I want to have the code feature-complete first (because I need it for some other stuff), only then will I try to optimize the code and build in additional algorithms. But if you find a bug in the factorization now please tell me! I also recently updated the todo list on the GiNaC website. Well, factorization will be done, but for the other two new points I have done nothing yet and there is always the chance that I will do nothing in the future. If you have interest in working on them, just do so (but maybe tell me in advance). In case you are planning to add other new features, I could announce them on the todo list as well. @Alexei: if you want to repair that GCD and is_polynomial business go ahead! That would be great! Regards, Jens

2 12

[SCM] GiNaC -- a C++ library for symbolic computations branch, master, updated. release_1-4-0-77-g55af760
by git＠ginac.de 27 Aug '08

27 Aug '08

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "GiNaC -- a C++ library for symbolic computations". The branch, master has been updated via 55af76071727bd6e2fd540ad58eac26dd961f9c9 (commit) via a79a813e7249f793859d1d3b443d1931dbab94b6 (commit) via 6ff92476a4dcd32f9a0c6f59c95c74812ea86fef (commit) via 55e50ea26252dff7432bdce8b010f9fbd058e907 (commit) via 1d09676118944532e6100c93383d1659b1253a60 (commit) via adb1dbd383ae6e5a999b5f8cba72a5c2bfd50c11 (commit) via c77689e7ac8d8f4dbca0f337b6e9acf2419010ff (commit) via 77b6a0304a48d6a306deeda18177680f16025b33 (commit) via 7d7131d3af3de5425b7fe80b1f587740294371bc (commit) via b65fcd7481a401ec23c284c91f6f4e883e967676 (commit) via 8474043d373a19e04008f476fa9b77e45734b604 (commit) via ff09c4f8103f53fe3b7a51eb3c33eff2e5a243f0 (commit) via b4b302fe5d54720b58a23568a4270e04ee1ca216 (commit) from 1222eac51cee964961d2aad889dc4ceccb144a36 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 55af76071727bd6e2fd540ad58eac26dd961f9c9 Author: Jens Vollinga <jensv(a)nikhef.nl> Date: Wed Aug 27 17:25:16 2008 +0200 Added internal code for multivariate factorization. commit a79a813e7249f793859d1d3b443d1931dbab94b6 Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Mon Aug 25 16:57:38 2008 +0400 gcd_pf_pow_pow: deobfuscate a little bit (no functional changes). Use if (foo) return bar(); return baz(); instead of if (foo) { return bar(); } else { return baz(); } This makes the code a little bit more readable. commit 6ff92476a4dcd32f9a0c6f59c95c74812ea86fef Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Mon Aug 25 16:57:02 2008 +0400 gcd_pf_pow: get rid of duplicate code. This function (which helps gcd() handle partially factored expressions) contained two copies of the same code. Remove one redundant copy. commit 55e50ea26252dff7432bdce8b010f9fbd058e907 Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Mon Aug 25 16:56:04 2008 +0400 introduce gcd_pf_pow_pow: gcd helper to handle partially factored expressions. gcd_pf_pow_pow handles the case when both arguments of gcd() are powers. N.B. the actual code moved from gcd_pf_pow, no functional changes. commit 1d09676118944532e6100c93383d1659b1253a60 Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Mon Aug 25 16:55:42 2008 +0400 gcd_pf_{pow, mul}: don't check if the arguments are polynomials. These functions gets called only from gcd(), which does this check on its own. commit adb1dbd383ae6e5a999b5f8cba72a5c2bfd50c11 Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Mon Aug 25 16:55:13 2008 +0400 gcd_pf_mul: get rid of duplicate code. This function (which helps gcd() handle partially factored expressions) contained two copies of the same code. Remove one redundant copy. commit c77689e7ac8d8f4dbca0f337b6e9acf2419010ff Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Mon Aug 25 16:54:46 2008 +0400 gcd(): allow user to override (some of) heuristics. GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 5: * gcd(): don't use heuristic GCD algorithm if gcd_options::no_heur_gcd flag is set. * gcd(): don't handle partially factored expressions in a special way if gcd_options::no_part_factored flag is set. commit 77b6a0304a48d6a306deeda18177680f16025b33 Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Mon Aug 25 16:54:10 2008 +0400 refactor gcd() a little bit (no functional changes). GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 4: refactor gcd() a little bit, so subsequent patch(es) won't be so big and ugly. commit 7d7131d3af3de5425b7fe80b1f587740294371bc Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Mon Aug 25 16:53:47 2008 +0400 introduce gcd_pf_mul: gcd helper to handle partially factored expressions. GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 3: Move the code handling products from gcd() into a separate function. This is *really* only code move, everything else should be considered a bug. commit b65fcd7481a401ec23c284c91f6f4e883e967676 Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Mon Aug 25 16:53:07 2008 +0400 introduce gcd_pf_pow: gcd helper to handle partially factored expressions. GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 2: Move the code handling powers from gcd() into a separate function. This is *really* only code move, everything else should be considered a bug. commit 8474043d373a19e04008f476fa9b77e45734b604 Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Mon Aug 25 16:52:45 2008 +0400 gcd: allow user to override (some of) heuristics. GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 1: * add new (optional) argument to gcd() to control its behaviour. * introduce gcd_options structure. N.B. No actual code changes so far, the actual handling of newly introduced options is the subject of further patches. commit ff09c4f8103f53fe3b7a51eb3c33eff2e5a243f0 Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Sun Aug 24 15:25:24 2008 +0400 inifscn_nstdsums: make functions handling Li/G transformations reentrant. Explicitly pass the dummy symbols instead of using a global variable. As a side effect, the code is more clear now (that's a bit subjective, but anyway). commit b4b302fe5d54720b58a23568a4270e04ee1ca216 Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Fri Aug 22 02:42:17 2008 +0400 Bail out if both autogen and autogen'erated file(s) are missing. This makes the error message(s) more helpful. ----------------------------------------------------------------------- Summary of changes: ginac/Makefile.am | 1 + ginac/factor.cpp | 657 +++++++++++++++++++++++++++++++++++++++++++- ginac/inifcns_nstdsums.cpp | 128 +++++---- ginac/normal.cpp | 306 ++++++++++----------- ginac/normal.h | 25 ++- 5 files changed, 897 insertions(+), 220 deletions(-) hooks/post-receive -- GiNaC -- a C++ library for symbolic computations

1 0

[PATCH 10/10] gcd_pf_pow_pow: deobfuscate a little bit (no functional changes).
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

Use if (foo) return bar(); return baz(); instead of if (foo) { return bar(); } else { return baz(); } This makes the code a little bit more readable. --- ginac/normal.cpp | 40 ++++++++++++++++++++-------------------- 1 files changed, 20 insertions(+), 20 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 6392e3f..09773d3 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1647,8 +1647,9 @@ static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb) const ex& exp_a = a.op(1); ex pb = b.op(0); const ex& exp_b = b.op(1); + + // a = p^n, b = p^m, gcd = p^min(n, m) if (p.is_equal(pb)) { - // a = p^n, b = p^m, gcd = p^min(n, m) if (exp_a < exp_b) { if (ca) *ca = _ex1; @@ -1662,31 +1663,30 @@ static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb) *cb = _ex1; return power(p, exp_b); } - } else { - ex p_co, pb_co; - ex p_gcd = gcd(p, pb, &p_co, &pb_co, false); - if (p_gcd.is_equal(_ex1)) { - // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==> - // gcd(a,b) = 1 + } + + ex p_co, pb_co; + ex p_gcd = gcd(p, pb, &p_co, &pb_co, false); + // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==> gcd(a,b) = 1 + if (p_gcd.is_equal(_ex1)) { if (ca) *ca = a; if (cb) *cb = b; return _ex1; // XXX: do I need to check for p_gcd = -1? - } else { - // there are common factors: - // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==> - // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m - if (exp_a < exp_b) { - return power(p_gcd, exp_a)* - gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false); - } else { - return power(p_gcd, exp_b)* - gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false); - } - } // p_gcd.is_equal(_ex1) - } // p.is_equal(pb) + } + + // there are common factors: + // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==> + // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m + if (exp_a < exp_b) { + ex pg = gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false); + return power(p_gcd, exp_a)*pg; + } else { + ex pg = gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false); + return power(p_gcd, exp_b)*pg; + } } static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb) -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH 09/10] gcd_pf_pow: get rid of duplicate code.
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

This function (which helps gcd() handle partially factored expressions) contained two copies of the same code. Remove one redundant copy. --- ginac/normal.cpp | 80 ++++++++++++++++++----------------------------------- 1 files changed, 27 insertions(+), 53 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index d5319b9..6392e3f 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1691,65 +1691,39 @@ static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb) static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb) { - if (is_exactly_a<power>(a)) { - ex p = a.op(0); - const ex& exp_a = a.op(1); - if (is_exactly_a<power>(b)) - return gcd_pf_pow_pow(a, b, ca, cb); - else { - if (p.is_equal(b)) { - // a = p^n, b = p, gcd = p - if (ca) - *ca = power(p, a.op(1) - 1); - if (cb) - *cb = _ex1; - return p; - } + if (is_exactly_a<power>(a) && is_exactly_a<power>(b)) + return gcd_pf_pow_pow(a, b, ca, cb); - ex p_co, bpart_co; - ex p_gcd = gcd(p, b, &p_co, &bpart_co, false); + if (is_exactly_a<power>(b) && (! is_exactly_a<power>(a))) + return gcd_pf_pow(b, a, cb, ca); - if (p_gcd.is_equal(_ex1)) { - // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1 - if (ca) - *ca = a; - if (cb) - *cb = b; - return _ex1; - } else { - // a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x)) - return p_gcd*gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false); - } - } // is_exactly_a<power>(b) + GINAC_ASSERT(is_exactly_a<power>(a)); - } else if (is_exactly_a<power>(b)) { - ex p = b.op(0); - if (p.is_equal(a)) { - // a = p, b = p^n, gcd = p - if (ca) - *ca = _ex1; - if (cb) - *cb = power(p, b.op(1) - 1); - return p; - } + ex p = a.op(0); + const ex& exp_a = a.op(1); + if (p.is_equal(b)) { + // a = p^n, b = p, gcd = p + if (ca) + *ca = power(p, a.op(1) - 1); + if (cb) + *cb = _ex1; + return p; + } - ex p_co, apart_co; - const ex& exp_b(b.op(1)); - ex p_gcd = gcd(a, p, &apart_co, &p_co, false); - if (p_gcd.is_equal(_ex1)) { - // b=p(x)^n, gcd(a, p) = 1 ==> gcd(a, b) == 1 - if (ca) - *ca = a; - if (cb) - *cb = b; - return _ex1; - } else { - // there are common factors: - // a(x) = g(x) A(x), b(x) = g(x)^n B(x)^n ==> gcd = g(x) gcd(g(x)^(n-1) A(x)^n, B(x)) + ex p_co, bpart_co; + ex p_gcd = gcd(p, b, &p_co, &bpart_co, false); - return p_gcd*gcd(apart_co, power(p_gcd, exp_b-1)*power(p_co, exp_b), ca, cb, false); - } // p_gcd.is_equal(_ex1) + // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1 + if (p_gcd.is_equal(_ex1)) { + if (ca) + *ca = a; + if (cb) + *cb = b; + return _ex1; } + // a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x)) + ex rg = gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false); + return p_gcd*rg; } static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb) -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH 08/10] introduce gcd_pf_pow_pow: gcd helper to handle partially factored expressions.
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

gcd_pf_pow_pow handles the case when both arguments of gcd() are powers. N.B. the actual code moved from gcd_pf_pow, no functional changes. --- ginac/normal.cpp | 98 +++++++++++++++++++++++++++++------------------------- 1 files changed, 53 insertions(+), 45 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index e97a7ce..d5319b9 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1639,56 +1639,64 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio return g; } +// gcd helper to handle partially factored polynomials (to avoid expanding +// large expressions). Both arguments should be powers. +static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb) +{ + ex p = a.op(0); + const ex& exp_a = a.op(1); + ex pb = b.op(0); + const ex& exp_b = b.op(1); + if (p.is_equal(pb)) { + // a = p^n, b = p^m, gcd = p^min(n, m) + if (exp_a < exp_b) { + if (ca) + *ca = _ex1; + if (cb) + *cb = power(p, exp_b - exp_a); + return power(p, exp_a); + } else { + if (ca) + *ca = power(p, exp_a - exp_b); + if (cb) + *cb = _ex1; + return power(p, exp_b); + } + } else { + ex p_co, pb_co; + ex p_gcd = gcd(p, pb, &p_co, &pb_co, false); + if (p_gcd.is_equal(_ex1)) { + // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==> + // gcd(a,b) = 1 + if (ca) + *ca = a; + if (cb) + *cb = b; + return _ex1; + // XXX: do I need to check for p_gcd = -1? + } else { + // there are common factors: + // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==> + // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m + if (exp_a < exp_b) { + return power(p_gcd, exp_a)* + gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false); + } else { + return power(p_gcd, exp_b)* + gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false); + } + } // p_gcd.is_equal(_ex1) + } // p.is_equal(pb) +} + static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb) { if (is_exactly_a<power>(a)) { ex p = a.op(0); const ex& exp_a = a.op(1); - if (is_exactly_a<power>(b)) { - ex pb = b.op(0); - const ex& exp_b = b.op(1); - if (p.is_equal(pb)) { - // a = p^n, b = p^m, gcd = p^min(n, m) - if (exp_a < exp_b) { - if (ca) - *ca = _ex1; - if (cb) - *cb = power(p, exp_b - exp_a); - return power(p, exp_a); - } else { - if (ca) - *ca = power(p, exp_a - exp_b); - if (cb) - *cb = _ex1; - return power(p, exp_b); - } - } else { - ex p_co, pb_co; - ex p_gcd = gcd(p, pb, &p_co, &pb_co, false); - if (p_gcd.is_equal(_ex1)) { - // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==> - // gcd(a,b) = 1 - if (ca) - *ca = a; - if (cb) - *cb = b; - return _ex1; - // XXX: do I need to check for p_gcd = -1? - } else { - // there are common factors: - // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==> - // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m - if (exp_a < exp_b) { - return power(p_gcd, exp_a)* - gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false); - } else { - return power(p_gcd, exp_b)* - gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false); - } - } // p_gcd.is_equal(_ex1) - } // p.is_equal(pb) - - } else { + if (is_exactly_a<power>(b)) + return gcd_pf_pow_pow(a, b, ca, cb); + else { if (p.is_equal(b)) { // a = p^n, b = p, gcd = p if (ca) -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH 07/10] gcd_pf_{pow, mul}: don't check if the arguments are polynomials.
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

These functions gets called only from gcd(), which does this check on its own. --- ginac/normal.cpp | 20 ++++++++++---------- 1 files changed, 10 insertions(+), 10 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 2bb3a43..e97a7ce 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1417,11 +1417,11 @@ static bool heur_gcd(ex& res, const ex& a, const ex& b, ex *ca, ex *cb, // gcd helper to handle partially factored polynomials (to avoid expanding // large expressions). At least one of the arguments should be a power. -static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args); +static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb); // gcd helper to handle partially factored polynomials (to avoid expanding // large expressions). At least one of the arguments should be a product. -static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args); +static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb); /** Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X) * and b(X) in Z[X]. Optionally also compute the cofactors of a and b, @@ -1467,10 +1467,10 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio // Partially factored cases (to avoid expanding large expressions) if (!(options & gcd_options::no_part_factored)) { if (is_exactly_a<mul>(a) || is_exactly_a<mul>(b)) - return gcd_pf_mul(a, b, ca, cb, check_args); + return gcd_pf_mul(a, b, ca, cb); #if FAST_COMPARE if (is_exactly_a<power>(a) || is_exactly_a<power>(b)) - return gcd_pf_pow(a, b, ca, cb, check_args); + return gcd_pf_pow(a, b, ca, cb); #endif } @@ -1639,7 +1639,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio return g; } -static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) +static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb) { if (is_exactly_a<power>(a)) { ex p = a.op(0); @@ -1664,7 +1664,7 @@ static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) } } else { ex p_co, pb_co; - ex p_gcd = gcd(p, pb, &p_co, &pb_co, check_args); + ex p_gcd = gcd(p, pb, &p_co, &pb_co, false); if (p_gcd.is_equal(_ex1)) { // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==> // gcd(a,b) = 1 @@ -1744,14 +1744,14 @@ static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) } } -static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) +static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb) { if (is_exactly_a<mul>(a) && is_exactly_a<mul>(b) && (b.nops() > a.nops())) - return gcd_pf_mul(b, a, cb, ca, check_args); + return gcd_pf_mul(b, a, cb, ca); if (is_exactly_a<mul>(b) && (!is_exactly_a<mul>(a))) - return gcd_pf_mul(b, a, cb, ca, check_args); + return gcd_pf_mul(b, a, cb, ca); GINAC_ASSERT(is_exactly_a<mul>(a)); size_t num = a.nops(); @@ -1760,7 +1760,7 @@ static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) ex part_b = b; for (size_t i=0; i<num; i++) { ex part_ca, part_cb; - g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, check_args)); + g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, false)); acc_ca.push_back(part_ca); part_b = part_cb; } -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH 06/10] gcd_pf_mul: get rid of duplicate code.
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

This function (which helps gcd() handle partially factored expressions) contained two copies of the same code. Remove one redundant copy. --- ginac/normal.cpp | 62 ++++++++++++++++++++---------------------------------- 1 files changed, 23 insertions(+), 39 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 5d044fc..2bb3a43 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1746,45 +1746,29 @@ static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) { - if (is_exactly_a<mul>(a)) { - if (is_exactly_a<mul>(b) && b.nops() > a.nops()) - goto factored_b; -factored_a: - size_t num = a.nops(); - exvector g; g.reserve(num); - exvector acc_ca; acc_ca.reserve(num); - ex part_b = b; - for (size_t i=0; i<num; i++) { - ex part_ca, part_cb; - g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, check_args)); - acc_ca.push_back(part_ca); - part_b = part_cb; - } - if (ca) - *ca = (new mul(acc_ca))->setflag(status_flags::dynallocated); - if (cb) - *cb = part_b; - return (new mul(g))->setflag(status_flags::dynallocated); - } else if (is_exactly_a<mul>(b)) { - if (is_exactly_a<mul>(a) && a.nops() > b.nops()) - goto factored_a; -factored_b: - size_t num = b.nops(); - exvector g; g.reserve(num); - exvector acc_cb; acc_cb.reserve(num); - ex part_a = a; - for (size_t i=0; i<num; i++) { - ex part_ca, part_cb; - g.push_back(gcd(part_a, b.op(i), &part_ca, &part_cb, check_args)); - acc_cb.push_back(part_cb); - part_a = part_ca; - } - if (ca) - *ca = part_a; - if (cb) - *cb = (new mul(acc_cb))->setflag(status_flags::dynallocated); - return (new mul(g))->setflag(status_flags::dynallocated); - } + if (is_exactly_a<mul>(a) && is_exactly_a<mul>(b) + && (b.nops() > a.nops())) + return gcd_pf_mul(b, a, cb, ca, check_args); + + if (is_exactly_a<mul>(b) && (!is_exactly_a<mul>(a))) + return gcd_pf_mul(b, a, cb, ca, check_args); + + GINAC_ASSERT(is_exactly_a<mul>(a)); + size_t num = a.nops(); + exvector g; g.reserve(num); + exvector acc_ca; acc_ca.reserve(num); + ex part_b = b; + for (size_t i=0; i<num; i++) { + ex part_ca, part_cb; + g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, check_args)); + acc_ca.push_back(part_ca); + part_b = part_cb; + } + if (ca) + *ca = (new mul(acc_ca))->setflag(status_flags::dynallocated); + if (cb) + *cb = part_b; + return (new mul(g))->setflag(status_flags::dynallocated); } /** Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X]. -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH 05/10] gcd(): allow user to override (some of) heuristics.
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 5: * gcd(): don't use heuristic GCD algorithm if gcd_options::no_heur_gcd flag is set. * gcd(): don't handle partially factored expressions in a special way if gcd_options::no_part_factored flag is set. --- ginac/normal.cpp | 40 ++++++++++++++++++++++------------------ 1 files changed, 22 insertions(+), 18 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 0227f4e..5d044fc 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1465,12 +1465,14 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio } // Partially factored cases (to avoid expanding large expressions) - if (is_exactly_a<mul>(a) || is_exactly_a<mul>(b)) - return gcd_pf_mul(a, b, ca, cb, check_args); + if (!(options & gcd_options::no_part_factored)) { + if (is_exactly_a<mul>(a) || is_exactly_a<mul>(b)) + return gcd_pf_mul(a, b, ca, cb, check_args); #if FAST_COMPARE - if (is_exactly_a<power>(a) || is_exactly_a<power>(b)) - return gcd_pf_pow(a, b, ca, cb, check_args); + if (is_exactly_a<power>(a) || is_exactly_a<power>(b)) + return gcd_pf_pow(a, b, ca, cb, check_args); #endif + } // Some trivial cases ex aex = a.expand(), bex = b.expand(); @@ -1601,23 +1603,25 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio // Try heuristic algorithm first, fall back to PRS if that failed ex g; - bool found = heur_gcd(g, aex, bex, ca, cb, var); - if (found) { - // heur_gcd have already computed cofactors... - if (g.is_equal(_ex1)) { - // ... but we want to keep them factored if possible. - if (ca) - *ca = a; - if (cb) - *cb = b; + if (!(options & gcd_options::no_heur_gcd)) { + bool found = heur_gcd(g, aex, bex, ca, cb, var); + if (found) { + // heur_gcd have already computed cofactors... + if (g.is_equal(_ex1)) { + // ... but we want to keep them factored if possible. + if (ca) + *ca = a; + if (cb) + *cb = b; + } + return g; } - return g; - } #if STATISTICS - else { - heur_gcd_failed++; - } + else { + heur_gcd_failed++; + } #endif + } g = sr_gcd(aex, bex, var); if (g.is_equal(_ex1)) { -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH 04/10] refactor gcd() a little bit (no functional changes).
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 4: refactor gcd() a little bit, so subsequent patch(es) won't be so big and ugly. --- ginac/normal.cpp | 41 ++++++++++++++++++++++------------------- 1 files changed, 22 insertions(+), 19 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 0cb9100..0227f4e 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1602,33 +1602,36 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio // Try heuristic algorithm first, fall back to PRS if that failed ex g; bool found = heur_gcd(g, aex, bex, ca, cb, var); - if (!found) { -#if STATISTICS - heur_gcd_failed++; -#endif - g = sr_gcd(aex, bex, var); - if (g.is_equal(_ex1)) { - // Keep cofactors factored if possible - if (ca) - *ca = a; - if (cb) - *cb = b; - } else { - if (ca) - divide(aex, g, *ca, false); - if (cb) - divide(bex, g, *cb, false); - } - } else { + if (found) { + // heur_gcd have already computed cofactors... if (g.is_equal(_ex1)) { - // Keep cofactors factored if possible + // ... but we want to keep them factored if possible. if (ca) *ca = a; if (cb) *cb = b; } + return g; } +#if STATISTICS + else { + heur_gcd_failed++; + } +#endif + g = sr_gcd(aex, bex, var); + if (g.is_equal(_ex1)) { + // Keep cofactors factored if possible + if (ca) + *ca = a; + if (cb) + *cb = b; + } else { + if (ca) + divide(aex, g, *ca, false); + if (cb) + divide(bex, g, *cb, false); + } return g; } -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH 03/10] introduce gcd_pf_mul: gcd helper to handle partially factored expressions.
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 3: Move the code handling products from gcd() into a separate function. This is *really* only code move, everything else should be considered a bug. --- ginac/normal.cpp | 89 +++++++++++++++++++++++++++++------------------------ 1 files changed, 49 insertions(+), 40 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 0af5aad..0cb9100 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1419,6 +1419,10 @@ static bool heur_gcd(ex& res, const ex& a, const ex& b, ex *ca, ex *cb, // large expressions). At least one of the arguments should be a power. static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args); +// gcd helper to handle partially factored polynomials (to avoid expanding +// large expressions). At least one of the arguments should be a product. +static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args); + /** Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X) * and b(X) in Z[X]. Optionally also compute the cofactors of a and b, * defined by a = ca * gcd(a, b) and b = cb * gcd(a, b). @@ -1461,46 +1465,8 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio } // Partially factored cases (to avoid expanding large expressions) - if (is_exactly_a<mul>(a)) { - if (is_exactly_a<mul>(b) && b.nops() > a.nops()) - goto factored_b; -factored_a: - size_t num = a.nops(); - exvector g; g.reserve(num); - exvector acc_ca; acc_ca.reserve(num); - ex part_b = b; - for (size_t i=0; i<num; i++) { - ex part_ca, part_cb; - g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, check_args)); - acc_ca.push_back(part_ca); - part_b = part_cb; - } - if (ca) - *ca = (new mul(acc_ca))->setflag(status_flags::dynallocated); - if (cb) - *cb = part_b; - return (new mul(g))->setflag(status_flags::dynallocated); - } else if (is_exactly_a<mul>(b)) { - if (is_exactly_a<mul>(a) && a.nops() > b.nops()) - goto factored_a; -factored_b: - size_t num = b.nops(); - exvector g; g.reserve(num); - exvector acc_cb; acc_cb.reserve(num); - ex part_a = a; - for (size_t i=0; i<num; i++) { - ex part_ca, part_cb; - g.push_back(gcd(part_a, b.op(i), &part_ca, &part_cb, check_args)); - acc_cb.push_back(part_cb); - part_a = part_ca; - } - if (ca) - *ca = part_a; - if (cb) - *cb = (new mul(acc_cb))->setflag(status_flags::dynallocated); - return (new mul(g))->setflag(status_flags::dynallocated); - } - + if (is_exactly_a<mul>(a) || is_exactly_a<mul>(b)) + return gcd_pf_mul(a, b, ca, cb, check_args); #if FAST_COMPARE if (is_exactly_a<power>(a) || is_exactly_a<power>(b)) return gcd_pf_pow(a, b, ca, cb, check_args); @@ -1771,6 +1737,49 @@ static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) } } +static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) +{ + if (is_exactly_a<mul>(a)) { + if (is_exactly_a<mul>(b) && b.nops() > a.nops()) + goto factored_b; +factored_a: + size_t num = a.nops(); + exvector g; g.reserve(num); + exvector acc_ca; acc_ca.reserve(num); + ex part_b = b; + for (size_t i=0; i<num; i++) { + ex part_ca, part_cb; + g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, check_args)); + acc_ca.push_back(part_ca); + part_b = part_cb; + } + if (ca) + *ca = (new mul(acc_ca))->setflag(status_flags::dynallocated); + if (cb) + *cb = part_b; + return (new mul(g))->setflag(status_flags::dynallocated); + } else if (is_exactly_a<mul>(b)) { + if (is_exactly_a<mul>(a) && a.nops() > b.nops()) + goto factored_a; +factored_b: + size_t num = b.nops(); + exvector g; g.reserve(num); + exvector acc_cb; acc_cb.reserve(num); + ex part_a = a; + for (size_t i=0; i<num; i++) { + ex part_ca, part_cb; + g.push_back(gcd(part_a, b.op(i), &part_ca, &part_cb, check_args)); + acc_cb.push_back(part_cb); + part_a = part_ca; + } + if (ca) + *ca = part_a; + if (cb) + *cb = (new mul(acc_cb))->setflag(status_flags::dynallocated); + return (new mul(g))->setflag(status_flags::dynallocated); + } +} + /** Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X]. * * @param a first multivariate polynomial -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0