18 files changed, 339 insertions, 196 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 5e6e649..8709a00 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -132,17 +132,22 @@ mat a f =
         let m g = do
             g (fi (rows a)) (fi (cols a)) p
         f m
--- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.
---
--- @\> flatten ('ident' 3)
--- 9 |> [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@
+
+{- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.
+
+>>> flatten (ident 3)
+fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]
+
+-}
 flatten :: Element t => Matrix t -> Vector t
 flatten = xdat . cmat
 
+{-
 type Mt t s = Int -> Int -> Ptr t -> s
--- not yet admitted by my haddock version
--- infixr 6 ::>
--- type t ::> s = Mt t s
+
+infixr 6 ::>
+type t ::> s = Mt t s
+-}
 
 -- | the inverse of 'Data.Packed.Matrix.fromLists'
 toLists :: (Element t) => Matrix t -> [[t]]
@@ -207,11 +212,11 @@ createMatrix ord r c = do
 {- | Creates a matrix from a vector by grouping the elements in rows with the desired number of columns. (GNU-Octave groups by columns. To do it you can define @reshapeF r = trans . reshape r@
 where r is the desired number of rows.)
 
-@\> reshape 4 ('fromList' [1..12])
+>>> reshape 4 (fromList [1..12])
 (3><4)
  [ 1.0,  2.0,  3.0,  4.0
  , 5.0,  6.0,  7.0,  8.0
- , 9.0, 10.0, 11.0, 12.0 ]@
+ , 9.0, 10.0, 11.0, 12.0 ]
 
 -}
 reshape :: Storable t => Int -> Vector t -> Matrix t
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index f0bd9aa..415c972 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -113,16 +113,20 @@ inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r
 
 {- | extracts the Vector elements to a list
 
-@> toList (linspace 5 (1,10))
-[1.0,3.25,5.5,7.75,10.0]@
+>>> toList (linspace 5 (1,10))
+[1.0,3.25,5.5,7.75,10.0]
 
 -}
 toList :: Storable a => Vector a -> [a]
 toList v = safeRead v $ peekArray (dim v)
 
-{- | An alternative to 'fromList' with explicit dimension. The input
+{- | Create a vector from a list of elements and explicit dimension. The input
      list is explicitly truncated if it is too long, so it may safely
      be used, for instance, with infinite lists.
+
+>>> 5 |> [1..]
+fromList [1.0,2.0,3.0,4.0,5.0]
+
 -}
 (|>) :: (Storable a) => Int -> [a] -> Vector a
 infixl 9 |>
@@ -159,8 +163,8 @@ at v n
 
 {- | takes a number of consecutive elements from a Vector
 
-@> subVector 2 3 (fromList [1..10])
-3 |> [3.0,4.0,5.0]@
+>>> subVector 2 3 (fromList [1..10])
+fromList [3.0,4.0,5.0]
 
 -}
 subVector :: Storable t => Int       -- ^ index of the starting element
@@ -172,8 +176,8 @@ subVector = Vector.slice
 
 {- | Reads a vector position:
 
-@> fromList [0..9] \@\> 7
-7.0@
+>>> fromList [0..9] @> 7
+7.0
 
 -}
 (@>) :: Storable t => Vector t -> Int -> t
@@ -183,8 +187,8 @@ infixl 9 @>
 
 {- | concatenate a list of vectors
 
-@> vjoin [fromList [1..5], constant 1 3]
-8 |> [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]@
+>>> vjoin [fromList [1..5::Double], konst 1 3]
+fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]
 
 -}
 vjoin :: Storable t => [Vector t] -> Vector t
@@ -205,8 +209,8 @@ vjoin as = unsafePerformIO $ do
 
 {- | Extract consecutive subvectors of the given sizes.
 
-@> takesV [3,4] (linspace 10 (1,10))
-[3 |> [1.0,2.0,3.0],4 |> [4.0,5.0,6.0,7.0]]@
+>>> takesV [3,4] (linspace 10 (1,10::Double))
+[fromList [1.0,2.0,3.0],fromList [4.0,5.0,6.0,7.0]]
 
 -}
 takesV :: Storable t => [Int] -> Vector t -> [Vector t]
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 5365c1c..f72bd15 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -116,13 +116,10 @@ Single row-column components are automatically expanded to match the
 corresponding common row and column:
 
 @
-> let disp = putStr . dispf 2
-> let vector xs = fromList xs :: Vector Double
-> let diagl = diag . vector
-> let rowm = asRow . vector
-
-> disp $ fromBlocks [[ident 5, 7, rowm[10,20]], [3, diagl[1,2,3], 0]]
+disp = putStr . dispf 2
+@
 
+>>> disp $ fromBlocks [[ident 5, 7, row[10,20]], [3, diagl[1,2,3], 0]]
 8x10
 1  0  0  0  0  7  7  7  10  20
 0  1  0  0  0  7  7  7  10  20
@@ -132,7 +129,6 @@ corresponding common row and column:
 3  3  3  3  3  1  0  0   0   0
 3  3  3  3  3  0  2  0   0   0
 3  3  3  3  3  0  0  3   0   0
-@
 
 -}
 fromBlocks :: Element t => [[Matrix t]] -> Matrix t
@@ -163,7 +159,19 @@ adaptBlocks ms = ms' where
 
 --------------------------------------------------------------------------------
 
--- | create a block diagonal matrix
+{- | create a block diagonal matrix
+
+>>>  disp 2 $ diagBlock [konst 1 (2,2), konst 2 (3,5), col [5,7]]
+7x8
+1  1  0  0  0  0  0  0
+1  1  0  0  0  0  0  0
+0  0  2  2  2  2  2  0
+0  0  2  2  2  2  2  0
+0  0  2  2  2  2  2  0
+0  0  0  0  0  0  0  5
+0  0  0  0  0  0  0  7
+
+-}
 diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t
 diagBlock ms = fromBlocks $ zipWith f ms [0..]
   where
