diff --git a/CHANGELOG b/CHANGELOG
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -1,3 +1,16 @@
+0.18.0.0
+--------
+
+    * Many new functions and instances in the Static module
+    
+    * meanCov and gaussianSample use Herm type
+
+    * thinQR, thinRQ
+    
+    * compactSVDTol
+
+    * unitary changed to normalize, also admits Vector (Complex Double)
+  
 0.17.0.0
 --------
 
diff --git a/THANKS.md b/THANKS.md
--- a/THANKS.md
+++ b/THANKS.md
@@ -192,7 +192,8 @@
 
 - Matt Peddie wrote the interfaces to the interpolation and simulated annealing modules.
 
-- "maxc01" solved uninstallability in FreeBSD and improved urandom
+- "maxc01" solved uninstallability in FreeBSD, improved urandom, and fixed a Windows
+  link error using rand_s.
 
 - "ntfrgl" added {take,drop}Last{Rows,Columns} and odeSolveVWith with generalized step control function
    and fixed link errors related to mod/mod_l.
@@ -203,4 +204,27 @@
 
 - Vassil Keremidchiev fixed the cabal options for OpenBlas, fixed several installation
   issues, and added support for stack-based build.
+
+- Greg Nwosu fixed arm compilation
+
+- Patrik Jansson changed meanCov and gaussianSample to use Herm type. Fixed stack.yaml.
+
+- Justin Le added NFData instances for Static types, added mapping and outer product
+  methods to Domain, and many other functions to the Static module.
+
+- Sidharth Kapur added Normed and numeric instances for several Static types,
+fixed the CPP issue in cabal files, and made many other contributions.
+
+- Matt Renaud improved the documentation.
+
+- Joshua Moerman fixed cabal/stack flags for windows.
+
+- Francesco Mazzoli, Niklas Hambüchen, Patrick Chilton, and Andras Slemmer
+  discovered a serious and subtle bug in the wrapper helpers causing memory corruption.
+  Andras Slemmer fixed the bug. Thank you all.
+
+- Kevin Slagle implemented thinQR and thinRQ, much faster than the original qr,
+  and added compactSVDTol.
+
+- "fedeinthemix" suggested a better name and a more general type for unitary.
 
diff --git a/hmatrix.cabal b/hmatrix.cabal
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.17.0.2
+Version:            0.18.0.0
 License:            BSD3
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -55,7 +55,6 @@
                         Internal.Devel
                         Internal.Vectorized
                         Internal.Matrix
-                        Internal.Foreign
                         Internal.ST
                         Internal.IO
                         Internal.Element
@@ -79,8 +78,7 @@
                         src/Internal/C/vector-aux.c
 
 
-    extensions:         ForeignFunctionInterface,
-                        CPP
+    extensions:         ForeignFunctionInterface
 
     ghc-options:        -Wall
                         -fno-warn-missing-signatures
@@ -94,7 +92,6 @@
     if arch(i386)
         cc-options:     -msse2
 
-    cpp-options:        -DBINARY
 
     if os(OSX)
         if flag(openblas)
@@ -124,7 +121,7 @@
 
     if os(windows)
         if flag(openblas)
-            extra-libraries:    libopenblas
+            extra-libraries:    libopenblas, libgcc_s_seh-1, libgfortran-3, libquadmath-0
         else
             extra-libraries:    blas lapack
 
diff --git a/src/Internal/Algorithms.hs b/src/Internal/Algorithms.hs
--- a/src/Internal/Algorithms.hs
+++ b/src/Internal/Algorithms.hs
@@ -30,6 +30,7 @@
 import Internal.Numeric
 import Data.List(foldl1')
 import qualified Data.Array as A
+import qualified Data.Vector.Storable as Vector
 import Internal.ST
 import Internal.Vectorized(range)
 import Control.DeepSeq
@@ -291,9 +292,13 @@
 
 -}
 compactSVD :: Field t  => Matrix t -> (Matrix t, Vector Double, Matrix t)
-compactSVD m = (u', subVector 0 d s, v') where
+compactSVD = compactSVDTol 1
+
+-- | @compactSVDTol r@ is similar to 'compactSVD' (for which @r=1@), but uses tolerance @tol=r*g*eps*(max rows cols)@ to distinguish nonzero singular values, where @g@ is the greatest singular value. If @g<r*eps@, then only one singular value is returned.
+compactSVDTol :: Field t  => Double -> Matrix t -> (Matrix t, Vector Double, Matrix t)
+compactSVDTol r m = (u', subVector 0 d s, v') where
     (u,s,v) = thinSVD m
-    d = rankSVD (1*eps) m s `max` 1
+    d = rankSVD (r*eps) m s `max` 1
     u' = takeColumns d u
     v' = takeColumns d v
 
@@ -475,9 +480,14 @@
 -- | QR factorization.
 --
 -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
