diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Daniel YU (c) 2018
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Daniel YU nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,7 @@
+# tensors
+
+[![Hackage](https://img.shields.io/badge/hackage-v0.1.0-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)
+
+
+Type level tensors in Haskell.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import           Distribution.Simple
+main = defaultMain
diff --git a/src/Data/Tensor.hs b/src/Data/Tensor.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Tensor.hs
@@ -0,0 +1,133 @@
+-- |
+-- Module:      Data.Tensor
+-- Copyright:   (c) 2018 Daniel YU
+-- License:     BSD3
+-- Maintainer:  Daniel YU <leptonyu@gmail.com>
+-- Stability:   experimental
+-- Portability: portable
+--
+-- Tensor In Haskell
+--
+-- In ghci
+--
+-- > λ> :set -XDataKinds
+-- > λ> :set -XOverloadedLists
+-- > λ> import Data.Tensor
+-- > λ> a = identity :: Tensor '[3,3] Int
+-- > λ> a
+-- > [[1,0,0],
+-- > [0,1,0],
+-- > [0,0,1]]
+-- > λ> b = [1..9] :: Tensor '[3,3] Int
+-- > λ> b
+-- > [[1,2,3],
+-- > [4,5,6],
+-- > [7,8,9]]
+-- > λ> a + b
+-- > [[2,2,3],
+-- > [4,6,6],
+-- > [7,8,10]]
+-- > λ> a - b
+-- > [[0,-2,-3],
+-- > [-4,-4,-6],
+-- > [-7,-8,-8]]
+-- > λ> a * b
+-- > [[1,0,0],
+-- > [0,5,0],
+-- > [0,0,9]]
+-- > λ> a `dot` b
+-- > [[1,2,3],
+-- > [4,5,6],
+-- > [7,8,9]]
+-- > λ> :t a `dyad` b
+-- > a `dyad` b :: Tensor '[3, 3, 3, 3] Int
+-- > λ> contraction a (i0,i1)
+-- > 3
+-- > λ> :t contraction a (i0,i1)
+-- > contraction a (i0,i1) :: Tensor '[] Int
+-- > λ> select a (i0,i0)
+-- > [1,0,0]
+-- > λ> select a (i0,i1)
+-- > [0,1,0]
+-- > λ> select a (i0,i2)
+-- > [0,0,1]
+-- > λ> c = 1 :: Tensor '[3,3] Int
+-- > λ> c
+-- > [[1,1,1],
+-- > [1,1,1],
+-- > [1,1,1]]
+-- > λ> d = [1..4] :: Tensor '[2,2] Int
+-- > λ> d
+-- > [[1,2],
+-- > [3,4]]
+-- > λ> transpose d
+-- > [[1,3],
+-- > [2,4]]
+
+module Data.Tensor(
+  -- * Tensor Definition
+    Tensor
+  , identity
+  , Scalar
+  , Vector
+  , Matrix
+  , SimpleTensor
+  -- ** Tensor Index
+  , TensorIndex
+  , Index
+  -- * Tensor Dimension
+  , TensorRank
+  , shape
+  , rank
+  -- * Tensor Operation
+  -- ** Reshape Tensor
+  , reshape
+  -- ** Clone Tensor
+  , clone
+  -- ** Transpose Tensor
+  , Transpose
+  , transpose
+  -- ** Dyadic Tensor
+  , dyad'
+  , dyad
+  -- ** Tensor Product
+  , DotTensor
+  , dot
+  -- ** Contraction Tensor
+  , ContractionCheck
+  , Contraction
+  , TensorDim
+  , DropIndex
+  , contraction
+  -- ** Tensor Selection
+  , (!)
+  , CheckDim
+  , CheckSelect
+  , Select
+  , select
+  , CheckSlice
+  , Slice
+  , slice
+  , expand
+  -- * Matrix Operation
+  , det
+  , lu
+  , det'
+  -- * Helper
+  , runTensor
+  , i0
+  , i1
+  , i2
+  , i3
+  , i4
+  , i5
+  , i6
+  , i7
+  , i8
+  , i9
+  ) where
+
+import           Data.Tensor.Index
+import           Data.Tensor.Matrix
+import           Data.Tensor.Tensor
+import           Data.Tensor.Type
diff --git a/src/Data/Tensor/Index.hs b/src/Data/Tensor/Index.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Tensor/Index.hs
@@ -0,0 +1,39 @@
+{-# LANGUAGE AllowAmbiguousTypes       #-}
+{-# LANGUAGE DataKinds                 #-}
+{-# LANGUAGE ExistentialQuantification #-}
+{-# LANGUAGE FlexibleInstances         #-}
+{-# LANGUAGE KindSignatures            #-}
+{-# LANGUAGE PolyKinds                 #-}
+{-# LANGUAGE ScopedTypeVariables       #-}
+{-# LANGUAGE TypeFamilies              #-}
+{-# LANGUAGE TypeInType                #-}
+{-# LANGUAGE TypeSynonymInstances      #-}
+{-# LANGUAGE UndecidableInstances      #-}
+
+module Data.Tensor.Index where
+
+import           Data.Proxy
+import           Data.Singletons
+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)
+
+instance forall s. SingI s => Bounded (TensorIndex s) where
+  minBound = toEnum 0
+  maxBound = let s = natsVal (Proxy :: Proxy s) in  toEnum (product s - 1)
+
+instance forall s. SingI s =>  Enum (TensorIndex s) where
+  toEnum i   = let s = natsVal (Proxy :: Proxy s) in TensorIndex $ viToti s i
+  fromEnum (TensorIndex i) = let s = natsVal (Proxy :: Proxy s) in tiTovi s i
+
+instance forall s. SingI s => IsList (TensorIndex s) where
+  type Item (TensorIndex s) = Int
+  fromList v =
+    let s = natsVal (Proxy :: Proxy s)
+    in if length v /= length s then error "length not match"
+        else if or (zipWith (\i n-> i <0 || i >= n) v s) then error "index overflow"
+          else TensorIndex v
+  toList (TensorIndex v) = v
diff --git a/src/Data/Tensor/Matrix.hs b/src/Data/Tensor/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Tensor/Matrix.hs
@@ -0,0 +1,92 @@
+{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}
+{-# LANGUAGE DataKinds           #-}
+{-# LANGUAGE FlexibleContexts    #-}
+{-# LANGUAGE PolyKinds           #-}
+{-# LANGUAGE RankNTypes          #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeFamilies        #-}
+module Data.Tensor.Matrix where
+
+import           Data.List          (foldl')
+import           Data.Proxy
+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
+dotM = dot
+
+trace :: (KnownNat a, Num n) => SimpleMatrix a n -> n
+trace t = let (Tensor f) = contraction t (i0,i1) 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
+-- > [[1,0,0],
+-- > [2,1,0],
+-- > [3,-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)
+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)
+  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])
+    gi a fs [x,y]
+      | x > a && y == a = - (fs [x,y])
+      | x == y = fs [a,a]
+      | otherwise = 0
+    gj a fs [x,y]
+      | x > a && y == a = fs [x,y]
+      | x == y = fs [a,a]
+      | otherwise = 0
+
+det' :: forall a n . (KnownNat a, Integral n) => SimpleMatrix a n -> n
+det' t =
+  let (l,u,m) = lu t
+      s = shape t
+      r = length s
+  in (go s r l * go s r u) `div` (m ^ (r+1))
+  where
+    go s' r' (Tensor f) = let fs = f s' in product $ fmap (\i -> fs [i,i]) ([0..r' - 1] :: [Int])
+
+-- | <https://en.wikipedia.org/wiki/Determinant Determinant> of n x n matrix
+--
+-- > λ> a = [2,0,1,3,0,0,2,2,1,2,1,34,3,2,34,4] :: Tensor '[4,4] Int
+-- > λ> a
+-- > [[2,0,1,3],
+-- > [0,0,2,2],
+-- > [1,2,1,34],
+-- > [3,2,34,4]]
+-- > λ> det a
+-- > 520
+--
+-- This implementation is not so fast, it can calculate 8 x 8 in 1 second with all the num none zero on my computer.
+-- It should be faster if more zero in the matrix.
+det :: forall a n. (KnownNat a, Num n, Eq n) => SimpleMatrix a n -> n
+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])
+    {-# 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'
diff --git a/src/Data/Tensor/Tensor.hs b/src/Data/Tensor/Tensor.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Tensor/Tensor.hs
@@ -0,0 +1,408 @@
+{-# 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  #-}
+
+module Data.Tensor.Tensor where
+
+import           Data.List                    (group, intercalate)
+import           Data.Proxy
+import           Data.Singletons
+import qualified Data.Singletons.Prelude      as N
+import qualified Data.Singletons.Prelude.List as N
+import           Data.Tensor.Index
+import           Data.Tensor.Type
+import qualified Data.Vector                  as V
+import           GHC.Exts                     (IsList (..))
+import           GHC.TypeLits
+
+-----------------------
+-- Tensor
+-----------------------
+
+-- | Definition of <https://en.wikipedia.org/wiki/Tensor Tensor>.
+-- `s` means shape of tensor.
+--
+-- > identity :: Tensor '[3,3] Int
+newtype Tensor (s :: [Nat]) n = Tensor { getValue :: Index -> Index -> n }
+
+-- | <https://en.wikipedia.org/wiki/Scalarr_(mathematics) Scalar> is rank 0 of tensor
+type Scalar n  = Tensor '[] n
+
+-- | <https://en.wikipedia.org/wiki/Vector_(mathematics_and_physics) Vector> is rank 1 of tensor
+type Vector s n = Tensor '[s] n
+
+-- | <https://en.wikipedia.org/wiki/Matrix_(mathematics) Matrix> is rank 2 of tensor
+type Matrix a b n = Tensor '[a,b] n
+
+-- | Simple Tensor is rank `r` tensor, has `n^r` dimension in total.
+--
+-- > SimpleTensor 2 3 Int == Matrix 3 3 Int == Tensor '[3,3] Int
+-- > 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 => Functor (Tensor s) where
+  fmap f (Tensor t) = Tensor (\s i -> f (t s i))
+
+instance SingI s => Applicative (Tensor s) where
+  pure n = Tensor $ \_ _ -> n
+  Tensor f <*> Tensor t = Tensor $ \s i -> f s i (t s i)
+
+instance SingI s => Foldable (Tensor s) where
+  foldMap f t = foldMap (f.(t !)) ([minBound..maxBound] :: [TensorIndex s])
+
+instance (SingI s, Show n) => Show (Tensor s n) where
+  show (Tensor f) = let s = natsVal (Proxy :: Proxy s) in go 0 [] s (f s)
+    where
+      {-# INLINE go #-}
+      go :: Int -> [Int] -> [Int] -> (Index -> n) -> String
+      go _ i []     fs = show $ fs (reverse i)
+      go z i [n]    fs = g2 n z "," $ fmap (\x -> show (fs $ reverse (x:i))) [0..n-1]
+      go z i (n:ns) fs = g2 n z ",\n" $ fmap (\x -> go (z+1) (x:i) ns fs) [0..n-1]
+      {-# 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]
+        | n > 9 = take 8 xs ++ [ "..", last xs]
+        | otherwise = xs
+
+-----------------------
+-- Tensor as Num
+-----------------------
+instance (SingI s, Num n) => Num (Tensor s n) where
+  (+) = zipWithTensor (+)
+  (*) = zipWithTensor (*)
+  abs = fmap abs
+  signum = fmap signum
+  negate = fmap negate
+  fromInteger = pure . fromInteger
+
+instance (SingI s, Fractional n) => Fractional (Tensor s n) where
+  fromRational = pure . fromRational
+  (/) = zipWithTensor (/)
+
+instance (SingI s, Floating n) => Floating (Tensor s n) where
+  pi      = pure pi
+  exp     = fmap exp
+  log     = fmap log
+  sqrt    = fmap sqrt
+  logBase = error "undefined"
+  sin     = fmap sin
+  cos     = fmap cos
+  tan     = fmap tan
+  asin    = fmap asin
+  acos    = fmap acos
+  atan    = fmap atan
+  sinh    = fmap sinh
+  cosh    = fmap cosh
+  tanh    = fmap tanh
+  asinh   = fmap asinh
+  acosh   = fmap acosh
+  atanh   = fmap atanh
+
+
+{-# INLINE generateTensor #-}
+generateTensor :: SingI s => (Index -> n) -> Proxy s -> Tensor s n
+generateTensor fn p =
+  let s  = natsVal p
+      ps = product s
+  in if ps == 0 then pure (fn [0]) else Tensor $ \_ -> fn
+
+{-# INLINE transformTensor #-}
+transformTensor
+  :: forall s s' n. SingI s
+  => (([Int], [Int]) -> [Int] -> [Int])
+  -> Tensor s  n
+  -> Tensor s' n
+transformTensor go (Tensor f) = let s = natsVal (Proxy :: Proxy s) in Tensor $ \s' i' -> f s (go (i',s') s)
+
+-- | Clone tensor to a new `V.Vector` based tensor
+clone :: SingI s => Tensor s n -> Tensor s n
+clone t =
+  let s = shape t
+      v = V.generate (product s) (\i -> t ! toEnum i)
+  in Tensor $ \_ i -> v V.! tiTovi s i
+
+{-# INLINE zipWithTensor #-}
+zipWithTensor :: SingI s => (n -> n -> n) -> Tensor s n -> Tensor s n -> Tensor s n
+zipWithTensor f t1 t2 = generateTensor (\i -> f (t1 ! TensorIndex i) (t2 ! TensorIndex i)) Proxy
+
+instance SingI s => IsList (Tensor s n) where
+  type Item (Tensor s n) = n
+  fromList v =
+    let s = natsVal (Proxy :: Proxy s)
+        l = product s
+    in if l /= length v
+      then error "length not match"
+      else let vv = V.fromList v in Tensor $ \s' i -> vv V.! tiTovi s' i
+  toList  t = let n = rank t - 1 in fmap (\i -> t ! toEnum i) [0..n]
+
+-----------------------
+-- Tensor Shape
+-----------------------
+-- | Shape of Tensor, is a list of integers, uniquely determine the shape of tensor.
+shape :: forall s n. SingI s => Tensor s n -> [Int]
+shape _ = natsVal (Proxy :: Proxy s)
+
+-- | Rank of Tensor
+rank :: SingI s => Tensor s n -> Int
+rank = length . shape
+
+-----------------------
+-- Tensor Operation
+-----------------------
+-- | Get value from tensor by index
+(!) :: SingI s => Tensor s n -> TensorIndex s -> n
+(!) t (TensorIndex i) = getValue t (shape t) i
+
+-- | Reshape a tensor to another tensor, with total dimensions are equal.
+reshape :: (N.Product s ~ N.Product s', SingI s) => Tensor s n -> Tensor s' n
+reshape = transformTensor go
+  where
+    {-# INLINE go #-}
+    go (i',s') s = viToti s $ tiTovi s' i'
+
+type Transpose (a :: [Nat]) = N.Reverse a
+
+-- | <https://en.wikipedia.org/wiki/Transpose Transpose> tensor completely
+--
+-- > λ> a = [1..9] :: Tensor '[3,3] Int
+-- > λ> a
+-- > [[1,2,3],
+-- > [4,5,6],
+-- > [7,8,9]]
+-- > λ> transpose a
+-- > [[1,4,7],
+-- > [2,5,8],
+-- > [3,6,9]]
+transpose :: SingI a => Tensor a n -> Tensor (Transpose a) n
+transpose  = transformTensor go
+  where
+    {-# INLINE go #-}
+    go (i',_) _ = reverse i'
+
+-- | 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
+
+dyad'
+  :: ( r ~ (N.++) s t
+     , SingI s
+     , SingI t
+     , SingI r)
+  => (n -> m -> o)
+  -> Tensor s n
+  -> Tensor t m
+  -> Tensor r o
+dyad' f t1 t2 =
+  let l = rank t1
+  in generateTensor (\i -> let (ti1,ti2) = splitAt l i in f (t1 ! TensorIndex ti1) (t2 ! TensorIndex ti2)) Proxy
+
+-- | <https://en.wikipedia.org/wiki/Dyadics Dyadic Tensor>
+--
+-- > λ> a = [1..4] :: Tensor '[2,2] Int
+-- > λ> a
+-- > [[1,2],
+-- > [3,4]]
+-- > λ> :t a `dyad` a
+-- > a `dyad` a :: Tensor '[2, 2, 2, 2] Int
+-- > λ> a `dyad` a
+-- > [[[[1,2],
+-- > [3,4]],
+-- > [[2,4],
+-- > [6,8]]],
+-- > [[[3,6],
+-- > [9,12]],
+-- > [[4,8],
+-- > [12,16]]]]
+dyad
+  :: ( r ~ (N.++) s t
+     , SingI s
+     , SingI t
+     , SingI r
+     , Num n)
+  => Tensor s n -> Tensor t n -> Tensor r n
+dyad = dyad' (*)
+
+
+type DotTensor s1 s2 = (N.++) (N.Init s1) (N.Tail s2)
+
+-- | Tensor Product
+--
+-- > λ> a = [1..4] :: Tensor '[2,2] Int
+-- > λ> a
+-- > [[1,2],
+-- > [3,4]]
+-- > λ> a `dot` a
+-- > [[7,10],
+-- > [15,22]]
+--
+-- > dot a b == contraction (dyad a b) (rank a - 1, rank a)
+--
+-- For rank 2 tensor, it is just matrix product.
+dot
+  :: ( N.Last s ~ N.Head s'
+     , SingI (DotTensor s s')
+     , SingI s
+     , SingI s'
+     , Num n)
+  => Tensor s n
+  -> Tensor s' n
+  -> Tensor (DotTensor s s') n
+dot t1 t2 =
+  let s1 = shape t1
+      n  = last s1
+      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
+
+
+type ContractionCheck 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 DropIndex (s :: [Nat]) (i :: Nat) = (N.++) (N.Fst (N.SplitAt i s)) (N.Tail (N.Snd (N.SplitAt i s)))
+
+-- | Contraction Tensor
+--
+-- > λ> a = [1..16] :: Tensor '[4,4] Int
+-- > λ> a
+-- > [[1,2,3,4],
+-- > [5,6,7,8],
+-- > [9,10,11,12],
+-- > [13,14,15,16]]
+-- > λ> contraction a (i0,i1)
+-- > 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
+     , s' ~ Contraction s x y
+     , TensorDim s x ~ TensorDim s y
+     , KnownNat x
+     , KnownNat y
+     , SingI s
+     , SingI s'
+     , KnownNat  (TensorDim s x)
+     , Num n)
+  => Tensor s  n
+  -> (Proxy x, Proxy y)
+  -> Tensor s' n
+contraction t@(Tensor f) (px, py) =
+  let x  = toNat px
+      y  = toNat py
+      n  = toNat (Proxy :: Proxy (TensorDim s x))
+      s  = shape t
+  in generateTensor (go x (y-x-1) n (f s) ) Proxy
+  where
+    {-# INLINE go #-}
+    go a b n fs i =
+      let (r1,rt) = splitAt a i
+          (r3,r4) = splitAt b rt
+      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)
+-- > [1,0,0,0]
+-- > λ> select a (i0,i1)
+-- > [0,1,0,0]
+select
+  :: ( CheckSelect dim i s ~ 'True
+     , s' ~ Select dim s
+     , SingI s
+     , KnownNat dim
+     , KnownNat i)
+  => Tensor s n
+  -> (Proxy dim, Proxy i)
+  -> Tensor s' n
+select t (pd, pid) =
+  let dim = toNat pd
+      ind = toNat pid
+  in transformTensor (go dim ind) t
+  where
+    {-# INLINE go #-}
+    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)]
+
+-- | Slice tensor
+--
+-- > λ> a = identity :: Tensor '[4,4] Int
+-- > λ> a
+-- > [[1,0,0,0],
+-- > [0,1,0,0],
+-- > [0,0,1,0],
+-- > [0,0,0,1]]
+-- > λ> slice a (i0,(i1,i3))
+-- > [[0,1,0,0],
+-- > [0,0,1,0]]
+-- > λ> slice a (i1,(i1,i3))
+-- > [[0,0],
+-- > [1,0],
+-- > [0,1],
+-- > [0,0]]
+slice
+  :: ( CheckSlice dim from to s ~ 'True
+     , s' ~ Slice dim from to s
+     , KnownNat dim
+     , KnownNat from
+     , KnownNat ((N.-) to from)
+     , SingI s)
+  => Tensor s n
+  -> (Proxy dim, (Proxy from, Proxy to))
+  -> Tensor s' n
+slice t (pd, (pa,_)) =
+  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
+-- > λ> a
+-- > [[1,0],
+-- > [0,1]]
+-- > λ> expand a :: Tensor '[4,4] Int
+-- > [[1,0,1,0],
+-- > [0,1,0,1],
+-- > [1,0,1,0],
+-- > [0,1,0,1]]
+expand
+  :: (TensorRank s ~ TensorRank s'
+     , SingI s)
+  => Tensor s n
+  -> Tensor s' n
+expand = transformTensor go
+  where
+    {-# INLINE go #-}
+    go (i',_) = zipWith mod i'
+
+-- | Convert tensor to untyped function, for internal usage.
+runTensor :: SingI s => Tensor s n -> [Int] -> n
+runTensor t@(Tensor f) = f (shape t)
diff --git a/src/Data/Tensor/Type.hs b/src/Data/Tensor/Type.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Tensor/Type.hs
@@ -0,0 +1,56 @@
+{-# 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')
+import           Data.Proxy
+import           Data.Singletons
+import           Data.Singletons.Prelude.List
+import           GHC.TypeLits
+import           Unsafe.Coerce
+
+type Index = [Int]
+
+toNat :: KnownNat s => Proxy s -> Int
+toNat = unsafeCoerce . natVal
+
+natsVal :: forall (s::[Nat]). SingI s => Proxy s -> Index
+natsVal _ = case (sing :: Sing s) of
+  SNil         -> []
+  (SCons x xs) -> unsafeCoerce <$> (fromSing x: fromSing xs)
+
+viToti :: Index -> Int -> Index
+viToti s i = snd $ foldl' go (i,[]) (reverse s)
+  where
+    {-# INLINE go #-}
+    go (i',x) n = let (d,r) = divMod i' n in (d,r:x)
+
+tiTovi :: Index -> Index -> Int
+tiTovi = go 0
+  where
+    {-# INLINE go #-}
+    go i (n:ns) (ind:inds) = go (i * n + ind) ns inds
+    go i _ _               = i
+
+-----------------------
+-- Tensor Type Index
+-----------------------
+i0  = Proxy :: Proxy 0
+i1  = Proxy :: Proxy 1
+i2  = Proxy :: Proxy 2
+i3  = Proxy :: Proxy 3
+i4  = Proxy :: Proxy 4
+i5  = Proxy :: Proxy 5
+i6  = Proxy :: Proxy 6
+i7  = Proxy :: Proxy 7
+i8  = Proxy :: Proxy 8
+i9  = Proxy :: Proxy 9
diff --git a/tensors.cabal b/tensors.cabal
new file mode 100644
--- /dev/null
+++ b/tensors.cabal
@@ -0,0 +1,60 @@
+-- This file has been generated from package.yaml by hpack version 0.28.2.
+--
+-- see: https://github.com/sol/hpack
+--
+-- hash: e3f073098b5158a03b88e4fc4de1a4df407ac0ba2001b330f35933f570819334
+
+name:           tensors
+version:        0.1.0
+synopsis:       Tensor in Haskell
+description:    Tensor use type level programming in haskell.
+category:       Library
+homepage:       https://github.com/leptonyu/tensors#readme
+author:         Daniel YU
+maintainer:     Daniel YU <leptonyu@gmail.com>
+copyright:      (c) 2018 Daniel YU
+license:        BSD3
+license-file:   LICENSE
+build-type:     Simple
+cabal-version:  >= 1.10
+extra-source-files:
+    README.md
+
+library
+  exposed-modules:
+      Data.Tensor
+  other-modules:
+      Data.Tensor.Type
+      Data.Tensor.Index
+      Data.Tensor.Tensor
+      Data.Tensor.Matrix
+  hs-source-dirs:
+      src
+  ghc-options: -Wall -Wcompat -Wincomplete-record-updates -Wincomplete-uni-patterns -Wredundant-constraints -fno-warn-orphans -fno-warn-missing-signatures
+  build-depends:
+      base >=4.7 && <5
+    , singletons
+    , vector
+  default-language: Haskell2010
+
+test-suite spec
+  type: exitcode-stdio-1.0
+  main-is: Spec.hs
+  other-modules:
+      Data.Tensor
+      Data.Tensor.Index
+      Data.Tensor.Matrix
+      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
+  build-depends:
+      QuickCheck
+    , base >=4.7 && <5
+    , hspec ==2.*
+    , singletons
+    , vector
+  default-language: Haskell2010
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -0,0 +1,24 @@
+{-# LANGUAGE FlexibleContexts    #-}
+{-# LANGUAGE OverloadedStrings   #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Main where
+
+import           Data.Tensor.Type
+import           Test.Hspec
+import           Test.QuickCheck
+
+main = hspec spec
+
+
+spec :: Spec
+spec = do
+  describe "Data.Tensor" specTensor
+
+
+specTensor = do
+  context "viToti" $ do
+    it "quickCheck" $ property $
+      \s0 -> let s = take 5 $ fmap (`mod` 10) s0
+                 n = [0..product s - 1]
+            in fmap (tiTovi s . viToti s) n == n