@@ -186,12 +194,13 @@ fliprl m = fromColumns . reverse . toColumns $ m
 
 {- | creates a rectangular diagonal matrix:
 
-@> diagRect 7 (fromList [10,20,30]) 4 5 :: Matrix Double
+>>> diagRect 7 (fromList [10,20,30]) 4 5 :: Matrix Double
 (4><5)
  [ 10.0,  7.0,  7.0, 7.0, 7.0
  ,  7.0, 20.0,  7.0, 7.0, 7.0
  ,  7.0,  7.0, 30.0, 7.0, 7.0
- ,  7.0,  7.0,  7.0, 7.0, 7.0 ]@
+ ,  7.0,  7.0,  7.0, 7.0, 7.0 ]
+
 -}
 diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t
 diagRect z v r c = ST.runSTMatrix $ do
@@ -208,10 +217,10 @@ takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (co
 
 {- | An easy way to create a matrix:
 
-@\> (2><3)[1..6]
+>>> (2><3)[2,4,7,-3,11,0]
 (2><3)
- [ 1.0, 2.0, 3.0
- , 4.0, 5.0, 6.0 ]@
+ [  2.0,  4.0, 7.0
+ , -3.0, 11.0, 0.0 ]
 
 This is the format produced by the instances of Show (Matrix a), which
 can also be used for input.
@@ -219,12 +228,11 @@ can also be used for input.
 The input list is explicitly truncated, so that it can
 safely be used with lists that are too long (like infinite lists).
 
-Example:
-
-@\> (2><3)[1..]
+>>> (2><3)[1..]
 (2><3)
  [ 1.0, 2.0, 3.0
- , 4.0, 5.0, 6.0 ]@
+ , 4.0, 5.0, 6.0 ]
+
 
 -}
 (><) :: (Storable a) => Int -> Int -> [a] -> Matrix a
@@ -253,20 +261,35 @@ dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt
 
 {- | Creates a 'Matrix' from a list of lists (considered as rows).
 
-@\> fromLists [[1,2],[3,4],[5,6]]
+>>> fromLists [[1,2],[3,4],[5,6]]
 (3><2)
  [ 1.0, 2.0
  , 3.0, 4.0
- , 5.0, 6.0 ]@
+ , 5.0, 6.0 ]
+
 -}
 fromLists :: Element t => [[t]] -> Matrix t
 fromLists = fromRows . map fromList
 
 -- | creates a 1-row matrix from a vector
+--
+-- >>> asRow (fromList [1..5])
+--  (1><5)
+--   [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
+--
 asRow :: Storable a => Vector a -> Matrix a
 asRow v = reshape (dim v) v
 
 -- | creates a 1-column matrix from a vector
+--
+-- >>> asColumn (fromList [1..5])
+-- (5><1)
+--  [ 1.0
+--  , 2.0
+--  , 3.0
+--  , 4.0
+--  , 5.0 ]
+--
 asColumn :: Storable a => Vector a -> Matrix a
 asColumn v = reshape 1 v
 
@@ -307,12 +330,12 @@ extractRows l m = fromRows $ extract (toRows m) l
 
 {- | creates matrix by repetition of a matrix a given number of rows and columns
 
-@> repmat (ident 2) 2 3 :: Matrix Double
+>>> repmat (ident 2) 2 3
 (4><6)
  [ 1.0, 0.0, 1.0, 0.0, 1.0, 0.0
  , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0
  , 1.0, 0.0, 1.0, 0.0, 1.0, 0.0
- , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ]@
+ , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ]
 
 -}
 repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t
