diff --git a/README b/README
--- a/README
+++ b/README
@@ -33,20 +33,30 @@
 
 INSTALLATION --------------------------------------
 
-$ runhaskell Setup.lhs configure --prefix=$HOME --user
-$ runhaskell Setup.lhs build
-$ runhaskell Setup.lhs haddock
-$ runhaskell Setup.lhs install
+Automatic (using cabal-install and HackageDB):
 
+    $ cabal install hmatrix
+
+Manual:
+
+    Install storable-complex from HackageDB and then
+
+    $ runhaskell Setup.lhs configure --prefix=$HOME --user
+    $ runhaskell Setup.lhs build
+    $ runhaskell Setup.lhs haddock
+    $ runhaskell Setup.lhs install
+
 See below for installation on Windows.
 
 TESTS ---------------------------------------------
 
-$ runhaskell examples/tests
+$ ghci
+GHCi, version 6.8.2: http://www.haskell.org/ghc/  :? for help
+Prelude> Numeric.LinearAlgebra.Tests.runTests 20
 
 Additional tests with big matrices (taking a few minutes):
 
-$ runhaskell examples/tests --big
+$ runhaskell examples/experiments bigtests
 
 EXAMPLES ------------------------------------------------------
 
diff --git a/examples/experiments/bigtests.hs b/examples/experiments/bigtests.hs
new file mode 100644
--- /dev/null
+++ b/examples/experiments/bigtests.hs
@@ -0,0 +1,60 @@
+module Main where
+
+import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra.Tests
+import System.Random(randomRs,mkStdGen)
+import Test.HUnit hiding (test)
+import System(getArgs)
+
+
+pseudorandomR seed (n,m) = reshape m $ fromList $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed
+
+pseudorandomC seed (n,m) = toComplex (pseudorandomR seed (n,m), pseudorandomR (seed+1) (n,m))
+
+bigmat = m + trans m
+    where m = pseudorandomR 18 (1000,1000) :: Matrix Double
+bigmatc = mc + ctrans mc
+    where mc = pseudorandomC 19 (1000,1000) :: Matrix (Complex Double)
+
+utest str b = TestCase $ assertBool str b
+
+feye n = flipud (ident n) :: Matrix Double
+
+infixl 4 |~|
+a |~| b = dist a b < 10^^(-10)
+
+dist a b = r
+    where norm = pnorm Infinity
+          na = norm a
+          nb = norm b
+          nab = norm (a-b)
+          mx = max na nb
+          mn = min na nb
+          r = if mn < eps
+                then mx
+                else nab/mx
+
+square m = rows m == cols m
+
+unitary m = square m && m <> ctrans m |~| ident (rows m)
+
+eigProp m = complex m <> v |~| v <> diag s
+    where (s, v) = eig m
+
+eigSHProp m = m <> v |~| v <> real (diag s)
+              && unitary v
+              && m |~| v <> real (diag s) <> ctrans v
+    where (s, v) = eigSH m
+
+bigtests = do
+    putStrLn "--------- big matrices -----"
+    runTestTT $ TestList
+     [ utest "eigS" $ eigSHProp bigmat
+     , utest "eigH" $ eigSHProp bigmatc
+     , utest "eigR" $ eigProp   bigmat
+     , utest "eigC" $ eigProp   bigmatc
+     , utest "det"  $ det (feye 1000) == 1 && det (feye 1002) == -1
+     ]
+    return ()
+
+main = bigtests
diff --git a/examples/pca1.hs b/examples/pca1.hs
--- a/examples/pca1.hs
+++ b/examples/pca1.hs
@@ -26,7 +26,7 @@
     decode x = x <> vp + m
     m = mean dataSet
     c = cov dataSet
-    (_,v) = eigSH c
+    (_,v) = eigSH' c
     vp = takeRows n (trans v)
 
 norm = pnorm PNorm2
