2 files changed, 51 insertions, 13 deletions

diff --git a/examples/speed.hs b/examples/speed.hs
index 865e3e1..22d7220 100644
--- a/examples/speed.hs
+++ b/examples/speed.hs
@@ -1,21 +1,42 @@
 {- speed tests
 
-$ ghc --make -O speed
-In my machine:
-$ ./speed 5 100000 1
-(3><3)
- [ -1.7877285611885504e-2, 0.0,    -0.9998401885597121
- ,                    0.0, 1.0,                    0.0
- ,     0.9998401885597168, 0.0, -1.7877285611891697e-2 ]
-0.29 CPU seconds
+GNU-Octave (see speed.m in this folder):
 
-GNU-Octave:
 ./speed.m
   -0.017877255967426   0.000000000000000  -0.999840189089781
    0.000000000000000   1.000000000000000   0.000000000000000
    0.999840189089763   0.000000000000000  -0.017877255967417
 9.69 seconds
 
+Mathematica:
+
+rot[a_]:={{ Cos[a], 0, Sin[a]},
+          { 0,      1,      0},
+          { -Sin[a],0,Cos[a]}}//N
+
+test := Timing[
+    n = 100000;
+    angles = Range[0.,n-1]/(n-1);
+    Fold[Dot,IdentityMatrix[3],rot/@angles] // MatrixForm
+]
+
+2.08013 Second
+     {{\(-0.017877255967432837`\), "0.`", \(-0.9998401890898042`\)},
+      {"0.`", "1.`", "0.`"},
+      {"0.9998401890898042`", "0.`", \(-0.017877255967432837`\)}}
+
+$ ghc --make -O speed
+
+$ ./speed 5 100000 1
+(3><3)
+ [ -1.7877255967425523e-2, 0.0,    -0.9998401890897632
+ ,                    0.0, 1.0,                    0.0
+ ,      0.999840189089781, 0.0, -1.7877255967416586e-2 ]
+0.33 CPU seconds
+
+cos 50000 = -0.0178772559665563
+sin 50000 = -0.999840189089790
+
 -}
 
 import Numeric.LinearAlgebra
@@ -31,7 +52,7 @@ timing act = do
     t0 <- getCPUTime
     act
     t1 <- getCPUTime
-    printf "%.2f CPU seconds\n" $ (fromIntegral ((t1 - t0) `div` (10^10)) / 100 :: Double)
+    printf "%.2f CPU seconds\n" $ (fromIntegral ((t1 - t0) `div` (10^10)) / 100 :: Double) :: IO ()
 
 op a b = trans $ (trans a) <> (trans b)
 
@@ -52,9 +73,8 @@ rot a = (3><3) [ c,0,s
     where c = cos a
           s = sin a
 
-
-fun n r = foldl1' (<>) (map r [0,delta..1]) where delta = 1 /(fromIntegral n-1)
-
+fun n r = foldl1' (<>) (map r angles)
+    where angles = toList $ linspace n (0,1)
 
 main = do
     args <- getArgs
diff --git a/examples/tests.hs b/examples/tests.hs
index 7b85a9b..90437bb 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -13,6 +13,7 @@ import Test.QuickCheck hiding (test)
 import Test.HUnit hiding ((~:),test)
 import System.Random(randomRs,mkStdGen)
 import System.Info
+import Data.List(foldl1')
 
 type RM = Matrix Double
 type CM = Matrix (Complex Double)
@@ -310,6 +311,20 @@ mulF a b = multiply' ColumnMajor a b
 
 ---------------------------------------------------------------------
 
+rot :: Double -> Matrix Double
+rot a = (3><3) [ c,0,s
+               , 0,1,0
+               ,-s,0,c ]
+    where c = cos a
+          s = sin a
+
+fun n = foldl1' (<>) (map rot angles)
+    where angles = toList $ linspace n (0,1)
+
+rotTest = fun (10^5) :~12~: rot 5E4
+
+---------------------------------------------------------------------
+
 tests = do
     putStrLn "--------- internal -----"
     quickCheck ((\m -> m == trans m).sym :: Sym Double -> Bool)
@@ -322,6 +337,9 @@ tests = do
     quickCheck $ \m -> m == asFortran (m :: CM)
     quickCheck $ \m -> m == (asC . asFortran) (m :: RM)
     quickCheck $ \m -> m == (asC . asFortran) (m :: CM)
+    runTestTT $ TestList
+     [ test "1E5 rots" rotTest
+     ]
     putStrLn "--------- multiply ----"
     quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: RM)
     quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: CM)