@@ -376,13 +399,14 @@ mk c g = \k -> g (divMod k c)
 
 {- |
 
-@ghci> mapMatrixWithIndexM_ (\\(i,j) v -> printf \"m[%.0f,%.0f] = %.f\\n\" i j v :: IO()) ((2><3)[1 :: Double ..])
+>>> mapMatrixWithIndexM_ (\(i,j) v -> printf "m[%d,%d] = %.f\n" i j v :: IO()) ((2><3)[1 :: Double ..])
 m[0,0] = 1
 m[0,1] = 2
 m[0,2] = 3
 m[1,0] = 4
 m[1,1] = 5
-m[1,2] = 6@
+m[1,2] = 6
+
 -}
 mapMatrixWithIndexM_
   :: (Element a, Num a, Monad m) =>
@@ -393,11 +417,12 @@ mapMatrixWithIndexM_ g m = mapVectorWithIndexM_ (mk c g) . flatten $ m
 
 {- |
 
-@ghci> mapMatrixWithIndexM (\\(i,j) v -> Just $ 100*v + 10*i + j) (ident 3:: Matrix Double)
+>>> mapMatrixWithIndexM (\(i,j) v -> Just $ 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
 Just (3><3)
  [ 100.0,   1.0,   2.0
  ,  10.0, 111.0,  12.0
- ,  20.0,  21.0, 122.0 ]@
+ ,  20.0,  21.0, 122.0 ]
+
 -}
 mapMatrixWithIndexM
   :: (Element a, Storable b, Monad m) =>
@@ -407,11 +432,13 @@ mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . fla
       c = cols m
 
 {- |
-@ghci> mapMatrixWithIndex (\\(i,j) v -> 100*v + 10*i + j) (ident 3:: Matrix Double)
+
+>>> mapMatrixWithIndex (\\(i,j) v -> 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
 (3><3)
  [ 100.0,   1.0,   2.0
  ,  10.0, 111.0,  12.0
- ,  20.0,  21.0, 122.0 ]@
+ ,  20.0,  21.0, 122.0 ]
+
  -}
 mapMatrixWithIndex
   :: (Element a, Storable b) =>
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs
index d2fd351..5339c7e 100644
--- a/lib/Numeric/Container.hs
+++ b/lib/Numeric/Container.hs
@@ -85,8 +85,8 @@ constant = constantD-- about 2x faster
 
 {- | Creates a real vector containing a range of values:
 
-@\> linspace 5 (-3,7)
-5 |> [-3.0,-0.5,2.0,4.5,7.0]@
+>>> linspace 5 (-3,7)
+fromList [-3.0,-0.5,2.0,4.5,7.0]@
 
 Logarithmic spacing can be defined as follows:
 
@@ -105,7 +105,32 @@ cdot u v = udot (conj u) v
 class Contraction a b c | a b -> c, a c -> b, b c -> a
   where
     infixl 7 <>
-    -- | matrix-matrix product, matrix-vector product, unconjugated dot product
+    {- | matrix-matrix product, matrix-vector product, unconjugated dot product
+
+>>> let a = (3><4) [1..] :: Matrix Double
+
+>>> a
+(3><4)
+ [ 1.0,  2.0,  3.0,  4.0
+ , 5.0,  6.0,  7.0,  8.0
+ , 9.0, 10.0, 11.0, 12.0 ]
+
+>>> disp 2 (a <> trans a)
+3x3
+ 30   70  110
+ 70  174  278
+110  278  446
+
+>>> a <> fromList [1,0,2,-1::Double]
+fromList [3.0,11.0,19.0]
+
+>>> fromList [1,2,3::Double] <> a
+fromList [38.0,44.0,50.0,56.0]
+
+>>> fromList [1,i] <> fromList[2*i+1,3]
+1.0 :+ 5.0
+
+-}
     (<>) :: a -> b -> c
 
 instance Product t => Contraction (Vector t) (Vector t) t where
@@ -135,9 +160,14 @@ instance LSDiv Matrix Matrix where
 
 --------------------------------------------------------
 
--- | dot product : @u · v = 'cdot' u v@
---
--- unicode 0x00b7, Alt-Gr .
+{- | dot product : @u · v = 'cdot' u v@
+
+ unicode 0x00b7, Alt-Gr .
+
+>>> fromList [1,i] · fromList[2*i+1,3]
+1.0 :+ (-1.0)
+
+-}
 (·) :: (Container Vector t, Product t) => Vector t -> Vector t -> t
 infixl 7 ·
 u · v = cdot u v
@@ -147,6 +177,16 @@ u · v = cdot u v
 -- bidirectional type inference
 class Konst e d c | d -> c, c -> d
   where
+    -- |
+    -- >>> konst 7 3 :: Vector Float
+    -- fromList [7.0,7.0,7.0]
+    --
+    -- >>> konst i (3::Int,4::Int)
+    -- (3><4)
+    --  [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
+    --  , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
+    --  , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ]
+    --
     konst :: e -> d -> c e
 
 instance Container Vector e => Konst e Int Vector
@@ -161,6 +201,19 @@ instance Container Vector e => Konst e (Int,Int) Matrix
 
 class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f
   where
+    -- |
+    -- >>> build 5 (**2) :: Vector Double
+    -- fromList [0.0,1.0,4.0,9.0,16.0]
+    --
+    -- Hilbert matrix of order N:
+    --
+    -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double
+    -- >>> putStr . dispf 2 $ hilb 3
+    -- 3x3
+    -- 1.00  0.50  0.33
+    -- 0.50  0.33  0.25
+    -- 0.33  0.25  0.20
+    --
     build :: d -> f -> c e
 
 instance Container Vector e => Build Int (e -> e) Vector e
diff --git a/lib/Numeric/ContainerBoot.hs b/lib/Numeric/ContainerBoot.hs
index d61633e..a333489 100644
--- a/lib/Numeric/ContainerBoot.hs
+++ b/lib/Numeric/ContainerBoot.hs
@@ -66,6 +66,11 @@ type instance ArgOf Matrix a = a -> a -> a
 -- | Basic element-by-element functions for numeric containers
 class (Complexable c, Fractional e, Element e) => Container c e where
     -- | create a structure with a single element
+    --
+    -- >>> let v = fromList [1..3::Double]
+    -- >>> v / scalar (norm2 v)
+    -- fromList [0.2672612419124244,0.5345224838248488,0.8017837257372732]
+    --
     scalar      :: e -> c e
     -- | complex conjugate
     conj        :: c e -> c e
@@ -114,8 +119,9 @@ class (Complexable c, Fractional e, Element e) => Container c e where
 
     -- | A more efficient implementation of @cmap (\\x -> if x>0 then 1 else 0)@
     --
-    -- @> step $ linspace 5 (-1,1::Double)
-    -- 5 |> [0.0,0.0,0.0,1.0,1.0]@
+    -- >>> step $ linspace 5 (-1,1::Double)
+    -- 5 |> [0.0,0.0,0.0,1.0,1.0]
+    --
     
     step :: RealElement e => c e -> c e
 
@@ -123,11 +129,12 @@ class (Complexable c, Fractional e, Element e) => Container c e where
     --
     -- Arguments with any dimension = 1 are automatically expanded: 
     --
-    -- @> cond ((1>\<4)[1..]) ((3>\<1)[1..]) 0 100 ((3>\<4)[1..]) :: Matrix Double
+    -- >>> cond ((1><4)[1..]) ((3><1)[1..]) 0 100 ((3><4)[1..]) :: Matrix Double
     -- (3><4)
     -- [ 100.0,   2.0,   3.0,  4.0
     -- ,   0.0, 100.0,   7.0,  8.0
-    -- ,   0.0,   0.0, 100.0, 12.0 ]@
+    -- ,   0.0,   0.0, 100.0, 12.0 ]
+    --
     
     cond :: RealElement e 
          => c e -- ^ a
@@ -139,16 +146,22 @@ class (Complexable c, Fractional e, Element e) => Container c e where
 
     -- | Find index of elements which satisfy a predicate
     --
-    -- @> find (>0) (ident 3 :: Matrix Double)
-    -- [(0,0),(1,1),(2,2)]@
+    -- >>> find (>0) (ident 3 :: Matrix Double)
+    -- [(0,0),(1,1),(2,2)]
+    --
 
     find :: (e -> Bool) -> c e -> [IndexOf c]
 
     -- | Create a structure from an association list
     --
-    -- @> assoc 5 0 [(2,7),(1,3)] :: Vector Double
-    -- 5 |> [0.0,3.0,7.0,0.0,0.0]@
-    
+    -- >>> assoc 5 0 [(3,7),(1,4)] :: Vector Double
+    -- fromList [0.0,4.0,0.0,7.0,0.0]
+    --
+    -- >>> assoc (2,3) 0 [((0,2),7),((1,0),2*i-3)] :: Matrix (Complex Double)
+    -- (2><3)
+    --  [    0.0 :+ 0.0, 0.0 :+ 0.0, 7.0 :+ 0.0
+    --  , (-3.0) :+ 2.0, 0.0 :+ 0.0, 0.0 :+ 0.0 ]
+    --
     assoc :: IndexOf c        -- ^ size
           -> e                -- ^ default value
           -> [(IndexOf c, e)] -- ^ association list
@@ -156,13 +169,19 @@ class (Complexable c, Fractional e, Element e) => Container c e where
 
     -- | Modify a structure using an update function
     --
-    -- @> accum (ident 5) (+) [((1,1),5),((0,3),3)] :: Matrix Double
+    -- >>> accum (ident 5) (+) [((1,1),5),((0,3),3)] :: Matrix Double
     -- (5><5)
     --  [ 1.0, 0.0, 0.0, 3.0, 0.0
     --  , 0.0, 6.0, 0.0, 0.0, 0.0
     --  , 0.0, 0.0, 1.0, 0.0, 0.0
     --  , 0.0, 0.0, 0.0, 1.0, 0.0
-    --  , 0.0, 0.0, 0.0, 0.0, 1.0 ]@
+    --  , 0.0, 0.0, 0.0, 0.0, 1.0 ]
+    --
+    -- computation of histogram:
+    --
+    -- >>> accum (konst 0 7) (+) (map (flip (,) 1) [4,5,4,1,5,2,5]) :: Vector Double
+    -- fromList [0.0,1.0,1.0,0.0,2.0,3.0,0.0]
+    --
     
     accum :: c e              -- ^ initial structure
           -> (e -> e -> e)    -- ^ update function
