[Axiom-developer] Algebra speedups

[Waldek Hebisch]

Below you will find a patch (applied to wh-sandbox 603) which
significantly reduces time to run the testsuite (on my machine
from 20 minutes down to 12 minutes).

Some remarks:
- the patch somewhat reduces abstraction (for example uses Lisp
  EQ to quickly confirm equality).  IMHO this is necessary
  (and in case of this patch moderate) price for speed.
- surprisingly, manipulations on operators and kernels take a lot
  of time, one could get more speed replacing linear search in
  SortedCache by smarter method (say binary search) but that would
  require changing data structures.  Also, a lot of time is spent
  searching for weight property -- this raises the question if
  we should have a weight field in the operator representation
  (currently accessing weight requires searching list of
  operator properties).
- this speedup is mainly fruit of sbcl port.  More precisely,
  the hot spots were identified using sbcl profiler.

I belive that there are many other places in Axiom where one
can easily gain significant speedups.  If you want to try compile
sbcl based Axiom.  After starting Axiom do:

)lisp (require :sb-sprof)

)lisp (sb-sprof:start-profiling)

Your computation

)lisp (sb-sprof:stop-profiling)

)lisp (sb-sprof:report)

Notes:  the report is typically rather long, so use something like
script program to capture it.  Also, computation between 5-50 
seconds give best results (shorter are less accurate, longer may
overflow memory when generationg report).

And now the patch:

diff -ru por/wh-20070530/src/algebra/elemntry.spad.pamphlet 
wh-20070530/src/algebra/elemntry.spad.pamphlet
--- por/wh-20070530/src/algebra/elemntry.spad.pamphlet  Fri Jun  1 21:17:44 2007
+++ wh-20070530/src/algebra/elemntry.spad.pamphlet      Sun Jun 10 03:34:26 2007
@@ -623,18 +623,29 @@
       zero? x => 1
       is?(x, oplog) => first argument kernel x
       x < 0 and empty? variables x => inv iexp(-x)
-      h  := inv(2::F)
-      i  := iisqrt1()
-      s2 := h * iisqrt2()
-      s3 := h * iisqrt3()
-      u  := specialTrigs(x / i, [[1,false],[-1,false], [i,false], [-i,false],
-            [h + i * s3,false], [-h + i * s3, false], [-h - i * s3, false],
-             [h - i * s3, false], [s2 + i * s2, false], [-s2 + i * s2, false],
-              [-s2 - i * s2, false], [s2 - i * s2, false], [s3 + i * h, false],
-               [-s3 + i * h, false], [-s3 - i * h, false], [s3 - i * h, 
false]])
-      u case F => u :: F
+      R has RetractableTo Z =>
+        i  := iisqrt1()
+        xi := x / i
+        y := xi / pi()
+        -- this test saves us a lot of effort in the common case
+        -- when no trigonometic simplifiaction is possible
+        retractIfCan(y)@Union(Fraction Z, "failed") case "failed" =>
+          kernel(opexp, x)
+        h  := inv(2::F)
+        s2 := h * iisqrt2()
+        s3 := h * iisqrt3()
+        u  := specialTrigs(xi, [[1,false],[-1,false], [i,false],
+                 [-i,false], [h + i * s3,false], [-h + i * s3, false],
+                 [-h - i * s3, false], [h - i * s3, false],
+                 [s2 + i * s2, false], [-s2 + i * s2, false],
+                 [-s2 - i * s2, false], [s2 - i * s2, false],
+                 [s3 + i * h, false], [-s3 + i * h, false],
+                 [-s3 - i * h, false], [s3 - i * h, false]])
+        u case F => u :: F
+        kernel(opexp, x)
       kernel(opexp, x)
 
+
 -- THIS DETERMINES WHEN TO PERFORM THE log exp f -> f SIMPLIFICATION
 -- CURRENT BEHAVIOR:
 --     IF R IS COMPLEX(S) THEN ONLY ELEMENTS WHICH ARE RETRACTABLE TO R
diff -ru por/wh-20070530/src/algebra/float.spad.pamphlet 
wh-20070530/src/algebra/float.spad.pamphlet
--- por/wh-20070530/src/algebra/float.spad.pamphlet     Fri Jun  1 21:17:44 2007
+++ wh-20070530/src/algebra/float.spad.pamphlet Sun Jun 10 03:37:25 2007
@@ -636,6 +636,8 @@
    shift(x:%,n:I) == [x.mantissa,x.exponent+n]
 
    x = y ==
+      x.exponent = y.exponent =>
+          x.mantissa = y.mantissa
       order x = order y and sign x = sign y and zero? (x - y)
    x < y ==
       y.mantissa = 0 => x.mantissa < 0
diff -ru por/wh-20070530/src/algebra/gdpoly.spad.pamphlet 
wh-20070530/src/algebra/gdpoly.spad.pamphlet
--- por/wh-20070530/src/algebra/gdpoly.spad.pamphlet    Fri Jun  1 21:17:44 2007
+++ wh-20070530/src/algebra/gdpoly.spad.pamphlet        Sun Jun 10 03:38:56 2007
@@ -87,7 +87,18 @@
         ground?(p) => leadingCoefficient p
         "failed"
 
-      degree(p: %,v: OV) == degree(univariate(p,v))
+      degree(p: %,v: OV) ==
+         -- degree(univariate(p,v))
+         zero?(p) => 0
+         res : NonNegativeInteger := 0
+         locv := lookup v
+         while not empty? p repeat
+             t := first p
+             j := t.k.locv
+             if j > res then res := j
+             p := rest p
+         res
+
       minimumDegree(p: %,v: OV) == minimumDegree(univariate(p,v))
       differentiate(p: %,v: OV) ==
             multivariate(differentiate(univariate(p,v)),v)
diff -ru por/wh-20070530/src/algebra/op.spad.pamphlet 
wh-20070530/src/algebra/op.spad.pamphlet
--- por/wh-20070530/src/algebra/op.spad.pamphlet        Fri Jun  1 21:17:44 2007
+++ wh-20070530/src/algebra/op.spad.pamphlet    Sat Jun  9 19:53:20 2007
@@ -156,6 +156,7 @@
 -- property EQUAL? contains a function f: (BOP, BOP) -> Boolean
 -- such that f(o1, o2) is true iff o1 = o2
     op1 = op2 ==
+      (EQ$Lisp)(op1, op2) => true
       name(op1) ^= name(op2) => false
       op1.narg ^= op2.narg => false
       brace(keys properties op1)^=$Set(String) brace(keys properties op2) => 
false
diff -ru por/wh-20070530/src/algebra/vector.spad.pamphlet 
wh-20070530/src/algebra/vector.spad.pamphlet
--- por/wh-20070530/src/algebra/vector.spad.pamphlet    Fri Jun  1 21:17:44 2007
+++ wh-20070530/src/algebra/vector.spad.pamphlet        Sun Jun 10 03:40:22 2007
@@ -344,7 +344,11 @@
         same?: % -> Boolean
         same? z == every?(#1 = z(minIndex z), z)
  
-        x = y == _and/[qelt(x,i)$Rep = qelt(y,i)$Rep for i in 1..dim]
+        x = y ==
+          for i in 1..dim repeat
+            ^(qelt(x,i)$Rep = qelt(y,i)$Rep) => return false
+          true
+          -- _and/[qelt(x,i)$Rep = qelt(y,i)$Rep for i in 1..dim]
  
         retract(z:%):R ==
           same? z => z(minIndex z)

  
-- 
                              Waldek Hebisch
address@hidden