diff --git a/hmatrix-gsl-stats.cabal b/hmatrix-gsl-stats.cabal
--- a/hmatrix-gsl-stats.cabal
+++ b/hmatrix-gsl-stats.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix-gsl-stats
-Version:            0.4.1.4
+Version:            0.4.1.5
 License:            BSD3
 License-file:       LICENSE
 Copyright:          (c) A.V.H. McPhail 2010, 2011, 2013, 2015, 2016
@@ -37,7 +37,7 @@
     Build-Depends:      base >= 4 && < 5, binary,
                         vector,
                         storable-complex,
-                        hmatrix >= 0.17 && < 0.18
+                        hmatrix >= 0.18
 			
     Default-Extensions: ForeignFunctionInterface
 
diff --git a/lib/Numeric/GSL/Distribution/Continuous.hs b/lib/Numeric/GSL/Distribution/Continuous.hs
--- a/lib/Numeric/GSL/Distribution/Continuous.hs
+++ b/lib/Numeric/GSL/Distribution/Continuous.hs
@@ -1,8 +1,9 @@
 {-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Distribution.Continuous
--- Copyright   :  (c) A. V. H. McPhail 2010, 2015
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2015, 2016
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -60,7 +61,7 @@
 
 -----------------------------------------------------------------------------
 
-infixl 1 #
+infixr 1 #
 a # b = applyRaw a b
 {-# INLINE (#) #-}
 
@@ -128,7 +129,7 @@
             -> Vector Double   -- ^ result
 random_0p_v  d s l = unsafePerformIO $ do
    r <- createVector l
-   (distribution_random_zero_param_v (fromIntegral s) (fromei d)) # r #| "random_0p_v"
+   (r # id) (distribution_random_zero_param_v (fromIntegral s) (fromei d)) #| "random_0p_v"
    return r
 
 -- | probability of a variate take a value outside the argument
@@ -179,7 +180,7 @@
             -> Vector Double   -- ^ result
 random_1p_v d s p l = unsafePerformIO $ do
    r <- createVector l
-   (distribution_random_one_param_v (fromIntegral s) (fromei d) p) # r #| "random_1p_v"
+   (r # id) (distribution_random_one_param_v (fromIntegral s) (fromei d) p) #| "random_1p_v"
    return r
 
 -- | probability of a variate take a value outside the argument
@@ -229,7 +230,7 @@
             -> Vector Double   -- ^ result
 random_2p_v d s p1 p2 l = unsafePerformIO $ do
    r <- createVector l
-   (distribution_random_two_param_v (fromIntegral s) (fromei d) p1 p2) # r #| "random_2p_v"
+   (r # id) (distribution_random_two_param_v (fromIntegral s) (fromei d) p1 p2) #| "random_2p_v"
    return r
 
 -- | probability of a variate take a value outside the argument
@@ -295,7 +296,7 @@
             -> Vector Double   -- ^ result
 random_3p_v d s p1 p2 p3 l = unsafePerformIO $ do
    r <- createVector l
-   (distribution_random_three_param_v (fromIntegral s) (fromei d) p1 p2 p3) # r #| "random_3p_v"
+   (r # id) (distribution_random_three_param_v (fromIntegral s) (fromei d) p1 p2 p3) #| "random_3p_v"
    return r
 
 -- | probability of a variate take a value outside the argument
@@ -325,7 +326,7 @@
           -> Vector Double   -- ^ result
 random_mp d s p = unsafePerformIO $ do
                   r <- createVector $ size p
-                  (distribution_random_multi_param (fromIntegral s) (fromei d)) # p # r #| "random_mp"
+                  (p # r # id) (distribution_random_multi_param (fromIntegral s) (fromei d)) #| "random_mp"
                   return r
 
 -- | draw a sample from a multi parameter distribution
@@ -336,7 +337,7 @@
 random_mp_s (RNG rng) d p = do
   let l = size p
   r <- createVector l
-  withForeignPtr rng $ \rg -> (distribution_random_multi_param_s rg (fromei d)) # p # r #| "random_mp_s"
+  withForeignPtr rng $ \rg -> (p # r # id) (distribution_random_multi_param_s rg (fromei d)) #| "random_mp_s"
   return r
   
 -- | probability of a variate take a value outside the argument
@@ -351,7 +352,7 @@
     where density_only f' d' p' q' = if f' /= Density
                                               then error "distribution has no CDF"
                                               else alloca $ \r -> do
-                                                                  (distribution_dist_multi_param (fromei f') (fromei d') r) # p' # q' #| "density_mp"
+                                                                  (p' # q' # id) (distribution_dist_multi_param (fromei f') (fromei d') r) #| "density_mp"
                                                                   r' <- peek r
                                                                   return r'
 
@@ -402,7 +403,7 @@
 random_biv_v d s p1 p2 p3 l = unsafePerformIO $ do
    r1 <- createVector l
    r2 <- createVector l
-   (distribution_random_bivariate_v (fromIntegral s) (fromei d) p1 p2 p3) # r1 # r2 #| "random_biv_v"
+   (r1 # r2 # id) (distribution_random_bivariate_v (fromIntegral s) (fromei d) p1 p2 p3) #| "random_biv_v"
    return (r1,r2)
 
 -- | probability of a variate take a value outside the argument
@@ -433,7 +434,7 @@
                  -> Vector Double -- result
 spherical_vector s vs = unsafePerformIO $ do
                         r <- createVector vs
-                        (distribution_spherical_vector (fromIntegral s)) # r #| "spherical_vector"
+                        (r # id) (distribution_spherical_vector (fromIntegral s)) #| "spherical_vector"
                         return r
 
 foreign import ccall "distribution-aux.h spherical_vector" distribution_spherical_vector :: CInt -> CInt -> Ptr Double -> IO CInt
diff --git a/lib/Numeric/GSL/Distribution/Discrete.hs b/lib/Numeric/GSL/Distribution/Discrete.hs
--- a/lib/Numeric/GSL/Distribution/Discrete.hs
+++ b/lib/Numeric/GSL/Distribution/Discrete.hs
@@ -2,7 +2,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Distribution.Discrete
--- Copyright   :  (c) A. V. H. McPhail 2010, 2015
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2015, 2016
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -56,7 +56,7 @@
 
 -----------------------------------------------------------------------------
 
-infixl 1 #
+infixr 1 #
 a # b = applyRaw a b
 {-# INLINE (#) #-}
 
@@ -116,7 +116,7 @@
             -> Vector Word32   -- ^ result
 random_1p_v d s p l = unsafePerformIO $ do
    r <- createVector l
-   (distribution_discrete_one_param_v (fromIntegral s) (fromei d) p) # r #| "random_1p_v"
+   (r # id) (distribution_discrete_one_param_v (fromIntegral s) (fromei d) p) #| "random_1p_v"
    return (map fromIntegral r)
 
 -- | probability of a variate take a value outside the argument
@@ -179,7 +179,7 @@
             -> Vector Word32   -- ^ result
 random_2p_v d s p1 p2 l = unsafePerformIO $ do
    r <- createVector l
-   (distribution_discrete_two_param_v (fromIntegral s) (fromei d) p1 (fromIntegral p2)) # r #| "random_2p_v"
+   (r # id) (distribution_discrete_two_param_v (fromIntegral s) (fromei d) p1 (fromIntegral p2)) #| "random_2p_v"
    return (map fromIntegral r)
 
 -- | probability of a variate take a value outside the argument
@@ -239,7 +239,7 @@
             -> Vector Word32   -- ^ result
 random_3p_v d s p1 p2 p3 l = unsafePerformIO $ do
    r <- createVector l
-   (distribution_discrete_three_param_v (fromIntegral s) (fromei d) (fromIntegral p1) (fromIntegral p2) (fromIntegral p3)) # r #| "random_3p_v"
+   (r # id) (distribution_discrete_three_param_v (fromIntegral s) (fromei d) (fromIntegral p1) (fromIntegral p2) (fromIntegral p3)) #| "random_3p_v"
    return (map fromIntegral r)
 
 -- | probability of a variate take a value outside the argument
@@ -272,7 +272,7 @@
           -> Vector Word32   -- ^ result
 random_mp d s t p = unsafePerformIO $ do
                     r <- createVector $ size p
-                    (distribution_discrete_multi_param (fromIntegral s) (fromei d) (fromIntegral t)) # p # r #| "random_mp"
+                    (p # r # id) (distribution_discrete_multi_param (fromIntegral s) (fromei d) (fromIntegral t)) #| "random_mp"
                     return $ map (\x -> (fromIntegral x) :: Word32) r
 
 -- | draw a sample from a multi-parameter distribution
@@ -283,7 +283,7 @@
             -> IO (Vector Word32) -- ^ result
 random_mp_s (RNG rng) d t p = withForeignPtr rng $ \rg -> do
   r <- createVector $ size p
-  (distribution_discrete_multi_param_s rg (fromei d) (fromIntegral t)) # p # r #| "random_mp_s"
+  (p # r # id) (distribution_discrete_multi_param_s rg (fromei d) (fromIntegral t)) #| "random_mp_s"
   return $ map (\x -> (fromIntegral x) :: Word32) r
 
 -- | probability of a variate take a value outside the argument
@@ -296,11 +296,11 @@
                      case d of
                             Multinomial -> density_only f d p q
     where density_only f' d' p' q' = if f' /= Density
-                                              then error "distribution has no CDF"
-                                              else alloca $ \r -> do
-                                                                  (distribution_dist_multi_param (fromei f') (fromei d') r) # p' # (map (\x -> (fromIntegral x) :: CUInt) q') #| "density_mp"
-                                                                  r' <- peek r
-                                                                  return r'
+                                     then error "distribution has no CDF"
+                                     else alloca $ \r -> do
+                                                         (p' # (map (\x -> (fromIntegral x) :: CUInt) q') # id) (distribution_dist_multi_param (fromei f') (fromei d') r) #| "density_mp"
+                                                         r' <- peek r
+                                                         return r'
 
 foreign import ccall "distribution-aux.h discrete_mp" distribution_discrete_multi_param :: CInt -> CInt -> CUInt -> CInt -> Ptr Double -> CInt -> Ptr CUInt -> IO CInt
 foreign import ccall "distribution-aux.h discrete_mp_s" distribution_discrete_multi_param_s :: RNGHandle -> CInt -> CUInt -> CInt -> Ptr Double -> CInt -> Ptr CUInt -> IO CInt
diff --git a/lib/Numeric/GSL/Fitting/Linear.hs b/lib/Numeric/GSL/Fitting/Linear.hs
--- a/lib/Numeric/GSL/Fitting/Linear.hs
+++ b/lib/Numeric/GSL/Fitting/Linear.hs
@@ -1,7 +1,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Fitting.Linear
--- Copyright   :  (c) A. V. H. McPhail 2010
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2016
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -45,7 +45,7 @@
 
 -----------------------------------------------------------------------------
 
-infixl 1 #
+infixr 1 #
 a # b = applyRaw a b
 {-# INLINE (#) #-}
 
@@ -62,7 +62,7 @@
                        alloca $ \cov00 -> 
                            alloca $ \cov01 -> 
                                alloca $ \cov11 -> do
-                                                  (fitting_linear c0 c1 chi_sq cov00 cov01 cov11) # x # y #| "linear"
+                                                  (x # y # id) (fitting_linear c0 c1 chi_sq cov00 cov01 cov11) #| "linear"
                                                   c0' <- peek c0
                                                   c1' <- peek c1
                                                   cov00' <- peek cov00
@@ -90,7 +90,7 @@
                        alloca $ \cov00 -> 
                            alloca $ \cov01 -> 
                                alloca $ \cov11 -> do
-                                                  (fitting_linear_w c0 c1 chi_sq cov00 cov01 cov11) # x # w # y #| "linear_w"
+                                                  (x # w # y # id) (fitting_linear_w c0 c1 chi_sq cov00 cov01 cov11) #| "linear_w"
                                                   c0' <- peek c0
                                                   c1' <- peek c1
                                                   cov00' <- peek cov00
@@ -139,7 +139,7 @@
                cov <- createMatrix RowMajor p p
                c <- createVector p
                alloca$ \chi_sq -> do
-                   (fitting_multifit chi_sq) # (cmat x) # y # c # cov #| "multifit"
+                   (cmat x # y # c # cov # id) (fitting_multifit chi_sq) #| "multifit"
                    chi_sq' <- peek chi_sq
                    return (c,cov,chi_sq')
 
@@ -159,7 +159,7 @@
                    cov <- createMatrix RowMajor p p
                    c <- createVector p
                    alloca$ \chi_sq -> do
-                       (fitting_multifit_w chi_sq) # (cmat x) # w # y # c # cov #| "multifit"
+                       (cmat x # w # y # c # cov # id) (fitting_multifit_w chi_sq) #| "multifit"
                        chi_sq' <- peek chi_sq
                        return (c,cov,chi_sq')
 
@@ -177,7 +177,7 @@
 multifit_est x c cov = unsafePerformIO $ do
                        alloca $ \y ->
                            alloca $ \e -> do
-                               (fitting_multifit_est y e) # x # c # cov #| "multifit_estimate"
+                               (x # c # cov # id) (fitting_multifit_est y e) #| "multifit_estimate"
                                y' <- peek y
                                e' <- peek e
                                return (y',e')
diff --git a/lib/Numeric/GSL/Histogram.hs b/lib/Numeric/GSL/Histogram.hs
--- a/lib/Numeric/GSL/Histogram.hs
+++ b/lib/Numeric/GSL/Histogram.hs
@@ -3,7 +3,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Histogram
--- Copyright   :  (c) A. V. H. McPhail 2010
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2016
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -77,8 +77,8 @@
 -----------------------------------------------------------------------------
 
 data Hist
-infixl 1 #
-a # b = applyRaw a b
+infixr 1 #
+a # b = apply a b
 {-# INLINE (#) #-}
 
 -----------------------------------------------------------------------------
@@ -136,7 +136,7 @@
                let sz = fromIntegral $ size v - 1
                h <- histogram_new sz
                h' <- newForeignPtr histogram_free h
-               (\d p -> withForeignPtr h' $ \f -> histogram_set_ranges f (fromIntegral d) p) # v #| "fromRanges"
+               (v # id) (\d p -> withForeignPtr h' $ \f -> histogram_set_ranges f (fromIntegral d) p) #| "fromRanges"
                return $ H (fromIntegral sz) h'
 
 -- | create a histogram with n bins and lower and upper limits
@@ -198,7 +198,7 @@
 fromVectors :: Vector Double -> Vector Double -> Histogram
 fromVectors r b = unsafePerformIO $ do
                   h@(H _ h') <- fromRangesIO r
-                  (\rs r' bs b' -> withForeignPtr h' $ \h'' -> histogram_from_vectors h'' rs r' bs b') # r # b #| "fromVectors"
+                  (r # b # id) (\rs r' bs b' -> withForeignPtr h' $ \h'' -> histogram_from_vectors h'' rs r' bs b') #| "fromVectors"
                   return h
                   
 foreign import ccall "gsl-histogram.h from_vectors" histogram_from_vectors :: HistHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
@@ -208,7 +208,7 @@
 toVectors (H b h) = unsafePerformIO $ do
                     rs <- createVector (b+1)
                     bs <- createVector b
-                    (\s1 p1 s2 p2 -> withForeignPtr h $ \f -> histogram_to_vectors f s1 p1 s2 p2) # rs # bs #| "toVectors"
+                    (rs # bs # id) (\s1 p1 s2 p2 -> withForeignPtr h $ \f -> histogram_to_vectors f s1 p1 s2 p2) #| "toVectors"
                     return (rs,bs)
 
 foreign import ccall "gsl-histogram.h to_vectors" histogram_to_vectors :: HistHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
@@ -265,7 +265,7 @@
 -- | add 1.0 to the correct bin for each element of the vector, fails silently if the value is outside the range
 incrementVectorIO :: Histogram -> Vector Double -> IO ()
 incrementVectorIO (H _ h) v = do
-                             (\s p -> withForeignPtr h (\f -> histogram_increment_vector f s p)) # v #| "incrementVector"
+                             (v # id) (\s p -> withForeignPtr h (\f -> histogram_increment_vector f s p)) #| "incrementVector"
                              return ()
 
 -- | adds the weight (second Double) to the bin appropriate for the value (first Double)
@@ -281,7 +281,7 @@
 -- | add the weight (second vector) to the correct bin for each element of the first vector, fails silently if the value is outside the range
 accumulateVectorIO :: Histogram -> Vector Double -> Vector Double -> IO ()
 accumulateVectorIO (H _ h) v w = do
-                                (\s1 p1 s2 p2 -> withForeignPtr h (\f -> histogram_accumulate_vector f s1 p1 s2 p2)) # v # w #| "accumulateVector"
+                                (v # w # id) (\s1 p1 s2 p2 -> withForeignPtr h (\f -> histogram_accumulate_vector f s1 p1 s2 p2)) #| "accumulateVector"
                                 return ()
 
 -- | returns the contents of the i-th bin
diff --git a/lib/Numeric/GSL/Histogram2D.hs b/lib/Numeric/GSL/Histogram2D.hs
--- a/lib/Numeric/GSL/Histogram2D.hs
+++ b/lib/Numeric/GSL/Histogram2D.hs
@@ -3,7 +3,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Histogram2D
--- Copyright   :  (c) A. V. H. McPhail 2010
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2016
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -83,7 +83,7 @@
 
 -----------------------------------------------------------------------------
 
-infixl 1 #
+infixr 1 #
 a # b = applyRaw a b
 {-# INLINE (#) #-}
 
@@ -160,7 +160,7 @@
                    let sy = fromIntegral $ size w - 1
                    h <- histogram2d_new sx sy
                    h' <- newForeignPtr histogram2d_free h
-                   (\xs xp ys yp -> withForeignPtr h' (\f -> histogram2d_set_ranges f xp xs yp ys)) # v # w #| "fromRanges"
+                   (v # w # id) (\xs xp ys yp -> withForeignPtr h' (\f -> histogram2d_set_ranges f xp xs yp ys)) #| "fromRanges"
                    return $ H (fromIntegral sx) (fromIntegral sy) h'
 
 -- | create a histogram with n bins and lower and upper limits
@@ -221,7 +221,7 @@
                     rx <- createVector (bx+1)
                     ry <- createVector (by+1)
                     bs <- createMatrix RowMajor bx by
-                    (\s1 p1 s2 p2 sx sy p -> withForeignPtr h $ \f -> histogram_to_matrix f s1 p1 s2 p2 sx sy p) # rx # ry # bs #| "toMatrix"
+                    (rx # ry # bs # id) (\s1 p1 s2 p2 sx sy p -> withForeignPtr h $ \f -> histogram_to_matrix f s1 p1 s2 p2 sx sy p) #| "toMatrix"
                     return (rx,ry,bs)
 
 foreign import ccall "histogram-aux.h to_matrix" histogram_to_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt
@@ -243,7 +243,7 @@
            -> Histogram2D              -- ^result
 fromMatrix x y w = unsafePerformIO $ do
                    h@(H _ _ h') <- fromRangesIO x y
-                   (\xs x' ys y' rs cs b -> withForeignPtr h' $ \h'' -> histogram_from_matrix h'' xs x' ys y' rs cs b) # x # y # (cmat w) #| "fromMatrix"
+                   (x # y # (cmat w) # id) (\xs x' ys y' rs cs b -> withForeignPtr h' $ \h'' -> histogram_from_matrix h'' xs x' ys y' rs cs b) #| "fromMatrix"
                    return h
 
 foreign import ccall "histogram-aux.h from_matrix" histogram_from_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt
@@ -301,7 +301,7 @@
 -- | add 1.0 to the correct bin for each element of the vector pair, fails silently if the value is outside the range
 incrementVectorIO :: Histogram2D -> Vector Double -> Vector Double -> IO ()
 incrementVectorIO (H _ _ h) x y = do
-                                (\xs xp ys yp -> withForeignPtr h (\f -> histogram2d_increment_matrix f xs xp ys yp)) # x # y #| "incrementVector"
+                                (x # y # id) (\xs xp ys yp -> withForeignPtr h (\f -> histogram2d_increment_matrix f xs xp ys yp)) #| "incrementVector"
                                 return ()
 
 -- | add 1.0  to the correct bin for each element of the list, fails silently if the value is outside the range
@@ -317,7 +317,7 @@
 -- | add the weight (third) to the correct bin for each vector pair element, fails silently if the value is outside the range
 accumulateVectorIO :: Histogram2D -> Vector Double -> Vector Double -> Vector Double -> IO ()
 accumulateVectorIO (H _ _ h) x y w = do
-                                (\xs xp ys yp ws wp -> withForeignPtr h (\f -> histogram2d_accumulate_matrix f xs xp ys yp ws wp)) # x # y # w #| "accumulateVector"
+                                (x # y # w # id) (\xs xp ys yp ws wp -> withForeignPtr h (\f -> histogram2d_accumulate_matrix f xs xp ys yp ws wp)) #| "accumulateVector"
                                 return ()
 
 -- | add the weight (snd) to the correct bin for each (fst) element of the list, fails silently if the value is outside the range
@@ -422,7 +422,7 @@
 count :: Histogram2D -> (Vector Double, Vector Double) -> Vector Double
 count (H _ _ h) (x,y) = unsafePerformIO $ do
                r <- createVector $ size x
-               (\xs' x' ys' y' rs' r' -> withForeignPtr h $ \h' -> histogram2d_count h' xs' x' ys' y' rs' r') # x # y # r #| "histogram2d_count"
+               (x # y # r # id) (\xs' x' ys' y' rs' r' -> withForeignPtr h $ \h' -> histogram2d_count h' xs' x' ys' y' rs' r') #| "histogram2d_count"
                return r
 
 foreign import ccall "histogram-aux.h hist2d_count" histogram2d_count :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
@@ -435,7 +435,7 @@
 countPaired :: Histogram2D -> Vector (Double,Double) -> Vector Double
 countPaired (H _ _ h) x = unsafePerformIO $ do
                r <- createVector $ length x
-               (\xs' x' rs' r' -> withForeignPtr h $ \h' -> histogram2d_count_pair h' xs' x' rs' r') # x # r #| "histogram2d_count_pair"
+               (x # r # id) (\xs' x' rs' r' -> withForeignPtr h $ \h' -> histogram2d_count_pair h' xs' x' rs' r') #| "histogram2d_count_pair"
                return r
 
 foreign import ccall "histogram-aux.h hist2d_count_pair" histogram2d_count_pair :: Hist2DHandle -> CInt -> Ptr (Double,Double) -> CInt -> Ptr Double -> IO CInt
diff --git a/lib/Numeric/GSL/Permutation.hs b/lib/Numeric/GSL/Permutation.hs
--- a/lib/Numeric/GSL/Permutation.hs
+++ b/lib/Numeric/GSL/Permutation.hs
@@ -2,7 +2,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Permutation
--- Copyright   :  (c) A. V. H. McPhail 2010
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2016
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -61,8 +61,8 @@
 
 -----------------------------------------------------------------------------
 
-infixl 1 #
-a # b = applyRaw a b
+infixr 1 #
+a # b = apply a b
 {-# INLINE (#) #-}
 
 -----------------------------------------------------------------------------
@@ -234,14 +234,14 @@
 permute :: Permutation -> Vector Double -> Vector Double
 permute (P n p) v = unsafePerformIO $ do
                     r <- createVector n
-                    (\vs vp rs rp -> withForeignPtr p $ \p' -> permutation_permute p' vs vp rs rp) # v # r #| "permute"
+                    (v # r # id) (\vs vp rs rp -> withForeignPtr p $ \p' -> permutation_permute p' vs vp rs rp) #| "permute"
                     return r
 
 -- | apply the inverse permutation to a vector
 inverse_permute :: Permutation -> Vector Double -> Vector Double
 inverse_permute (P n p) v = unsafePerformIO $ do
                     r <- createVector n
-                    (\vs vp rs rp -> withForeignPtr p $ \p' -> permutation_permute_inverse p' vs vp rs rp) # v # r #| "permute"
+                    (v # r #id) (\vs vp rs rp -> withForeignPtr p $ \p' -> permutation_permute_inverse p' vs vp rs rp) #| "permute"
                     return r
 
 -- | multiply two permutations, P = PA * PB
diff --git a/lib/Numeric/GSL/Sort.hs b/lib/Numeric/GSL/Sort.hs
--- a/lib/Numeric/GSL/Sort.hs
+++ b/lib/Numeric/GSL/Sort.hs
@@ -1,7 +1,8 @@
+{-# LANGUAGE FlexibleContexts #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Sort
--- Copyright   :  (c) A. V. H. McPhail 2010, 2015
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2015, 2016
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -27,15 +28,16 @@
 
 import System.IO.Unsafe(unsafePerformIO)
 
-infixl 1 #
-a # b = applyRaw a b
+infixr 1 #
+a # b = apply a b
 {-# INLINE (#) #-}
 
 -- | sort the elements of a vector into ascending order
 sort :: Vector Double -> Vector Double
 sort v = unsafePerformIO $ do
-         r <- createVector (size v)
-         sort_sort # v # r #| "sort"
-         return r
+  r <- createVector (size v)
+  (v # r # id) sort_sort #| "sort"
+  return r
 
 foreign import ccall "sort-aux.h sort" sort_sort :: CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
+
diff --git a/lib/Numeric/GSL/Statistics.hs b/lib/Numeric/GSL/Statistics.hs
--- a/lib/Numeric/GSL/Statistics.hs
+++ b/lib/Numeric/GSL/Statistics.hs
@@ -1,7 +1,8 @@
+{-# LANGUAGE FlexibleContexts #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Statistics
--- Copyright   :  (c) A. V. H. McPhail 2010
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2016
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -15,7 +16,7 @@
 -----------------------------------------------------------------------------
 
 module Numeric.GSL.Statistics (
-                               mean
+{-                               mean
                               , variance,variance_m,variance_pm
                               , stddev,stddev_m,stddev_pm
                               , tot_sumsq,tot_sumsq_m
@@ -36,7 +37,7 @@
                               , correlation
                               --
                               , median, quantile
-                ) where
+-}                ) where
 
 -----------------------------------------------------------------------------
 
@@ -55,8 +56,8 @@
 
 -----------------------------------------------------------------------------
 
-infixl 1 #
-a # b = applyRaw a b
+infixr 1 #
+a # b = apply a b
 {-# INLINE (#) #-}
 
 -----------------------------------------------------------------------------
@@ -65,15 +66,17 @@
 
 -----------------------------------------------------------------------------
 
+getD1 :: (Ptr Double -> CInt -> Ptr Double -> IO CInt) -> String -> Vector Double -> Double
 getD1 f s v = unsafePerformIO $ do
                 alloca $ \r -> do
-                   (f r) # v #| s
+                   (v # id) (f r) #| s
                    r' <- peek r
                    return r'
 
+getD2 :: (Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt) -> String -> Vector Double -> Vector Double -> Double
 getD2 f s v w = unsafePerformIO $ do
                    alloca $ \r -> do
-                      (f r) # v # w #| s
+                      (v # w # id) (f r) #| s
                       r' <- peek r
                       return r'
 
@@ -312,4 +315,5 @@
 foreign import ccall "statistics-aux.h quantile" statistics_quantile :: Double -> PD -> CInt -> PD -> IO CInt
 
 -----------------------------------------------------------------------------
+