@@ -382,11 +401,12 @@ vXm v m = flatten $ (asRow v) `mXm` m
 
 {- | Outer product of two vectors.
 
-@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
+>>> fromList [1,2,3] `outer` fromList [5,2,3]
 (3><3)
  [  5.0, 2.0, 3.0
  , 10.0, 4.0, 6.0
- , 15.0, 6.0, 9.0 ]@
+ , 15.0, 6.0, 9.0 ]
+
 -}
 outer :: (Product t) => Vector t -> Vector t -> Matrix t
 outer u v = asColumn u `multiply` asRow v
@@ -402,7 +422,7 @@ m2=(4><3)
  ,  7.0,  8.0,  9.0
  , 10.0, 11.0, 12.0 ]@
 
-@\> kronecker m1 m2
+>>> kronecker m1 m2
 (8><9)
  [  1.0,  2.0,  3.0,   2.0,   4.0,   6.0,  0.0,  0.0,  0.0
  ,  4.0,  5.0,  6.0,   8.0,  10.0,  12.0,  0.0,  0.0,  0.0
@@ -411,7 +431,8 @@ m2=(4><3)
  ,  0.0,  0.0,  0.0,  -1.0,  -2.0,  -3.0,  3.0,  6.0,  9.0
  ,  0.0,  0.0,  0.0,  -4.0,  -5.0,  -6.0, 12.0, 15.0, 18.0
  ,  0.0,  0.0,  0.0,  -7.0,  -8.0,  -9.0, 21.0, 24.0, 27.0
- ,  0.0,  0.0,  0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
+ ,  0.0,  0.0,  0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]
+
 -}
 kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t
 kronecker a b = fromBlocks
diff --git a/lib/Numeric/GSL/Differentiation.hs b/lib/Numeric/GSL/Differentiation.hs
index d2a332c..93c5007 100644
--- a/lib/Numeric/GSL/Differentiation.hs
+++ b/lib/Numeric/GSL/Differentiation.hs
@@ -55,11 +55,11 @@ foreign import ccall safe "gsl-aux.h deriv"
 
 {- | Adaptive central difference algorithm, /gsl_deriv_central/. For example:
 
-> > let deriv = derivCentral 0.01 
-> > deriv sin (pi/4)
->(0.7071067812000676,1.0600063101654055e-10)
-> > cos (pi/4)
->0.7071067811865476 
+>>> let deriv = derivCentral 0.01
+>>> deriv sin (pi/4)
+(0.7071067812000676,1.0600063101654055e-10)
+>>> cos (pi/4)
+0.7071067811865476
 
 -}
 derivCentral :: Double                  -- ^ initial step size
diff --git a/lib/Numeric/GSL/Fitting.hs b/lib/Numeric/GSL/Fitting.hs
index 6343b76..c4f3a91 100644
--- a/lib/Numeric/GSL/Fitting.hs
+++ b/lib/Numeric/GSL/Fitting.hs
@@ -13,7 +13,8 @@ Nonlinear Least-Squares Fitting
 
 The example program in the GSL manual (see examples/fitting.hs):
 
-@dat = [
+@
+dat = [
  ([0.0],([6.0133918608118675],0.1)),
  ([1.0],([5.5153769909966535],0.1)),
  ([2.0],([5.261094606015287],0.1)),
@@ -25,8 +26,9 @@ expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]
 expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]
 
 (sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
+@
 
-\> path
+>>> path
 (6><5)
  [ 1.0,  76.45780563978782, 1.6465931240727802, 1.8147715267618197e-2, 0.6465931240727797
  , 2.0, 37.683816318260355,  2.858760367632973,  8.092094813253975e-2, 1.4479636296208662
@@ -34,10 +36,10 @@ expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]
  , 4.0,  5.630494933603935,  5.021755718065913,   0.10287787128056883, 1.0338835440862608
  , 5.0,  5.443976278682909,  5.045204331329302,   0.10405523433131504,  1.019416067207375
  , 6.0, 5.4439736648994685,  5.045357818922331,   0.10404905846029407, 1.0192487112786812 ]