+-- Note: the current implementation is very slow for large matrices. 'thinQR' is much faster.
 qr :: Field t => Matrix t -> (Matrix t, Matrix t)
 qr = {-# SCC "qr" #-} unpackQR . qr'
 
+-- | A version of 'qr' which returns only the @min (rows m) (cols m)@ columns of @q@ and rows of @r@.
+thinQR :: Field t => Matrix t -> (Matrix t, Matrix t)
+thinQR = {-# SCC "thinQR" #-} thinUnpackQR . qr'
+
 -- | Compute the QR decomposition of a matrix in compact form.
 qrRaw :: Field t => Matrix t -> QR t
 qrRaw m = QR x v
@@ -494,9 +504,17 @@
 -- | RQ factorization.
 --
 -- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.
+-- Note: the current implementation is very slow for large matrices. 'thinRQ' is much faster.
 rq :: Field t => Matrix t -> (Matrix t, Matrix t)
-rq m =  {-# SCC "rq" #-} (r,q) where
-    (q',r') = qr $ trans $ rev1 m
+rq = {-# SCC "rq" #-} rqFromQR qr
+
+-- | A version of 'rq' which returns only the @min (rows m) (cols m)@ columns of @r@ and rows of @q@.
+thinRQ :: Field t => Matrix t -> (Matrix t, Matrix t)
+thinRQ = {-# SCC "thinQR" #-} rqFromQR thinQR
+
+rqFromQR :: Field t => (Matrix t -> (Matrix t, Matrix t)) -> Matrix t -> (Matrix t, Matrix t)
+rqFromQR qr0 m = (r,q) where
+    (q',r') = qr0 $ trans $ rev1 m
     r = rev2 (trans r')
     q = rev2 (trans q')
     rev1 = flipud . fliprl
@@ -723,6 +741,12 @@
           vs = zipWith zh [1..mn] cs
           hs = zipWith haussholder (toList tau) vs
           q = foldl1' mXm hs
+
+thinUnpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
+thinUnpackQR (pq, tau) = (q, r)
+  where mn = uncurry min $ size pq
+        q = qrgr mn $ QR pq tau
+        r = fromRows $ zipWith (\i v -> Vector.replicate i 0 Vector.++ Vector.drop i v) [0..mn-1] (toRows pq)
 
 unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t)
 unpackHess hf m
diff --git a/src/Internal/C/vector-aux.c b/src/Internal/C/vector-aux.c
--- a/src/Internal/C/vector-aux.c
+++ b/src/Internal/C/vector-aux.c
@@ -945,6 +945,8 @@
 /* Windows use thread-safe random
    See: http://stackoverflow.com/questions/143108/is-windows-rand-s-thread-safe
 */
+#if defined (__APPLE__) || (__FreeBSD__)
+
 inline double urandom() {
     /* the probalility of matching will be theoretically p^3(in fact, it is not)
        p is matching probalility of random().
@@ -959,6 +961,22 @@
     return (double)nrand48(state) / (double)max_random;
 }
 
+#else
+
+#define _CRT_RAND_S
+inline double urandom() {
+    unsigned int number;
+    errno_t err;
+    err = rand_s(&number);
+    if (err!=0) {
+        printf("something wrong\n");
+        return -1;
+    }
+    return (double)number / (double)UINT_MAX;
+}
+
+#endif
+
 double gaussrand(int *phase, double *pV1, double *pV2, double *pS)
 {
 	double V1=*pV1, V2=*pV2, S=*pS;
@@ -985,10 +1003,39 @@
 
 }
 
+#if defined(_WIN32) || defined(WIN32)
+
 int random_vector(unsigned int seed, int code, DVEC(r)) {
     int phase = 0;
     double V1,V2,S;
 
+    srand(seed);
+
+    int k;
+    switch (code) {
+      case 0: { // uniform
+        for (k=0; k<rn; k++) {
+            rp[k] = urandom();
+        }
+        OK
+      }
+      case 1: { // gaussian
+        for (k=0; k<rn; k++) {
+            rp[k] = gaussrand(&phase,&V1,&V2,&S);
+        }
+        OK
+      }
+
+      default: ERROR(BAD_CODE);
+    }
+}
+
+#else
+
+int random_vector(unsigned int seed, int code, DVEC(r)) {
+    int phase = 0;
+    double V1,V2,S;
+
     srandom(seed);
 
     int k;
@@ -1009,6 +1056,8 @@
       default: ERROR(BAD_CODE);
     }
 }
+
+#endif
 
 #else
 
diff --git a/src/Internal/Container.hs b/src/Internal/Container.hs
--- a/src/Internal/Container.hs
+++ b/src/Internal/Container.hs
@@ -28,7 +28,7 @@
 import Internal.Matrix
 import Internal.Element
 import Internal.Numeric
-import Internal.Algorithms(Field,linearSolveSVD)
+import Internal.Algorithms(Field,linearSolveSVD,Herm,mTm)
 
 ------------------------------------------------------------------
 
@@ -206,14 +206,14 @@
  , -1.0127225830525157e-2,     3.0373954915729318 ])
 
 -}
-meanCov :: Matrix Double -> (Vector Double, Matrix Double)
+meanCov :: Matrix Double -> (Vector Double, Herm Double)
 meanCov x = (med,cov) where
     r    = rows x
     k    = 1 / fromIntegral r
     med  = konst k r `vXm` x
     meds = konst 1 r `outer` med
     xc   = x `sub` meds
-    cov  = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc)
+    cov  = scale (recip (fromIntegral (r-1))) (mTm xc)
 
 --------------------------------------------------------------------------------
 
@@ -293,5 +293,3 @@
     | otherwise = error $ "out of range index in remap"
   where
     [i',j'] = conformMs [i,j]
-    
-
diff --git a/src/Internal/Devel.hs b/src/Internal/Devel.hs
--- a/src/Internal/Devel.hs
+++ b/src/Internal/Devel.hs
@@ -80,8 +80,8 @@
   where
     type Trans c b
     type TransRaw c b
-    apply      :: (Trans c b) -> c -> b
-    applyRaw   :: (TransRaw c b) -> c -> b
+    apply      :: c -> (b -> IO r) -> (Trans c b) -> IO r
+    applyRaw   :: c -> (b -> IO r) -> (TransRaw c b) -> IO r
     infixl 1 `apply`, `applyRaw`
 
 instance Storable t => TransArray (Vector t)
@@ -92,4 +92,3 @@
     {-# INLINE apply #-}
     applyRaw = avec
     {-# INLINE applyRaw #-}
-
diff --git a/src/Internal/Element.hs b/src/Internal/Element.hs
--- a/src/Internal/Element.hs
+++ b/src/Internal/Element.hs
@@ -3,7 +3,6 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE CPP #-}
 
 -----------------------------------------------------------------------------
 -- |
@@ -35,7 +34,6 @@
 
 -------------------------------------------------------------------
 
-#ifdef BINARY
 
 import Data.Binary
 
@@ -48,7 +46,6 @@
           v <- get
           return (reshape c v)
 
-#endif
 
 -------------------------------------------------------------------
 
diff --git a/src/Internal/Foreign.hs b/src/Internal/Foreign.hs
deleted file mode 100644
--- a/src/Internal/Foreign.hs
+++ /dev/null
@@ -1,102 +0,0 @@
-{-# LANGUAGE MagicHash, UnboxedTuples #-}
--- | FFI and hmatrix helpers.
---
--- Sample usage, to upload a perspective matrix to a shader.
---
--- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3) 
--- @
---
-
-module Internal.Foreign 
-    ( app
-    , appVector, appVectorLen
-    , appMatrix, appMatrixLen, appMatrixRaw, appMatrixRawLen
-    , unsafeMatrixToVector, unsafeMatrixToForeignPtr
-    ) where
-
-import Foreign.C.Types(CInt)
-import Internal.Vector
-import Internal.Matrix
-import qualified Data.Vector.Storable as S
-import Foreign (Ptr, ForeignPtr, Storable)
-import GHC.Base (IO(..), realWorld#)
-
-{-# INLINE unsafeInlinePerformIO #-}
--- | If we use unsafePerformIO, it may not get inlined, so in a function that returns IO (which are all safe uses of app* in this module), there would be
--- unecessary calls to unsafePerformIO or its internals.
-unsafeInlinePerformIO :: IO a -> a
-unsafeInlinePerformIO (IO f) = case f realWorld# of
-    (# _, x #) -> x
-
-{-# INLINE app #-}
--- | Only useful since it is left associated with a precedence of 1, unlike 'Prelude.$', which is right associative.
--- e.g.
---
--- @
--- someFunction
---     \`appMatrixLen\` m
---     \`appVectorLen\` v
---     \`app\` other
---     \`app\` arguments
---     \`app\` go here
--- @
---
--- One could also write:
---
--- @
--- (someFunction 
---     \`appMatrixLen\` m
---     \`appVectorLen\` v) 
---     other 
---     arguments 
---     (go here)
--- @
---
-app :: (a -> b) -> a -> b
-app f = f
-
-{-# INLINE appVector #-}
-appVector :: Storable a => (Ptr a -> b) -> Vector a -> b
-appVector f x = unsafeInlinePerformIO (S.unsafeWith x (return . f))
-
-{-# INLINE appVectorLen #-}
-appVectorLen :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b
-appVectorLen f x = unsafeInlinePerformIO (S.unsafeWith x (return . f (fromIntegral (S.length x))))
-
-{-# INLINE appMatrix #-}
-appMatrix :: Element a => (Ptr a -> b) -> Matrix a -> b
-appMatrix f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f))
-
-{-# INLINE appMatrixLen #-}
-appMatrixLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
-appMatrixLen f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f r c))
-  where
-    r = fromIntegral (rows x)
-    c = fromIntegral (cols x)
-
-{-# INLINE appMatrixRaw #-}
-appMatrixRaw :: Storable a => (Ptr a -> b) -> Matrix a -> b
-appMatrixRaw f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f))
-
-{-# INLINE appMatrixRawLen #-}
-appMatrixRawLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
-appMatrixRawLen f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f r c))
-  where
-    r = fromIntegral (rows x)
-    c = fromIntegral (cols x)
-
-infixl 1 `app`
-infixl 1 `appVector`
-infixl 1 `appMatrix`
-infixl 1 `appMatrixRaw`
-
-{-# INLINE unsafeMatrixToVector #-}
--- | This will disregard the order of the matrix, and simply return it as-is. 
--- If the order of the matrix is RowMajor, this function is identical to 'flatten'.
-unsafeMatrixToVector :: Matrix a -> Vector a
-unsafeMatrixToVector = xdat
-
-{-# INLINE unsafeMatrixToForeignPtr #-}
-unsafeMatrixToForeignPtr :: Storable a => Matrix a -> (ForeignPtr a, Int)
-unsafeMatrixToForeignPtr m = S.unsafeToForeignPtr0 (xdat m)
-
diff --git a/src/Internal/LAPACK.hs b/src/Internal/LAPACK.hs
--- a/src/Internal/LAPACK.hs
+++ b/src/Internal/LAPACK.hs
@@ -18,7 +18,7 @@
 
 import Internal.Devel
 import Internal.Vector
-import Internal.Matrix hiding ((#))
+import Internal.Matrix hiding ((#), (#!))
 import Internal.Conversion
 import Internal.Element
 import Foreign.Ptr(nullPtr)
@@ -28,10 +28,13 @@
 
 -----------------------------------------------------------------------------------
 
-infixl 1 #
+infixr 1 #
 a # b = apply a b
 {-# INLINE (#) #-}
 
+a #! b = a # b # id
+{-# INLINE (#!) #-}
+
 -----------------------------------------------------------------------------------
 
 type TMMM t = t ::> t ::> t ::> Ok
@@ -56,7 +59,7 @@
     when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++
                                        show (rows a,cols a) ++ " x " ++ show (rows b, cols b)
     s <- createMatrix ColumnMajor (rows a) (cols b)
-    f (isT a) (isT b) # (tt a) # (tt b) # s #| st
+    ((tt a) # (tt b) #! s) (f (isT a) (isT b)) #| st
     return s
 
 -- | Matrix product based on BLAS's /dgemm/.
@@ -80,7 +83,7 @@
     when (cols a /= rows b) $ error $
         "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b
     s <- createMatrix ColumnMajor (rows a) (cols b)
-    c_multiplyI m # a # b # s #|"c_multiplyI"
+    (a # b #! s) (c_multiplyI m) #|"c_multiplyI"
     return s
 
 multiplyL :: Z -> Matrix Z -> Matrix Z -> Matrix Z
@@ -88,7 +91,7 @@
     when (cols a /= rows b) $ error $
         "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b
     s <- createMatrix ColumnMajor (rows a) (cols b)
-    c_multiplyL m # a # b # s #|"c_multiplyL"
+    (a # b #! s) (c_multiplyL m) #|"c_multiplyL"
     return s
 
 -----------------------------------------------------------------------------
@@ -121,7 +124,7 @@
     u <- createMatrix ColumnMajor r r
     s <- createVector (min r c)
     v <- createMatrix ColumnMajor c c
-    f # a # u # s # v #| st
+    (a # u # s #! v) f #| st
     return (u,s,v)
   where
     r = rows x
@@ -149,7 +152,7 @@
     u <- createMatrix ColumnMajor r q
     s <- createVector q
     v <- createMatrix ColumnMajor q c
-    f # a # u # s # v #| st
+    (a # u # s #! v) f #| st
     return (u,s,v)
   where
     r = rows x
@@ -176,7 +179,7 @@
 svAux f st x = unsafePerformIO $ do
     a <- copy ColumnMajor x
     s <- createVector q
-    g # a # s #| st
+    (a #! s) g #| st
     return s
   where
     r = rows x
@@ -197,7 +200,7 @@
     a <- copy ColumnMajor x
     s <- createVector q
     v <- createMatrix ColumnMajor c c
-    g # a # s # v #| st
+    (a # s #! v) g #| st
     return (s,v)
   where
     r = rows x
@@ -218,7 +221,7 @@
     a <- copy ColumnMajor x
     u <- createMatrix ColumnMajor r r
     s <- createVector q
-    g # a # u # s #| st
+    (a # u #! s) g #| st
     return (u,s)
   where
     r = rows x
@@ -237,7 +240,7 @@
     a <- copy ColumnMajor m
     l <- createVector r
     v <- createMatrix ColumnMajor r r
-    g # a # l # v #| st
+    (a # l #! v) g #| st
     return (l,v)
   where
     r = rows m
@@ -252,7 +255,7 @@
 eigOnlyAux f st m = unsafePerformIO $ do
     a <- copy ColumnMajor m
     l <- createVector r
-    g # a # l #| st
+    (a #! l) g #| st
     return l
   where
     r = rows m
@@ -277,7 +280,7 @@
     a <- copy ColumnMajor m
     l <- createVector r
     v <- createMatrix ColumnMajor r r
-    g # a # l # v #| "eigR"
+    (a # l #! v) g #| "eigR"
     return (l,v)
   where
     r = rows m
@@ -305,7 +308,7 @@
 eigSHAux f st m = unsafePerformIO $ do
     l <- createVector r
     v <- copy ColumnMajor m
-    f # l # v #| st
+    (l #! v) f #| st
     return (l,v)
   where
     r = rows m
@@ -356,7 +359,7 @@
     | n1==n2 && n1==r = unsafePerformIO . g $ do
         a' <- copy ColumnMajor a
         s  <- copy ColumnMajor b
-        f # a' # s #| st
+        (a' #! s) f #| st
         return s
     | otherwise = error $ st ++ " of nonsquare matrix"
   where
@@ -387,7 +390,7 @@
 linearSolveSQAux2 g f st a b
     | n1==n2 && n1==r = unsafePerformIO . g $ do
         s <- copy ColumnMajor b
-        f # a # s #| st
+        (a #! s) f #| st
         return s
     | otherwise = error $ st ++ " of nonsquare matrix"
   where
@@ -415,7 +418,7 @@
         a' <- copy ColumnMajor a
         r  <- createMatrix ColumnMajor (max m n) nrhs
         setRect 0 0 b r
-        f # a' # r #| st
+        (a' #! r) f #| st
         return r
     | otherwise = error $ "different number of rows in linearSolve ("++st++")"
   where
@@ -458,7 +461,7 @@
 
 cholAux f st a = do
     r <- copy ColumnMajor a
-    f # r #| st
+    (r # id) f #| st
     return r
 
 -- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.
@@ -495,7 +498,7 @@
 qrAux f st a = unsafePerformIO $ do
     r <- copy ColumnMajor a
     tau <- createVector mn
-    f # tau # r #| st
+    (tau #! r) f #| st
     return (r,tau)
   where
     m = rows a
@@ -514,7 +517,7 @@
 
 qrgrAux f st n (a, tau) = unsafePerformIO $ do
     res <- copy ColumnMajor (subMatrix (0,0) (rows a,n) a)
-    f # (subVector 0 n tau') # res #| st
+    ((subVector 0 n tau') #! res) f #| st
     return res
   where
     tau' = vjoin [tau, constantD 0 n]
@@ -534,7 +537,7 @@
 hessAux f st a = unsafePerformIO $ do
     r <- copy ColumnMajor a
     tau <- createVector (mn-1)
-    f # tau # r #| st
+    (tau #! r) f #| st
     return (r,tau)
   where
     m = rows a
@@ -556,7 +559,7 @@
 schurAux f st a = unsafePerformIO $ do
     u <- createMatrix ColumnMajor n n
     s <- copy ColumnMajor a
-    f # u # s #| st
+    (u #! s) f #| st
     return (u,s)
   where
     n = rows a
@@ -576,7 +579,7 @@
 luAux f st a = unsafePerformIO $ do
     lu <- copy ColumnMajor a
     piv <- createVector (min n m)
-    f # piv # lu #| st
+    (piv #! lu) f #| st
     return (lu, map (pred.round) (toList piv))
   where
     n = rows a
@@ -598,7 +601,7 @@
 lusAux f st a piv b
     | n1==n2 && n2==n =unsafePerformIO $ do
          x <- copy ColumnMajor b
-         f # a # piv' # x #| st
+         (a # piv' #! x) f #| st
          return x
     | otherwise = error st
   where
@@ -622,7 +625,7 @@
 ldlAux f st a = unsafePerformIO $ do
     ldl <- copy ColumnMajor a
     piv <- createVector (rows a)
-    f # piv # ldl #| st
+    (piv #! ldl) f #| st
     return (ldl, map (pred.round) (toList piv))
 
 -----------------------------------------------------------------------------------
@@ -637,4 +640,3 @@
 -- | Solve a complex linear system from a precomputed LDL decomposition ('ldlC'), using LAPACK's /zsytrs/.
 ldlsC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
 ldlsC a piv b = lusAux zsytrs "ldlsC" (fmat a) piv b
-
diff --git a/src/Internal/Matrix.hs b/src/Internal/Matrix.hs
--- a/src/Internal/Matrix.hs
+++ b/src/Internal/Matrix.hs
@@ -5,10 +5,9 @@
 {-# LANGUAGE TypeOperators            #-}
 {-# LANGUAGE TypeFamilies             #-}
 {-# LANGUAGE ViewPatterns             #-}
+{-# LANGUAGE DeriveGeneric            #-}
 {-# LANGUAGE ConstrainedClassMethods  #-}
 
-
-
 -- |
 -- Module      :  Internal.Matrix
 -- Copyright   :  (c) Alberto Ruiz 2007-15
@@ -23,7 +22,7 @@
 
 import Internal.Vector
 import Internal.Devel
-import Internal.Vectorized hiding ((#))
+import Internal.Vectorized hiding ((#), (#!))
 import Foreign.Marshal.Alloc ( free )
 import Foreign.Marshal.Array(newArray)
 import Foreign.Ptr ( Ptr )
@@ -111,15 +110,15 @@
 
 -- C-Haskell matrix adapters
 {-# INLINE amatr #-}
-amatr :: Storable a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
-amatr f x = inlinePerformIO (unsafeWith (xdat x) (return . f r c))
+amatr :: Storable a => Matrix a -> (f -> IO r) -> (CInt -> CInt -> Ptr a -> f) -> IO r
+amatr x f g = unsafeWith (xdat x) (f . g r c)
   where
     r  = fi (rows x)
     c  = fi (cols x)
 
 {-# INLINE amat #-}
-amat :: Storable a => (CInt -> CInt -> CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
-amat f x = inlinePerformIO (unsafeWith (xdat x) (return . f r c sr sc))
+amat :: Storable a => Matrix a -> (f -> IO r) -> (CInt -> CInt -> CInt -> CInt -> Ptr a -> f) -> IO r
+amat x f g = unsafeWith (xdat x) (f . g r c sr sc)
   where
     r  = fi (rows x)
     c  = fi (cols x)
@@ -136,10 +135,13 @@
     applyRaw = amatr
     {-# INLINE applyRaw #-}
 
-infixl 1 #
+infixr 1 #
 a # b = apply a b
 {-# INLINE (#) #-}
 
+a #! b = a # b # id
+{-# INLINE (#!) #-}
+
 --------------------------------------------------------------------------------
 
 copy ord m = extractR ord m 0 (idxs[0,rows m-1]) 0 (idxs[0,cols m-1])
@@ -427,7 +429,8 @@
     let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr
         nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc
     r <- createMatrix ord nr nc
-    f moder modec # vr # vc # m # r  #|"extract"
+    (vr # vc # m #! r) (f moder modec)  #|"extract"
+
     return r
 
 type Extr x = CInt -> CInt -> CIdxs (CIdxs (OM x (OM x (IO CInt))))
@@ -441,7 +444,7 @@
 
 ---------------------------------------------------------------
 
-setRectAux f i j m r = f (fi i) (fi j) # m # r #|"setRect"
+setRectAux f i j m r = (m #! r) (f (fi i) (fi j)) #|"setRect"
 
 type SetRect x = I -> I -> x ::> x::> Ok
 
@@ -456,7 +459,7 @@
 
 sortG f v = unsafePerformIO $ do
     r <- createVector (dim v)
-    f # v # r #|"sortG"
+    (v #! r) f #|"sortG"
     return r
 
 sortIdxD = sortG c_sort_indexD
@@ -483,7 +486,7 @@
 
 compareG f u v = unsafePerformIO $ do
     r <- createVector (dim v)
-    f # u # v # r #|"compareG"
+    (u # v #! r) f #|"compareG"
     return r
 
 compareD = compareG c_compareD
@@ -500,7 +503,7 @@
 
 selectG f c u v w = unsafePerformIO $ do
     r <- createVector (dim v)
-    f # c # u # v # w # r #|"selectG"
+    (c # u # v # w #! r) f #|"selectG"
     return r
 
 selectD = selectG c_selectD
@@ -523,7 +526,7 @@
 
 remapG f i j m = unsafePerformIO $ do
     r <- createMatrix RowMajor (rows i) (cols i)
-    f # i # j # m # r #|"remapG"
+    (i # j # m #! r) f #|"remapG"
     return r
 
 remapD = remapG c_remapD
@@ -546,7 +549,7 @@
 
 rowOpAux f c x i1 i2 j1 j2 m = do
     px <- newArray [x]
-    f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2) # m #|"rowOp"
+    (m # id) (f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2)) #|"rowOp"
     free px
 
 type RowOp x = CInt -> Ptr x -> CInt -> CInt -> CInt -> CInt -> x ::> Ok
@@ -562,7 +565,7 @@
 
 --------------------------------------------------------------------------------
 
-gemmg f v m1 m2 m3 = f # v # m1 # m2 # m3 #|"gemmg"
+gemmg f v m1 m2 m3 = (v # m1 # m2 #! m3) f #|"gemmg"
 
 type Tgemm x = x :> x ::> x ::> x ::> Ok
 
@@ -590,10 +593,9 @@
 saveMatrix name format m = do
     cname   <- newCString name
     cformat <- newCString format
-    c_saveMatrix cname cformat # m #|"saveMatrix"
+    (m # id) (c_saveMatrix cname cformat) #|"saveMatrix"
     free cname
     free cformat
     return ()
 
 --------------------------------------------------------------------------------
-
diff --git a/src/Internal/Modular.hs b/src/Internal/Modular.hs
--- a/src/Internal/Modular.hs
+++ b/src/Internal/Modular.hs
@@ -72,7 +72,7 @@
   where
     compare a b = compare (unMod a) (unMod b)
 
-instance (Integral t, KnownNat m, Integral (Mod m t)) => Real (Mod m t)
+instance (Integral t, Real t, KnownNat m, Integral (Mod m t)) => Real (Mod m t)
   where
     toRational x = toInteger x % 1
 
@@ -114,7 +114,7 @@
   where
     show = show . unMod
 
-instance forall n t . (Integral t, KnownNat n) => Num (Mod n t)
+instance (Integral t, KnownNat n) => Num (Mod n t)
   where
     (+) = l2 (\m a b -> (a + b) `mod` (fromIntegral m))
     (*) = l2 (\m a b -> (a * b) `mod` (fromIntegral m))
@@ -159,11 +159,11 @@
         m' = fromIntegral . natVal $ (undefined :: Proxy m)
 
 
-instance forall m . KnownNat m => CTrans (Mod m I)
-instance forall m . KnownNat m => CTrans (Mod m Z)
+instance KnownNat m => CTrans (Mod m I)
+instance KnownNat m => CTrans (Mod m Z)
 
 
-instance forall m . KnownNat m => Container Vector (Mod m I)
+instance KnownNat m => Container Vector (Mod m I)
   where
     conj' = id
     size' = dim
@@ -203,7 +203,7 @@
     fromZ'   = vmod . fromZ'
     toZ'     = toZ' . f2i
 
-instance forall m . KnownNat m => Container Vector (Mod m Z)
+instance KnownNat m => Container Vector (Mod m Z)
   where
     conj' = id
     size' = dim
@@ -311,7 +311,7 @@
 
 lift2m f a b = liftMatrix vmod (f (f2iM a) (f2iM b))
 
-instance forall m . KnownNat m => Num (Vector (Mod m I))
+instance KnownNat m => Num (Vector (Mod m I))
   where
     (+) = lift2 (+)
     (*) = lift2 (*)
@@ -321,7 +321,7 @@
     negate = lift1 negate
     fromInteger x = fromInt (fromInteger x)
 
-instance forall m . KnownNat m => Num (Vector (Mod m Z))
+instance KnownNat m => Num (Vector (Mod m Z))
   where
     (+) = lift2 (+)
     (*) = lift2 (*)
diff --git a/src/Internal/Numeric.hs b/src/Internal/Numeric.hs
--- a/src/Internal/Numeric.hs
+++ b/src/Internal/Numeric.hs
@@ -1,4 +1,3 @@
-{-# LANGUAGE CPP #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
diff --git a/src/Internal/Random.hs b/src/Internal/Random.hs
--- a/src/Internal/Random.hs
+++ b/src/Internal/Random.hs
@@ -31,13 +31,13 @@
 gaussianSample :: Seed
                -> Int -- ^ number of rows
                -> Vector Double -- ^ mean vector
-               -> Matrix Double -- ^ covariance matrix
+               -> Herm Double   -- ^ covariance matrix
                -> Matrix Double -- ^ result
 gaussianSample seed n med cov = m where
     c = dim med
     meds = konst' 1 n `outer` med
     rs = reshape c $ randomVector seed Gaussian (c * n)
-    m = rs `mXm` cholSH cov `add` meds
+    m = rs `mXm` chol cov `add` meds
 
 -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
 -- uniform distribution.
diff --git a/src/Internal/Sparse.hs b/src/Internal/Sparse.hs
--- a/src/Internal/Sparse.hs
+++ b/src/Internal/Sparse.hs
@@ -144,13 +144,13 @@
 gmXv SparseR { gmCSR = CSR{..}, .. } v = unsafePerformIO $ do
     dim v /= nCols ~!~ printf "gmXv (CSR): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)
     r <- createVector nRows
-    c_smXv # csrVals # csrCols # csrRows # v # r #|"CSRXv"
+    (csrVals # csrCols # csrRows # v #! r) c_smXv #|"CSRXv"
     return r
 
 gmXv SparseC { gmCSC = CSC{..}, .. } v = unsafePerformIO $ do
     dim v /= nCols ~!~ printf "gmXv (CSC): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)
     r <- createVector nRows
-    c_smTXv # cscVals # cscRows # cscCols # v # r #|"CSCXv"
+    (cscVals # cscRows # cscCols # v #! r) c_smTXv #|"CSCXv"
     return r
 
 gmXv Diag{..} v
@@ -211,4 +211,3 @@
     tr (Diag v n m) = Diag v m n
     tr (Dense a n m) = Dense (tr a) m n
     tr' = tr
-
diff --git a/src/Internal/Static.hs b/src/Internal/Static.hs
--- a/src/Internal/Static.hs
+++ b/src/Internal/Static.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE CPP #-}
 #if __GLASGOW_HASKELL__ >= 708
 
 {-# LANGUAGE DataKinds #-}
@@ -11,6 +12,8 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE TypeOperators #-}
 {-# LANGUAGE ViewPatterns #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE DeriveGeneric #-}
 
 {- |
 Module      :  Internal.Static
@@ -28,18 +31,35 @@
 import Numeric.LinearAlgebra hiding (konst,size,R,C)
 import Internal.Vector as D hiding (R,C)
 import Internal.ST
+import Control.DeepSeq
 import Data.Proxy(Proxy)
 import Foreign.Storable(Storable)
 import Text.Printf
 
+import Data.Binary
+import GHC.Generics (Generic)
+import Data.Proxy (Proxy(..))
+
 --------------------------------------------------------------------------------
 
 type ℝ = Double
 type ℂ = Complex Double
 
 newtype Dim (n :: Nat) t = Dim t
-  deriving Show
+  deriving (Show, Generic)
 
+instance (KnownNat n, Binary a) => Binary (Dim n a) where
+  get = do
+    k <- get
+    let n = natVal (Proxy :: Proxy n)
+    if n == k
+      then Dim <$> get
+      else fail ("Expected dimension " ++ (show n) ++ ", but found dimension " ++ (show k))
+
+  put (Dim x) = do
+    put (natVal (Proxy :: Proxy n))
+    put x
+
 lift1F
   :: (c t -> c t)
   -> Dim n (c t) -> Dim n (c t)
@@ -50,18 +70,22 @@
   -> Dim n (c t) -> Dim n (c t) -> Dim n (c t)
 lift2F f (Dim u) (Dim v) = Dim (f u v)
 
+instance NFData t => NFData (Dim n t) where
+    rnf (Dim (force -> !_)) = ()
+
 --------------------------------------------------------------------------------
 
 newtype R n = R (Dim n (Vector ℝ))
-  deriving (Num,Fractional,Floating)
+  deriving (Num,Fractional,Floating,Generic,Binary)
 
 newtype C n = C (Dim n (Vector ℂ))
-  deriving (Num,Fractional,Floating)
+  deriving (Num,Fractional,Floating,Generic)
 
 newtype L m n = L (Dim m (Dim n (Matrix ℝ)))
-
-newtype M m n = M (Dim m (Dim n (Matrix  ℂ)))
+  deriving (Generic, Binary)
 
+newtype M m n = M (Dim m (Dim n (Matrix ℂ)))
+  deriving (Generic)
 
 mkR :: Vector ℝ -> R n
 mkR = R . Dim
@@ -75,6 +99,18 @@
 mkM :: Matrix ℂ -> M m n
 mkM x = M (Dim (Dim x))
 
+instance NFData (R n) where
+    rnf (R (force -> !_)) = ()
+
+instance NFData (C n) where
+    rnf (C (force -> !_)) = ()
+
+instance NFData (L n m) where
+    rnf (L (force -> !_)) = ()
+
+instance NFData (M n m) where
+    rnf (M (force -> !_)) = ()
+
 --------------------------------------------------------------------------------
 
 type V n t = Dim n (Vector t)
@@ -173,7 +209,7 @@
 singleM m = rows m == 1 && cols m == 1
 
 
-instance forall n. KnownNat n => Sized ℂ (C n) Vector
+instance KnownNat n => Sized ℂ (C n) Vector
   where
     size _ = fromIntegral . natVal $ (undefined :: Proxy n)
     konst x = mkC (LA.scalar x)
@@ -189,7 +225,7 @@
         r = mkC v :: C n
 
 
-instance forall n. KnownNat n => Sized ℝ (R n) Vector
+instance KnownNat n => Sized ℝ (R n) Vector
   where
     size _ = fromIntegral . natVal $ (undefined :: Proxy n)
     konst x = mkR (LA.scalar x)
@@ -206,7 +242,7 @@
 
 
 
-instance forall m n . (KnownNat m, KnownNat n) => Sized ℝ (L m n) Matrix
+instance (KnownNat m, KnownNat n) => Sized ℝ (L m n) Matrix
   where
     size _ = ((fromIntegral . natVal) (undefined :: Proxy m)
              ,(fromIntegral . natVal) (undefined :: Proxy n))
@@ -224,7 +260,7 @@
         r = mkL x :: L m n
 
 
-instance forall m n . (KnownNat m, KnownNat n) => Sized ℂ (M m n) Matrix
+instance (KnownNat m, KnownNat n) => Sized ℂ (M m n) Matrix
   where
     size _ = ((fromIntegral . natVal) (undefined :: Proxy m)
              ,(fromIntegral . natVal) (undefined :: Proxy n))
@@ -282,7 +318,7 @@
 
 --------------------------------------------------------------------------------
 
-instance forall n . KnownNat n => Show (R n)
+instance KnownNat n => Show (R n)
   where
     show s@(R (Dim v))
       | singleV v = "("++show (v!0)++" :: R "++show d++")"
@@ -290,7 +326,7 @@
       where
         d = size s
 
-instance forall n . KnownNat n => Show (C n)
+instance KnownNat n => Show (C n)
   where
     show s@(C (Dim v))
       | singleV v = "("++show (v!0)++" :: C "++show d++")"
@@ -298,7 +334,7 @@
       where
         d = size s
 
-instance forall m n . (KnownNat m, KnownNat n) => Show (L m n)
+instance (KnownNat m, KnownNat n) => Show (L m n)
   where
     show (isDiag -> Just (z,y,(m',n'))) = printf "(diag %s %s :: L %d %d)" (show z) (drop 9 $ show y) m' n'
     show s@(L (Dim (Dim x)))
@@ -307,7 +343,7 @@
       where
         (m',n') = size s
 
-instance forall m n . (KnownNat m, KnownNat n) => Show (M m n)
+instance (KnownNat m, KnownNat n) => Show (M m n)
   where
     show (isDiagC -> Just (z,y,(m',n'))) = printf "(diag %s %s :: M %d %d)" (show z) (drop 9 $ show y) m' n'
     show s@(M (Dim (Dim x)))
@@ -318,7 +354,7 @@
 
 --------------------------------------------------------------------------------
 
-instance forall n t . (Num (Vector t), Numeric t )=> Num (Dim n (Vector t))
+instance (Num (Vector t), Numeric t )=> Num (Dim n (Vector t))
   where
     (+) = lift2F (+)
     (*) = lift2F (*)
@@ -482,6 +518,18 @@
     (**)  = lift2MD (**)
     pi    = M pi
 
+instance Additive (R n) where
+    add = (+)
+
+instance Additive (C n) where
+    add = (+)
+
+instance (KnownNat m, KnownNat n) => Additive (L m n) where
+    add = (+)
+
+instance (KnownNat m, KnownNat n) => Additive (M m n) where
+    add = (+)
+
 --------------------------------------------------------------------------------
 
 
@@ -518,6 +566,17 @@
         putStr "C " >> putStr (tail . dropWhile (/='x') $ su)
 
 --------------------------------------------------------------------------------
+
+overMatL' :: (KnownNat m, KnownNat n)
+          => (LA.Matrix ℝ -> LA.Matrix ℝ) -> L m n -> L m n
+overMatL' f = mkL . f . unwrap
+{-# INLINE overMatL' #-}
+
+overMatM' :: (KnownNat m, KnownNat n)
+          => (LA.Matrix ℂ -> LA.Matrix ℂ) -> M m n -> M m n
+overMatM' f = mkM . f . unwrap
+{-# INLINE overMatM' #-}
+
 
 #else
 
diff --git a/src/Internal/Util.hs b/src/Internal/Util.hs
--- a/src/Internal/Util.hs
+++ b/src/Internal/Util.hs
@@ -41,7 +41,7 @@
     ℕ,ℤ,ℝ,ℂ,iC,
     Normed(..), norm_Frob, norm_nuclear,
     magnit,
-    unitary,
+    normalize,
     mt,
     (~!~),
     pairwiseD2,
@@ -255,6 +255,7 @@
 -- ^ 2-norm of real vector
 norm = pnorm PNorm2
 
+-- | p-norm for vectors, operator norm for matrices
 class Normed a
   where
     norm_0   :: a -> R
@@ -319,10 +320,11 @@
     norm_2 = norm_2 . double
     norm_Inf = norm_Inf . double
 
-
+-- | Frobenius norm (Schatten p-norm with p=2)
 norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R
 norm_Frob = norm_2 . flatten
 
+-- | Sum of singular values (Schatten p-norm with p=1)
 norm_nuclear :: Field t => Matrix t -> R
 norm_nuclear = sumElements . singularValues
 
@@ -341,8 +343,8 @@
 
 
 -- | Obtains a vector in the same direction with 2-norm=1
-unitary :: Vector Double -> Vector Double
-unitary v = v / scalar (norm v)
+normalize :: (Normed (Vector t), Num (Vector t), Field t) => Vector t -> Vector t
+normalize v = v / real (scalar (norm_2 v))
 
 
 -- | trans . inv
diff --git a/src/Internal/Vector.hs b/src/Internal/Vector.hs
--- a/src/Internal/Vector.hs
+++ b/src/Internal/Vector.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-}
+{-# LANGUAGE MagicHash, UnboxedTuples, BangPatterns, FlexibleContexts #-}
 {-# LANGUAGE TypeSynonymInstances #-}
 
 
@@ -39,12 +39,10 @@
 import qualified Data.Vector.Storable as Vector
 import Data.Vector.Storable(Vector, fromList, unsafeToForeignPtr, unsafeFromForeignPtr, unsafeWith)
 
-#ifdef BINARY
 import Data.Binary
 import Control.Monad(replicateM)
 import qualified Data.ByteString.Internal as BS
 import Data.Vector.Storable.Internal(updPtr)
-#endif
 
 type I = CInt
 type Z = Int64
@@ -64,13 +62,13 @@
 -- | Number of elements
 dim :: (Storable t) => Vector t -> Int
 dim = Vector.length
+{-# INLINE dim #-}
 
 
 -- C-Haskell vector adapter
 {-# INLINE avec #-}
-avec :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b
-avec f v = inlinePerformIO (unsafeWith v (return . f (fromIntegral (Vector.length v))))
-infixl 1 `avec`
+avec :: Storable a => Vector a -> (f -> IO r) -> ((CInt -> Ptr a -> f) -> IO r)
+avec v f g = unsafeWith v $ \ptr -> f (g (fromIntegral (Vector.length v)) ptr)
 
 -- allocates memory for a new vector
 createVector :: Storable a => Int -> IO (Vector a)
@@ -142,6 +140,7 @@
                         -> Vector t  -- ^ source
                         -> Vector t  -- ^ result
 subVector = Vector.slice
+{-# INLINE subVector #-}
 
 
 
@@ -201,7 +200,7 @@
 
 ---------------------------------------------------------------
 
--- | transforms a complex vector into a real vector with alternating real and imaginary parts 
+-- | transforms a complex vector into a real vector with alternating real and imaginary parts
 asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a
 asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n)
     where (fp,i,n) = unsafeToForeignPtr v
@@ -246,7 +245,7 @@
 {-# INLINE zipVectorWith #-}
 
 -- | unzipWith for Vectors
-unzipVectorWith :: (Storable (a,b), Storable c, Storable d) 
+unzipVectorWith :: (Storable (a,b), Storable c, Storable d)
                    => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d)
 unzipVectorWith f u = unsafePerformIO $ do
       let n = dim u
@@ -257,7 +256,7 @@
               unsafeWith w $ \pw -> do
                   let go (-1) = return ()
                       go !k   = do z <- peekElemOff pu k
-                                   let (x,y) = f z 
+                                   let (x,y) = f z
                                    pokeElemOff      pv k x
                                    pokeElemOff      pw k y
                                    go (k-1)
@@ -305,11 +304,11 @@
     return w
     where mapVectorM' w' !k !t
               | k == t               = do
-                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                        y <- f x
                                        return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
               | otherwise            = do
-                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                        y <- f x
                                        _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
                                        mapVectorM' w' (k+1) t
@@ -324,7 +323,7 @@
                                     x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                     f x
               | otherwise         = do
-                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                     _ <- f x
                                     mapVectorM' (k+1) t
 {-# INLINE mapVectorM_ #-}
@@ -338,11 +337,11 @@
     return w
     where mapVectorM' w' !k !t
               | k == t               = do
-                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                        y <- f k x
                                        return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
               | otherwise            = do
-                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                        y <- f k x
                                        _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
                                        mapVectorM' w' (k+1) t
@@ -357,7 +356,7 @@
                                     x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                     f k x
               | otherwise         = do
-                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                     _ <- f k x
                                     mapVectorM' (k+1) t
 {-# INLINE mapVectorWithIndexM_ #-}
@@ -380,7 +379,6 @@
 --------------------------------------------------------------------------------
 
 
-#ifdef BINARY
 
 -- a 64K cache, with a Double taking 13 bytes in Bytestring,
 -- implies a chunk size of 5041
@@ -432,7 +430,6 @@
 
     -- get = fmap bs2v get
 
-#endif
 
 
 -------------------------------------------------------------------
@@ -458,4 +455,3 @@
 unzipVector = unzipVectorWith id
 
 -------------------------------------------------------------------
-
diff --git a/src/Internal/Vectorized.hs b/src/Internal/Vectorized.hs
--- a/src/Internal/Vectorized.hs
+++ b/src/Internal/Vectorized.hs
@@ -27,10 +27,13 @@
 import System.IO.Unsafe(unsafePerformIO)
 import Control.Monad(when)
 
-infixl 1 #
+infixr 1 #
 a # b = applyRaw a b
 {-# INLINE (#) #-}
 
+a #! b = a # b # id
+{-# INLINE (#!) #-}
+
 fromei x = fromIntegral (fromEnum x) :: CInt
 
 data FunCodeV = Sin
@@ -103,7 +106,7 @@
 
 sumg f x = unsafePerformIO $ do
     r <- createVector 1
-    f # x # r #| "sum"
+    (x #! r) f #| "sum"
     return $ r @> 0
 
 type TVV t = t :> t :> Ok
@@ -139,7 +142,7 @@
 
 prodg f x = unsafePerformIO $ do
     r <- createVector 1
-    f # x # r #| "prod"
+    (x #! r) f #| "prod"
     return $ r @> 0
 
 
@@ -154,24 +157,24 @@
 
 toScalarAux fun code v = unsafePerformIO $ do
     r <- createVector 1
-    fun (fromei code) # v # r #|"toScalarAux"
+    (v #! r) (fun (fromei code)) #|"toScalarAux"
     return (r @> 0)
 
 vectorMapAux fun code v = unsafePerformIO $ do
     r <- createVector (dim v)
-    fun (fromei code) # v # r #|"vectorMapAux"
+    (v #! r) (fun (fromei code)) #|"vectorMapAux"
     return r
 
 vectorMapValAux fun code val v = unsafePerformIO $ do
     r <- createVector (dim v)
     pval <- newArray [val]
-    fun (fromei code) pval # v # r #|"vectorMapValAux"
+    (v #! r) (fun (fromei code) pval) #|"vectorMapValAux"
     free pval
     return r
 
 vectorZipAux fun code u v = unsafePerformIO $ do
     r <- createVector (dim u)
-    fun (fromei code) # u # v # r #|"vectorZipAux"
+    (u # v #! r) (fun (fromei code)) #|"vectorZipAux"
     return r
 
 ---------------------------------------------------------------------
@@ -368,7 +371,7 @@
              -> Vector Double
 randomVector seed dist n = unsafePerformIO $ do
     r <- createVector n
-    c_random_vector (fi seed) ((fi.fromEnum) dist) # r #|"randomVector"
+    (r # id) (c_random_vector (fi seed) ((fi.fromEnum) dist)) #|"randomVector"
     return r
 
 foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> Double :> Ok
@@ -377,7 +380,7 @@
 
 roundVector v = unsafePerformIO $ do
     r <- createVector (dim v)
-    c_round_vector # v # r #|"roundVector"
+    (v #! r) c_round_vector #|"roundVector"
     return r
 
 foreign import ccall unsafe "round_vector" c_round_vector :: TVV Double
@@ -391,7 +394,7 @@
 range :: Int -> Vector I
 range n = unsafePerformIO $ do
     r <- createVector n
-    c_range_vector # r #|"range"
+    (r # id) c_range_vector #|"range"
     return r
 
 foreign import ccall unsafe "range_vector" c_range_vector :: CInt :> Ok
@@ -431,7 +434,7 @@
 
 tog f v = unsafePerformIO $ do
     r <- createVector (dim v)
-    f # v # r #|"tog"
+    (v #! r) f #|"tog"
     return r
 
 foreign import ccall unsafe "float2double" c_float2double :: Float :> Double :> Ok
@@ -450,7 +453,7 @@
 
 stepg f v = unsafePerformIO $ do
     r <- createVector (dim v)
-    f # v # r #|"step"
+    (v #! r) f #|"step"
     return r
 
 stepD :: Vector Double -> Vector Double
@@ -475,7 +478,7 @@
 
 conjugateAux fun x = unsafePerformIO $ do
     v <- createVector (dim x)
-    fun # x # v #|"conjugateAux"
+    (x #! v) fun #|"conjugateAux"
     return v
 
 conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
@@ -493,7 +496,7 @@
         let n = dim v
         r <- createVector n
         let f _ s _ d =  copyArray d s n >> return 0
-        f # v # r #|"cloneVector"
+        (v #! r) f #|"cloneVector"
         return r
 
 --------------------------------------------------------------------------------
@@ -501,7 +504,7 @@
 constantAux fun x n = unsafePerformIO $ do
     v <- createVector n
     px <- newArray [x]
-    fun px # v #|"constantAux"
+    (v # id) (fun px) #|"constantAux"
     free px
     return v
 
@@ -515,4 +518,3 @@
 foreign import ccall unsafe "constantL" cconstantL :: TConst Z
 
 ----------------------------------------------------------------------
-
diff --git a/src/Numeric/LinearAlgebra.hs b/src/Numeric/LinearAlgebra.hs
--- a/src/Numeric/LinearAlgebra.hs
+++ b/src/Numeric/LinearAlgebra.hs
@@ -22,7 +22,8 @@
 
     -- * Numeric classes
     -- |
-    -- The standard numeric classes are defined elementwise:
+    -- The standard numeric classes are defined elementwise (commonly referred to
+    -- as the Hadamard product or the Schur product):
     --
     -- >>>  vector [1,2,3] * vector [3,0,-2]
     -- fromList [3.0,0.0,-6.0]
@@ -116,6 +117,7 @@
     svd,
     thinSVD,
     compactSVD,
+    compactSVDTol,
     singularValues,
     leftSV, rightSV,
 
@@ -125,7 +127,7 @@
     geigSH,
 
     -- * QR
-    qr, rq, qrRaw, qrgr,
+    qr, thinQR, rq, thinRQ, qrRaw, qrgr,
 
     -- * Cholesky
     chol, mbChol,
@@ -152,7 +154,7 @@
     Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
 
     -- * Misc
-    meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, magnit,
+    meanCov, rowOuters, pairwiseD2, normalize, peps, relativeError, magnit,
     haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv,
     iC, sym, mTm, trustSym, unSym,
     -- * Auxiliary classes
diff --git a/src/Numeric/LinearAlgebra/Devel.hs b/src/Numeric/LinearAlgebra/Devel.hs
--- a/src/Numeric/LinearAlgebra/Devel.hs
+++ b/src/Numeric/LinearAlgebra/Devel.hs
@@ -12,13 +12,6 @@
 --------------------------------------------------------------------------------
 
 module Numeric.LinearAlgebra.Devel(
-    -- * FFI helpers
-    -- | Sample usage, to upload a perspective matrix to a shader.
-    --
-    -- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3)
-    -- @
-    module Internal.Foreign,
-
     -- * FFI tools
     -- | See @examples/devel@ in the repository.
     
@@ -66,7 +59,6 @@
 
 ) where
 
-import Internal.Foreign
 import Internal.Devel
 import Internal.ST
 import Internal.Vector
diff --git a/src/Numeric/LinearAlgebra/Static.hs b/src/Numeric/LinearAlgebra/Static.hs
--- a/src/Numeric/LinearAlgebra/Static.hs
+++ b/src/Numeric/LinearAlgebra/Static.hs
@@ -42,19 +42,24 @@
     blockAt,
     matrix,
     -- * Complex
-    C, M, Her, her, 𝑖,
+    ℂ, C, M, Her, her, 𝑖,
     -- * Products
     (<>),(#>),(<.>),
     -- * Linear Systems
     linSolve, (<\>),
     -- * Factorizations
     svd, withCompactSVD, svdTall, svdFlat, Eigen(..),
-    withNullspace, qr, chol,
+    withNullspace, withOrth, qr, chol,
+    -- * Norms
+    Normed(..),
+    -- * Random arrays
+    Seed, RandDist(..),
+    randomVector, rand, randn, gaussianSample, uniformSample,
     -- * Misc
-    mean,
+    mean, meanCov,
     Disp(..), Domain(..),
-    withVector, withMatrix,
-    toRows, toColumns,
+    withVector, withMatrix, exactLength, exactDims,
+    toRows, toColumns, withRows, withColumns,
     Sized(..), Diag(..), Sym, sym, mTm, unSym, (<·>)
 ) where
 
@@ -65,11 +70,16 @@
     row,col,vector,matrix,linspace,toRows,toColumns,
     (<\>),fromList,takeDiag,svd,eig,eigSH,
     eigenvalues,eigenvaluesSH,build,
-    qr,size,dot,chol,range,R,C,sym,mTm,unSym)
+    qr,size,dot,chol,range,R,C,sym,mTm,unSym,
+    randomVector,rand,randn,gaussianSample,uniformSample,meanCov)
 import qualified Numeric.LinearAlgebra as LA
-import Data.Proxy(Proxy)
+import qualified Numeric.LinearAlgebra.Devel as LA
+import Data.Proxy(Proxy(..))
 import Internal.Static
 import Control.Arrow((***))
+import Text.Printf
+import Data.Type.Equality ((:~:)(Refl))
+import qualified Data.Bifunctor as BF (first)
 
 ud1 :: R n -> Vector ℝ
 ud1 (R (Dim v)) = v
@@ -218,19 +228,14 @@
     takeDiag :: m -> d
 
 
-instance forall n . (KnownNat n) => Diag (L n n) (R n)
-  where
-    takeDiag m = mkR (LA.takeDiag (extract m))
-
-
-instance forall m n . (KnownNat m, KnownNat n, m <= n+1) => Diag (L m n) (R m)
+instance KnownNat n => Diag (L n n) (R n)
   where
-    takeDiag m = mkR (LA.takeDiag (extract m))
+    takeDiag x = mkR (LA.takeDiag (extract x))
 
 
-instance forall m n . (KnownNat m, KnownNat n, n <= m+1) => Diag (L m n) (R n)
+instance KnownNat n => Diag (M n n) (C n)
   where
-    takeDiag m = mkR (LA.takeDiag (extract m))
+    takeDiag x = mkC (LA.takeDiag (extract x))
 
 
 --------------------------------------------------------------------------------
@@ -289,7 +294,21 @@
 her m = Her $ (m + LA.tr m)/2
 
 
+instance (KnownNat n) => Disp (Sym n)
+  where
+    disp n (Sym x) = do
+        let a = extract x
+        let su = LA.dispf n a
+        printf "Sym %d" (cols a) >> putStr (dropWhile (/='\n') $ su)
 
+instance (KnownNat n) => Disp (Her n)
+  where
+    disp n (Her x) = do
+        let a = extract x
+        let su = LA.dispcf n a
+        printf "Her %d" (cols a) >> putStr (dropWhile (/='\n') $ su)
+
+
 instance KnownNat n => Eigen (Sym n) (R n) (L n n)
   where
     eigenvalues (Sym (extract -> m)) =  mkR . LA.eigenvaluesSH . LA.trustSym $ m
@@ -319,6 +338,15 @@
        Nothing -> error "static/dynamic mismatch"
        Just (SomeNat (_ :: Proxy k)) -> f (mkL a :: L n k)
 
+withOrth
+    :: forall m n z . (KnownNat m, KnownNat n)
+    => L m n
+    -> (forall k. (KnownNat k) => L n k -> z)
+    -> z
+withOrth (LA.orth . extract -> a) f =
+    case someNatVal $ fromIntegral $ cols a of
+       Nothing -> error "static/dynamic mismatch"
+       Just (SomeNat (_ :: Proxy k)) -> f (mkL a :: L n k)
 
 withCompactSVD
     :: forall m n z . (KnownNat m, KnownNat n)
@@ -365,11 +393,31 @@
 toRows :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R n]
 toRows (LA.toRows . extract -> vs) = map mkR vs
 
+withRows
+    :: forall n z . KnownNat n
+    => [R n]
+    -> (forall m . KnownNat m => L m n -> z)
+    -> z
+withRows (LA.fromRows . map extract -> m) f =
+    case someNatVal $ fromIntegral $ LA.rows m of
+       Nothing -> error "static/dynamic mismatch"
+       Just (SomeNat (_ :: Proxy m)) -> f (mkL m :: L m n)
 
 toColumns :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R m]
 toColumns (LA.toColumns . extract -> vs) = map mkR vs
 
+withColumns
+    :: forall m z . KnownNat m
+    => [R m]
+    -> (forall n . KnownNat n => L m n -> z)
+    -> z
+withColumns (LA.fromColumns . map extract -> m) f =
+    case someNatVal $ fromIntegral $ LA.cols m of
+       Nothing -> error "static/dynamic mismatch"
+       Just (SomeNat (_ :: Proxy n)) -> f (mkL m :: L m n)
 
+
+
 --------------------------------------------------------------------------------
 
 build
@@ -392,6 +440,15 @@
        Nothing -> error "static/dynamic mismatch"
        Just (SomeNat (_ :: Proxy m)) -> f (mkR v :: R m)
 
+-- | Useful for constraining two dependently typed vectors to match each
+-- other in length when they are unknown at compile-time.
+exactLength
+    :: forall n m . (KnownNat n, KnownNat m)
+    => R m
+    -> Maybe (R n)
+exactLength v = do
+    Refl <- sameNat (Proxy :: Proxy n) (Proxy :: Proxy m)
+    return $ mkR (unwrap v)
 
 withMatrix
     :: forall z
@@ -407,6 +464,64 @@
                Just (SomeNat (_ :: Proxy n)) ->
                   f (mkL a :: L m n)
 
+-- | Useful for constraining two dependently typed matrices to match each
+-- other in dimensions when they are unknown at compile-time.
+exactDims
+    :: forall n m j k . (KnownNat n, KnownNat m, KnownNat j, KnownNat k)
+    => L m n
+    -> Maybe (L j k)
+exactDims m = do
+    Refl <- sameNat (Proxy :: Proxy m) (Proxy :: Proxy j)
+    Refl <- sameNat (Proxy :: Proxy n) (Proxy :: Proxy k)
+    return $ mkL (unwrap m)
+
+randomVector
+    :: forall n . KnownNat n
+    => Seed
+    -> RandDist
+    -> R n
+randomVector s d = mkR (LA.randomVector s d
+                          (fromInteger (natVal (Proxy :: Proxy n)))
+                       )
+
+rand
+    :: forall m n . (KnownNat m, KnownNat n)
+    => IO (L m n)
+rand = mkL <$> LA.rand (fromInteger (natVal (Proxy :: Proxy m)))
+                       (fromInteger (natVal (Proxy :: Proxy n)))
+
+randn
+    :: forall m n . (KnownNat m, KnownNat n)
+    => IO (L m n)
+randn = mkL <$> LA.randn (fromInteger (natVal (Proxy :: Proxy m)))
+                         (fromInteger (natVal (Proxy :: Proxy n)))
+
+gaussianSample
+    :: forall m n . (KnownNat m, KnownNat n)
+    => Seed
+    -> R n
+    -> Sym n
+    -> L m n
+gaussianSample s (extract -> mu) (Sym (extract -> sigma)) =
+    mkL $ LA.gaussianSample s (fromInteger (natVal (Proxy :: Proxy m)))
+                            mu (LA.trustSym sigma)
+
+uniformSample
+    :: forall m n . (KnownNat m, KnownNat n)
+    => Seed
+    -> R n    -- ^ minimums of each row
+    -> R n    -- ^ maximums of each row
+    -> L m n
+uniformSample s (extract -> mins) (extract -> maxs) =
+    mkL $ LA.uniformSample s (fromInteger (natVal (Proxy :: Proxy m)))
+                           (zip (LA.toList mins) (LA.toList maxs))
+
+meanCov
+    :: forall m n . (KnownNat m, KnownNat n, 1 <= m)
+    => L m n
+    -> (R n, Sym n)
+meanCov (extract -> vs) = mkR *** (Sym . mkL . LA.unSym) $ LA.meanCov vs
+
 --------------------------------------------------------------------------------
 
 class Domain field vec mat | mat -> vec field, vec -> mat field, field -> mat vec
@@ -416,6 +531,15 @@
     dot :: forall n . (KnownNat n) => vec n -> vec n -> field
     cross :: vec 3 -> vec 3 -> vec 3
     diagR ::  forall m n k . (KnownNat m, KnownNat n, KnownNat k) => field -> vec k -> mat m n
+    dvmap :: forall n. KnownNat n => (field -> field) -> vec n -> vec n
+    dmmap :: forall n m. (KnownNat m, KnownNat n) => (field -> field) -> mat n m -> mat n m
+    outer :: forall n m. (KnownNat m, KnownNat n) => vec n -> vec m -> mat n m
+    zipWithVector :: forall n. KnownNat n => (field -> field -> field) -> vec n -> vec n -> vec n
+    det :: forall n. KnownNat n => mat n n -> field
+    invlndet :: forall n. KnownNat n => mat n n -> (mat n n, (field, field))
+    expm :: forall n. KnownNat n => mat n n -> mat n n
+    sqrtm :: forall n. KnownNat n => mat n n -> mat n n
+    inv :: forall n. KnownNat n => mat n n -> mat n n
 
 
 instance Domain ℝ R L
@@ -425,6 +549,15 @@
     dot = dotR
     cross = crossR
     diagR = diagRectR
+    dvmap = mapR
+    dmmap = mapL
+    outer = outerR
+    zipWithVector = zipWithR
+    det = detL
+    invlndet = invlndetL
+    expm = expmL
+    sqrtm = sqrtmL
+    inv = invL
 
 instance Domain ℂ C M
   where
@@ -433,6 +566,15 @@
     dot = dotC
     cross = crossC
     diagR = diagRectC
+    dvmap = mapC
+    dmmap = mapM'
+    outer = outerC
+    zipWithVector = zipWithC
+    det = detM
+    invlndet = invlndetM
+    expm = expmM
+    sqrtm = sqrtmM
+    inv = invM
 
 --------------------------------------------------------------------------------
 
@@ -472,6 +614,33 @@
     z2 = x!2*y!0-x!0*y!2
     z3 = x!0*y!1-x!1*y!0
 
+outerR :: (KnownNat m, KnownNat n) => R n -> R m -> L n m
+outerR (extract -> x) (extract -> y) = mkL (LA.outer x y)
+
+mapR :: KnownNat n => (ℝ -> ℝ) -> R n -> R n
+mapR f (unwrap -> v) = mkR (LA.cmap f v)
+
+zipWithR :: KnownNat n => (ℝ -> ℝ -> ℝ) -> R n -> R n -> R n
+zipWithR f (extract -> x) (extract -> y) = mkR (LA.zipVectorWith f x y)
+
+mapL :: (KnownNat n, KnownNat m) => (ℝ -> ℝ) -> L n m -> L n m
+mapL f = overMatL' (LA.cmap f)
+
+detL :: KnownNat n => Sq n -> ℝ
+detL = LA.det . unwrap
+
+invlndetL :: KnownNat n => Sq n -> (L n n, (ℝ, ℝ))
+invlndetL = BF.first mkL . LA.invlndet . unwrap
+
+expmL :: KnownNat n => Sq n -> Sq n
+expmL = overMatL' LA.expm
+
+sqrtmL :: KnownNat n => Sq n -> Sq n
+sqrtmL = overMatL' LA.sqrtm
+
+invL :: KnownNat n => Sq n -> Sq n
+invL = overMatL' LA.inv
+
 --------------------------------------------------------------------------------
 
 mulC :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => M m k -> M k n -> M m n
@@ -510,6 +679,33 @@
     z2 = x!2*y!0-x!0*y!2
     z3 = x!0*y!1-x!1*y!0
 
+outerC :: (KnownNat m, KnownNat n) => C n -> C m -> M n m
+outerC (extract -> x) (extract -> y) = mkM (LA.outer x y)
+
+mapC :: KnownNat n => (ℂ -> ℂ) -> C n -> C n
+mapC f (unwrap -> v) = mkC (LA.cmap f v)
+
+zipWithC :: KnownNat n => (ℂ -> ℂ -> ℂ) -> C n -> C n -> C n
+zipWithC f (extract -> x) (extract -> y) = mkC (LA.zipVectorWith f x y)
+
+mapM' :: (KnownNat n, KnownNat m) => (ℂ -> ℂ) -> M n m -> M n m
+mapM' f = overMatM' (LA.cmap f)
+
+detM :: KnownNat n => M n n -> ℂ
+detM = LA.det . unwrap
+
+invlndetM :: KnownNat n => M n n -> (M n n, (ℂ, ℂ))
+invlndetM = BF.first mkM . LA.invlndet . unwrap
+
+expmM :: KnownNat n => M n n -> M n n
+expmM = overMatM' LA.expm
+
+sqrtmM :: KnownNat n => M n n -> M n n
+sqrtmM = overMatM' LA.sqrtm
+
+invM :: KnownNat n => M n n -> M n n
+invM = overMatM' LA.inv
+
 --------------------------------------------------------------------------------
 
 diagRectR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℝ -> R k -> L m n
@@ -616,3 +812,67 @@
   where
     checkT _ = test
 
+--------------------------------------------------------------------------------
+
+instance KnownNat n => Normed (R n)
+  where
+    norm_0 v = norm_0 (extract v)
+    norm_1 v = norm_1 (extract v)
+    norm_2 v = norm_2 (extract v)
+    norm_Inf v = norm_Inf (extract v)
+
+instance (KnownNat m, KnownNat n) => Normed (L m n)
+  where
+    norm_0 m = norm_0 (extract m)
+    norm_1 m = norm_1 (extract m)
+    norm_2 m = norm_2 (extract m)
+    norm_Inf m = norm_Inf (extract m)
+
+mkSym f = Sym . f . unSym
+mkSym2 f x y = Sym (f (unSym x) (unSym y))
+
+instance KnownNat n =>  Num (Sym n)
+  where
+    (+) = mkSym2 (+)
+    (*) = mkSym2 (*)
+    (-) = mkSym2 (-)
+    abs = mkSym abs
+    signum = mkSym signum
+    negate = mkSym negate
+    fromInteger = Sym . fromInteger
+
+instance KnownNat n => Fractional (Sym n)
+  where
+    fromRational = Sym . fromRational
+    (/) = mkSym2 (/)
+
+instance KnownNat n => Floating (Sym n)
+  where
+    sin   = mkSym sin
+    cos   = mkSym cos
+    tan   = mkSym tan
+    asin  = mkSym asin
+    acos  = mkSym acos
+    atan  = mkSym atan
+    sinh  = mkSym sinh
+    cosh  = mkSym cosh
+    tanh  = mkSym tanh
+    asinh = mkSym asinh
+    acosh = mkSym acosh
+    atanh = mkSym atanh
+    exp   = mkSym exp
+    log   = mkSym log
+    sqrt  = mkSym sqrt
+    (**)  = mkSym2 (**)
+    pi    = Sym pi
+
+instance KnownNat n => Additive (Sym n) where
+    add = (+)
+
+instance KnownNat n => Transposable (Sym n) (Sym n) where
+    tr  = id
+    tr' = id
+
+instance KnownNat n => Transposable (Her n) (Her n) where
+    tr          = id
+    tr' (Her m) = Her (tr' m)