From a2771aa1669604a1537e67b51ad364c8a6e011c7 Mon Sep 17 00:00:00 2001
From: Joe Groff <arcata@gmail.com>
Date: Wed, 30 Sep 2009 12:58:32 -0500
Subject: [PATCH] pit math.matrices and math.matrices.simd against each other
 in calculating matrix exponentials

---
 .../matrix-exponential-scalar.factor          | 24 +++++++++++++++++++
 .../matrix-exponential-simd.factor            | 18 ++++++++++++++
 extra/math/matrices/simd/simd.factor          | 15 +++++++++++-
 3 files changed, 56 insertions(+), 1 deletion(-)
 create mode 100644 extra/benchmark/matrix-exponential-scalar/matrix-exponential-scalar.factor
 create mode 100644 extra/benchmark/matrix-exponential-simd/matrix-exponential-simd.factor

diff --git a/extra/benchmark/matrix-exponential-scalar/matrix-exponential-scalar.factor b/extra/benchmark/matrix-exponential-scalar/matrix-exponential-scalar.factor
new file mode 100644
index 0000000000..de4bf1ffe7
--- /dev/null
+++ b/extra/benchmark/matrix-exponential-scalar/matrix-exponential-scalar.factor
@@ -0,0 +1,24 @@
+USING: locals math math.combinatorics math.matrices
+prettyprint sequences typed ;
+IN: benchmark.matrix-exponential-scalar
+
+:: e^m ( m iterations -- e^m )
+    {
+        { 0.0 0.0 0.0 0.0 }
+        { 0.0 0.0 0.0 0.0 }
+        { 0.0 0.0 0.0 0.0 }
+        { 0.0 0.0 0.0 0.0 }
+    }
+    iterations iota [| i |
+        m i m^n i factorial >float m/n m+
+    ] each ;
+
+:: matrix-e ( -- )
+    f :> result!
+    4 identity-matrix :> i4
+    10000 [
+        i4 20 e^m result!
+    ] times
+    result . ;
+
+MAIN: matrix-e
diff --git a/extra/benchmark/matrix-exponential-simd/matrix-exponential-simd.factor b/extra/benchmark/matrix-exponential-simd/matrix-exponential-simd.factor
new file mode 100644
index 0000000000..a23b3f2843
--- /dev/null
+++ b/extra/benchmark/matrix-exponential-simd/matrix-exponential-simd.factor
@@ -0,0 +1,18 @@
+USING: locals math math.combinatorics math.matrices.simd
+prettyprint sequences typed ;
+IN: benchmark.matrix-exponential-simd
+
+TYPED:: e^m4 ( m: matrix4 iterations: fixnum -- e^m: matrix4 )
+    zero-matrix4
+    iterations iota [| i |
+        m i m4^n i factorial >float m4/n m4+
+    ] each ;
+
+:: matrix-e ( -- )
+    f :> result!
+    10000 [
+        identity-matrix4 20 e^m4 result!
+    ] times
+    result . ;
+
+MAIN: matrix-e
diff --git a/extra/math/matrices/simd/simd.factor b/extra/math/matrices/simd/simd.factor
index 014cd86265..0c4c3e1866 100644
--- a/extra/math/matrices/simd/simd.factor
+++ b/extra/math/matrices/simd/simd.factor
@@ -1,6 +1,6 @@
 ! (c)Joe Groff bsd license
 USING: accessors classes.struct generalizations kernel locals
-math math.functions math.matrices.simd math.vectors
+math math.combinatorics math.functions math.matrices.simd math.vectors
 math.vectors.simd sequences sequences.private specialized-arrays
 typed ;
 QUALIFIED-WITH: alien.c-types c
@@ -105,6 +105,19 @@ CONSTANT: identity-matrix4
         }
     }
 
+CONSTANT: zero-matrix4
+    S{ matrix4 f
+        float-4-array{
+            float-4{ 0.0 0.0 0.0 0.0 }
+            float-4{ 0.0 0.0 0.0 0.0 }
+            float-4{ 0.0 0.0 0.0 0.0 }
+            float-4{ 0.0 0.0 0.0 0.0 }
+        }
+    }
+
+TYPED:: m4^n ( m: matrix4 n: fixnum -- m^n: matrix4 )
+    identity-matrix4 n [ m m4. ] times ;
+
 TYPED:: scale-matrix4 ( factors: float-4 -- matrix: matrix4 )
     matrix4 (struct) :> c