-\> sol
+>>> sol
 [(5.045357818922331,6.027976702418132e-2),
 (0.10404905846029407,3.157045047172834e-3),
-(1.0192487112786812,3.782067731353722e-2)]@
+(1.0192487112786812,3.782067731353722e-2)]
 
 -}
 -----------------------------------------------------------------------------
diff --git a/lib/Numeric/GSL/Fourier.hs b/lib/Numeric/GSL/Fourier.hs
index 4ef19b3..86aedd6 100644
--- a/lib/Numeric/GSL/Fourier.hs
+++ b/lib/Numeric/GSL/Fourier.hs
@@ -35,8 +35,8 @@ foreign import ccall unsafe "gsl-aux.h fft" c_fft ::  CInt -> TCVCV
 
 {- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave.
 
-@> fft ('fromList' [1,2,3,4])
-vector (4) [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]@
+>>> fft (fromList [1,2,3,4])
+fromList [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]
 
 -}
 fft :: Vector (Complex Double) -> Vector (Complex Double)
diff --git a/lib/Numeric/GSL/Integration.hs b/lib/Numeric/GSL/Integration.hs
index 7c1cdb9..5f0a415 100644
--- a/lib/Numeric/GSL/Integration.hs
+++ b/lib/Numeric/GSL/Integration.hs
@@ -35,16 +35,16 @@ eps = 1e-12
 
 {- | conversion of Haskell functions into function pointers that can be used in the C side
 -}
-foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) 
+foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double))
 
 --------------------------------------------------------------------
 {- | Numerical integration using /gsl_integration_qags/ (adaptive integration with singularities). For example:
 
-@\> let quad = integrateQAGS 1E-9 1000 
-\> let f a x = x**(-0.5) * log (a*x)
-\> quad (f 1) 0 1
-(-3.999999999999974,4.871658632055187e-13)@
- 
+>>> let quad = integrateQAGS 1E-9 1000
+>>> let f a x = x**(-0.5) * log (a*x)
+>>> quad (f 1) 0 1
+(-3.999999999999974,4.871658632055187e-13)
+
 -}
 
 integrateQAGS :: Double               -- ^ precision (e.g. 1E-9)
@@ -56,7 +56,7 @@ integrateQAGS :: Double               -- ^ precision (e.g. 1E-9)
 integrateQAGS prec n f a b = unsafePerformIO $ do
     r <- malloc
     e <- malloc
-    fp <- mkfun (\x _ -> f x) 
+    fp <- mkfun (\x _ -> f x)
     c_integrate_qags fp a b eps prec (fromIntegral n) r e // check "integrate_qags"
     vr <- peek r
     ve <- peek e
@@ -73,10 +73,10 @@ foreign import ccall safe "integrate_qags" c_integrate_qags
 -----------------------------------------------------------------
 {- | Numerical integration using /gsl_integration_qng/ (useful for fast integration of smooth functions). For example:
 
-@\> let quad = integrateQNG 1E-6 
-\> quad (\\x -> 4\/(1+x*x)) 0 1 
-(3.141592653589793,3.487868498008632e-14)@
- 
+>>> let quad = integrateQNG 1E-6
+>>> quad (\x -> 4/(1+x*x)) 0 1
+(3.141592653589793,3.487868498008632e-14)
+
 -}
 
 integrateQNG :: Double               -- ^ precision (e.g. 1E-9)
@@ -87,7 +87,7 @@ integrateQNG :: Double               -- ^ precision (e.g. 1E-9)
 integrateQNG prec f a b = unsafePerformIO $ do
     r <- malloc
     e <- malloc
-    fp <- mkfun (\x _ -> f x) 
+    fp <- mkfun (\x _ -> f x)
     c_integrate_qng fp a b eps prec r e  // check "integrate_qng"
     vr <- peek r
     ve <- peek e
@@ -102,14 +102,14 @@ foreign import ccall safe "integrate_qng" c_integrate_qng
     -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt
 
 --------------------------------------------------------------------
-{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS). 
+{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS).
 For example:
 
-@\> let quad = integrateQAGI 1E-9 1000 
-\> let f a x = exp(-a * x^2)
-\> quad (f 0.5) 
-(2.5066282746310002,6.229215880648858e-11)@
- 
+>>> let quad = integrateQAGI 1E-9 1000
+>>> let f a x = exp(-a * x^2)
+>>> quad (f 0.5)
+(2.5066282746310002,6.229215880648858e-11)
+
 -}
 
 integrateQAGI :: Double               -- ^ precision (e.g. 1E-9)
@@ -119,7 +119,7 @@ integrateQAGI :: Double               -- ^ precision (e.g. 1E-9)
 integrateQAGI prec n f = unsafePerformIO $ do
     r <- malloc
     e <- malloc
-    fp <- mkfun (\x _ -> f x) 
+    fp <- mkfun (\x _ -> f x)
     c_integrate_qagi fp eps prec (fromIntegral n) r e // check "integrate_qagi"
     vr <- peek r
     ve <- peek e
@@ -134,14 +134,14 @@ foreign import ccall safe "integrate_qagi" c_integrate_qagi
     -> CInt -> Ptr Double -> Ptr Double -> IO CInt
 
 --------------------------------------------------------------------
-{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf). 
+{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf).
 For example:
 
-@\> let quad = integrateQAGIU 1E-9 1000 
-\> let f a x = exp(-a * x^2)
-\> quad (f 0.5) 0
-(1.2533141373155001,3.114607940324429e-11)@
- 
+>>> let quad = integrateQAGIU 1E-9 1000
+>>> let f a x = exp(-a * x^2)
+>>> quad (f 0.5) 0
+(1.2533141373155001,3.114607940324429e-11)
+
 -}
 
 integrateQAGIU :: Double               -- ^ precision (e.g. 1E-9)
