[GiNaC-devel] [PATCH] {add, power}::eval(): convert terms into primitive polynomials when possible.

[Sheplyakov Alexei]

Hello!

On Thu, Jul 12, 2007 at 11:23:02PM +0200, Richard B. Kreckel wrote:
> 
> Alexei Sheplyakov wrote:
> >The following transformations are done automagically now:
> >(-x+y)*(x-y) -> -(x-y)^2
> >(3*x+y)/(6*x+2*y)^2 -> 1/4*(y+3*x)^(-1)
> >
> >Please note: exam_series13() from check/exam_pseries.cpp fails due to
> >this patch, because the expression 1/x*(-(1+x)/(1-x)) + (1+x)/x/(1-x)
> >gets evaluated to zero.
> 
> This patch appears to come with a slight performance hit. 

Indeed, mul::eval() has more work to be done.

advantages:
- subsequent operations, in particular gcd, are faster.
- no more ugly stuff like (-x+y)*(x-y).

disadvantages:
- subsequent operations, in particualr expand, may be slower.
- there is some overhead in the trivial case, e.g. (x*y*z).eval().

> On my machine, 
> 14x14 Vandermonde matrix determinants degrades from 19s to 32s,
> Lewis-Wester test O1 from 7.5s to 11s,

On my box results are 14.847s versus 16.949s 

> and Lewis-Wester test A from 0.055s to 0.22s.

Ouch. I do observe some slowdown (from 0.292s to 0.3s) too, but it is
not *that* bad.

> (I'm seriously surprised by the latter since it only 
> involves integers, but it's reproducible.)

a/b is actually mul(a, power(b, _ex_1)), and it has seq.size() >= 2.
My patch adds empty loop over 2 elements, so there should be some slowdown.

> Looking at the code it appears to me like mul::eval recurses now. 

mul::eval calls evalchildren, so it is recursive even without my patch.

> Wouldn't it be faster to collect the coefficients in one sweep? 

I've tried this (see the patch [1]), it has almost no effect.

> Well, I recognize that most terms will be trivial in the following
> recursions, but still.

My patch causes overhead even in the trivial case, e.g. when the
polynomial is already unit normal. Consider e.g. 

e = (x+1)*(2*x^2+x+1)*...*(100*x^100+99*x^99+...+1);

mul::eval will do 100 of useless checks (each of them is relatively
expansive on its own). Any ideas how to avoid this? E.g. what is 
the best way to check if polynomial is already primitive (unit normal)?

> Also, you write that "in some very unlucky event it can even loop 
> forever". I don't quite understand. 

The main variable(s) is (are) random (the first term of the add object).
During re-evaluation expression ordering might change, so the recursion
is not guaranteed to terminate.

> Do you have an example for this?

No, this is purely theoretical.

[1]

[PATCH] {mul,power}::eval(): remove extra recursion (divide by integer in place).

---
 ginac/mul.cpp    |   42 +++++++++++++++++++++++++++++++-----------
 ginac/power.cpp  |   16 +++++++++++++---
 ginac/compiler.h |   12 ++++++++++++
 3 files changed, 56 insertions(+), 14 deletions(-)

diff --git a/ginac/mul.cpp b/ginac/mul.cpp
index 4ceeb79..f9f7be5 100644
--- a/ginac/mul.cpp
+++ b/ginac/mul.cpp
@@ -34,6 +34,7 @@
 #include "lst.h"
 #include "archive.h"
 #include "utils.h"
+#include "compiler.h"
 
 namespace GiNaC {
 
@@ -440,14 +441,18 @@ ex mul::eval(int level) const
 		                ex_to<numeric>(addref.overall_coeff).
 		                mul_dyn(ex_to<numeric>(overall_coeff))))
 		      ->setflag(status_flags::dynallocated | status_flags::evaluated);
