pit math.matrices and math.matrices.simd against each other in calculating matrix exponentials

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

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

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