diff --git a/examples/tests.hs b/examples/tests.hs
deleted file mode 100644
--- a/examples/tests.hs
+++ /dev/null
@@ -1,506 +0,0 @@
-{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
-
-module Main where
-
-import Numeric.GSL hiding (sin,cos,exp,choose)
-import Numeric.LinearAlgebra
-import Numeric.LinearAlgebra.LAPACK
-import qualified Numeric.GSL.Matrix as GSL
-import Test.QuickCheck hiding (test)
-import Test.HUnit hiding ((~:),test)
-import System.Random(randomRs,mkStdGen)
-import System.Info
-import Data.List(foldl1', transpose)
-import System(getArgs)
-import Debug.Trace(trace)
-
-debug x = trace (show x) x
-
-type RM = Matrix Double
-type CM = Matrix (Complex Double)
-
--- relative error
-dist :: (Normed t, Num t) => t -> t -> Double
-dist a b = r
-    where norm = pnorm Infinity
-          na = norm a
-          nb = norm b
-          nab = norm (a-b)
-          mx = max na nb
-          mn = min na nb
-          r = if mn < eps
-                then mx
-                else nab/mx
-
-infixl 4 |~|
-a |~| b = a :~10~: b
-
-data Aprox a = (:~) a Int
-
-(~:) :: (Normed a, Num a) => Aprox a -> a -> Bool
-a :~n~: b = dist a b < 10^^(-n)
-
-
-maxdim = 10
-
-instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where
-    arbitrary = do
-        r <- arbitrary
-        i <- arbitrary
-        return (r:+i)
-    coarbitrary = undefined
-
-instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where 
-   arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)
-                  m <- choose (1,maxdim)
-                  n <- choose (1,maxdim)
-                  l <- vector (m*n)
-                  ctype <- arbitrary
-                  let h = if ctype then (m><n) else (m>|<n)
-                  trMode <- arbitrary
-                  let tr = if trMode then trans else id
-                  return $ tr (h l)
-   coarbitrary = undefined
-
-data PairM a = PairM (Matrix a) (Matrix a) deriving Show
-instance (Num a, Element a, Arbitrary a) => Arbitrary (PairM a) where
-    arbitrary = do
-        a <- choose (1,maxdim)
-        b <- choose (1,maxdim)
-        c <- choose (1,maxdim)
-        l1 <- vector (a*b)
-        l2 <- vector (b*c)
-        return $ PairM ((a><b) (map fromIntegral (l1::[Int]))) ((b><c) (map fromIntegral (l2::[Int])))
-        --return $ PairM ((a><b) l1) ((b><c) l2)
-    coarbitrary = undefined
-
-data SqM a = SqM (Matrix a) deriving Show
-sqm (SqM a) = a
-instance (Element a, Arbitrary a) => Arbitrary (SqM a) where
-    arbitrary = do
-        n <- choose (1,maxdim)
-        l <- vector (n*n)
-        return $ SqM $ (n><n) l
-    coarbitrary = undefined
-
-data Sym a = Sym (Matrix a) deriving Show
-sym (Sym a) = a
-instance (Linear Vector a, Arbitrary a) => Arbitrary (Sym a) where
-    arbitrary = do
-        SqM m <- arbitrary
-        return $ Sym (m + trans m)
-    coarbitrary = undefined
-
-data Her = Her (Matrix (Complex Double)) deriving Show
-her (Her a) = a
-instance {-(Field a, Arbitrary a, Num a) =>-} Arbitrary Her where
-    arbitrary = do
-        SqM m <- arbitrary
-        return $ Her (m + ctrans m)
-    coarbitrary = undefined
-
-data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show
-instance (Num a, Field a, Arbitrary a) => Arbitrary (PairSM a) where
-    arbitrary = do
-        a <- choose (1,maxdim)
-        c <- choose (1,maxdim)
-        l1 <- vector (a*a)
-        l2 <- vector (a*c)
-        return $ PairSM ((a><a) (map fromIntegral (l1::[Int]))) ((a><c) (map fromIntegral (l2::[Int])))
-        --return $ PairSM ((a><a) l1) ((a><c) l2)
-    coarbitrary = undefined
-
-instance (Field a, Arbitrary a) => Arbitrary (Vector a) where 
-   arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)
-                  m <- choose (1,maxdim^2)
-                  l <- vector m
-                  return $ fromList l
-   coarbitrary = undefined
-
-data PairV a = PairV (Vector a) (Vector a)
-instance (Field a, Arbitrary a) => Arbitrary (PairV a) where 
-   arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)
-                  m <- choose (1,maxdim^2)
-                  l1 <- vector m
-                  l2 <- vector m
-                  return $ PairV (fromList l1) (fromList l2)
-   coarbitrary = undefined
-
-----------------------------------------------------------------------
-
-test str b = TestCase $ assertBool str b
-
-----------------------------------------------------------------------
-
-pseudorandomR seed (n,m) = reshape m $ fromList $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed 
-
-pseudorandomC seed (n,m) = toComplex (pseudorandomR seed (n,m), pseudorandomR (seed+1) (n,m))
-
-bigmat = m + trans m :: RM
-    where m = pseudorandomR 18 (1000,1000)
-bigmatc = mc + ctrans mc ::CM
-    where mc = pseudorandomC 19 (1000,1000)
-
-----------------------------------------------------------------------
-
-
-m = (3><3)
- [ 1, 2, 3
- , 4, 5, 7
- , 2, 8, 4 :: Double
- ]
-
-mc = (3><3)
- [ 1, 2, 3
- , 4, 5, 7
- , 2, 8, i
- ]
-
-
-mr = (3><4)
- [ 1, 2, 3, 4,
-   2, 4, 6, 8,
-   1, 1, 1, 2:: Double
- ]
-
-mrc = (3><4)
- [ 1, 2, 3, 4,
-   2, 4, 6, 8,
-   i, i, i, 2
- ]
-
-a = (3><4)
- [ 1, 0, 0, 0
- , 0, 2, 0, 0
- , 0, 0, 0, 0 :: Double
- ]
-
-b = (3><4)
- [ 1, 0, 0, 0
- , 0, 2, 3, 0
- , 0, 0, 4, 0 :: Double
- ]
-
-ac = (2><3) [1 .. 6::Double]
-bc = (3><4) [7 .. 18::Double]
-
-af = (2>|<3) [1,4,2,5,3,6::Double]
-bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double]
-
--------------------------------------------------------
-
-feye n = flipud (ident n) :: Matrix Double
-
-
-luTest1 m = m |~| p <> l <> u
-    where (l,u,p,_) = lu m
-
-detTest1 = det m == 26
-        && det mc == 38 :+ (-3)
-        && det (feye 2) == -1
-
-detTest2 m = s d1 |~| s d2
-    where d1 = det m
-          d2 = det' m * det q
-          det' m = product $ toList $ takeDiag r
-          (q,r) = qr m
-          s x = fromList [x]
-
-invTest m = degenerate m || m <> inv m |~| ident (rows m)
-
-pinvTest m =  m <> p <> m |~| m
-           && p <> m <> p |~| p
-           && hermitian (m<>p)
-           && hermitian (p<>m)
-    where p = pinv m
-
-square m = rows m == cols m
-
-unitary m = square m && m <> ctrans m |~| ident (rows m)
-
-hermitian m = m |~| ctrans m
-
-upperTriang m = rows m == 1 || down == z
-    where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
-          z = constant 0 (dim down)
-
-upperHessenberg m = rows m < 3 || down == z
-    where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m))
-          z = constant 0 (dim down)
-
-svdTest svd m = u <> real d <> trans v |~| m
-          && unitary u && unitary v
-    where (u,d,v) = full svd m
-
-svdTest' svd m = m |~| 0 || u <> real (diag s) <> trans v |~| m
-    where (u,s,v) = economy svd m
-
-eigTest m = complex m <> v |~| v <> diag s
-    where (s, v) = eig m
-
-eigTestSH m = m <> v |~| v <> real (diag s)
-              && unitary v
-              && m |~| v <> real (diag s) <> ctrans v
-    where (s, v) = eigSH m
-
-zeros (r,c) = reshape c (constant 0 (r*c))
-
-ones (r,c) = zeros (r,c) + 1
-
-degenerate m = rank m < min (rows m) (cols m) 
-
-prec = 1E-15
-
-singular m = s1 < prec || s2/s1 < prec
-    where (_,ss,_) = svd m
-          s = toList ss
-          s1 = maximum s
-          s2 = minimum s
-
-nullspaceTest m = null nl || m <> n |~| zeros (r,c) -- 0
-    where nl = nullspacePrec 1 m
-          n = fromColumns nl
-          r = rows m
-          c = cols m - rank m
-
---------------------------------------------------------------------
-
-polyEval cs x = foldr (\c ac->ac*x+c) 0 cs
-
-polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p))
-
-
-polySolveTest = test "polySolve" (polySolveTest' [1,2,3,4])
-
----------------------------------------------------------------------
-
-quad f a b = fst $ integrateQAGS 1E-9 100 f a b
-
--- A multiple integral can be easily defined using partial application
-quad2 f a b g1 g2 = quad h a b
-    where h x = quad (f x) (g1 x) (g2 x)
-
-volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) 
-                        0 r (const 0) (\x->sqrt (r*r-x*x))
-
-epsTol = 1E-8::Double
-
-integrateTest = test "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < epsTol)
-
----------------------------------------------------------------------
-
-besselTest = test "bessel_J0_e" ( abs (r-expected) < e )
-    where (r,e) = bessel_J0_e 5.0
-          expected = -0.17759677131433830434739701
-
-exponentialTest = test "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 )
-    where (v,e,err) = exp_e10_e 30.0
-          expected = exp 30.0
-
-gammaTest = test "gamma" (gamma 5 == 24.0)
-
----------------------------------------------------------------------
-
-cholRTest = chol ((2><2) [1,2,2,9::Double]) == (2><2) [1,2,0,2.23606797749979]
-cholCTest = chol ((2><2) [1,2,2,9::Complex Double]) == (2><2) [1,2,0,2.23606797749979]
-
----------------------------------------------------------------------
-
-qrTest qr m = q <> r |~| m && unitary q && upperTriang r
-    where (q,r) = qr m
-
----------------------------------------------------------------------
-
-hessTest m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h
-    where (p,h) = hess m
-
----------------------------------------------------------------------
-
-schurTest1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s
-    where (u,s) = schur m
-
-schurTest2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme
-    where (u,s) = schur m
-
----------------------------------------------------------------------
-
-nd1 = (3><3) [ 1/2, 1/4, 1/4
-             , 0/1, 1/2, 1/4
-             , 1/2, 1/4, 1/2 :: Double]
-
-nd2 = (2><2) [1, 0, 1, 1:: Complex Double]
-
-expmTest1 = expm nd1 :~14~: (3><3)
- [ 1.762110887278176
- , 0.478085470590435
- , 0.478085470590435
- , 0.104719410945666
- , 1.709751181805343
- , 0.425725765117601
- , 0.851451530235203
- , 0.530445176063267
- , 1.814470592751009 ]
-
-expmTest2 = expm nd2 :~15~: (2><2)
- [ 2.718281828459045
- , 0.000000000000000
- , 2.718281828459045
- , 2.718281828459045 ]
-
-expmTestDiag m = expm (logm m) |~| complex m
-    where logm m = matFunc Prelude.log m
-
-
-
----------------------------------------------------------------------
-
-asFortran m = (rows m >|< cols m) $ toList (flatten $ trans  m)
-asC m = (rows m >< cols m) $ toList (flatten m)
-
-mulC a b = a <> b
-mulF a b = trans $ trans b <> trans a
-
--------------------------------------------------------------------------
-
-multiplyG a b = reshape (cols b) $ fromList $ concat $ multiplyL (toLists a) (toLists b)
-    where multiplyL a b = [[dotL x y | y <- transpose b] | x <- a]
-          dotL a b = sum (zipWith (*) a b)
-
-r >|< c = f where
-    f l | dim v == r*c = reshapeF r v
-        | otherwise    = error "(>|<)"
-        where v = fromList l
-    reshapeF r = trans . reshape r
-
----------------------------------------------------------------------
-
-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
-    setErrorHandlerOff
-    putStrLn "--------- internal -----"
-    quickCheck ((\m -> m == trans m).sym :: Sym Double -> Bool)
-    quickCheck ((\m -> m == trans m).sym :: Sym (Complex Double) -> Bool)
-    quickCheck $ \l -> null l || (toList . fromList) l == (l :: [Double])
-    quickCheck $ \l -> null l || (toList . fromList) l == (l :: [Complex Double])
-    quickCheck $ \m -> m == asC (m :: RM)
-    quickCheck $ \m -> m == asC (m :: CM)
-    quickCheck $ \m -> m == asFortran (m :: RM)
-    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)
-    quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == trans (mulF (trans m2) (trans m1 :: RM))
-    quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == trans (mulF (trans m2) (trans m1 :: CM))
-    quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: RM)
-    quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: CM)
-    putStrLn "--------- lu ---------"
-    quickCheck (luTest1 :: RM->Bool)
-    quickCheck (luTest1 :: CM->Bool)
-    quickCheck (detTest2 . sqm  :: SqM Double -> Bool)
-    quickCheck (detTest2 . sqm  :: SqM (Complex Double) -> Bool)
-    runTestTT $ TestList
-     [ test "det1" detTest1
-     ]
-    putStrLn "--------- svd ---------"
-    quickCheck (svdTest svdR)
-    quickCheck (svdTest svdRdd)
-    quickCheck (svdTest svdC)
-    quickCheck (svdTest' svdR)
-    quickCheck (svdTest' svdRdd)
-    quickCheck (svdTest' svdC)
-    quickCheck (svdTest' GSL.svdg)
-    putStrLn "--------- eig ---------"
-    quickCheck (eigTest  . sqm :: SqM Double -> Bool)
-    quickCheck (eigTest  . sqm :: SqM (Complex Double) -> Bool)
-    quickCheck (eigTestSH . sym :: Sym Double -> Bool)
-    quickCheck (eigTestSH . her :: Her -> Bool)
-    putStrLn "--------- inv ------"
-    quickCheck (invTest . sqm   :: SqM Double -> Bool)
-    quickCheck (invTest . sqm   :: SqM (Complex Double) -> Bool)
-    putStrLn "--------- pinv ------"
-    quickCheck (pinvTest ::RM->Bool)
-    if os == "mingw32"
-        then putStrLn "complex pinvTest skipped in this OS"
-        else quickCheck (pinvTest ::CM->Bool)
-    putStrLn "--------- chol ------"
-    runTestTT $ TestList
-     [ test "cholR" cholRTest
-     , test "cholC" cholRTest
-     ]
-    putStrLn "--------- qr ---------"
-    quickCheck (qrTest GSL.qr)
-    quickCheck (qrTest (GSL.unpackQR . GSL.qrPacked))
-    quickCheck (qrTest (    unpackQR . GSL.qrPacked))
-    quickCheck (qrTest qr ::RM->Bool)
-    quickCheck (qrTest qr ::CM->Bool)
-    putStrLn "--------- hess --------"
-    quickCheck (hessTest . sqm ::SqM Double->Bool)
-    quickCheck (hessTest . sqm ::SqM (Complex Double) -> Bool)
-    putStrLn "--------- schur --------"
-    quickCheck (schurTest2 . sqm ::SqM Double->Bool)
-    if os == "mingw32"
-        then putStrLn "complex schur skipped in this OS"
-        else quickCheck (schurTest1 . sqm ::SqM (Complex Double) -> Bool)
-    putStrLn "--------- expm --------"
-    runTestTT $ TestList
-     [ test "expmd" (expmTestDiag $ (2><2) [1,2,3,5 :: Double])
-     , test "expm1" (expmTest1)
-     , test "expm2" (expmTest2)
-     ]
-    putStrLn "--------- nullspace ------"
-    quickCheck (nullspaceTest :: RM -> Bool)
-    quickCheck (nullspaceTest :: CM -> Bool)
-    putStrLn "--------- vector operations ------"
-    quickCheck $ (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM))
-    quickCheck $ (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM))
-    quickCheck $ (\u -> cos u * tan u |~| sin (u::RM))
-    quickCheck $ (\u -> (cos u * tan u) :~6~: sin (u::CM))
-    runTestTT $ TestList
-     [ test "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)
-     , test "arith2" $ (((1+i) .* ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( (140*i-51).*1 :: CM)
-     , test "arith3" $ exp (i.*ones(10,10)*pi) + 1 |~| 0
-     , test "<\\>"   $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3]
-     ]
-    putStrLn "--------- GSL ------"
-    quickCheck $ \v -> ifft (fft v) |~| v
-    runTestTT $ TestList
-     [ gammaTest
-     , besselTest
-     , exponentialTest
-     , integrateTest
-     , polySolveTest
-     ]
-
-bigtests = do
-    putStrLn "--------- big matrices -----"
-    runTestTT $ TestList
-     [ test "eigS" $ eigTestSH bigmat
-     , test "eigH" $ eigTestSH bigmatc
-     , test "eigR" $ eigTest   bigmat
-     , test "eigC" $ eigTest   bigmatc
-     , test "det"  $ det (feye 1000) == 1 && det (feye 1002) == -1
-     ]
-
-main = do
-    args <- getArgs
-    if "--big" `elem` args
-        then bigtests
-        else tests
diff --git a/hmatrix.cabal b/hmatrix.cabal
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.2.0.0
+Version:            0.2.1.0
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -13,26 +13,25 @@
                     .
                     More information: <http://alberrto.googlepages.com/gslhaskell>
 Category:           Numerical, Math
-tested-with:        GHC ==6.6.1, GHC ==6.8.1, GHC ==6.8.2
+tested-with:        GHC ==6.8.2
 
 cabal-version:      >=1.2
+build-type:         Simple
 
 flag splitBase
     description:    Choose the new smaller, split-up base package.
 
 library
     if flag(splitBase)
-      build-depends:    base >= 3, array
+      build-depends:    base >= 3, array, QuickCheck, HUnit, storable-complex
     else
-      build-depends:    base < 3
-
-    ghc-options:        -O
+      build-depends:    base < 3, QuickCheck, HUnit, storable-complex
 
     if os(windows)
-      ghc-options:      -fasm
+      ghc-options:
     else
       if arch(x86_64)
-         ghc-options:   -fasm
+         ghc-options:
       else
          ghc-options:   -fvia-C
 
@@ -84,12 +83,15 @@
                         Numeric.LinearAlgebra.Instances,
                         Numeric.LinearAlgebra.Interface,
                         Numeric.LinearAlgebra.Algorithms,
-                        Graphics.Plot
+                        Graphics.Plot,
+                        Numeric.LinearAlgebra.Tests
     other-modules:      Data.Packed.Internal,
                         Data.Packed.Internal.Common,
                         Data.Packed.Internal.Vector,
                         Data.Packed.Internal.Matrix,
-                        Numeric.GSL.Special.Internal
+                        Numeric.GSL.Special.Internal,
+                        Numeric.LinearAlgebra.Tests.Instances
+                        Numeric.LinearAlgebra.Tests.Properties
     C-sources:          lib/Data/Packed/Internal/auxi.c,
                         lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,
                         lib/Numeric/GSL/gsl-aux.c
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -22,16 +22,8 @@
 import Debug.Trace
 import Foreign.C.String(peekCString)
 import Foreign.C.Types
+import Foreign.Storable.Complex
 
-----------------------------------------------------------------------
-instance (Storable a, RealFloat a) => Storable (Complex a) where    --
-    alignment x = alignment (realPart x)                            --
-    sizeOf x    = 2 * sizeOf (realPart x)                           --
-    peek p = do                                                     --
-        [re,im] <- peekArray 2 (castPtr p)                          --
-        return (re :+ im)                                           --
-    poke p (a :+ b) = pokeArray (castPtr p) [a,b]                   --
-----------------------------------------------------------------------
 
 -- | @debug x = trace (show x) x@
 debug :: (Show a) => a -> a
diff --git a/lib/Numeric/GSL/Special/Airy.hs b/lib/Numeric/GSL/Special/Airy.hs
--- a/lib/Numeric/GSL/Special/Airy.hs
+++ b/lib/Numeric/GSL/Special/Airy.hs
@@ -38,8 +38,8 @@
 , airy_zero_Ai_deriv_e
 , airy_zero_Ai_deriv
 , airy_zero_Bi_deriv_e
-, airy_zero_Bi_deriv,
-Precision(..)
+, airy_zero_Bi_deriv
+, Precision(..)
 ) where
 
 import Foreign(Ptr)
diff --git a/lib/Numeric/GSL/Special/Debye.hs b/lib/Numeric/GSL/Special/Debye.hs
--- a/lib/Numeric/GSL/Special/Debye.hs
+++ b/lib/Numeric/GSL/Special/Debye.hs
@@ -23,10 +23,10 @@
 , debye_3
 , debye_4_e
 , debye_4
---, debye_5_e
---, debye_5
---, debye_6_e
---, debye_6
+, debye_5_e
+, debye_5
+, debye_6_e
+, debye_6
 ) where
 
 import Foreign(Ptr)
@@ -92,27 +92,27 @@
 -- | wrapper for int gsl_sf_debye_5_e(double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_debye_5_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
---debye_5_e :: Double -> (Double,Double)
---debye_5_e x = createSFR "debye_5_e" $ gsl_sf_debye_5_e x
---foreign import ccall "debye.h gsl_sf_debye_5_e" gsl_sf_debye_5_e :: Double -> Ptr () -> IO CInt
+debye_5_e :: Double -> (Double,Double)
+debye_5_e x = createSFR "debye_5_e" $ gsl_sf_debye_5_e x
+foreign import ccall "debye.h gsl_sf_debye_5_e" gsl_sf_debye_5_e :: Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_debye_5(double x);
 --
 --   <http://www.google.com/search?q=gsl_sf_debye_5&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
---debye_5 :: Double -> Double
---debye_5 = gsl_sf_debye_5
---foreign import ccall "debye.h gsl_sf_debye_5" gsl_sf_debye_5 :: Double -> Double
+debye_5 :: Double -> Double
+debye_5 = gsl_sf_debye_5
+foreign import ccall "debye.h gsl_sf_debye_5" gsl_sf_debye_5 :: Double -> Double
 
 -- | wrapper for int gsl_sf_debye_6_e(double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_debye_6_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
---debye_6_e :: Double -> (Double,Double)
---debye_6_e x = createSFR "debye_6_e" $ gsl_sf_debye_6_e x
---foreign import ccall "debye.h gsl_sf_debye_6_e" gsl_sf_debye_6_e :: Double -> Ptr () -> IO CInt
+debye_6_e :: Double -> (Double,Double)
+debye_6_e x = createSFR "debye_6_e" $ gsl_sf_debye_6_e x
+foreign import ccall "debye.h gsl_sf_debye_6_e" gsl_sf_debye_6_e :: Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_debye_6(double x);
 --
 --   <http://www.google.com/search?q=gsl_sf_debye_6&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
---debye_6 :: Double -> Double
---debye_6 = gsl_sf_debye_6
---foreign import ccall "debye.h gsl_sf_debye_6" gsl_sf_debye_6 :: Double -> Double
+debye_6 :: Double -> Double
+debye_6 = gsl_sf_debye_6
+foreign import ccall "debye.h gsl_sf_debye_6" gsl_sf_debye_6 :: Double -> Double
diff --git a/lib/Numeric/GSL/Special/Ellint.hs b/lib/Numeric/GSL/Special/Ellint.hs
--- a/lib/Numeric/GSL/Special/Ellint.hs
+++ b/lib/Numeric/GSL/Special/Ellint.hs
@@ -19,10 +19,10 @@
 , ellint_Kcomp
 , ellint_Ecomp_e
 , ellint_Ecomp
---, ellint_Pcomp_e
---, ellint_Pcomp
---, ellint_Dcomp_e
---, ellint_Dcomp
+, ellint_Pcomp_e
+, ellint_Pcomp
+, ellint_Dcomp_e
+, ellint_Dcomp
 , ellint_F_e
 , ellint_F
 , ellint_E_e
@@ -76,30 +76,30 @@
 -- | wrapper for int gsl_sf_ellint_Pcomp_e(double k,double n,gsl_mode_t mode,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_ellint_Pcomp_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
---ellint_Pcomp_e :: Double -> Double -> Precision -> (Double,Double)
---ellint_Pcomp_e k n mode = createSFR "ellint_Pcomp_e" $ gsl_sf_ellint_Pcomp_e k n  (precCode mode)
---foreign import ccall "ellint.h gsl_sf_ellint_Pcomp_e" gsl_sf_ellint_Pcomp_e :: Double -> Double -> Gsl_mode_t -> Ptr () -> IO CInt
+ellint_Pcomp_e :: Double -> Double -> Precision -> (Double,Double)
+ellint_Pcomp_e k n mode = createSFR "ellint_Pcomp_e" $ gsl_sf_ellint_Pcomp_e k n  (precCode mode)
+foreign import ccall "ellint.h gsl_sf_ellint_Pcomp_e" gsl_sf_ellint_Pcomp_e :: Double -> Double -> Gsl_mode_t -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_ellint_Pcomp(double k,double n,gsl_mode_t mode);
 --
 --   <http://www.google.com/search?q=gsl_sf_ellint_Pcomp&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
---ellint_Pcomp :: Double -> Double -> Precision -> Double
---ellint_Pcomp k n mode = gsl_sf_ellint_Pcomp k n  (precCode mode)
---foreign import ccall "ellint.h gsl_sf_ellint_Pcomp" gsl_sf_ellint_Pcomp :: Double -> Double -> Gsl_mode_t -> Double
+ellint_Pcomp :: Double -> Double -> Precision -> Double
+ellint_Pcomp k n mode = gsl_sf_ellint_Pcomp k n  (precCode mode)
+foreign import ccall "ellint.h gsl_sf_ellint_Pcomp" gsl_sf_ellint_Pcomp :: Double -> Double -> Gsl_mode_t -> Double
 
 -- | wrapper for int gsl_sf_ellint_Dcomp_e(double k,gsl_mode_t mode,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_ellint_Dcomp_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
---ellint_Dcomp_e :: Double -> Precision -> (Double,Double)
---ellint_Dcomp_e k mode = createSFR "ellint_Dcomp_e" $ gsl_sf_ellint_Dcomp_e k  (precCode mode)
---foreign import ccall "ellint.h gsl_sf_ellint_Dcomp_e" gsl_sf_ellint_Dcomp_e :: Double -> Gsl_mode_t -> Ptr () -> IO CInt
+ellint_Dcomp_e :: Double -> Precision -> (Double,Double)
+ellint_Dcomp_e k mode = createSFR "ellint_Dcomp_e" $ gsl_sf_ellint_Dcomp_e k  (precCode mode)
+foreign import ccall "ellint.h gsl_sf_ellint_Dcomp_e" gsl_sf_ellint_Dcomp_e :: Double -> Gsl_mode_t -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_ellint_Dcomp(double k,gsl_mode_t mode);
 --
 --   <http://www.google.com/search?q=gsl_sf_ellint_Dcomp&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
---ellint_Dcomp :: Double -> Precision -> Double
---ellint_Dcomp k mode = gsl_sf_ellint_Dcomp k  (precCode mode)
---foreign import ccall "ellint.h gsl_sf_ellint_Dcomp" gsl_sf_ellint_Dcomp :: Double -> Gsl_mode_t -> Double
+ellint_Dcomp :: Double -> Precision -> Double
+ellint_Dcomp k mode = gsl_sf_ellint_Dcomp k  (precCode mode)
+foreign import ccall "ellint.h gsl_sf_ellint_Dcomp" gsl_sf_ellint_Dcomp :: Double -> Gsl_mode_t -> Double
 
 -- | wrapper for int gsl_sf_ellint_F_e(double phi,double k,gsl_mode_t mode,gsl_sf_result* result);
 --
diff --git a/lib/Numeric/GSL/Special/Legendre.hs b/lib/Numeric/GSL/Special/Legendre.hs
--- a/lib/Numeric/GSL/Special/Legendre.hs
+++ b/lib/Numeric/GSL/Special/Legendre.hs
@@ -15,24 +15,42 @@
 ------------------------------------------------------------
 
 module Numeric.GSL.Special.Legendre(
-  legendre_Pl
+  legendre_Pl_e
+, legendre_Pl
+, legendre_P1_e
+, legendre_P2_e
+, legendre_P3_e
 , legendre_P1
 , legendre_P2
 , legendre_P3
+, legendre_Q0_e
 , legendre_Q0
+, legendre_Q1_e
 , legendre_Q1
+, legendre_Ql_e
 , legendre_Ql
+, legendre_Plm_e
 , legendre_Plm
+, legendre_sphPlm_e
 , legendre_sphPlm
 , legendre_array_size
+, conicalP_half_e
 , conicalP_half
+, conicalP_mhalf_e
 , conicalP_mhalf
+, conicalP_0_e
 , conicalP_0
+, conicalP_1_e
 , conicalP_1
+, conicalP_sph_reg_e
 , conicalP_sph_reg
+, conicalP_cyl_reg_e
 , conicalP_cyl_reg
+, legendre_H3d_0_e
 , legendre_H3d_0
+, legendre_H3d_1_e
 , legendre_H3d_1
+, legendre_H3d_e
 , legendre_H3d
 ) where
 
@@ -40,12 +58,12 @@
 import Foreign.C.Types(CInt)
 import Numeric.GSL.Special.Internal
 
--- | wrapper for int gsl_sf_legendre_Pl_e(int l,double x,double* result);
+-- | wrapper for int gsl_sf_legendre_Pl_e(int l,double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_legendre_Pl_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-legendre_Pl_e :: CInt -> Double -> Ptr Double -> CInt
-legendre_Pl_e = gsl_sf_legendre_Pl_e
-foreign import ccall "legendre.h gsl_sf_legendre_Pl_e" gsl_sf_legendre_Pl_e :: CInt -> Double -> Ptr Double -> CInt
+legendre_Pl_e :: CInt -> Double -> (Double,Double)
+legendre_Pl_e l x = createSFR "legendre_Pl_e" $ gsl_sf_legendre_Pl_e l x
+foreign import ccall "legendre.h gsl_sf_legendre_Pl_e" gsl_sf_legendre_Pl_e :: CInt -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_legendre_Pl(int l,double x);
 --
@@ -68,26 +86,26 @@
 legendre_Pl_deriv_array = gsl_sf_legendre_Pl_deriv_array
 foreign import ccall "legendre.h gsl_sf_legendre_Pl_deriv_array" gsl_sf_legendre_Pl_deriv_array :: CInt -> Double -> Ptr Double -> Ptr Double -> CInt
 
--- | wrapper for int gsl_sf_legendre_P1_e(double x,double* result);
+-- | wrapper for int gsl_sf_legendre_P1_e(double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_legendre_P1_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-legendre_P1_e :: Double -> Ptr Double -> CInt
-legendre_P1_e = gsl_sf_legendre_P1_e
-foreign import ccall "legendre.h gsl_sf_legendre_P1_e" gsl_sf_legendre_P1_e :: Double -> Ptr Double -> CInt
+legendre_P1_e :: Double -> (Double,Double)
+legendre_P1_e x = createSFR "legendre_P1_e" $ gsl_sf_legendre_P1_e x
+foreign import ccall "legendre.h gsl_sf_legendre_P1_e" gsl_sf_legendre_P1_e :: Double -> Ptr () -> IO CInt
 
--- | wrapper for int gsl_sf_legendre_P2_e(double x,double* result);
+-- | wrapper for int gsl_sf_legendre_P2_e(double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_legendre_P2_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-legendre_P2_e :: Double -> Ptr Double -> CInt
-legendre_P2_e = gsl_sf_legendre_P2_e
-foreign import ccall "legendre.h gsl_sf_legendre_P2_e" gsl_sf_legendre_P2_e :: Double -> Ptr Double -> CInt
+legendre_P2_e :: Double -> (Double,Double)
+legendre_P2_e x = createSFR "legendre_P2_e" $ gsl_sf_legendre_P2_e x
+foreign import ccall "legendre.h gsl_sf_legendre_P2_e" gsl_sf_legendre_P2_e :: Double -> Ptr () -> IO CInt
 
--- | wrapper for int gsl_sf_legendre_P3_e(double x,double* result);
+-- | wrapper for int gsl_sf_legendre_P3_e(double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_legendre_P3_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-legendre_P3_e :: Double -> Ptr Double -> CInt
-legendre_P3_e = gsl_sf_legendre_P3_e
-foreign import ccall "legendre.h gsl_sf_legendre_P3_e" gsl_sf_legendre_P3_e :: Double -> Ptr Double -> CInt
+legendre_P3_e :: Double -> (Double,Double)
+legendre_P3_e x = createSFR "legendre_P3_e" $ gsl_sf_legendre_P3_e x
+foreign import ccall "legendre.h gsl_sf_legendre_P3_e" gsl_sf_legendre_P3_e :: Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_legendre_P1(double x);
 --
@@ -110,12 +128,12 @@
 legendre_P3 = gsl_sf_legendre_P3
 foreign import ccall "legendre.h gsl_sf_legendre_P3" gsl_sf_legendre_P3 :: Double -> Double
 
--- | wrapper for int gsl_sf_legendre_Q0_e(double x,double* result);
+-- | wrapper for int gsl_sf_legendre_Q0_e(double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_legendre_Q0_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-legendre_Q0_e :: Double -> Ptr Double -> CInt
-legendre_Q0_e = gsl_sf_legendre_Q0_e
-foreign import ccall "legendre.h gsl_sf_legendre_Q0_e" gsl_sf_legendre_Q0_e :: Double -> Ptr Double -> CInt
+legendre_Q0_e :: Double -> (Double,Double)
+legendre_Q0_e x = createSFR "legendre_Q0_e" $ gsl_sf_legendre_Q0_e x
+foreign import ccall "legendre.h gsl_sf_legendre_Q0_e" gsl_sf_legendre_Q0_e :: Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_legendre_Q0(double x);
 --
@@ -124,12 +142,12 @@
 legendre_Q0 = gsl_sf_legendre_Q0
 foreign import ccall "legendre.h gsl_sf_legendre_Q0" gsl_sf_legendre_Q0 :: Double -> Double
 
--- | wrapper for int gsl_sf_legendre_Q1_e(double x,double* result);
+-- | wrapper for int gsl_sf_legendre_Q1_e(double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_legendre_Q1_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-legendre_Q1_e :: Double -> Ptr Double -> CInt
-legendre_Q1_e = gsl_sf_legendre_Q1_e
-foreign import ccall "legendre.h gsl_sf_legendre_Q1_e" gsl_sf_legendre_Q1_e :: Double -> Ptr Double -> CInt
+legendre_Q1_e :: Double -> (Double,Double)
+legendre_Q1_e x = createSFR "legendre_Q1_e" $ gsl_sf_legendre_Q1_e x
+foreign import ccall "legendre.h gsl_sf_legendre_Q1_e" gsl_sf_legendre_Q1_e :: Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_legendre_Q1(double x);
 --
@@ -138,12 +156,12 @@
 legendre_Q1 = gsl_sf_legendre_Q1
 foreign import ccall "legendre.h gsl_sf_legendre_Q1" gsl_sf_legendre_Q1 :: Double -> Double
 
--- | wrapper for int gsl_sf_legendre_Ql_e(int l,double x,double* result);
+-- | wrapper for int gsl_sf_legendre_Ql_e(int l,double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_legendre_Ql_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-legendre_Ql_e :: CInt -> Double -> Ptr Double -> CInt
-legendre_Ql_e = gsl_sf_legendre_Ql_e
-foreign import ccall "legendre.h gsl_sf_legendre_Ql_e" gsl_sf_legendre_Ql_e :: CInt -> Double -> Ptr Double -> CInt
+legendre_Ql_e :: CInt -> Double -> (Double,Double)
+legendre_Ql_e l x = createSFR "legendre_Ql_e" $ gsl_sf_legendre_Ql_e l x
+foreign import ccall "legendre.h gsl_sf_legendre_Ql_e" gsl_sf_legendre_Ql_e :: CInt -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_legendre_Ql(int l,double x);
 --
@@ -152,12 +170,12 @@
 legendre_Ql = gsl_sf_legendre_Ql
 foreign import ccall "legendre.h gsl_sf_legendre_Ql" gsl_sf_legendre_Ql :: CInt -> Double -> Double
 
--- | wrapper for int gsl_sf_legendre_Plm_e(int l,int m,double x,double* result);
+-- | wrapper for int gsl_sf_legendre_Plm_e(int l,int m,double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_legendre_Plm_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-legendre_Plm_e :: CInt -> CInt -> Double -> Ptr Double -> CInt
-legendre_Plm_e = gsl_sf_legendre_Plm_e
-foreign import ccall "legendre.h gsl_sf_legendre_Plm_e" gsl_sf_legendre_Plm_e :: CInt -> CInt -> Double -> Ptr Double -> CInt
+legendre_Plm_e :: CInt -> CInt -> Double -> (Double,Double)
+legendre_Plm_e l m x = createSFR "legendre_Plm_e" $ gsl_sf_legendre_Plm_e l m x
+foreign import ccall "legendre.h gsl_sf_legendre_Plm_e" gsl_sf_legendre_Plm_e :: CInt -> CInt -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_legendre_Plm(int l,int m,double x);
 --
@@ -180,12 +198,12 @@
 legendre_Plm_deriv_array = gsl_sf_legendre_Plm_deriv_array
 foreign import ccall "legendre.h gsl_sf_legendre_Plm_deriv_array" gsl_sf_legendre_Plm_deriv_array :: CInt -> CInt -> Double -> Ptr Double -> Ptr Double -> CInt
 
--- | wrapper for int gsl_sf_legendre_sphPlm_e(int l,int m,double x,double* result);
+-- | wrapper for int gsl_sf_legendre_sphPlm_e(int l,int m,double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_legendre_sphPlm_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-legendre_sphPlm_e :: CInt -> CInt -> Double -> Ptr Double -> CInt
-legendre_sphPlm_e = gsl_sf_legendre_sphPlm_e
-foreign import ccall "legendre.h gsl_sf_legendre_sphPlm_e" gsl_sf_legendre_sphPlm_e :: CInt -> CInt -> Double -> Ptr Double -> CInt
+legendre_sphPlm_e :: CInt -> CInt -> Double -> (Double,Double)
+legendre_sphPlm_e l m x = createSFR "legendre_sphPlm_e" $ gsl_sf_legendre_sphPlm_e l m x
+foreign import ccall "legendre.h gsl_sf_legendre_sphPlm_e" gsl_sf_legendre_sphPlm_e :: CInt -> CInt -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_legendre_sphPlm(int l,int m,double x);
 --
@@ -215,12 +233,12 @@
 legendre_array_size = gsl_sf_legendre_array_size
 foreign import ccall "legendre.h gsl_sf_legendre_array_size" gsl_sf_legendre_array_size :: CInt -> CInt -> CInt
 
--- | wrapper for int gsl_sf_conicalP_half_e(double lambda,double x,double* result);
+-- | wrapper for int gsl_sf_conicalP_half_e(double lambda,double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_conicalP_half_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-conicalP_half_e :: Double -> Double -> Ptr Double -> CInt
-conicalP_half_e = gsl_sf_conicalP_half_e
-foreign import ccall "legendre.h gsl_sf_conicalP_half_e" gsl_sf_conicalP_half_e :: Double -> Double -> Ptr Double -> CInt
+conicalP_half_e :: Double -> Double -> (Double,Double)
+conicalP_half_e lambda x = createSFR "conicalP_half_e" $ gsl_sf_conicalP_half_e lambda x
+foreign import ccall "legendre.h gsl_sf_conicalP_half_e" gsl_sf_conicalP_half_e :: Double -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_conicalP_half(double lambda,double x);
 --
@@ -229,12 +247,12 @@
 conicalP_half = gsl_sf_conicalP_half
 foreign import ccall "legendre.h gsl_sf_conicalP_half" gsl_sf_conicalP_half :: Double -> Double -> Double
 
--- | wrapper for int gsl_sf_conicalP_mhalf_e(double lambda,double x,double* result);
+-- | wrapper for int gsl_sf_conicalP_mhalf_e(double lambda,double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_conicalP_mhalf_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-conicalP_mhalf_e :: Double -> Double -> Ptr Double -> CInt
-conicalP_mhalf_e = gsl_sf_conicalP_mhalf_e
-foreign import ccall "legendre.h gsl_sf_conicalP_mhalf_e" gsl_sf_conicalP_mhalf_e :: Double -> Double -> Ptr Double -> CInt
+conicalP_mhalf_e :: Double -> Double -> (Double,Double)
+conicalP_mhalf_e lambda x = createSFR "conicalP_mhalf_e" $ gsl_sf_conicalP_mhalf_e lambda x
+foreign import ccall "legendre.h gsl_sf_conicalP_mhalf_e" gsl_sf_conicalP_mhalf_e :: Double -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_conicalP_mhalf(double lambda,double x);
 --
@@ -243,12 +261,12 @@
 conicalP_mhalf = gsl_sf_conicalP_mhalf
 foreign import ccall "legendre.h gsl_sf_conicalP_mhalf" gsl_sf_conicalP_mhalf :: Double -> Double -> Double
 
--- | wrapper for int gsl_sf_conicalP_0_e(double lambda,double x,double* result);
+-- | wrapper for int gsl_sf_conicalP_0_e(double lambda,double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_conicalP_0_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-conicalP_0_e :: Double -> Double -> Ptr Double -> CInt
-conicalP_0_e = gsl_sf_conicalP_0_e
-foreign import ccall "legendre.h gsl_sf_conicalP_0_e" gsl_sf_conicalP_0_e :: Double -> Double -> Ptr Double -> CInt
+conicalP_0_e :: Double -> Double -> (Double,Double)
+conicalP_0_e lambda x = createSFR "conicalP_0_e" $ gsl_sf_conicalP_0_e lambda x
+foreign import ccall "legendre.h gsl_sf_conicalP_0_e" gsl_sf_conicalP_0_e :: Double -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_conicalP_0(double lambda,double x);
 --
@@ -257,12 +275,12 @@
 conicalP_0 = gsl_sf_conicalP_0
 foreign import ccall "legendre.h gsl_sf_conicalP_0" gsl_sf_conicalP_0 :: Double -> Double -> Double
 
--- | wrapper for int gsl_sf_conicalP_1_e(double lambda,double x,double* result);
+-- | wrapper for int gsl_sf_conicalP_1_e(double lambda,double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_conicalP_1_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-conicalP_1_e :: Double -> Double -> Ptr Double -> CInt
-conicalP_1_e = gsl_sf_conicalP_1_e
-foreign import ccall "legendre.h gsl_sf_conicalP_1_e" gsl_sf_conicalP_1_e :: Double -> Double -> Ptr Double -> CInt
+conicalP_1_e :: Double -> Double -> (Double,Double)
+conicalP_1_e lambda x = createSFR "conicalP_1_e" $ gsl_sf_conicalP_1_e lambda x
+foreign import ccall "legendre.h gsl_sf_conicalP_1_e" gsl_sf_conicalP_1_e :: Double -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_conicalP_1(double lambda,double x);
 --
@@ -271,12 +289,12 @@
 conicalP_1 = gsl_sf_conicalP_1
 foreign import ccall "legendre.h gsl_sf_conicalP_1" gsl_sf_conicalP_1 :: Double -> Double -> Double
 
--- | wrapper for int gsl_sf_conicalP_sph_reg_e(int l,double lambda,double x,double* result);
+-- | wrapper for int gsl_sf_conicalP_sph_reg_e(int l,double lambda,double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_conicalP_sph_reg_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-conicalP_sph_reg_e :: CInt -> Double -> Double -> Ptr Double -> CInt
-conicalP_sph_reg_e = gsl_sf_conicalP_sph_reg_e
-foreign import ccall "legendre.h gsl_sf_conicalP_sph_reg_e" gsl_sf_conicalP_sph_reg_e :: CInt -> Double -> Double -> Ptr Double -> CInt
+conicalP_sph_reg_e :: CInt -> Double -> Double -> (Double,Double)
+conicalP_sph_reg_e l lambda x = createSFR "conicalP_sph_reg_e" $ gsl_sf_conicalP_sph_reg_e l lambda x
+foreign import ccall "legendre.h gsl_sf_conicalP_sph_reg_e" gsl_sf_conicalP_sph_reg_e :: CInt -> Double -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_conicalP_sph_reg(int l,double lambda,double x);
 --
@@ -285,12 +303,12 @@
 conicalP_sph_reg = gsl_sf_conicalP_sph_reg
 foreign import ccall "legendre.h gsl_sf_conicalP_sph_reg" gsl_sf_conicalP_sph_reg :: CInt -> Double -> Double -> Double
 
--- | wrapper for int gsl_sf_conicalP_cyl_reg_e(int m,double lambda,double x,double* result);
+-- | wrapper for int gsl_sf_conicalP_cyl_reg_e(int m,double lambda,double x,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_conicalP_cyl_reg_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-conicalP_cyl_reg_e :: CInt -> Double -> Double -> Ptr Double -> CInt
-conicalP_cyl_reg_e = gsl_sf_conicalP_cyl_reg_e
-foreign import ccall "legendre.h gsl_sf_conicalP_cyl_reg_e" gsl_sf_conicalP_cyl_reg_e :: CInt -> Double -> Double -> Ptr Double -> CInt
+conicalP_cyl_reg_e :: CInt -> Double -> Double -> (Double,Double)
+conicalP_cyl_reg_e m lambda x = createSFR "conicalP_cyl_reg_e" $ gsl_sf_conicalP_cyl_reg_e m lambda x
+foreign import ccall "legendre.h gsl_sf_conicalP_cyl_reg_e" gsl_sf_conicalP_cyl_reg_e :: CInt -> Double -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_conicalP_cyl_reg(int m,double lambda,double x);
 --
@@ -299,12 +317,12 @@
 conicalP_cyl_reg = gsl_sf_conicalP_cyl_reg
 foreign import ccall "legendre.h gsl_sf_conicalP_cyl_reg" gsl_sf_conicalP_cyl_reg :: CInt -> Double -> Double -> Double
 
--- | wrapper for int gsl_sf_legendre_H3d_0_e(double lambda,double eta,double* result);
+-- | wrapper for int gsl_sf_legendre_H3d_0_e(double lambda,double eta,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_legendre_H3d_0_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-legendre_H3d_0_e :: Double -> Double -> Ptr Double -> CInt
-legendre_H3d_0_e = gsl_sf_legendre_H3d_0_e
-foreign import ccall "legendre.h gsl_sf_legendre_H3d_0_e" gsl_sf_legendre_H3d_0_e :: Double -> Double -> Ptr Double -> CInt
+legendre_H3d_0_e :: Double -> Double -> (Double,Double)
+legendre_H3d_0_e lambda eta = createSFR "legendre_H3d_0_e" $ gsl_sf_legendre_H3d_0_e lambda eta
+foreign import ccall "legendre.h gsl_sf_legendre_H3d_0_e" gsl_sf_legendre_H3d_0_e :: Double -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_legendre_H3d_0(double lambda,double eta);
 --
@@ -313,12 +331,12 @@
 legendre_H3d_0 = gsl_sf_legendre_H3d_0
 foreign import ccall "legendre.h gsl_sf_legendre_H3d_0" gsl_sf_legendre_H3d_0 :: Double -> Double -> Double
 
--- | wrapper for int gsl_sf_legendre_H3d_1_e(double lambda,double eta,double* result);
+-- | wrapper for int gsl_sf_legendre_H3d_1_e(double lambda,double eta,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_legendre_H3d_1_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-legendre_H3d_1_e :: Double -> Double -> Ptr Double -> CInt
-legendre_H3d_1_e = gsl_sf_legendre_H3d_1_e
-foreign import ccall "legendre.h gsl_sf_legendre_H3d_1_e" gsl_sf_legendre_H3d_1_e :: Double -> Double -> Ptr Double -> CInt
+legendre_H3d_1_e :: Double -> Double -> (Double,Double)
+legendre_H3d_1_e lambda eta = createSFR "legendre_H3d_1_e" $ gsl_sf_legendre_H3d_1_e lambda eta
+foreign import ccall "legendre.h gsl_sf_legendre_H3d_1_e" gsl_sf_legendre_H3d_1_e :: Double -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_legendre_H3d_1(double lambda,double eta);
 --
@@ -327,12 +345,12 @@
 legendre_H3d_1 = gsl_sf_legendre_H3d_1
 foreign import ccall "legendre.h gsl_sf_legendre_H3d_1" gsl_sf_legendre_H3d_1 :: Double -> Double -> Double
 
--- | wrapper for int gsl_sf_legendre_H3d_e(int l,double lambda,double eta,double* result);
+-- | wrapper for int gsl_sf_legendre_H3d_e(int l,double lambda,double eta,gsl_sf_result* result);
 --
 --   <http://www.google.com/search?q=gsl_sf_legendre_H3d_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>
-legendre_H3d_e :: CInt -> Double -> Double -> Ptr Double -> CInt
-legendre_H3d_e = gsl_sf_legendre_H3d_e
-foreign import ccall "legendre.h gsl_sf_legendre_H3d_e" gsl_sf_legendre_H3d_e :: CInt -> Double -> Double -> Ptr Double -> CInt
+legendre_H3d_e :: CInt -> Double -> Double -> (Double,Double)
+legendre_H3d_e l lambda eta = createSFR "legendre_H3d_e" $ gsl_sf_legendre_H3d_e l lambda eta
+foreign import ccall "legendre.h gsl_sf_legendre_H3d_e" gsl_sf_legendre_H3d_e :: CInt -> Double -> Double -> Ptr () -> IO CInt
 
 -- | wrapper for double gsl_sf_legendre_H3d(int l,double lambda,double eta);
 --
diff --git a/lib/Numeric/GSL/Special/auto.hs b/lib/Numeric/GSL/Special/auto.hs
--- a/lib/Numeric/GSL/Special/auto.hs
+++ b/lib/Numeric/GSL/Special/auto.hs
@@ -1,3 +1,5 @@
+#!/usr/bin/env runhaskell
+
 -- automatic generation of wrappers for simple GSL special functions
 
 import Text.ParserCombinators.Parsec
diff --git a/lib/Numeric/GSL/Special/autoall.sh b/lib/Numeric/GSL/Special/autoall.sh
--- a/lib/Numeric/GSL/Special/autoall.sh
+++ b/lib/Numeric/GSL/Special/autoall.sh
@@ -1,41 +1,40 @@
 #!/bin/bash
 
-#Airy.hs
-#  include Precision (..) in the export list
-#Exp.hs
-#  remove extern inline definition, qualify name
-#Coupling
-#  remove deprecated INCORRECT
-#Trig.hs
-#  qualify names
-#Legendre.hs
-#  remove extern inline
-#Log.hs
-#  remove extern inline, qualify name
+function rep {
+    ./replace.hs "$1" "$2" < $3 > /tmp/tmp-rep
+    cp /tmp/tmp-rep $3
+}
 
-#runhaskell auto airy
-runhaskell auto bessel
-runhaskell auto clausen
-runhaskell auto coulomb
-#runhaskell auto coupling
-runhaskell auto dawson
-runhaskell auto debye
-runhaskell auto dilog
-runhaskell auto elementary
-runhaskell auto ellint
-runhaskell auto erf
-#runhaskell auto exp
-runhaskell auto expint
-runhaskell auto fermi_dirac
-runhaskell auto gamma
-runhaskell auto gegenbauer
-runhaskell auto hyperg
-runhaskell auto laguerre
-runhaskell auto lambert
-#runhaskell auto legendre legendre.h
-#runhaskell auto log
-runhaskell auto pow_int
-runhaskell auto psi
-runhaskell auto synchrotron
-#runhaskell auto trig
-runhaskell auto zeta
+
+./auto.hs airy
+rep ') where' ', Precision(..)\n) where' Airy.hs
+./auto.hs bessel
+./auto.hs clausen
+./auto.hs coulomb
+./auto.hs coupling
+rep ', coupling_6j_INCORRECT_e\n, coupling_6j_INCORRECT\n' '' Coupling.hs
+./auto.hs dawson
+./auto.hs debye
+./auto.hs dilog
+./auto.hs elementary
+./auto.hs ellint
+./auto.hs erf
+./auto.hs exp
+rep ', exp\n' ', Numeric.GSL.Special.Exp.exp\n' Exp.hs
+./auto.hs expint
+./auto.hs fermi_dirac
+./auto.hs gamma
+./auto.hs gegenbauer
+./auto.hs hyperg
+./auto.hs laguerre
+./auto.hs lambert
+./auto.hs legendre gsl_sf_legendre.h
+./auto.hs log
+rep ', log\n' ', Numeric.GSL.Special.Log.log\n' Log.hs
+./auto.hs pow_int
+./auto.hs psi
+./auto.hs synchrotron
+./auto.hs trig
+rep ', sin\n' ', Numeric.GSL.Special.Trig.sin\n' Trig.hs
+rep ', cos\n' ', Numeric.GSL.Special.Trig.cos\n' Trig.hs
+./auto.hs zeta
diff --git a/lib/Numeric/GSL/Special/gsl_sf_exp.h b/lib/Numeric/GSL/Special/gsl_sf_exp.h
deleted file mode 100644
--- a/lib/Numeric/GSL/Special/gsl_sf_exp.h
+++ /dev/null
@@ -1,146 +0,0 @@
-/* specfunc/gsl_sf_exp.h
- * 
- * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
- * 
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or (at
- * your option) any later version.
- * 
- * This program is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * General Public License for more details.
- * 
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- */
-
-/* Author:  G. Jungman */
-
-#ifndef __GSL_SF_EXP_H__
-#define __GSL_SF_EXP_H__
-
-#include <gsl/gsl_sf_result.h>
-#include <gsl/gsl_precision.h>
-
-#undef __BEGIN_DECLS
-#undef __END_DECLS
-#ifdef __cplusplus
-# define __BEGIN_DECLS extern "C" {
-# define __END_DECLS }
-#else
-# define __BEGIN_DECLS /* empty */
-# define __END_DECLS /* empty */
-#endif
-
-__BEGIN_DECLS
-
-/* Provide an exp() function with GSL semantics,
- * i.e. with proper error checking, etc.
- *
- * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
- */
-int gsl_sf_exp_e(const double x, gsl_sf_result * result);
-double gsl_sf_exp(const double x);
-
-
-/* Exp(x)
- *
- * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
- */
-int gsl_sf_exp_e10_e(const double x, gsl_sf_result_e10 * result);
-
-
-/* Exponentiate and multiply by a given factor:  y * Exp(x)
- *
- * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
- */
-int gsl_sf_exp_mult_e(const double x, const double y, gsl_sf_result * result);
-double gsl_sf_exp_mult(const double x, const double y);
-
-
-/* Exponentiate and multiply by a given factor:  y * Exp(x)
- *
- * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
- */
-int gsl_sf_exp_mult_e10_e(const double x, const double y, gsl_sf_result_e10 * result);
-
-
-/* exp(x)-1
- *
- * exceptions: GSL_EOVRFLW
- */
-int gsl_sf_expm1_e(const double x, gsl_sf_result * result);
-double gsl_sf_expm1(const double x);
-
-
-/* (exp(x)-1)/x = 1 + x/2 + x^2/(2*3) + x^3/(2*3*4) + ...
- *
- * exceptions: GSL_EOVRFLW
- */
-int gsl_sf_exprel_e(const double x, gsl_sf_result * result);
-double gsl_sf_exprel(const double x);
-
-
-/* 2(exp(x)-1-x)/x^2 = 1 + x/3 + x^2/(3*4) + x^3/(3*4*5) + ...
- *
- * exceptions: GSL_EOVRFLW
- */
-int gsl_sf_exprel_2_e(double x, gsl_sf_result * result);
-double gsl_sf_exprel_2(const double x);
-
-
-/* Similarly for the N-th generalization of
- * the above. The so-called N-relative exponential
- *
- * exprel_N(x) = N!/x^N (exp(x) - Sum[x^k/k!, {k,0,N-1}])
- *             = 1 + x/(N+1) + x^2/((N+1)(N+2)) + ...
- *             = 1F1(1,1+N,x)
- */
-int gsl_sf_exprel_n_e(const int n, const double x, gsl_sf_result * result);
-double gsl_sf_exprel_n(const int n, const double x);
-
-
-/* Exponentiate a quantity with an associated error.
- */
-int gsl_sf_exp_err_e(const double x, const double dx, gsl_sf_result * result);
-
-/* Exponentiate a quantity with an associated error.
- */
-int gsl_sf_exp_err_e10_e(const double x, const double dx, gsl_sf_result_e10 * result);
-
-
-/* Exponentiate and multiply by a given factor:  y * Exp(x),
- * for quantities with associated errors.
- *
- * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
- */
-int gsl_sf_exp_mult_err_e(const double x, const double dx, const double y, const double dy, gsl_sf_result * result);
-
-
-/* Exponentiate and multiply by a given factor:  y * Exp(x),
- * for quantities with associated errors.
- *
- * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
- */
-int gsl_sf_exp_mult_err_e10_e(const double x, const double dx, const double y, const double dy, gsl_sf_result_e10 * result);
-
-__END_DECLS
-
-
-#ifdef HAVE_INLINE
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_errno.h>
-
-__BEGIN_DECLS
-
-
-
-__END_DECLS
-
-#endif /* HAVE_INLINE */
-
-
-#endif /* __GSL_SF_EXP_H__ */
diff --git a/lib/Numeric/GSL/Special/gsl_sf_legendre.h b/lib/Numeric/GSL/Special/gsl_sf_legendre.h
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/GSL/Special/gsl_sf_legendre.h
@@ -0,0 +1,319 @@
+/* specfunc/gsl_sf_legendre.h
+ * 
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2004 Gerard Jungman
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/* Author:  G. Jungman */
+
+#ifndef __GSL_SF_LEGENDRE_H__
+#define __GSL_SF_LEGENDRE_H__
+
+#include <gsl/gsl_sf_result.h>
+
+#undef __BEGIN_DECLS
+#undef __END_DECLS
+#ifdef __cplusplus
+# define __BEGIN_DECLS extern "C" {
+# define __END_DECLS }
+#else
+# define __BEGIN_DECLS /* empty */
+# define __END_DECLS /* empty */
+#endif
+
+__BEGIN_DECLS
+
+
+/* P_l(x)   l >= 0; |x| <= 1
+ *
+ * exceptions: GSL_EDOM
+ */
+int     gsl_sf_legendre_Pl_e(const int l, const double x, gsl_sf_result * result);
+double  gsl_sf_legendre_Pl(const int l, const double x);
+
+
+/* P_l(x) for l=0,...,lmax; |x| <= 1
+ *
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_legendre_Pl_array(
+  const int lmax, const double x,
+  double * result_array
+  );
+
+
+/* P_l(x) and P_l'(x) for l=0,...,lmax; |x| <= 1
+ *
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_legendre_Pl_deriv_array(
+  const int lmax, const double x,
+  double * result_array,
+  double * result_deriv_array
+  );
+
+
+/* P_l(x), l=1,2,3
+ *
+ * exceptions: none
+ */
+int gsl_sf_legendre_P1_e(double x, gsl_sf_result * result);
+int gsl_sf_legendre_P2_e(double x, gsl_sf_result * result);
+int gsl_sf_legendre_P3_e(double x, gsl_sf_result * result);
+double gsl_sf_legendre_P1(const double x);
+double gsl_sf_legendre_P2(const double x);
+double gsl_sf_legendre_P3(const double x);
+
+
+/* Q_0(x), x > -1, x != 1
+ *
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_legendre_Q0_e(const double x, gsl_sf_result * result);
+double gsl_sf_legendre_Q0(const double x);
+
+
+/* Q_1(x), x > -1, x != 1
+ *
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_legendre_Q1_e(const double x, gsl_sf_result * result);
+double gsl_sf_legendre_Q1(const double x);
+
+
+/* Q_l(x), x > -1, x != 1, l >= 0
+ *
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_legendre_Ql_e(const int l, const double x, gsl_sf_result * result);
+double gsl_sf_legendre_Ql(const int l, const double x);
+
+
+/* P_l^m(x)  m >= 0; l >= m; |x| <= 1.0
+ *
+ * Note that this function grows combinatorially with l.
+ * Therefore we can easily generate an overflow for l larger
+ * than about 150.
+ *
+ * There is no trouble for small m, but when m and l are both large,
+ * then there will be trouble. Rather than allow overflows, these
+ * functions refuse to calculate when they can sense that l and m are
+ * too big.
+ *
+ * If you really want to calculate a spherical harmonic, then DO NOT
+ * use this. Instead use legendre_sphPlm() below, which  uses a similar
+ * recursion, but with the normalized functions.
+ *
+ * exceptions: GSL_EDOM, GSL_EOVRFLW
+ */
+int     gsl_sf_legendre_Plm_e(const int l, const int m, const double x, gsl_sf_result * result);
+double  gsl_sf_legendre_Plm(const int l, const int m, const double x);
+
+
+/* P_l^m(x)  m >= 0; l >= m; |x| <= 1.0
+ * l=|m|,...,lmax
+ *
+ * exceptions: GSL_EDOM, GSL_EOVRFLW
+ */
+int gsl_sf_legendre_Plm_array(
+  const int lmax, const int m, const double x,
+  double * result_array
+  );
+
+
+/* P_l^m(x)  and d(P_l^m(x))/dx;  m >= 0; lmax >= m; |x| <= 1.0
+ * l=|m|,...,lmax
+ *
+ * exceptions: GSL_EDOM, GSL_EOVRFLW
+ */
+int gsl_sf_legendre_Plm_deriv_array(
+  const int lmax, const int m, const double x,
+  double * result_array,
+  double * result_deriv_array
+  );
+
+
+/* P_l^m(x), normalized properly for use in spherical harmonics
+ * m >= 0; l >= m; |x| <= 1.0
+ *
+ * There is no overflow problem, as there is for the
+ * standard normalization of P_l^m(x).
+ *
+ * Specifically, it returns:
+ *
+ *        sqrt((2l+1)/(4pi)) sqrt((l-m)!/(l+m)!) P_l^m(x)
+ *
+ * exceptions: GSL_EDOM
+ */
+int     gsl_sf_legendre_sphPlm_e(const int l, int m, const double x, gsl_sf_result * result);
+double  gsl_sf_legendre_sphPlm(const int l, const int m, const double x);
+
+
+/* sphPlm(l,m,x) values
+ * m >= 0; l >= m; |x| <= 1.0
+ * l=|m|,...,lmax
+ *
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_legendre_sphPlm_array(
+  const int lmax, int m, const double x,
+  double * result_array
+  );
+
+
+/* sphPlm(l,m,x) and d(sphPlm(l,m,x))/dx values
+ * m >= 0; l >= m; |x| <= 1.0
+ * l=|m|,...,lmax
+ *
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_legendre_sphPlm_deriv_array(
+  const int lmax, const int m, const double x,
+  double * result_array,
+  double * result_deriv_array
+  );
+
+
+
+/* size of result_array[] needed for the array versions of Plm
+ * (lmax - m + 1)
+ */
+int gsl_sf_legendre_array_size(const int lmax, const int m);
+
+
+/* Irregular Spherical Conical Function
+ * P^{1/2}_{-1/2 + I lambda}(x)
+ *
+ * x > -1.0
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_conicalP_half_e(const double lambda, const double x, gsl_sf_result * result);
+double gsl_sf_conicalP_half(const double lambda, const double x);
+
+
+/* Regular Spherical Conical Function
+ * P^{-1/2}_{-1/2 + I lambda}(x)
+ *
+ * x > -1.0
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_conicalP_mhalf_e(const double lambda, const double x, gsl_sf_result * result);
+double gsl_sf_conicalP_mhalf(const double lambda, const double x);
+
+
+/* Conical Function
+ * P^{0}_{-1/2 + I lambda}(x)
+ *
+ * x > -1.0
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_conicalP_0_e(const double lambda, const double x, gsl_sf_result * result);
+double gsl_sf_conicalP_0(const double lambda, const double x);
+
+
+/* Conical Function
+ * P^{1}_{-1/2 + I lambda}(x)
+ *
+ * x > -1.0
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_conicalP_1_e(const double lambda, const double x, gsl_sf_result * result);
+double gsl_sf_conicalP_1(const double lambda, const double x);
+
+
+/* Regular Spherical Conical Function
+ * P^{-1/2-l}_{-1/2 + I lambda}(x)
+ *
+ * x > -1.0, l >= -1
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_conicalP_sph_reg_e(const int l, const double lambda, const double x, gsl_sf_result * result);
+double gsl_sf_conicalP_sph_reg(const int l, const double lambda, const double x);
+
+
+/* Regular Cylindrical Conical Function
+ * P^{-m}_{-1/2 + I lambda}(x)
+ *
+ * x > -1.0, m >= -1
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_conicalP_cyl_reg_e(const int m, const double lambda, const double x, gsl_sf_result * result);
+double gsl_sf_conicalP_cyl_reg(const int m, const double lambda, const double x);
+
+
+/* The following spherical functions are specializations
+ * of Legendre functions which give the regular eigenfunctions
+ * of the Laplacian on a 3-dimensional hyperbolic space.
+ * Of particular interest is the flat limit, which is
+ * Flat-Lim := {lambda->Inf, eta->0, lambda*eta fixed}.
+ */
+  
+/* Zeroth radial eigenfunction of the Laplacian on the
+ * 3-dimensional hyperbolic space.
+ *
+ * legendre_H3d_0(lambda,eta) := sin(lambda*eta)/(lambda*sinh(eta))
+ * 
+ * Normalization:
+ * Flat-Lim legendre_H3d_0(lambda,eta) = j_0(lambda*eta)
+ *
+ * eta >= 0.0
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_legendre_H3d_0_e(const double lambda, const double eta, gsl_sf_result * result);
+double gsl_sf_legendre_H3d_0(const double lambda, const double eta);
+
+
+/* First radial eigenfunction of the Laplacian on the
+ * 3-dimensional hyperbolic space.
+ *
+ * legendre_H3d_1(lambda,eta) :=
+ *    1/sqrt(lambda^2 + 1) sin(lam eta)/(lam sinh(eta))
+ *    (coth(eta) - lambda cot(lambda*eta))
+ * 
+ * Normalization:
+ * Flat-Lim legendre_H3d_1(lambda,eta) = j_1(lambda*eta)
+ *
+ * eta >= 0.0
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_legendre_H3d_1_e(const double lambda, const double eta, gsl_sf_result * result);
+double gsl_sf_legendre_H3d_1(const double lambda, const double eta);
+
+
+/* l'th radial eigenfunction of the Laplacian on the
+ * 3-dimensional hyperbolic space.
+ *
+ * Normalization:
+ * Flat-Lim legendre_H3d_l(l,lambda,eta) = j_l(lambda*eta)
+ *
+ * eta >= 0.0, l >= 0
+ * exceptions: GSL_EDOM
+ */
+int gsl_sf_legendre_H3d_e(const int l, const double lambda, const double eta, gsl_sf_result * result);
+double gsl_sf_legendre_H3d(const int l, const double lambda, const double eta);
+
+
+/* Array of H3d(ell),  0 <= ell <= lmax
+ */
+int gsl_sf_legendre_H3d_array(const int lmax, const double lambda, const double eta, double * result_array);
+
+
+
+
+
+__END_DECLS
+
+#endif /* __GSL_SF_LEGENDRE_H__ */
diff --git a/lib/Numeric/GSL/Special/gsl_sf_log.h b/lib/Numeric/GSL/Special/gsl_sf_log.h
deleted file mode 100644
--- a/lib/Numeric/GSL/Special/gsl_sf_log.h
+++ /dev/null
@@ -1,90 +0,0 @@
-/* specfunc/gsl_sf_log.h
- * 
- * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
- * 
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or (at
- * your option) any later version.
- * 
- * This program is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * General Public License for more details.
- * 
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- */
-
-/* Author:  G. Jungman */
-
-#ifndef __GSL_SF_LOG_H__
-#define __GSL_SF_LOG_H__
-
-#include <gsl/gsl_sf_result.h>
-
-#undef __BEGIN_DECLS
-#undef __END_DECLS
-#ifdef __cplusplus
-# define __BEGIN_DECLS extern "C" {
-# define __END_DECLS }
-#else
-# define __BEGIN_DECLS /* empty */
-# define __END_DECLS /* empty */
-#endif
-
-__BEGIN_DECLS
-
-
-/* Provide a logarithm function with GSL semantics.
- *
- * exceptions: GSL_EDOM
- */
-int gsl_sf_log_e(const double x, gsl_sf_result * result);
-double gsl_sf_log(const double x);
-
-
-/* Log(|x|)
- *
- * exceptions: GSL_EDOM
- */
-int gsl_sf_log_abs_e(const double x, gsl_sf_result * result);
-double gsl_sf_log_abs(const double x);
-
-
-/* Complex Logarithm
- *   exp(lnr + I theta) = zr + I zi
- * Returns argument in [-pi,pi].
- *
- * exceptions: GSL_EDOM
- */
-int gsl_sf_complex_log_e(const double zr, const double zi, gsl_sf_result * lnr, gsl_sf_result * theta);
-
-
-/* Log(1 + x)
- *
- * exceptions: GSL_EDOM
- */
-int gsl_sf_log_1plusx_e(const double x, gsl_sf_result * result);
-double gsl_sf_log_1plusx(const double x);
-
-
-/* Log(1 + x) - x
- *
- * exceptions: GSL_EDOM
- */
-int gsl_sf_log_1plusx_mx_e(const double x, gsl_sf_result * result);
-double gsl_sf_log_1plusx_mx(const double x);
-
-
-#ifdef HAVE_INLINE
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_errno.h>
-
-#endif /* HAVE_INLINE */
-
-
-__END_DECLS
-
-#endif /* __GSL_SF_LOG_H__ */
diff --git a/lib/Numeric/GSL/Special/replace.hs b/lib/Numeric/GSL/Special/replace.hs
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/GSL/Special/replace.hs
@@ -0,0 +1,14 @@
+#!/usr/bin/env runhaskell
+
+import Data.List(isPrefixOf)
+import System(getArgs)
+
+rep (c,r) [] = []
+rep (c,r) f@(x:xs) 
+  | c `isPrefixOf` f = r ++ rep (c,r) (drop (length c) f)
+  | otherwise        = x:(rep (c,r) xs)
+
+main = do
+    args <- getArgs
+    let [p',r'] = map (rep ("\\n","\n")) args
+    interact $ rep (p',r')
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -69,7 +69,7 @@
 import Data.Array
 
 -- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
-class (Normed (Matrix t), Linear Matrix t) => Field t where
+class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
     -- | Singular value decomposition using lapack's dgesvd or zgesvd.
     svd         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
     luPacked    :: Matrix t -> (Matrix t, [Int])
@@ -142,7 +142,7 @@
 
 -- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf.
 --
--- If @c = chol m@ then @m == c \<> ctrans c@.
+-- If @c = chol m@ then @m == ctrans c \<> c@.
 chol :: Field t => Matrix t ->  Matrix t
 chol m | m `equal` ctrans m = cholSH m
        | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -1,4 +1,4 @@
-{-# OPTIONS_GHC -fglasgow-exts #-}
+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Linear
@@ -60,18 +60,7 @@
     divide = vectorZipC Div
     equal u v = dim u == dim v && vectorMax (liftVector magnitude (sub u v)) == 0.0
 
-instance Linear Matrix Double where
-    scale x = liftMatrix (scale x)
-    scaleRecip x = liftMatrix (scaleRecip x)
-    addConstant x = liftMatrix (addConstant x)
-    add = liftMatrix2 add
-    sub = liftMatrix2 sub
-    mul = liftMatrix2 mul
-    divide = liftMatrix2 divide
-    equal a b = cols a == cols b && flatten a `equal` flatten b
-
-
-instance Linear Matrix (Complex Double) where
+instance (Linear Vector a, Container Matrix a) => (Linear Matrix a) where
     scale x = liftMatrix (scale x)
     scaleRecip x = liftMatrix (scaleRecip x)
     addConstant x = liftMatrix (addConstant x)
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -0,0 +1,197 @@
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.LinearAlgebra.Tests
+Copyright   :  (c) Alberto Ruiz 2007
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  portable
+
+Some tests.
+
+-}
+
+module Numeric.LinearAlgebra.Tests(
+--  module Numeric.LinearAlgebra.Tests.Instances,
+--  module Numeric.LinearAlgebra.Tests.Properties,
+  qCheck, runTests
+--, runBigTests
+) where
+
+import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra.Tests.Instances
+import Numeric.LinearAlgebra.Tests.Properties
+import Test.QuickCheck
+import Test.HUnit hiding ((~:),test)
+import System.Info
+import Data.List(foldl1')
+import Numeric.GSL hiding (sin,cos,exp,choose)
+
+qCheck n = check defaultConfig {configSize = const n}
+
+utest str b = TestCase $ assertBool str b
+
+feye n = flipud (ident n) :: Matrix Double
+
+detTest1 = det m == 26
+        && det mc == 38 :+ (-3)
+        && det (feye 2) == -1
+    where
+        m = (3><3) 
+            [ 1, 2, 3
+            , 4, 5, 7
+            , 2, 8, 4 :: Double
+            ]
+        mc = (3><3)
+            [ 1, 2, 3
+            , 4, 5, 7
+            , 2, 8, i
+            ]
+
+--------------------------------------------------------------------
+
+polyEval cs x = foldr (\c ac->ac*x+c) 0 cs
+
+polySolveProp p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p))
+
+---------------------------------------------------------------------
+
+quad f a b = fst $ integrateQAGS 1E-9 100 f a b
+
+-- A multiple integral can be easily defined using partial application
+quad2 f a b g1 g2 = quad h a b
+    where h x = quad (f x) (g1 x) (g2 x)
+
+volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) 
+                        0 r (const 0) (\x->sqrt (r*r-x*x))
+
+---------------------------------------------------------------------
+
+besselTest = utest "bessel_J0_e" ( abs (r-expected) < e )
+    where (r,e) = bessel_J0_e 5.0
+          expected = -0.17759677131433830434739701
+
+exponentialTest = utest "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 )
+    where (v,e,err) = exp_e10_e 30.0
+          expected = exp 30.0
+
+---------------------------------------------------------------------
+
+nd1 = (3><3) [ 1/2, 1/4, 1/4
+             , 0/1, 1/2, 1/4
+             , 1/2, 1/4, 1/2 :: Double]
+
+nd2 = (2><2) [1, 0, 1, 1:: Complex Double]
+
+expmTest1 = expm nd1 :~14~: (3><3)
+ [ 1.762110887278176
+ , 0.478085470590435
+ , 0.478085470590435
+ , 0.104719410945666
+ , 1.709751181805343
+ , 0.425725765117601
+ , 0.851451530235203
+ , 0.530445176063267
+ , 1.814470592751009 ]
+
+expmTest2 = expm nd2 :~15~: (2><2)
+ [ 2.718281828459045
+ , 0.000000000000000
+ , 2.718281828459045
+ , 2.718281828459045 ]
+
+
+---------------------------------------------------------------------
+
+rot :: Double -> Matrix Double
+rot a = (3><3) [ c,0,s
+               , 0,1,0
+               ,-s,0,c ]
+    where c = cos a
+          s = sin a
+
+rotTest = fun (10^5) :~12~: rot 5E4
+    where fun n = foldl1' (<>) (map rot angles)
+              where angles = toList $ linspace n (0,1)
+
+
+-- | All tests must pass with a maximum dimension of about 20
+--  (some tests may fail with bigger sizes due to precision loss).
+runTests :: Int  -- ^ maximum dimension
+         -> IO ()
+runTests n = do
+    setErrorHandlerOff
+    let test p = qCheck n p
+    putStrLn "------ lu"
+    test (luProp    . rM)
+    test (luProp    . cM)
+    putStrLn "------ inv"
+    test (invProp   . rSqWC)
+    test (invProp   . cSqWC)
+    putStrLn "------ pinv"
+    test (pinvProp  . rM)
+    if os == "mingw32"
+        then putStrLn "complex pinvTest skipped in this OS"
+        else test (pinvProp  . cM)
+    putStrLn "------ det"
+    test (detProp   . rSqWC)
+    test (detProp   . cSqWC)
+    putStrLn "------ svd"
+    test (svdProp1  . rM)
+    test (svdProp1  . cM)
+    test (svdProp2  . rM)
+    test (svdProp2  . cM)
+    putStrLn "------ eig"
+    test (eigSHProp . rHer)
+    test (eigSHProp . cHer)
+    test (eigProp   . rSq)
+    test (eigProp   . cSq)
+    putStrLn "------ nullSpace"
+    test (nullspaceProp . rM)
+    test (nullspaceProp . cM)
+    putStrLn "------ qr"
+    test (qrProp     . rM)
+    test (qrProp     . cM)
+    putStrLn "------ hess"
+    test (hessProp   . rSq)
+    test (hessProp   . cSq)
+    putStrLn "------ schur"
+    test (schurProp2 . rSq)
+    if os == "mingw32"
+        then putStrLn "complex schur skipped in this OS"
+        else test (schurProp1 . cSq)
+    putStrLn "------ chol"
+    test (cholProp   . rPosDef)
+    test (cholProp   . cPosDef)
+    putStrLn "------ expm"
+    test (expmDiagProp . rSqWC)
+    test (expmDiagProp . cSqWC)
+    putStrLn "------ fft"
+    test (\v -> ifft (fft v) |~| v)
+    putStrLn "------ vector operations"
+    test (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM))
+    test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM))
+    test (\u -> cos u * tan u |~| sin (u::RM))
+    test (\u -> (cos u * tan u) |~| sin (u::CM))
+    putStrLn "------ some unit tests"
+    runTestTT $ TestList
+        [ utest "1E5 rots" rotTest
+        , utest "det1" detTest1
+        , utest "expm1" (expmTest1)
+        , utest "expm2" (expmTest2)
+        , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)
+        , utest "arith2" $ (((1+i) .* ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( (140*i-51).*1 :: CM)
+        , utest "arith3" $ exp (i.*ones(10,10)*pi) + 1 |~| 0
+        , utest "<\\>"   $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3]
+        , utest "gamma" (gamma 5 == 24.0)
+        , besselTest
+        , exponentialTest
+        , utest "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < 1E-8)
+        , utest "polySolve" (polySolveProp [1,2,3,4])
+        ]
+    return ()
+
+-- -- | Some additional tests on big matrices. They take a few minutes.
+-- runBigTests :: IO ()
+-- runBigTests = undefined
diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -0,0 +1,142 @@
+{-# OPTIONS #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.LinearAlgebra.Tests.Instances
+Copyright   :  (c) Alberto Ruiz 2008
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  portable
+
+Arbitrary instances for vectors, matrices.
+
+-}
+
+module Numeric.LinearAlgebra.Tests.Instances(
+    Sq(..),     rSq,cSq,
+    Rot(..),    rRot,cRot,
+    Her(..),    rHer,cHer,
+    WC(..),     rWC,cWC,
+    SqWC(..),   rSqWC, cSqWC,
+    PosDef(..), rPosDef, cPosDef,
+    RM,CM, rM,cM
+) where
+
+import Numeric.LinearAlgebra
+import Test.QuickCheck
+import Control.Monad(replicateM)
+
+instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where
+    arbitrary = do
+        r <- arbitrary
+        i <- arbitrary
+        return (r:+i)
+    coarbitrary = undefined
+
+chooseDim = sized $ \m -> choose (1,max 1 m)
+
+instance (Field a, Arbitrary a) => Arbitrary (Vector a) where 
+   arbitrary = do m <- chooseDim
+                  l <- vector m
+                  return $ fromList l
+   coarbitrary = undefined
+
+instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where 
+    arbitrary = do
+        m <- chooseDim
+        n <- chooseDim
+        l <- vector (m*n)
+        return $ (m><n) l
+    coarbitrary = undefined
+
+-- a square matrix
+newtype (Sq a) = Sq (Matrix a) deriving Show
+instance (Element a, Arbitrary a) => Arbitrary (Sq a) where
+    arbitrary = do
+        n <- chooseDim
+        l <- vector (n*n)
+        return $ Sq $ (n><n) l
+    coarbitrary = undefined
+
+-- a unitary matrix
+newtype (Rot a) = Rot (Matrix a) deriving Show
+instance (Field a, Arbitrary a) => Arbitrary (Rot a) where
+    arbitrary = do
+        Sq m <- arbitrary
+        let (q,_) = qr m
+        return (Rot q)
+    coarbitrary = undefined
+
+-- a complex hermitian or real symmetric matrix
+newtype (Her a) = Her (Matrix a) deriving Show
+instance (Field a, Arbitrary a) => Arbitrary (Her a) where
+    arbitrary = do
+        Sq m <- arbitrary
+        let m' = m/2
+        return $ Her (m' + ctrans m')
+    coarbitrary = undefined
+
+-- a well-conditioned general matrix (the singular values are between 1 and 100)
+newtype (WC a) = WC (Matrix a) deriving Show
+instance (Field a, Arbitrary a) => Arbitrary (WC a) where
+    arbitrary = do
+        m <- arbitrary
+        let (u,_,v) = svd m
+            r = rows m
+            c = cols m
+            n = min r c
+        sv <- replicateM n (choose (1,100))
+        let s = diagRect (fromList sv) r c
+        return $ WC (u <> real s <> trans v)
+    coarbitrary = undefined
+
+-- a well-conditioned square matrix (the singular values are between 1 and 100)
+newtype (SqWC a) = SqWC (Matrix a) deriving Show
+instance (Field a, Arbitrary a) => Arbitrary (SqWC a) where
+    arbitrary = do
+        Sq m <- arbitrary
+        let (u,_,v) = svd m
+            n = rows m
+        sv <- replicateM n (choose (1,100))
+        let s = diag (fromList sv)
+        return $ SqWC (u <> real s <> trans v)
+    coarbitrary = undefined
+
+-- a positive definite square matrix (the eigenvalues are between 0 and 100)
+newtype (PosDef a) = PosDef (Matrix a) deriving Show
+instance (Field a, Arbitrary a) => Arbitrary (PosDef a) where
+    arbitrary = do
+        Her m <- arbitrary
+        let (_,v) = eigSH m
+            n = rows m
+        l <- replicateM n (choose (0,100))
+        let s = diag (fromList l)
+            p = v <> real s <> ctrans v
+        return $ PosDef (0.5 .* p + 0.5 .* ctrans p)
+    coarbitrary = undefined
+
+type RM = Matrix Double
+type CM = Matrix (Complex Double)
+
+rM m = m :: RM
+cM m = m :: CM
+
+rHer (Her m) = m :: RM
+cHer (Her m) = m :: CM
+
+rRot (Rot m) = m :: RM
+cRot (Rot m) = m :: CM
+
+rSq  (Sq m)  = m :: RM
+cSq  (Sq m)  = m :: CM
+
+rWC (WC m) = m :: RM
+cWC (WC m) = m :: CM
+
+rSqWC (SqWC m) = m :: RM
+cSqWC (SqWC m) = m :: CM
+
+rPosDef (PosDef m) = m :: RM
+cPosDef (PosDef m) = m :: CM
+
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -0,0 +1,149 @@
+{-# OPTIONS #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.LinearAlgebra.Tests.Properties
+Copyright   :  (c) Alberto Ruiz 2008
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  portable
+
+Testing properties.
+
+-}
+
+module Numeric.LinearAlgebra.Tests.Properties (
+    dist, (|~|), (~:), Aprox((:~)),
+    zeros, ones,
+    square,
+    unitary,
+    hermitian,
+    wellCond,
+    positiveDefinite,
+    upperTriang,
+    upperHessenberg,
+    luProp,
+    invProp,
+    pinvProp,
+    detProp,
+    nullspaceProp,
+    svdProp1, svdProp2,
+    eigProp, eigSHProp,
+    qrProp,
+    hessProp,
+    schurProp1, schurProp2,
+    cholProp,
+    expmDiagProp
+) where
+
+import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra.Tests.Instances(Sq(..),Her(..),Rot(..))
+import Test.QuickCheck
+
+-- relative error
+dist :: (Normed t, Num t) => t -> t -> Double
+dist a b = r
+    where norm = pnorm Infinity
+          na = norm a
+          nb = norm b
+          nab = norm (a-b)
+          mx = max na nb
+          mn = min na nb
+          r = if mn < eps
+                then mx
+                else nab/mx
+
+infixl 4 |~|
+a |~| b = a :~10~: b
+--a |~| b = dist a b < 10^^(-10)
+
+data Aprox a = (:~) a Int
+(~:) :: (Normed a, Num a) => Aprox a -> a -> Bool
+a :~n~: b = dist a b < 10^^(-n)
+
+------------------------------------------------------
+
+square m = rows m == cols m
+
+unitary m = square m && m <> ctrans m |~| ident (rows m)
+
+hermitian m = square m && m |~| ctrans m
+
+wellCond m = rcond m > 1/100
+
+positiveDefinite m = minimum (toList e) > 0
+    where (e,v) = eigSH m
+
+upperTriang m = rows m == 1 || down == z
+    where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
+          z = constant 0 (dim down)
+
+upperHessenberg m = rows m < 3 || down == z
+    where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m))
+          z = constant 0 (dim down)
+
+zeros (r,c) = reshape c (constant 0 (r*c))
+
+ones (r,c) = zeros (r,c) + 1
+
+-----------------------------------------------------
+
+luProp m = m |~| p <> l <> u && det p == s
+    where (l,u,p,s) = lu m
+
+invProp m = m <> inv m |~| ident (rows m)
+
+pinvProp m =  m <> p <> m |~| m
+           && p <> m <> p |~| p
+           && hermitian (m<>p)
+           && hermitian (p<>m)
+    where p = pinv m
+
+detProp m = s d1 |~| s d2
+    where d1 = det m
+          d2 = det' m * det q
+          det' m = product $ toList $ takeDiag r
+          (q,r) = qr m
+          s x = fromList [x]
+
+nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c))
+    where nl = nullspacePrec 1 m
+          n = fromColumns nl
+          r = rows m
+          c = cols m - rank m
+
+svdProp1 m = u <> real d <> trans v |~| m
+          && unitary u && unitary v
+    where (u,d,v) = full svd m
+
+svdProp2 m = (m |~| 0) `trivial` ((m |~| 0) || u <> real (diag s) <> trans v |~| m)
+    where (u,s,v) = economy svd m
+
+eigProp m = complex m <> v |~| v <> diag s
+    where (s, v) = eig m
+
+eigSHProp m = m <> v |~| v <> real (diag s)
+              && unitary v
+              && m |~| v <> real (diag s) <> ctrans v
+    where (s, v) = eigSH m
+
+qrProp m = q <> r |~| m && unitary q && upperTriang r
+    where (q,r) = qr m
+
+hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h
+    where (p,h) = hess m
+
+schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s
+    where (u,s) = schur m
+
+schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme
+    where (u,s) = schur m
+
+cholProp m = m |~| ctrans c <> c && upperTriang c
+    where c = chol m
+          pos = positiveDefinite m
+
+expmDiagProp m = expm (logm m) :~ 7 ~: complex m
+    where logm m = matFunc log m
+