diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # tensors
 
-[![Hackage](https://img.shields.io/badge/hackage-v0.1.0-orange.svg)](https://hackage.haskell.org/package/tensors)
+[![Hackage](https://img.shields.io/badge/hackage-v0.1.1-orange.svg)](https://hackage.haskell.org/package/tensors)
 [![Build Status](https://travis-ci.org/leptonyu/tensors.svg?branch=master)](https://travis-ci.org/leptonyu/tensors)
 
 
diff --git a/src/Data/Tensor.hs b/src/Data/Tensor.hs
--- a/src/Data/Tensor.hs
+++ b/src/Data/Tensor.hs
@@ -77,6 +77,7 @@
   , Index
   -- * Tensor Dimension
   , TensorRank
+  , Shape
   , shape
   , rank
   -- * Tensor Operation
@@ -87,6 +88,9 @@
   -- ** Transpose Tensor
   , Transpose
   , transpose
+  , Swapaxes
+  , CheckSwapaxes
+  , swapaxes
   -- ** Dyadic Tensor
   , dyad'
   , dyad
@@ -94,7 +98,7 @@
   , DotTensor
   , dot
   -- ** Contraction Tensor
-  , ContractionCheck
+  , CheckContraction
   , Contraction
   , TensorDim
   , DropIndex
@@ -108,11 +112,41 @@
   , CheckSlice
   , Slice
   , slice
+  -- ** Tensor Manipulation
   , expand
+  , CheckConcatenate
+  , Concatenate
+  , concatenate
+  , CheckInsert
+  , Insert
+  , insert
+  , append
+  -- ** Tensor Space
+  , linspace
+  , geospace
+  , logspace
+  , CheckGrid
+  , grid
+  , meshgrid2
+  , meshgrid3
+  , meshgrid4
+  , meshgrid5
+  , meshgrid6
+  , meshgrid7
+  , meshgrid8
+  , meshgrid9
   -- * Matrix Operation
+  , SimpleMatrix
   , det
   , lu
   , det'
+  , dotM
+  , trace
+  , diag
+  -- * Statistical Functions
+  , average
+  , var
+  , std
   -- * Helper
   , runTensor
   , i0
@@ -129,5 +163,7 @@
 
 import           Data.Tensor.Index
 import           Data.Tensor.Matrix
+import           Data.Tensor.Space
+import           Data.Tensor.Statistics
 import           Data.Tensor.Tensor
 import           Data.Tensor.Type
diff --git a/src/Data/Tensor/Index.hs b/src/Data/Tensor/Index.hs
--- a/src/Data/Tensor/Index.hs
+++ b/src/Data/Tensor/Index.hs
@@ -1,25 +1,16 @@
-{-# LANGUAGE AllowAmbiguousTypes       #-}
-{-# LANGUAGE DataKinds                 #-}
-{-# LANGUAGE ExistentialQuantification #-}
-{-# LANGUAGE FlexibleInstances         #-}
-{-# LANGUAGE KindSignatures            #-}
-{-# LANGUAGE PolyKinds                 #-}
-{-# LANGUAGE ScopedTypeVariables       #-}
-{-# LANGUAGE TypeFamilies              #-}
-{-# LANGUAGE TypeInType                #-}
-{-# LANGUAGE TypeSynonymInstances      #-}
-{-# LANGUAGE UndecidableInstances      #-}
-
+{-# LANGUAGE UndecidableInstances #-}
 module Data.Tensor.Index where
 
 import           Data.Proxy
+import           Data.Reflection
 import           Data.Singletons
+import qualified Data.Singletons.Prelude.List as N
 import           Data.Tensor.Type
 import           GHC.Exts
 import           GHC.TypeLits
 
 -- | Tensor Index, used to locate each point of tensor
-newtype TensorIndex (shape :: [Nat]) = TensorIndex [Int] deriving (Eq,Show,Ord)
+newtype TensorIndex (shape :: [Nat]) = TensorIndex Index deriving (Eq,Show,Ord)
 
 instance forall s. SingI s => Bounded (TensorIndex s) where
   minBound = toEnum 0
@@ -37,3 +28,21 @@
         else if or (zipWith (\i n-> i <0 || i >= n) v s) then error "index overflow"
           else TensorIndex v
   toList (TensorIndex v) = v
+
+
+-- | Tensor rank.
+type TensorRank (s :: [Nat]) = N.Length s
+
+-- type TensorRankConstraint s i = N.And '[ (N.>=) i 0, (N.<) i (TensorRank s)]
+data TensorRankIndex (shape :: [Nat]) = forall (i :: Nat). KnownNat i => TensorRankIndex (Proxy i)
+
+instance SingI s => Show (TensorRankIndex s) where
+  show = show . fromEnum
+
+instance forall s. (SingI s, KnownNat (TensorRank s - 1)) => Bounded (TensorRankIndex s) where
+  minBound = TensorRankIndex i0
+  maxBound = TensorRankIndex (Proxy :: Proxy (TensorRank s - 1))
+
+instance forall (s::[Nat]). (SingI s) => Enum (TensorRankIndex s) where
+  toEnum i =  reifyNat (toInteger i) TensorRankIndex
+  fromEnum (TensorRankIndex p) = fromInteger $ natVal p
diff --git a/src/Data/Tensor/Matrix.hs b/src/Data/Tensor/Matrix.hs
--- a/src/Data/Tensor/Matrix.hs
+++ b/src/Data/Tensor/Matrix.hs
@@ -1,70 +1,91 @@
 {-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}
-{-# LANGUAGE DataKinds           #-}
-{-# LANGUAGE FlexibleContexts    #-}
-{-# LANGUAGE PolyKinds           #-}
-{-# LANGUAGE RankNTypes          #-}
-{-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE TypeFamilies        #-}
+{-# OPTIONS_GHC -fno-warn-incomplete-uni-patterns #-}
 module Data.Tensor.Matrix where
 
 import           Data.List          (foldl')
 import           Data.Proxy
+import           Data.Ratio
 import           Data.Tensor.Tensor
 import           Data.Tensor.Type
 import           GHC.TypeLits
 
 type SimpleMatrix a n = Matrix a a n
 
-dotM :: (KnownNat a, Num n) => SimpleMatrix a n -> SimpleMatrix a n -> SimpleMatrix a n
+-- | <https://en.wikipedia.org/wiki/Matrix_multiplication Matrix multiplication>
+dotM :: (KnownNat a, Num n, Eq n) => SimpleMatrix a n -> SimpleMatrix a n -> SimpleMatrix a n
 dotM = dot
 
+diag :: SimpleMatrix a n -> Vector a n
+diag (Tensor t) = Tensor $ \[s] [i] -> t [s,s] [i,i]
+
+-- | <https://en.wikipedia.org/wiki/Trace_(linear_algebra) Matrix trace>
 trace :: (KnownNat a, Num n) => SimpleMatrix a n -> n
-trace t = let (Tensor f) = contraction t (i0,i1) in f [] []
+trace t = let (Tensor f) = contraction (i0,i1) t in f [] []
 
 -- | <https://en.wikipedia.org/wiki/LU_decomposition LU decomposition> of n x n matrix
 --
 -- > λ> a = [1,2,3,2,5,7,3,5,3]:: Tensor '[3,3] Int
--- > λ> (l,u,m) = lu a
+-- > λ> (l,u,p) = lu a
 -- > λ> l
--- > [[1,0,0],
--- > [2,1,0],
--- > [3,-1,1]]
+-- > [[1 % 1,0 % 1,0 % 1],
+-- > [2 % 1,1 % 1,0 % 1],
+-- > [3 % 1,(-1) % 1,1 % 1]]
 -- > λ> u
--- > [[1,2,3],
--- > [0,1,1],
--- > [0,0,-5]]
--- > λ> m
--- > 1
-lu :: forall a n . (KnownNat a, Integral n) => SimpleMatrix a n -> (SimpleMatrix a n, SimpleMatrix a n, n)
+-- > [[1 % 1,2 % 1,3 % 1],
+-- > [0 % 1,1 % 1,1 % 1],
+-- > [0 % 1,0 % 1,(-5) % 1]]
+-- > λ> p
+-- > [[1,0,0],
+-- > [0,1,0],
+-- > [0,0,1]]
+lu :: forall a n . (KnownNat a, Integral n)
+   => SimpleMatrix a n
+   -> (SimpleMatrix a (Ratio n), SimpleMatrix a (Ratio n), SimpleMatrix a n)
 lu t =
   let a  = toNat (Proxy :: Proxy a)
-      (l,u,_,m) = foldl' go (identity, t, [a,a], 1) ([0..a-1] :: [Int])
-      p  = minimum (fmap (abs.(`gcd` m)) l)  `min` minimum (fmap (abs.(`gcd` m)) u)
-      g  = (`div` (p * signum m))
-  in (fmap g l,fmap g u, g m)
+      (l,u,p,_) = foldl' go (identity, fmap (% 1) t, identity, [a,a]) ([0..a-1] :: [Int])
+  in (l, u, p)
   where
-    go :: (SimpleMatrix a n, SimpleMatrix a n, [Int], n) -> Int -> (SimpleMatrix a n, SimpleMatrix a n,[Int],n)
-    go (l,u@(Tensor f),s,m) i =
-      let li = Tensor $ \_ -> gi i (f s)
-          lj = Tensor $ \_ -> gj i (f s)
-      in (l `dotM` lj, li `dotM` u, s, m * f s [i,i])
+    {-# INLINE go #-}
+    go (l,u@(Tensor f),p@(Tensor fp),s) i =
+      let ii = f s [i,i]
+      in if ii == 0 then
+          let is = filter (\j -> f s [i,j] /= 0) [i+1..head s-1]
+          in if null is then (l,u,p,s)
+            else let j  = head is
+                     u' = swap i j (f  s)
+                     p' = swap i j (fp s)
+                 in go (l,u',p',s) i
+        else
+          let ij = denominator ii % numerator ii
+              li = Tensor $ \_ i' -> ij * gi i (f s) i'
+              lj = Tensor $ \_ i' -> ij * gj i (f s) i'
+          in (l `dotM` lj, li `dotM` u, p, s)
+    {-# INLINE gi #-}
     gi a fs [x,y]
       | x > a && y == a = - (fs [x,y])
       | x == y = fs [a,a]
       | otherwise = 0
+    {-# INLINE gj #-}
     gj a fs [x,y]
       | x > a && y == a = fs [x,y]
       | x == y = fs [a,a]
       | otherwise = 0
+    {-# INLINE swap #-}
+    swap a b g = Tensor $ \_ [x,y] -> if y == a then g [x,b] else if y == b then g [x,a] else g [x,y]
 
+-- | determiant using `lu` decomposition
 det' :: forall a n . (KnownNat a, Integral n) => SimpleMatrix a n -> n
 det' t =
-  let (l,u,m) = lu t
+  let (l,u,p) = lu t
       s = shape t
-      r = length s
-  in (go s r l * go s r u) `div` (m ^ (r+1))
+      r = head s
+      v = go s r u
+      w = det p
+  in if v == 0 then 0 else go s r l * v * w
   where
-    go s' r' (Tensor f) = let fs = f s' in product $ fmap (\i -> fs [i,i]) ([0..r' - 1] :: [Int])
+    {-# INLINE go #-}
+    go s' r' (Tensor f) = let fs = f s' in numerator $ product $ fmap (\i -> fs [i,i]) ([0..r' - 1] :: [Int])
 
 -- | <https://en.wikipedia.org/wiki/Determinant Determinant> of n x n matrix
 --
@@ -83,10 +104,9 @@
 det = let n = toNat (Proxy :: Proxy a) in go n . runTensor
   where
     {-# INLINE go #-}
-    go :: Int -> ([Int] -> n) -> n
     go 1 f = f [0,0]
-    go n f = sum $ zipWith (g2 f n) ([0.. n-1] :: [Int]) (cycle [1, -1])
+    go n f = sum $ fmap (g2 f n) ([0.. n-1] :: [Int])
     {-# INLINE g2 #-}
-    g2 f n i sign = case f [0,i] of
-      0 -> 0
-      v -> let f' [x,y] = if y >= i then f [x+1,y +1] else f [x+1,y] in  sign * v * go (n-1) f'
+    g2 f n i =
+      let f' [x,y] = if y >= i then f [x+1,y +1] else f [x+1,y]
+      in f [0,i] `mult` (if even i then go (n-1) f' else - (go (n-1) f'))
diff --git a/src/Data/Tensor/Space.hs b/src/Data/Tensor/Space.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Tensor/Space.hs
@@ -0,0 +1,122 @@
+module Data.Tensor.Space where
+
+import           Data.Proxy
+import qualified Data.Singletons.Prelude      as N
+import qualified Data.Singletons.Prelude.List as N
+import           Data.Tensor.Tensor
+import           Data.Tensor.Type
+import qualified Data.Vector                  as V
+import           GHC.TypeLits
+
+-- | Return evenly spaced numbers over a specified interval.
+--
+-- > λ> linspace 1 10 :: Tensor '[10] Float
+-- > [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,..,10.0]
+linspace :: forall s n. (KnownNat s, Fractional n) => n -> n -> Vector s n
+linspace = space (\d s i -> s + d * fromInteger i) (\s e c -> (e - s) / fromInteger (c - 1))
+
+-- | Return numbers spaced evenly on a log scale (a geometric progression).
+--
+-- > λ> geospace 1 10 :: Vector 10 Float
+-- > [1.0,1.2915497,1.6681006,2.1544347,2.7825596,3.593814,4.641589,5.994843,..,10.000001]
+geospace :: forall s n. (KnownNat s, Floating n, Eq n) => n -> n -> Vector s n
+geospace start end
+  | start == 0 = error "divided by zero"
+  | otherwise  = space (\d s i -> s * (d ** fromInteger i)) (\s e c -> (e / s) ** (1 / fromInteger (c - 1))) start end
+
+-- | Return numbers spaced evenly on a log scale.
+--
+-- In linear space, the sequence starts at base ** start (base to the power of start) and ends with base ** stop
+--
+-- > λ> logspace 10 2 3 :: Vector 4 Float
+-- > [100.0,215.44347,464.15887,1000.0]
+logspace :: forall s n. (KnownNat s, Floating n, Eq n) => n -> n -> n -> Vector s n
+logspace base s e = geospace (base ** s) (base ** e)
+
+space :: forall s n. (KnownNat s) => (n -> n -> Integer -> n) -> (n -> n -> Integer -> n) -> n -> n -> Vector s n
+space = go (natVal (Proxy :: Proxy s))
+  where
+    go count g f start end
+      | count <= 1 = pure start
+      | otherwise =
+        let d  = f start end count
+            v  = V.generate (fromInteger count) (g d start . toInteger)
+        in Tensor $ \s i -> v V.! tiTovi s i
+
+
+type CheckGrid i a s = N.And '[ (N.>=) i 0, (N.<) i (N.Length s)]
+-- | Grid
+--
+-- > λ> a = [1..2] :: Vector 2 Int
+-- > λ> a
+-- > [1,2]
+-- > λ> grid i0 a :: Tensor '[2,2,2] Int
+-- > [[[1,1],
+-- > [1,1]],
+-- > [[2,2],
+-- > [2,2]]]
+-- > λ> grid i1 a :: Tensor '[2,2,2] Int
+-- > [[[1,1],
+-- > [2,2]],
+-- > [[1,1],
+-- > [2,2]]]
+-- > λ> grid i2 a :: Tensor '[2,2,2] Int
+-- > [[[1,2],
+-- > [1,2]],
+-- > [[1,2],
+-- > [1,2]]]
+grid
+  :: (CheckDim i s ~ 'True
+    , KnownNat i
+    , KnownNat (TensorDim s i))
+  => Proxy i
+  -> Vector (TensorDim s i) n
+  -> Tensor s n
+grid p v@(Tensor t) =
+  let i = toNat p
+      a = shape v
+  in Tensor $ \_ ind -> t a [ind !! i]
+
+-- | Generate 2 dimension grid coordinates by two vectors.
+meshgrid2 :: (s ~ '[a,b], KnownNat a, KnownNat b) => Vector a n -> Vector b n -> [Tensor s n]
+meshgrid2 a b = [grid i0 a, grid i1 b]
+
+-- | Generate 3 dimension grid coordinates by three vectors.
+meshgrid3 :: (s ~ '[a,b,c], KnownNat a, KnownNat b, KnownNat c) => Vector a n -> Vector b n -> Vector c n -> [Tensor s n]
+meshgrid3 a b c = [grid i0 a, grid i1 b, grid i2 c]
+
+-- | Generate 4 dimension grid coordinates by four vectors.
+meshgrid4
+  :: (s ~ '[a,b,c,d], KnownNat a, KnownNat b, KnownNat c, KnownNat d)
+  => Vector a n -> Vector b n -> Vector c n -> Vector d n -> [Tensor s n]
+meshgrid4 a b c d= [grid i0 a, grid i1 b, grid i2 c, grid i3 d]
+
+-- | Generate 5 dimension grid coordinates by five vectors.
+meshgrid5
+  :: (s ~ '[a,b,c,d,e], KnownNat a, KnownNat b, KnownNat c, KnownNat d, KnownNat e)
+  => Vector a n -> Vector b n -> Vector c n -> Vector d n -> Vector e n -> [Tensor s n]
+meshgrid5 a b c d e = [grid i0 a, grid i1 b, grid i2 c, grid i3 d, grid i4 e]
+
+-- | Generate 6 dimension grid coordinates by six vectors.
+meshgrid6
+  :: (s ~ '[a,b,c,d,e,f], KnownNat a, KnownNat b, KnownNat c, KnownNat d, KnownNat e, KnownNat f)
+  => Vector a n -> Vector b n -> Vector c n -> Vector d n -> Vector e n -> Vector f n -> [Tensor s n]
+meshgrid6 a b c d e f = [grid i0 a, grid i1 b, grid i2 c, grid i3 d, grid i4 e, grid i5 f]
+
+-- | Generate 7 dimension grid coordinates by seven vectors.
+meshgrid7
+  :: (s ~ '[a,b,c,d,e,f,g], KnownNat a, KnownNat b, KnownNat c, KnownNat d, KnownNat e, KnownNat f, KnownNat g)
+  => Vector a n -> Vector b n -> Vector c n -> Vector d n -> Vector e n -> Vector f n -> Vector g n -> [Tensor s n]
+meshgrid7 a b c d e f g = [grid i0 a, grid i1 b, grid i2 c, grid i3 d, grid i4 e, grid i5 f, grid i6 g]
+
+-- | Generate 8 dimension grid coordinates by eight vectors.
+meshgrid8
+  :: (s ~ '[a,b,c,d,e,f,g,h], KnownNat a, KnownNat b, KnownNat c, KnownNat d, KnownNat e, KnownNat f, KnownNat g, KnownNat h)
+  => Vector a n -> Vector b n -> Vector c n -> Vector d n -> Vector e n -> Vector f n -> Vector g n -> Vector h n -> [Tensor s n]
+meshgrid8 a b c d e f g h = [grid i0 a, grid i1 b, grid i2 c, grid i3 d, grid i4 e, grid i5 f, grid i6 g, grid i7 h]
+
+-- | Generate 9 dimension grid coordinates by nine vectors.
+meshgrid9
+  :: (s ~ '[a,b,c,d,e,f,g,h,i], KnownNat a, KnownNat b, KnownNat c, KnownNat d, KnownNat e, KnownNat f, KnownNat g, KnownNat h, KnownNat i)
+  => Vector a n -> Vector b n -> Vector c n -> Vector d n -> Vector e n -> Vector f n -> Vector g n -> Vector h n -> Vector i n -> [Tensor s n]
+meshgrid9 a b c d e f g h i = [grid i0 a, grid i1 b, grid i2 c, grid i3 d, grid i4 e, grid i5 f, grid i6 g, grid i7 h, grid i8 i]
diff --git a/src/Data/Tensor/Statistics.hs b/src/Data/Tensor/Statistics.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Tensor/Statistics.hs
@@ -0,0 +1,31 @@
+module Data.Tensor.Statistics where
+
+import           Data.Singletons
+import           Data.Tensor.Tensor
+
+-- | Average of tensor
+--
+-- > λ> average (identity :: Tensor '[3,3] Float)
+-- > 0.33333334
+average :: forall s n. (SingI s, Fractional n) => Tensor s n -> n
+average t =
+  let v = sum t
+      s = fromInteger $ toInteger $ product $ shape t
+  in v / s
+
+-- | Variance of tensor
+--
+-- > λ> var ([1,2,3,4] :: Vector 4 Double )
+-- > 1.25
+var :: forall s n. (SingI s, Fractional n) => Tensor s n -> n
+var t =
+  let m = average t
+      r = fmap (\v -> let x = v - m in x * x) t
+  in average r
+
+-- | Standard Deviation of tensor
+--
+-- > λ> std ([1,2,3,4] :: Vector 4 Double )
+-- > 1.118033988749895
+std :: forall s n. (SingI s, Floating n) => Tensor s n -> n
+std = sqrt . var
diff --git a/src/Data/Tensor/Tensor.hs b/src/Data/Tensor/Tensor.hs
--- a/src/Data/Tensor/Tensor.hs
+++ b/src/Data/Tensor/Tensor.hs
@@ -1,21 +1,9 @@
-{-# OPTIONS_GHC -fno-warn-redundant-constraints   #-}
 {-# OPTIONS_GHC -fno-warn-incomplete-uni-patterns #-}
-{-# LANGUAGE DataKinds             #-}
-{-# LANGUAGE FlexibleContexts      #-}
-{-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE KindSignatures        #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE PolyKinds             #-}
-{-# LANGUAGE ScopedTypeVariables   #-}
-{-# LANGUAGE Strict                #-}
-{-# LANGUAGE TypeFamilies          #-}
-{-# LANGUAGE TypeOperators         #-}
-{-# LANGUAGE TypeSynonymInstances  #-}
-{-# LANGUAGE UndecidableInstances  #-}
+{-# LANGUAGE UndecidableInstances #-}
 
 module Data.Tensor.Tensor where
 
-import           Data.List                    (group, intercalate)
+import           Data.List                    (intercalate)
 import           Data.Proxy
 import           Data.Singletons
 import qualified Data.Singletons.Prelude      as N
@@ -34,7 +22,7 @@
 -- `s` means shape of tensor.
 --
 -- > identity :: Tensor '[3,3] Int
-newtype Tensor (s :: [Nat]) n = Tensor { getValue :: Index -> Index -> n }
+newtype Tensor (s :: [Nat]) n = Tensor { getValue :: Shape -> Index -> n }
 
 -- | <https://en.wikipedia.org/wiki/Scalarr_(mathematics) Scalar> is rank 0 of tensor
 type Scalar n  = Tensor '[] n
@@ -51,7 +39,8 @@
 -- > SimpleTensor r 0 Int == Scalar Int
 type SimpleTensor (r :: Nat) (dim :: Nat) n = N.If ((N.==) dim 0) (Scalar n) (Tensor (N.Replicate r dim) n)
 
-type TensorRank (s :: [Nat]) = N.Length s
+instance (SingI s, Eq n) => Eq (Tensor s n) where
+  f == g = all (\i -> f ! i == g ! i ) ([minBound..maxBound] :: [TensorIndex s])
 
 instance SingI s => Functor (Tensor s) where
   fmap f (Tensor t) = Tensor (\s i -> f (t s i))
@@ -74,8 +63,7 @@
       {-# INLINE g2 #-}
       g2 n z sep xs = let x = g3 n z xs in "[" ++ intercalate sep x ++ "]"
       {-# INLINE g3 #-}
-      g3 n z xs
-        | z > 3 = take 3 xs ++ [ "..", last xs]
+      g3 n _ xs
         | n > 9 = take 8 xs ++ [ "..", last xs]
         | otherwise = xs
 
@@ -119,7 +107,7 @@
 generateTensor fn p =
   let s  = natsVal p
       ps = product s
-  in if ps == 0 then pure (fn [0]) else Tensor $ \_ -> fn
+  in if ps == 0 then pure (fn [0]) else Tensor $ const fn
 
 {-# INLINE transformTensor #-}
 transformTensor
@@ -194,9 +182,59 @@
     {-# INLINE go #-}
     go (i',_) _ = reverse i'
 
+type CheckSwapaxes i j s = N.And '[ (N.>=) i 0, (N.<) i j, (N.<) j (N.Length s)]
+type Swapaxes i j s = N.Concat '[N.Take i s, '[(N.!!) s j], N.Tail (N.Drop i (N.Take j s)) , '[(N.!!) s i], N.Tail (N.Drop j s)]
+
+-- | Swapaxes any rank
+--
+-- > λ> a = [1..24] :: Tensor '[2,3,4] Int
+-- > λ> a
+-- > [[[1,2,3,4],
+-- > [5,6,7,8],
+-- > [9,10,11,12]],
+-- > [[13,14,15,16],
+-- > [17,18,19,20],
+-- > [21,22,23,24]]]
+-- > λ> swapaxes i0 i1 a
+-- > [[[1,2,3,4],
+-- > [13,14,15,16]],
+-- > [[5,6,7,8],
+-- > [17,18,19,20]],
+-- > [[9,10,11,12],
+-- > [21,22,23,24]]]
+-- > λ> :t swapaxes i0 i1 a
+-- > swapaxes i0 i1 a :: Tensor '[3, 2, 4] Int
+-- > λ> :t swapaxes i1 i2 a
+-- > swapaxes i1 i2 a :: Tensor '[2, 4, 3] Int
+--
+-- In rank 2 tensor, `swapaxes` is just `transpose`
+--
+-- > transpose == swapaxes i0 i1
+swapaxes
+  :: (Swapaxes i j s ~ s'
+    , CheckSwapaxes i j s ~ 'True
+    , SingI s
+    , KnownNat i
+    , KnownNat j)
+  => Proxy i
+  -> Proxy j
+  -> Tensor s n
+  -> Tensor s' n
+swapaxes px pj =
+  let i = toNat px
+      j = toNat pj
+      go (s,_) _ = take i s ++ [s !! j] ++ tail (drop i (take j s)) ++ [s!!i] ++ tail (drop j s)
+  in transformTensor go
+
 -- | Unit tensor of shape s, if all the indices are equal then return 1, otherwise return 0.
 identity :: forall s n . (SingI s, Num n) => Tensor s n
-identity = generateTensor ((\i -> if i == 1 then 1 else 0) . length . group) Proxy
+identity = generateTensor go Proxy
+  where
+    go []  = 0
+    go [_] = 1
+    go (a:b:cs)
+      | a /= b = 0
+      | otherwise = go (b:cs)
 
 dyad'
   :: ( r ~ (N.++) s t
@@ -233,9 +271,10 @@
      , SingI s
      , SingI t
      , SingI r
-     , Num n)
+     , Num n
+     , Eq n)
   => Tensor s n -> Tensor t n -> Tensor r n
-dyad = dyad' (*)
+dyad = dyad' mult
 
 
 type DotTensor s1 s2 = (N.++) (N.Init s1) (N.Tail s2)
@@ -258,7 +297,8 @@
      , SingI (DotTensor s s')
      , SingI s
      , SingI s'
-     , Num n)
+     , Num n
+     , Eq n)
   => Tensor s n
   -> Tensor s' n
   -> Tensor (DotTensor s s') n
@@ -268,13 +308,12 @@
       b  = length s1 - 1
   in generateTensor (\i ->
         let (ti1,ti2) = splitAt b i
-        in sum $ fmap (\(x,y) -> (t1 ! TensorIndex x) * (t2 ! TensorIndex y)) [(ti1++[x],x:ti2)| x <- [0..n-1]]) Proxy
+        in sum $ fmap (\(x,y) -> (t1 ! TensorIndex x) `mult` (t2 ! TensorIndex y)) [(ti1++[x],x:ti2)| x <- [0..n-1]]) Proxy
 
 
-type ContractionCheck s x y = N.And '[(N.<) x y, (N.>=) x 0, (N.<) y (TensorRank s)]
+type CheckContraction s x y = N.And '[(N.<) x y, (N.>=) x 0, (N.<) y (TensorRank s)]
 type Contraction s x y = DropIndex (DropIndex s y) x
-type family TensorDim (s :: [Nat]) (i :: Nat) :: Nat where
-  TensorDim s i = (N.!!) s i
+type TensorDim s i = (N.!!) s i
 type DropIndex (s :: [Nat]) (i :: Nat) = (N.++) (N.Fst (N.SplitAt i s)) (N.Tail (N.Snd (N.SplitAt i s)))
 
 -- | Contraction Tensor
@@ -285,13 +324,13 @@
 -- > [5,6,7,8],
 -- > [9,10,11,12],
 -- > [13,14,15,16]]
--- > λ> contraction a (i0,i1)
+-- > λ> contraction (i0,i1) a
 -- > 34
 --
 -- In rank 2 tensor, contraction of tensor is just the <https://en.wikipedia.org/wiki/Trace_(linear_algebra) trace>.
 contraction
   :: forall x y s s' n.
-     ( ContractionCheck s x y ~ 'True
+     ( CheckContraction s x y ~ 'True
      , s' ~ Contraction s x y
      , TensorDim s x ~ TensorDim s y
      , KnownNat x
@@ -300,10 +339,10 @@
      , SingI s'
      , KnownNat  (TensorDim s x)
      , Num n)
-  => Tensor s  n
-  -> (Proxy x, Proxy y)
+  => (Proxy x, Proxy y)
+  -> Tensor s  n
   -> Tensor s' n
-contraction t@(Tensor f) (px, py) =
+contraction (px, py) t@(Tensor f) =
   let x  = toNat px
       y  = toNat py
       n  = toNat (Proxy :: Proxy (TensorDim s x))
@@ -317,17 +356,15 @@
       in sum $ fmap fs [r1 ++ (j:r3) ++ (j:r4) | j <- [0..n-1]]
 
 type CheckDim dim s = N.And '[(N.>=) dim 0, (N.<) dim (N.Length s)]
-
 type CheckSelect dim i s = N.And '[ CheckDim dim s , (N.>=) i 0, (N.<) i ((N.!!) s dim) ]
-
 type Select i s = (N.++) (N.Take i s) (N.Tail (N.Drop i s))
 
 -- | Select `i` indexing of tensor
 --
 -- > λ> a = identity :: Tensor '[4,4] Int
--- > λ> select a (i0,i0)
+-- > λ> select (i0,i0) a
 -- > [1,0,0,0]
--- > λ> select a (i0,i1)
+-- > λ> select (i0,i1) a
 -- > [0,1,0,0]
 select
   :: ( CheckSelect dim i s ~ 'True
@@ -335,10 +372,10 @@
      , SingI s
      , KnownNat dim
      , KnownNat i)
-  => Tensor s n
-  -> (Proxy dim, Proxy i)
+  => (Proxy dim, Proxy i)
+  -> Tensor s n
   -> Tensor s' n
-select t (pd, pid) =
+select (pd, pid) t=
   let dim = toNat pd
       ind = toNat pid
   in transformTensor (go dim ind) t
@@ -347,7 +384,7 @@
     go d i (i',_) _ = let (a,b) = splitAt d i' in a ++ (i:b)
 
 type CheckSlice dim from to s = N.And '[ CheckDim dim s, CheckSelect dim from s, (N.<) from to , (N.<=) to ((N.!!) s dim)]
-type Slice dim from to s = N.Concat '[N.Take dim s, '[(N.-) to from] , N.Tail (N.Drop dim s)]
+type Slice dim from to s = N.Concat '[N.Take dim s, '[to - from] , N.Tail (N.Drop dim s)]
 
 -- | Slice tensor
 --
@@ -357,10 +394,10 @@
 -- > [0,1,0,0],
 -- > [0,0,1,0],
 -- > [0,0,0,1]]
--- > λ> slice a (i0,(i1,i3))
+-- > λ> slice (i0,(i1,i3)) a
 -- > [[0,1,0,0],
 -- > [0,0,1,0]]
--- > λ> slice a (i1,(i1,i3))
+-- > λ> slice (i1,(i1,i3)) a
 -- > [[0,0],
 -- > [1,0],
 -- > [0,1],
@@ -370,18 +407,16 @@
      , s' ~ Slice dim from to s
      , KnownNat dim
      , KnownNat from
-     , KnownNat ((N.-) to from)
+     , KnownNat (to - from)
      , SingI s)
-  => Tensor s n
-  -> (Proxy dim, (Proxy from, Proxy to))
+  => (Proxy dim, (Proxy from, Proxy to))
+  -> Tensor s n
   -> Tensor s' n
-slice t (pd, (pa,_)) =
+slice (pd, (pa,_)) t =
   let d = toNat pd
       a = toNat pa
   in transformTensor (\(i',_) _ -> let (x,y:ys) = splitAt d i' in x ++ (y+a:ys)) t
 
-type CheckExpand s s' = N.And '[(N.==) (TensorRank s) (TensorRank s')]
-
 -- | Expand tensor
 --
 -- > λ> a = identity :: Tensor '[2,2] Int
@@ -403,6 +438,114 @@
     {-# INLINE go #-}
     go (i',_) = zipWith mod i'
 
+type CheckConcatenate i a b = N.And '[ (N.==) (N.Length a) (N.Length b), (N.>=) i 0, (N.<) i (N.Length a), (N.==) (Select i a) (Select i b) ]
+type Concatenate i a b = N.Concat '[N.Take i a, '[(N.+) (TensorDim a i) (TensorDim b i)], N.Tail (N.Drop i a)]
+
+-- | Join a sequence of arrays along an existing axis.
+--
+-- > λ> a = [1..4] :: Tensor '[2,2] Int
+-- > λ> a
+-- > [[1,2],
+-- > [3,4]]
+-- > λ> b = [1,1,1,1] :: Tensor '[2,2] Int
+-- > λ> b
+-- > [[1,1],
+-- > [1,1]]
+-- > λ> concentrate i0 a b
+-- > [[1,2],
+-- > [3,4],
+-- > [1,1],
+-- > [1,1]]
+-- > λ> concentrate i1 a b
+-- > [[1,2,1,1],
+-- > [3,4,1,1]]
+concatenate
+  :: (CheckConcatenate i a b ~ 'True
+    , Concatenate i a b ~ c
+    , SingI a
+    , SingI b
+    , KnownNat i)
+  => Proxy i
+  -> Tensor a n
+  -> Tensor b n
+  -> Tensor c n
+concatenate p ta@(Tensor a) tb@(Tensor b) =
+  let i  = toNat p
+      sa = shape ta
+      sb = shape tb
+      n  = sa !! i
+  in Tensor $ \_ ind -> let (ai,x:bi) = splitAt i ind in if x >= n then b sb (ai ++ (x-n):bi) else a sa ind
+
+type CheckInsert dim i a b = N.And '[ CheckDim dim b, (N.==) a (Select dim b), (N.>=) i 0, (N.<=) i (TensorDim b dim)]
+type Insert dim a b = N.Concat '[N.Take dim b, '[ TensorDim b dim + 1 ], N.Tail (N.Drop dim b)]
+
+-- | Insert tensor to higher level tensor
+--
+-- > λ> a = [1,2] :: Vector 2 Float
+-- > λ> b = a `dyad` a
+-- > λ> b
+-- > [[1.0,2.0],
+-- > [2.0,4.0]]
+-- > λ> :t b
+-- > b :: Tensor '[2, 2] Float
+-- > λ> c = [1..4] :: Tensor '[1,2,2] Float
+-- > λ> c
+-- > [[[1.0,2.0],
+-- > [3.0,4.0]]]
+-- > λ> d = insert i0 i0 b c
+-- > λ> :t d
+-- > d :: Tensor '[2, 2, 2] Float
+-- > λ> d
+-- > [[[1.0,2.0],
+-- > [2.0,4.0]],
+-- > [[1.0,2.0],
+-- > [3.0,4.0]]]
+-- > λ> insert i0 i1 b c
+-- > [[[1.0,2.0],
+-- > [3.0,4.0]],
+-- > [[1.0,2.0],
+-- > [2.0,4.0]]]
+insert
+  :: (CheckInsert dim i a b ~ 'True
+    , KnownNat i
+    , KnownNat dim
+    , SingI a
+    , SingI b)
+  => Proxy dim
+  -> Proxy i
+  -> Tensor a n
+  -> Tensor b n
+  -> Tensor (Insert dim a b) n
+insert pd px a@(Tensor ta) b@(Tensor tb) =
+  let d = toNat pd
+      i = toNat px
+      sa = shape a
+      sb = shape b
+  in Tensor $ \_ ci -> let (xs,n:ys) = splitAt d ci in if n == i then ta sa (xs++ys) else if n < i then tb sb ci else tb sb (xs ++ ((n-1):ys))
+
+-- | Append tensor to the end of some dimension of other tensor
+--
+-- > λ> a = [1,2] :: Vector 2 Float
+-- > λ> a
+-- > [1.0,2.0]
+-- > λ> b = 3 :: Tensor '[] Float
+-- > λ> b
+-- > 3.0
+-- > λ> append i0 b a
+-- > [1.0,2.0,3.0]
+append
+  :: forall dim a b n.
+    (CheckInsert dim (TensorDim b dim) a b ~ 'True
+    , KnownNat (TensorDim b dim)
+    , KnownNat dim
+    , SingI a
+    , SingI b)
+  => Proxy dim
+  -> Tensor a n
+  -> Tensor b n
+  -> Tensor (Insert dim a b) n
+append pd = insert pd (Proxy :: Proxy (TensorDim b dim))
+
 -- | Convert tensor to untyped function, for internal usage.
-runTensor :: SingI s => Tensor s n -> [Int] -> n
+runTensor :: SingI s => Tensor s n -> Index -> n
 runTensor t@(Tensor f) = f (shape t)
diff --git a/src/Data/Tensor/Type.hs b/src/Data/Tensor/Type.hs
--- a/src/Data/Tensor/Type.hs
+++ b/src/Data/Tensor/Type.hs
@@ -1,14 +1,3 @@
-{-# LANGUAGE DataKinds             #-}
-{-# LANGUAGE FlexibleContexts      #-}
-{-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE KindSignatures        #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE PolyKinds             #-}
-{-# LANGUAGE RankNTypes            #-}
-{-# LANGUAGE ScopedTypeVariables   #-}
-{-# LANGUAGE TypeFamilies          #-}
-{-# LANGUAGE TypeOperators         #-}
-
 module Data.Tensor.Type where
 
 import           Data.List                    (foldl')
@@ -18,12 +7,13 @@
 import           GHC.TypeLits
 import           Unsafe.Coerce
 
+type Shape = [Int]
 type Index = [Int]
 
 toNat :: KnownNat s => Proxy s -> Int
 toNat = unsafeCoerce . natVal
 
-natsVal :: forall (s::[Nat]). SingI s => Proxy s -> Index
+natsVal :: forall proxy (s::[Nat]). SingI s => proxy s -> Index
 natsVal _ = case (sing :: Sing s) of
   SNil         -> []
   (SCons x xs) -> unsafeCoerce <$> (fromSing x: fromSing xs)
@@ -40,6 +30,11 @@
     {-# INLINE go #-}
     go i (n:ns) (ind:inds) = go (i * n + ind) ns inds
     go i _ _               = i
+
+mult :: (Eq a, Num a) => a -> a -> a
+mult a b = case a of
+  0 -> 0
+  c -> c * b
 
 -----------------------
 -- Tensor Type Index
diff --git a/tensors.cabal b/tensors.cabal
--- a/tensors.cabal
+++ b/tensors.cabal
@@ -1,11 +1,13 @@
--- This file has been generated from package.yaml by hpack version 0.28.2.
+cabal-version: 1.12
+
+-- This file has been generated from package.yaml by hpack version 0.31.1.
 --
 -- see: https://github.com/sol/hpack
 --
--- hash: e3f073098b5158a03b88e4fc4de1a4df407ac0ba2001b330f35933f570819334
+-- hash: dbd547840dd8d9ad019059dbf09c36ab86d33914720ce19ab71a4fd3926ef044
 
 name:           tensors
-version:        0.1.0
+version:        0.1.1
 synopsis:       Tensor in Haskell
 description:    Tensor use type level programming in haskell.
 category:       Library
@@ -16,7 +18,6 @@
 license:        BSD3
 license-file:   LICENSE
 build-type:     Simple
-cabal-version:  >= 1.10
 extra-source-files:
     README.md
 
@@ -28,11 +29,15 @@
       Data.Tensor.Index
       Data.Tensor.Tensor
       Data.Tensor.Matrix
+      Data.Tensor.Space
+      Data.Tensor.Statistics
   hs-source-dirs:
       src
-  ghc-options: -Wall -Wcompat -Wincomplete-record-updates -Wincomplete-uni-patterns -Wredundant-constraints -fno-warn-orphans -fno-warn-missing-signatures
+  default-extensions: AllowAmbiguousTypes DataKinds ExistentialQuantification FlexibleContexts FlexibleInstances FunctionalDependencies KindSignatures MultiParamTypeClasses PolyKinds RankNTypes ScopedTypeVariables TypeFamilies TypeOperators TypeSynonymInstances
+  ghc-options: -Wall -Wcompat -Wincomplete-record-updates -Wincomplete-uni-patterns -fno-warn-orphans -fno-warn-missing-signatures
   build-depends:
       base >=4.7 && <5
+    , reflection
     , singletons
     , vector
   default-language: Haskell2010
@@ -44,17 +49,21 @@
       Data.Tensor
       Data.Tensor.Index
       Data.Tensor.Matrix
+      Data.Tensor.Space
+      Data.Tensor.Statistics
       Data.Tensor.Tensor
       Data.Tensor.Type
       Paths_tensors
   hs-source-dirs:
       test
       src
-  ghc-options: -Wall -Wcompat -Wincomplete-record-updates -Wincomplete-uni-patterns -Wredundant-constraints -fno-warn-orphans -fno-warn-missing-signatures
+  default-extensions: AllowAmbiguousTypes DataKinds ExistentialQuantification FlexibleContexts FlexibleInstances FunctionalDependencies KindSignatures MultiParamTypeClasses PolyKinds RankNTypes ScopedTypeVariables TypeFamilies TypeOperators TypeSynonymInstances
+  ghc-options: -Wall -Wcompat -Wincomplete-record-updates -Wincomplete-uni-patterns -fno-warn-orphans -fno-warn-missing-signatures
   build-depends:
       QuickCheck
     , base >=4.7 && <5
     , hspec ==2.*
+    , reflection
     , singletons
     , vector
   default-language: Haskell2010
diff --git a/test/Spec.hs b/test/Spec.hs
--- a/test/Spec.hs
+++ b/test/Spec.hs
@@ -1,24 +1,64 @@
+{-# LANGUAGE DataKinds           #-}
 {-# LANGUAGE FlexibleContexts    #-}
 {-# LANGUAGE OverloadedStrings   #-}
+{-# LANGUAGE PolyKinds           #-}
+{-# LANGUAGE RankNTypes          #-}
 {-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeFamilies        #-}
 
 module Main where
 
+import           Data.Proxy
+import           Data.Reflection
+import           Data.Singletons
+import           Data.Tensor
 import           Data.Tensor.Type
+import           GHC.Exts         (fromList)
+import           GHC.TypeLits
 import           Test.Hspec
 import           Test.QuickCheck
+import           Unsafe.Coerce
 
 main = hspec spec
 
-
 spec :: Spec
 spec = do
   describe "Data.Tensor" specTensor
 
+newtype MagicNats r = MagicNats (forall (n :: [Nat]). SingI n => Proxy n -> r)
+reifyNats :: forall r. [Int] -> (forall (n :: [Nat]). SingI n => Proxy n -> r) -> r
+reifyNats n k = unsafeCoerce (MagicNats k :: MagicNats r) n Proxy
 
+instance SingI a => Reifies (a::[Nat]) [Int] where
+  reflect = natsVal
+
+normalize :: [Int] -> [Int]
+normalize = take 5 . fmap (`mod` 10)
+
 specTensor = do
   context "viToti" $ do
     it "quickCheck" $ property $
-      \s0 -> let s = take 5 $ fmap (`mod` 10) s0
+      \s0 -> let s = normalize s0
                  n = [0..product s - 1]
             in fmap (tiTovi s . viToti s) n == n
+    it "index" $ property $
+      \i -> let i' = mod i 1000 + 1 :: Int
+                go :: forall (s :: Nat). KnownNat s => Int -> Proxy s -> Bool
+                go x _ = trace (identity :: Tensor '[s,s] Int) == x
+            in reifyNat (toInteger i') (go i')
+    it "index range" $ do
+      fmap fromEnum ([minBound..maxBound] :: [TensorIndex '[2,2,2,2,2]]) `shouldBe` [0..31]
+    it "det" $ property $
+      \s0 -> let s = (if null s0 then [1..16] else take 16 $ cycle s0) :: [Int]
+                 a = fromList s :: Tensor '[4,4] Int
+              in det a == det' a
+    it "identity" $ property $
+      \i s0 -> let n = mod i 20 + 1 :: Int
+                   m = n * n
+                   s = if null s0 then [1..m] else take m $ cycle s0
+                   go :: forall (n :: Nat). KnownNat n => [Int] -> Proxy n -> Bool
+                   go v _ =
+                     let a = fromList v :: Tensor '[n,n] Int
+                         idt = identity :: Tensor '[n,n] Int
+                     in a `dot` idt == idt `dot` a
+               in reifyNat (toInteger n) (go s)