[GiNaC-devel] [PATCH 1/3] divide: try to avoid expanding partially factored expressions.
Alexei Sheplyakov
varg at pc7135.jinr.ru
Fri Oct 6 13:38:37 CEST 2006
Dear all,
this ginsh code snippet:
collect_common_factors((x*y*a+x*y*b)^(3) + (x*z + x*y)^(2))
gives
x^2*(x*(3*a^2*y^3*b+3*a*y^3*b^2+y^3*b^3+a^3*y^3)+2*z*y+z^2+y^2)
Thus instead of collecting common factors from the terms of sums
(as manual says) collect_common_factors [sometimes] expands term[s]!
---
ginac/normal.cpp | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
1 files changed, 66 insertions(+), 0 deletions(-)
diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 7b90030..49f0bd7 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -614,6 +614,72 @@ #endif
if (!get_first_symbol(a, x) && !get_first_symbol(b, x))
throw(std::invalid_argument("invalid expression in divide()"));
+ // Try to avoid expanding partially factored expressions.
+ if (is_exactly_a<mul>(b)) {
+ // Divide sequentially by each term
+ ex rem_new, rem_old = a;
+ for (size_t i=0; i < b.nops(); i++) {
+ if (! divide(rem_old, b.op(i), rem_new, false))
+ return false;
+ rem_old = rem_new;
+ }
+ q = rem_new;
+ return true;
+ } else if (is_exactly_a<power>(b)) {
+ const ex& bb(b.op(0));
+ int exp_b = ex_to<numeric>(b.op(1)).to_int();
+ ex rem_new, rem_old = a;
+ for (int i=exp_b; i>0; i--) {
+ if (! divide(rem_old, bb, rem_new, false))
+ return false;
+ rem_old = rem_new;
+ }
+ q = rem_new;
+ return true;
+ }
+
+ if (is_exactly_a<mul>(a)) {
+ // Divide sequentially each term. If some term in a is divisible
+ // by b we are done... and if not, we can't really say anything.
+ size_t i;
+ ex rem_i;
+ bool divisible_p = false;
+ for (i=0; i < a.nops(); ++i) {
+ if (divide(a.op(i), b, rem_i, false)) {
+ divisible_p = true;
+ break;
+ }
+ }
+ if (divisible_p) {
+ exvector resv;
+ resv.reserve(a.nops());
+ for (size_t j=0; j < a.nops(); j++) {
+ if (j==i)
+ resv.push_back(rem_i);
+ else
+ resv.push_back(a.op(j));
+ }
+ q = (new mul(resv))->setflag(status_flags::dynallocated);
+ return true;
+ }
+ } else if (is_exactly_a<power>(a)) {
+ // The base itself might be divisible by b, in that case we don't
+ // need to expand a
+ const ex& ab(a.op(0));
+ int a_exp = ex_to<numeric>(a.op(1)).to_int();
+ ex rem_i;
+ if (divide(ab, b, rem_i, false)) {
+ q = rem_i*power(ab, a_exp - 1);
+ return true;
+ }
+ for (int i=2; i < a_exp; i++) {
+ if (divide(power(ab, i), b, rem_i, false)) {
+ q = rem_i*power(ab, a_exp - i);
+ return true;
+ }
+ } // ... so we *really* need to expand expression.
+ }
+
// Polynomial long division (recursive)
ex r = a.expand();
if (r.is_zero()) {
--
All science is either physics or stamp collecting.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 827 bytes
Desc: Digital signature
URL: <http://www.ginac.de/pipermail/ginac-devel/attachments/20061006/bb94913f/attachment.sig>
More information about the GiNaC-devel
mailing list