@@ -152,7 +152,7 @@ integrateQAGIU :: Double               -- ^ precision (e.g. 1E-9)
 integrateQAGIU prec n f a = unsafePerformIO $ do
     r <- malloc
     e <- malloc
-    fp <- mkfun (\x _ -> f x) 
+    fp <- mkfun (\x _ -> f x)
     c_integrate_qagiu fp a eps prec (fromIntegral n) r e // check "integrate_qagiu"
     vr <- peek r
     ve <- peek e
@@ -167,14 +167,14 @@ foreign import ccall safe "integrate_qagiu" c_integrate_qagiu
     -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
 
 --------------------------------------------------------------------
-{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b). 
+{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b).
 For example:
 
-@\> let quad = integrateQAGIL 1E-9 1000 
-\> let f a x = exp(-a * x^2)
-\> quad (f 0.5) 0 
-(1.2533141373155001,3.114607940324429e-11)@
- 
+>>> let quad = integrateQAGIL 1E-9 1000
+>>> let f a x = exp(-a * x^2)
+>>> quad (f 0.5) 0
+(1.2533141373155001,3.114607940324429e-11)
+
 -}
 
 integrateQAGIL :: Double               -- ^ precision (e.g. 1E-9)
@@ -185,7 +185,7 @@ integrateQAGIL :: Double               -- ^ precision (e.g. 1E-9)
 integrateQAGIL prec n f b = unsafePerformIO $ do
     r <- malloc
     e <- malloc
-    fp <- mkfun (\x _ -> f x) 
+    fp <- mkfun (\x _ -> f x)
     c_integrate_qagil fp b eps prec (fromIntegral n) r e // check "integrate_qagil"
     vr <- peek r
     ve <- peek e
@@ -212,10 +212,10 @@ routines in QUADPACK, yet fails less often for difficult integrands.@
 
 For example:
 
-@\> let quad = integrateCQUAD 1E-12 1000 
-\> let f a x = exp(-a * x^2)
-\> quad (f 0.5) 2 5
-(5.7025405463957006e-2,9.678874441303705e-16,95)@
+>>> let quad = integrateCQUAD 1E-12 1000
+>>> let f a x = exp(-a * x^2)
+>>> quad (f 0.5) 2 5
+(5.7025405463957006e-2,9.678874441303705e-16,95)
 
 Unlike other quadrature methods, integrateCQUAD also returns the
 number of function evaluations required.
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs
index 7ccca81..1879dab 100644
--- a/lib/Numeric/GSL/Minimization.hs
+++ b/lib/Numeric/GSL/Minimization.hs
@@ -16,15 +16,15 @@ Minimization of a multidimensional function using some of the algorithms describ
 The example in the GSL manual:
 
 @
-
 f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
 
 main = do
     let (s,p) = minimize NMSimplex2 1E-2 30 [1,1] f [5,7]
     print s
     print p
+@
 
-\> main
+>>> main
 [0.9920430849306288,1.9969168063253182]
  0.000  512.500  1.130  6.500  5.000
  1.000  290.625  1.409  5.250  4.000
@@ -33,7 +33,6 @@ main = do
  ...
 22.000   30.001  0.013  0.992  1.997
 23.000   30.001  0.008  0.992  1.997
-@
 
 The path to the solution can be graphically shown by means of:
 
diff --git a/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs
index ab037bd..9a29085 100644
--- a/lib/Numeric/GSL/ODE.hs
+++ b/lib/Numeric/GSL/ODE.hs
@@ -13,9 +13,10 @@ Solution of ordinary differential equation (ODE) initial value problems.
 
 A simple example:
 
-@import Numeric.GSL
+@
+import Numeric.GSL.ODE
 import Numeric.LinearAlgebra
-import Graphics.Plot
+import Numeric.LinearAlgebra.Util(mplot)
 
 xdot t [x,v] = [v, -0.95*x - 0.1*v]
 
@@ -23,7 +24,8 @@ ts = linspace 100 (0,20 :: Double)
 
 sol = odeSolve xdot [10,0] ts
 
-main = mplot (ts : toColumns sol)@
+main = mplot (ts : toColumns sol)
+@
 
 -}
 -----------------------------------------------------------------------------
diff --git a/lib/Numeric/GSL/Polynomials.hs b/lib/Numeric/GSL/Polynomials.hs
index 694c003..290c615 100644
--- a/lib/Numeric/GSL/Polynomials.hs
+++ b/lib/Numeric/GSL/Polynomials.hs
@@ -27,22 +27,22 @@ import System.IO.Unsafe (unsafePerformIO)
 import Foreign.C.Types (CInt(..))
 #endif
 