-  } else if (seq_size >= 2) {
+  } else if ((seq_size >= 2) && (! (flags & status_flags::expanded) )) {
     // Strip the content and the unit part from the each term. Thus
     // things like (-x+a)*(3*x-3*a) automagically turn into - 3*(x-a)^2
 
     epvector::const_iterator last = seq.end();
     epvector::const_iterator i = seq.begin();
+		epvector::const_iterator j = seq.begin();
+    std::auto_ptr<epvector> s(new epvector);
+		numeric oc = *_num1_p;
+		bool something_changed = false;
     while (i!=last) {
-      if (! (is_a<add>(i->rest) && i->coeff.is_equal(_ex1))) {
+      if (likely((! i->coeff.is_equal(_ex1)) || (! is_a<add>(i->rest)))) {
         // power::eval has such a rule, no need to handle powers here
         ++i;
         continue;
@@ -464,32 +469,47 @@ ex mul::eval(int level) const
       // very unlucky event it can even loop forever). Hopefully the main
       // variable will be the same for all terms in *this
       const bool unit_normal = lead_coeff.is_pos_integer();
-      if ((c == *_num1_p) && ((! canonicalizable) || unit_normal)) {
+      if (likely((c == *_num1_p) && ((! canonicalizable) || unit_normal))) {
         ++i;
         continue;
       }
 
-      std::auto_ptr<epvector> s(new epvector);
-      s->reserve(seq.size());
+			if (! something_changed) {
+				s->reserve(seq.size());
+				something_changed = true;
+			}
 
-      epvector::const_iterator j=seq.begin();
-      while (j!=i) {
+      while (j!=i && (j!=last)) {
         s->push_back(*j);
         ++j;
       }
 
       if (! unit_normal)
         c = c.mul(*_num_1_p);
-      const ex primitive = (i->rest)/c;
-      s->push_back(expair(primitive, _ex1));
+			
+			oc = oc.mul(c);
+
+      const add& addref = ex_to<add>(i->rest);
+      add* primitive = new add(addref);
+      primitive->setflag(status_flags::dynallocated);
+      primitive->clearflag(status_flags::hash_calculated);
+      primitive->overall_coeff = ex_to<numeric>(primitive->overall_coeff).div_dyn(c); 
+      for (epvector::iterator ai = primitive->seq.begin();
+                              ai != primitive->seq.end(); ++ai)
+        ai->coeff = ex_to<numeric>(ai->coeff).div_dyn(c);
+      
+      s->push_back(expair(*primitive, _ex1));
+
+      ++i;
       ++j;
-
+    }
+    if (something_changed) {
       while (j!=last) {
         s->push_back(*j);
         ++j;
       }
       return (new mul(s,
-                      ex_to<numeric>(overall_coeff).mul_dyn(c))
+                      ex_to<numeric>(overall_coeff).mul_dyn(oc))
                       )->setflag(status_flags::dynallocated);
     }
   }
diff --git a/ginac/power.cpp b/ginac/power.cpp
index 238e12e..0103398 100644
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -478,9 +478,19 @@ ex power::eval(int level) const
       const bool canonicalizable = lead_coeff.is_integer();
       const bool unit_normal = lead_coeff.is_pos_integer();
 
-      if (icont != *_num1_p)
-        return (new mul(power(ebasis/icont, *num_exponent),
-                        power(icont, *num_exponent)))->setflag(status_flags::dynallocated);
+      if (icont != *_num1_p) {
+        const add& addref = ex_to<add>(ebasis);
+        add* addp = new add(addref);
+        addp->setflag(status_flags::dynallocated);
+        addp->clearflag(status_flags::hash_calculated);
+        addp->overall_coeff = ex_to<numeric>(addp->overall_coeff).div_dyn(icont);
+        for (epvector::iterator i = addp->seq.begin(); i != addp->seq.end(); ++i)
+          i->coeff = ex_to<numeric>(i->coeff).div_dyn(icont);
+
+        return (new mul(power(*addp, *num_exponent),
+              icont.power(*num_exponent)))->setflag(status_flags::dynallocated);
+      }
+
 
       if (canonicalizable && (! unit_normal)) {
         if (num_exponent->is_even())
diff --git a/ginac/compiler.h b/ginac/compiler.h
new file mode 100644
index 0000000..66f6e92
--- /dev/null
+++ b/ginac/compiler.h
@@ -0,0 +1,12 @@
+#ifndef GINAC_COMPILER_DEP_HH
+#define GINAC_COMPILER_DEP_HH
+
+#ifdef __GNUC__
+#define unlikely(cond) __builtin_expect((cond), 0)
+#define likely(cond) __builtin_expect((cond), 1)
+#else
+#define unlikely(cond) (cond)
+#define likely(cond) (cond)
+#endif
+
+#endif /* GINAC_COMPILER_DEP_HH */
-- 
1.4.4.4


Best regards,
 Alexei

-- 
All science is either physics or stamp collecting.

-------------- 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/20070806/14ec0a4c/attachment.pgp