-{- | Solution of general polynomial equations, using /gsl_poly_complex_solve/. For example,
-     the three solutions of x^3 + 8 = 0
+{- | Solution of general polynomial equations, using /gsl_poly_complex_solve/. 
+
+For example, the three solutions of x^3 + 8 = 0
+
+>>> polySolve [8,0,0,1]
+[(-2.0) :+ 0.0,1.0 :+ 1.7320508075688776,1.0 :+ (-1.7320508075688776)]
 
-@\> polySolve [8,0,0,1]
-[(-1.9999999999999998) :+ 0.0,
- 1.0 :+ 1.732050807568877,
- 1.0 :+ (-1.732050807568877)]@
 
 The example in the GSL manual: To find the roots of x^5 -1 = 0:
 
-@\> polySolve [-1, 0, 0, 0, 0, 1]
-[(-0.8090169943749475) :+ 0.5877852522924731,
-(-0.8090169943749475) :+ (-0.5877852522924731),
-0.30901699437494734 :+ 0.9510565162951536,
-0.30901699437494734 :+ (-0.9510565162951536),
-1.0 :+ 0.0]@
+>>> polySolve [-1, 0, 0, 0, 0, 1]
+[(-0.8090169943749472) :+ 0.5877852522924731,
+(-0.8090169943749472) :+ (-0.5877852522924731),
+0.30901699437494756 :+ 0.9510565162951535,
+0.30901699437494756 :+ (-0.9510565162951535),
+1.0000000000000002 :+ 0.0]
 
 -}  
 polySolve :: [Double] -> [Complex Double]
diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs
index 6da15e5..9d561c4 100644
--- a/lib/Numeric/GSL/Root.hs
+++ b/lib/Numeric/GSL/Root.hs
@@ -13,33 +13,23 @@ Multidimensional root finding.
 
 The example in the GSL manual:
 
-@import Numeric.GSL
-import Numeric.LinearAlgebra(format)
-import Text.Printf(printf)
-
-rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
-
-disp = putStrLn . format \"  \" (printf \"%.3f\")
-
-main = do
-    let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
-    print sol
-    disp path
-
-\> main
+>>> let rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
+>>> let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
+>>> sol
 [1.0,1.0]
- 0.000  -10.000  -5.000  11.000  -1050.000
- 1.000   -3.976  24.827   4.976     90.203
+>>> disp 3 path
+11x5
+ 1.000  -10.000  -5.000  11.000  -1050.000
  2.000   -3.976  24.827   4.976     90.203
  3.000   -3.976  24.827   4.976     90.203
- 4.000   -1.274  -5.680   2.274    -73.018
+ 4.000   -3.976  24.827   4.976     90.203
  5.000   -1.274  -5.680   2.274    -73.018
- 6.000    0.249   0.298   0.751      2.359
+ 6.000   -1.274  -5.680   2.274    -73.018
  7.000    0.249   0.298   0.751      2.359
- 8.000    1.000   0.878  -0.000     -1.218
- 9.000    1.000   0.989  -0.000     -0.108
-10.000    1.000   1.000   0.000      0.000
-@
+ 8.000    0.249   0.298   0.751      2.359
+ 9.000    1.000   0.878  -0.000     -1.218
+10.000    1.000   0.989  -0.000     -0.108
+11.000    1.000   1.000   0.000      0.000
 
 -}
 -----------------------------------------------------------------------------
diff --git a/lib/Numeric/IO.hs b/lib/Numeric/IO.hs
index dacfa8b..57275ac 100644
--- a/lib/Numeric/IO.hs
+++ b/lib/Numeric/IO.hs
@@ -39,29 +39,27 @@ format sep f m = table sep . map (map f) . toLists $ m
 
 {- | Show a matrix with \"autoscaling\" and a given number of decimal places.
 
-@disp = putStr . disps 2
-
-\> disp $ 120 * (3><4) [1..]
+>>> putStr . disps 2 $ 120 * (3><4) [1..]
 3x4  E3
  0.12  0.24  0.36  0.48
  0.60  0.72  0.84  0.96
  1.08  1.20  1.32  1.44
-@
+
 -}
 disps :: Int -> Matrix Double -> String
 disps d x = sdims x ++ "  " ++ formatScaled d x
 
 {- | Show a matrix with a given number of decimal places.
 
-@disp = putStr . dispf 3
+>>> dispf 2 (1/3 + ident 3)
+"3x3\n1.33  0.33  0.33\n0.33  1.33  0.33\n0.33  0.33  1.33\n"
+
+>>> putStr . dispf 2 $ (3><4)[1,1.5..]
+3x4
+1.00  1.50  2.00  2.50
+3.00  3.50  4.00  4.50
+5.00  5.50  6.00  6.50
 
-\> disp (1/3 + ident 4)
-4x4
-1.333  0.333  0.333  0.333
-0.333  1.333  0.333  0.333
-0.333  0.333  1.333  0.333
-0.333  0.333  0.333  1.333
-@
 -}
 dispf :: Int -> Matrix Double -> String
 dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x
@@ -81,11 +79,9 @@ formatScaled dec t = "E"++show o++"\n" ++ ss
 
 {- | Show a vector using a function for showing matrices.
 
-@disp = putStr . vecdisp ('dispf' 2)
-
-\> disp ('linspace' 10 (0,1))
+>>> putStr . vecdisp (dispf 2) $ linspace 10 (0,1)
 10 |> 0.00  0.11  0.22  0.33  0.44  0.56  0.67  0.78  0.89  1.00
-@
+
 -}
 vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String
 vecdisp f v
@@ -94,7 +90,12 @@ vecdisp f v
     . f . trans . reshape 1
     $ v
 
--- | Tool to display matrices with latex syntax.
+{- | Tool to display matrices with latex syntax.
+
+>>>  latexFormat "bmatrix" (dispf 2 $ ident 2)
+"\\begin{bmatrix}\n1  &  0\n\\\\\n0  &  1\n\\end{bmatrix}"
+
+-}
 latexFormat :: String -- ^ type of braces: \"matrix\", \"bmatrix\", \"pmatrix\", etc.
             -> String -- ^ Formatted matrix, with elements separated by spaces and newlines
             -> String
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 7223cd9..7c1c032 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -371,19 +371,21 @@ pinv = pinvTol 1
 
 {- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value.
 
-@\> let m = 'fromLists' [[1,0,    0]
-                    ,[0,1,    0]
-                    ,[0,0,1e-10]]
-\  --
-\> 'pinv' m
+@
+m = (3><3) [ 1, 0,    0
+           , 0, 1,    0
+           , 0, 0, 1e-10] :: Matrix Double
+@
+
+>>> pinv m
 1. 0.           0.
 0. 1.           0.
 0. 0. 10000000000.
-\  --
-\> pinvTol 1E8 m
+
+>>> pinvTol 1E8 m
 1. 0. 0.
 0. 1. 0.
-0. 0. 1.@
+0. 0. 1.
 
 -}
 
@@ -598,10 +600,11 @@ It only works with invertible matrices that have a real solution. For diagonaliz
 @m = (2><2) [4,9
            ,0,4] :: Matrix Double@
 
-@\>sqrtm m
+>>> sqrtm m
 (2><2)
  [ 2.0, 2.25
- , 0.0,  2.0 ]@
+ , 0.0,  2.0 ]
+
 -}
 sqrtm ::  Field t => Matrix t -> Matrix t
 sqrtm = sqrtmInv
diff --git a/lib/Numeric/LinearAlgebra/Data.hs b/lib/Numeric/LinearAlgebra/Data.hs
index 6e5dcef..a3639d5 100644
--- a/lib/Numeric/LinearAlgebra/Data.hs
+++ b/lib/Numeric/LinearAlgebra/Data.hs
@@ -12,6 +12,8 @@ Stability   :  provisional
 
 module Numeric.LinearAlgebra.Data(
     -- * Vector
+    -- | 1D arrays are storable vectors from the vector package.
+    
     Vector, (|>), dim, (@>),
     
     -- * Matrix
@@ -49,13 +51,13 @@ module Numeric.LinearAlgebra.Data(
 
     -- * Products
     
-    (<>), (·), scale, outer, kronecker, cross,
+    (<>), (·), outer, kronecker, cross,
     sumElements, prodElements, absSum,
     optimiseMult,
     
     corr, conv, corrMin, corr2, conv2,
     
-    LSDiv(..),
+    (<\>),
     
     -- * Random arrays
 
@@ -80,7 +82,8 @@ module Numeric.LinearAlgebra.Data(
     module Data.Complex,
 
     -- * Misc
-    meanCov, arctan2,
+    scale, meanCov, arctan2,
+    rows, cols,
     separable,
     fromArray2D
 
diff --git a/lib/Numeric/LinearAlgebra/Data/Devel.hs b/lib/Numeric/LinearAlgebra/Data/Devel.hs
index 6358d1f..88c980c 100644
--- a/lib/Numeric/LinearAlgebra/Data/Devel.hs
+++ b/lib/Numeric/LinearAlgebra/Data/Devel.hs
@@ -51,13 +51,13 @@ module Numeric.LinearAlgebra.Data.Devel(
     liftMatrix, liftMatrix2, liftMatrix2Auto,
 
     -- * Misc
-    Element, Container
+    Element, Container, Product, Contraction, LSDiv
 ) where
 
 import Data.Packed.Foreign
 import Data.Packed.Development
 import Data.Packed.ST
-import Numeric.Container(Container)
+import Numeric.Container(Container,Contraction,LSDiv,Product)
 import Data.Packed
 
 
diff --git a/lib/Numeric/LinearAlgebra/Util.hs b/lib/Numeric/LinearAlgebra/Util.hs
index be01bc1..21b6188 100644
--- a/lib/Numeric/LinearAlgebra/Util.hs
+++ b/lib/Numeric/LinearAlgebra/Util.hs
@@ -13,7 +13,7 @@ Stability   :  provisional
 {-# OPTIONS_HADDOCK hide #-}
 
 module Numeric.LinearAlgebra.Util(
-    
+
     -- * Convenience functions
     size, disp,
     zeros, ones,
@@ -65,8 +65,16 @@ import Numeric.LinearAlgebra.Util.Convolution
 import Graphics.Plot
 
 
+{- | print a real matrix with given number of digits after the decimal point
+
+>>> disp 5 $ ident 2 / 3
+2x2
+0.33333  0.00000
+0.00000  0.33333
+
+-}
 disp :: Int -> Matrix Double -> IO ()
--- ^ show a matrix with given number of digits after the decimal point
+
 disp n = putStrLn . dispf n
 
 -- | pseudorandom matrix with uniform elements between 0 and 1
@@ -82,7 +90,16 @@ randm d r c = do
 rand :: Int -> Int -> IO (Matrix Double)
 rand = randm Uniform
 
--- | pseudorandom matrix with normal elements
+{- | pseudorandom matrix with normal elements
+
+>>> x <- randn 3 5
+>>> disp 3 x
+3x5
+0.386  -1.141   0.491  -0.510   1.512
+0.069  -0.919   1.022  -0.181   0.745
+0.313  -0.670  -0.097  -1.575  -0.583
+
+-}
 randn :: Int -> Int -> IO (Matrix Double)
 randn = randm Gaussian
 
@@ -115,9 +132,17 @@ infixl 3 &
 (&) :: Vector Double -> Vector Double -> Vector Double
 a & b = vjoin [a,b]
 
--- | horizontal concatenation of real matrices
---
--- (0x00a6 broken bar)
+{- | horizontal concatenation of real matrices
+
+ (0x00a6 broken bar)
+
+>>> ident 3 ¦ konst 7 (3,4)
+(3><7)
+ [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0
+ , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0
+ , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ]
+
+-}
 infixl 3 ¦
 (¦) :: Matrix Double -> Matrix Double -> Matrix Double
 a ¦ b = fromBlocks [[a,b]]
@@ -141,7 +166,15 @@ row = asRow . fromList
 col :: [Double] -> Matrix Double
 col = asColumn . fromList
 
--- | extract selected rows
+{- | extract selected rows
+
+>>> (20><4) [1..] ? [2,1,1]
+(3><4)
+ [ 9.0, 10.0, 11.0, 12.0
+ , 5.0,  6.0,  7.0,  8.0
+ , 5.0,  6.0,  7.0,  8.0 ]
+
+-}
 infixl 9 ?
 (?) :: Element t => Matrix t -> [Int] -> Matrix t
 (?) = flip extractRows
@@ -172,7 +205,7 @@ norm = pnorm PNorm2
 unitary :: Vector Double -> Vector Double
 unitary v = v / scalar (norm v)
 
--- | (rows &&& cols)
+-- | ('rows' &&& 'cols')
 size :: Matrix t -> (Int, Int)
 size m = (rows m, cols m)