diff --git a/egison.cabal b/egison.cabal
--- a/egison.cabal
+++ b/egison.cabal
@@ -1,5 +1,5 @@
 Name:                egison
-Version:             3.6.5
+Version:             3.7.0
 Synopsis:            Programming language with non-linear pattern-matching against non-free data
 Description:
   An interpreter for Egison, a **pattern-matching-oriented**, purely functional programming language.
@@ -59,7 +59,7 @@
 
 Extra-Source-Files:  benchmark/Benchmark.hs
 
-Data-files:          lib/core/*.egi lib/math/*.egi lib/math/common/*.egi lib/math/algebra/*.egi lib/math/analysis/*.egi
+Data-files:          lib/core/*.egi lib/math/*.egi lib/math/common/*.egi lib/math/algebra/*.egi lib/math/analysis/*.egi lib/math/geometry/*.egi
                      sample/*.egi sample/io/*.egi
                      elisp/egison-mode.el
 
diff --git a/hs-src/Language/Egison.hs b/hs-src/Language/Egison.hs
--- a/hs-src/Language/Egison.hs
+++ b/hs-src/Language/Egison.hs
@@ -118,6 +118,7 @@
   , "lib/math/algebra/vector.egi"
   , "lib/math/algebra/matrix.egi"
   , "lib/math/algebra/tensor.egi"
+  , "lib/math/geometry/differential-form.egi"
   , "lib/core/base.egi"
   , "lib/core/collection.egi"
   , "lib/core/assoc.egi"
diff --git a/hs-src/Language/Egison/Core.hs b/hs-src/Language/Egison/Core.hs
--- a/hs-src/Language/Egison/Core.hs
+++ b/hs-src/Language/Egison/Core.hs
@@ -263,7 +263,7 @@
     (UserIndexedData val' is') -> return $ Value $ UserIndexedData val' (is' ++ js)
     _ -> return $ Value $ UserIndexedData val js
 
-evalExpr env (IndexedExpr expr indices) = do
+evalExpr env (IndexedExpr False expr indices) = do
   tensor <- case expr of
               (VarExpr var) -> do
                 let mObjRef = refVar env (show (Var var (map f indices)))
@@ -276,6 +276,7 @@
                       Subscript n -> evalExprDeep env n >>= return . Subscript
                       SupSubscript n -> evalExprDeep env n >>= return . SupSubscript
               ) indices
+  
   ret <- case tensor of
       (Value (ScalarData (Div (Plus [(Term 1 [(Symbol id name [], 1)])]) (Plus [(Term 1 [])])))) -> do
         js2 <- mapM (\i -> case i of
@@ -308,6 +309,54 @@
   f (Subscript _) = Subscript ()
   f (SupSubscript _) = SupSubscript ()
 
+evalExpr env (SubrefsExpr expr jsExpr) = do
+  js <- evalExpr env jsExpr >>= collectionToList >>= return . (map Subscript)
+  tensor <- case expr of
+              (VarExpr var) -> do
+                let mObjRef = refVar env (show (Var var (take (length js) (repeat (Subscript ())))))
+                case mObjRef of
+                  (Just objRef) -> evalRef objRef
+                  Nothing -> evalExpr env expr
+              _ -> evalExpr env expr
+  ret <- case tensor of
+      (Value (ScalarData _)) -> do
+        return $ tensor
+      (Value (TensorData (Tensor ns xs is))) -> do
+        tref (is ++ js) (Tensor ns xs (is ++ js)) >>= toTensor >>= tContract' >>= fromTensor >>= return . Value
+      (Intermediate (ITensor (Tensor ns xs is))) -> do
+        tref (is ++ js) (Tensor ns xs (is ++ js)) >>= toTensor >>= tContract' >>= fromTensor
+      _ -> throwError $ NotImplemented "subrefs"
+  return ret
+ where
+  f :: Index a -> Index ()
+  f (Superscript _) = Superscript ()
+  f (Subscript _) = Subscript ()
+  f (SupSubscript _) = SupSubscript ()
+
+evalExpr env (SuprefsExpr expr jsExpr) = do
+  js <- evalExpr env jsExpr >>= collectionToList >>= return . (map Superscript)
+  tensor <- case expr of
+              (VarExpr var) -> do
+                let mObjRef = refVar env (show (Var var (take (length js) (repeat (Superscript ())))))
+                case mObjRef of
+                  (Just objRef) -> evalRef objRef
+                  Nothing -> evalExpr env expr
+              _ -> evalExpr env expr
+  ret <- case tensor of
+      (Value (ScalarData _)) -> do
+        return $ tensor
+      (Value (TensorData (Tensor ns xs is))) -> do
+        tref (is ++ js) (Tensor ns xs (is ++ js)) >>= toTensor >>= tContract' >>= fromTensor >>= return . Value
+      (Intermediate (ITensor (Tensor ns xs is))) -> do
+        tref (is ++ js) (Tensor ns xs (is ++ js)) >>= toTensor >>= tContract' >>= fromTensor
+      _ -> throwError $ NotImplemented "suprefs"
+  return ret
+ where
+  f :: Index a -> Index ()
+  f (Superscript _) = Superscript ()
+  f (Subscript _) = Subscript ()
+  f (SupSubscript _) = SupSubscript ()
+
 evalExpr env (LambdaExpr names expr) = do
   names' <- mapM (\name -> case name of
                              (TensorArg name') -> return name'
@@ -504,10 +553,39 @@
     _ -> applyFunc env func (Value (makeTuple args))
 
 evalExpr env (ApplyExpr func arg) = do
-  func <- evalExpr env func
+  func <- evalExpr env func >>= appendDFscripts 0
   arg <- evalExpr env arg
+--  arg <- evalExpr env arg >>= fromTupleWHNF
+--  let k = fromIntegral (length arg)
+--  arg <-  mapM (\(_,j) -> appendDFscripts 0 j) (zip [1..k] arg) >>= makeITuple
   case func of
     Value (TensorData t@(Tensor ns fs js)) -> do
+      tMap (\f -> applyFunc env (Value f) arg >>= evalWHNF) t >>= fromTensor >>= return . Value >>= removeDFscripts
+    Intermediate (ITensor t@(Tensor ns fs js)) -> do
+      tMap (\f -> applyFunc env f arg) t >>= fromTensor
+    Value (MemoizedFunc name ref hashRef env names body) -> do
+      indices <- evalWHNF arg
+      indices' <- mapM fromEgison $ fromTupleValue indices
+      hash <- liftIO $ readIORef hashRef
+      case HL.lookup indices' hash of
+        Just objRef -> do
+          evalRef objRef
+        Nothing -> do
+          whnf <- applyFunc env (Value (Func Nothing env names body)) arg
+          retRef <- newEvaluatedObjectRef whnf
+          hash <- liftIO $ readIORef hashRef
+          liftIO $ writeIORef hashRef (HL.insert indices' retRef hash)
+          writeObjectRef ref (Value (MemoizedFunc name ref hashRef env names body))
+          return whnf
+    _ -> applyFunc env func arg >>= removeDFscripts
+
+evalExpr env (WedgeApplyExpr func arg) = do
+  func <- evalExpr env func >>= appendDFscripts 0
+  arg <- evalExpr env arg >>= fromTupleWHNF
+  let k = fromIntegral (length arg)
+  arg <-  mapM (\(i,j) -> appendDFscripts i j) (zip [1..k] arg) >>= makeITuple
+  case func of
+    Value (TensorData t@(Tensor ns fs js)) -> do
       tMap (\f -> applyFunc env (Value f) arg >>= evalWHNF) t >>= fromTensor >>= return . Value
     Intermediate (ITensor t@(Tensor ns fs js)) -> do
       tMap (\f -> applyFunc env f arg) t >>= fromTensor
@@ -525,7 +603,7 @@
           liftIO $ writeIORef hashRef (HL.insert indices' retRef hash)
           writeObjectRef ref (Value (MemoizedFunc name ref hashRef env names body))
           return whnf
-    _ -> applyFunc env func arg
+    _ -> applyFunc env func arg >>= removeDFscripts
 
 evalExpr env (MemoizeExpr memoizeFrame expr) = do
   mapM (\(x, y, z) -> do x' <- evalExprDeep env x
@@ -1348,6 +1426,11 @@
 makeTuple [] = Tuple []
 makeTuple [x] = x
 makeTuple xs = Tuple xs
+
+makeITuple :: [WHNFData] -> EgisonM WHNFData
+makeITuple [] = return $ Intermediate (ITuple [])
+makeITuple [x] = return $ x
+makeITuple xs = mapM newEvaluatedObjectRef xs >>= (return . Intermediate . ITuple)
 
 --
 -- String
diff --git a/hs-src/Language/Egison/Desugar.hs b/hs-src/Language/Egison/Desugar.hs
--- a/hs-src/Language/Egison/Desugar.hs
+++ b/hs-src/Language/Egison/Desugar.hs
@@ -134,12 +134,18 @@
 
 desugar (ArrayRefExpr expr nums) =
   case nums of
-    (TupleExpr nums') -> desugar $ IndexedExpr expr (map Subscript nums')
-    _ -> desugar $ IndexedExpr expr [Subscript nums]
+    (TupleExpr nums') -> desugar $ IndexedExpr False expr (map Subscript nums')
+    _ -> desugar $ IndexedExpr False expr [Subscript nums]
   
-desugar (IndexedExpr expr indices) = 
-  IndexedExpr <$> desugar expr <*> (mapM desugarIndex indices)
+desugar (IndexedExpr b expr indices) = 
+  IndexedExpr b <$> desugar expr <*> (mapM desugarIndex indices)
 
+desugar (SubrefsExpr expr1 expr2) = 
+  SubrefsExpr <$> desugar expr1 <*> desugar expr2
+
+desugar (SuprefsExpr expr1 expr2) = 
+  SuprefsExpr <$> desugar expr1 <*> desugar expr2
+
 desugar (PowerExpr expr1 expr2) = do
   expr1' <- desugar expr1
   expr2' <- desugar expr2
@@ -361,6 +367,11 @@
 desugar (QuoteFunctionExpr expr) = do
   expr' <- desugar expr
   return $ QuoteFunctionExpr expr'
+
+desugar (WedgeExpr (ApplyExpr expr0 expr1)) = do
+  expr0' <- desugar expr0
+  expr1' <- desugar expr1
+  return $ WedgeApplyExpr expr0' expr1'
 
 desugar expr = return expr
 
diff --git a/hs-src/Language/Egison/Parser.hs b/hs-src/Language/Egison/Parser.hs
--- a/hs-src/Language/Egison/Parser.hs
+++ b/hs-src/Language/Egison/Parser.hs
@@ -177,7 +177,7 @@
 
 expr :: Parser EgisonExpr
 expr = P.lexeme lexer (do expr0 <- expr' <|> quoteExpr'
-                          expr1 <- option expr0 $ IndexedExpr expr0 <$> many1 (try (char '_' >> expr' >>= return . Subscript)
+                          expr1 <- option expr0 $ IndexedExpr False expr0 <$> many1 (try (char '_' >> expr' >>= return . Subscript)
                                                                            <|> try (char '~' >> expr' >>= return . Superscript)
                                                                            <|> try (string "~_" >> expr' >>= return . SupSubscript))
                           expr2 <- option expr1 $ UserIndexedExpr expr1 <$> many1 (try $ char '|' >> expr' >>= return . Userscript)
@@ -200,6 +200,7 @@
              <|> collectionExpr
 --             <|> quoteExpr
              <|> quoteFunctionExpr
+             <|> wedgeExpr
              <|> parens (ifExpr
                          <|> lambdaExpr
                          <|> memoizedLambdaExpr
@@ -241,6 +242,8 @@
                          <|> parExpr
                          <|> pseqExpr
                          <|> pmapExpr
+                         <|> subrefsExpr
+                         <|> suprefsExpr
                          )
              <?> "expression")
 
@@ -286,6 +289,9 @@
 quoteExpr :: Parser EgisonExpr
 quoteExpr = char '\'' >> QuoteExpr <$> expr
 
+wedgeExpr :: Parser EgisonExpr
+wedgeExpr = char '!' >> WedgeExpr <$> expr
+
 quoteFunctionExpr :: Parser EgisonExpr
 quoteFunctionExpr = char '`' >> QuoteFunctionExpr <$> expr
 
@@ -572,6 +578,12 @@
 pmapExpr :: Parser EgisonExpr
 pmapExpr = keywordPmap >> PmapExpr <$> expr <*> expr
 
+subrefsExpr :: Parser EgisonExpr
+subrefsExpr = keywordSubrefs >> SubrefsExpr <$> expr <*> expr
+
+suprefsExpr :: Parser EgisonExpr
+suprefsExpr = keywordSuprefs >> SuprefsExpr <$> expr <*> expr
+
 -- Patterns
 
 pattern :: Parser EgisonPattern
@@ -808,6 +820,8 @@
   , "par"
   , "pseq"
   , "pmap"
+  , "subrefs"
+  , "suprefs"
   , "something"
   , "undefined"]
   
@@ -887,6 +901,8 @@
 keywordPar                  = reserved "par"
 keywordPseq                 = reserved "pseq"
 keywordPmap                 = reserved "pmap"
+keywordSubrefs              = reserved "subrefs"
+keywordSuprefs              = reserved "suprefs"
 
 sign :: Num a => Parser (a -> a)
 sign = (char '-' >> return negate)
diff --git a/hs-src/Language/Egison/Primitives.hs b/hs-src/Language/Egison/Primitives.hs
--- a/hs-src/Language/Egison/Primitives.hs
+++ b/hs-src/Language/Egison/Primitives.hs
@@ -151,8 +151,8 @@
              , ("modulo",    integerBinaryOp mod)
              , ("quotient",   integerBinaryOp quot)
              , ("remainder", integerBinaryOp rem)
-             , ("abs", rationalUnaryOp abs)
-             , ("neg", rationalUnaryOp negate)
+             , ("b.abs", rationalUnaryOp abs)
+             , ("b.neg", rationalUnaryOp negate)
                
              , ("eq?",  eq)
              , ("lt?",  lt)
@@ -186,6 +186,7 @@
 
              , ("tensor-size", tensorSize')
              , ("tensor-to-list", tensorToList')
+             , ("df-order", dfOrder')
 
              , ("itof", integerToFloat)
              , ("rtof", rationalToFloat)
@@ -448,13 +449,19 @@
 tensorSize' = oneArg' $ tensorSize''
  where
   tensorSize'' (TensorData (Tensor ns _ _)) = return . Collection . Sq.fromList $ map toEgison ns
-  tensorSize'' _ = return . Collection $ Sq.fromList $ [toEgison (1 :: Integer)]
+  tensorSize'' _ = return . Collection $ Sq.fromList $ []
 
 tensorToList' :: PrimitiveFunc
 tensorToList' = oneArg' $ tensorToList''
  where
   tensorToList'' (TensorData (Tensor _ xs _)) = return . Collection . Sq.fromList $ V.toList xs
   tensorToList'' x = return . Collection $ Sq.fromList $ [x]
+
+dfOrder' :: PrimitiveFunc
+dfOrder' = oneArg' $ dfOrder''
+ where
+  dfOrder'' (TensorData (Tensor ns _ is)) = return (toEgison ((fromIntegral ((length ns) - (length is))) :: Integer))
+  dfOrder'' _ = return (toEgison (0 :: Integer))
 
 --
 -- Transform
diff --git a/hs-src/Language/Egison/Types.hs b/hs-src/Language/Egison/Types.hs
--- a/hs-src/Language/Egison/Types.hs
+++ b/hs-src/Language/Egison/Types.hs
@@ -45,6 +45,8 @@
     , tref
     , enumTensorIndices
     , tTranspose'
+    , appendDFscripts
+    , removeDFscripts
     , tMap
     , tMap2
     , tMapN
@@ -169,7 +171,7 @@
 import Data.HashMap.Strict (HashMap)
 import qualified Data.HashMap.Strict as HashMap
 
-import Data.List (intercalate, sort, sortBy, find, findIndex, splitAt, (\\), elem, delete, deleteBy, any)
+import Data.List (intercalate, sort, sortBy, find, findIndex, splitAt, (\\), elem, delete, deleteBy, any, partition)
 import Data.Text (Text)
 import qualified Data.Text as T
 
@@ -201,7 +203,9 @@
   | FloatExpr Double Double
   | VarExpr String
   | FreshVarExpr
-  | IndexedExpr EgisonExpr [Index EgisonExpr]
+  | IndexedExpr Bool EgisonExpr [Index EgisonExpr]
+  | SubrefsExpr EgisonExpr EgisonExpr
+  | SuprefsExpr EgisonExpr EgisonExpr
   | UserIndexedExpr EgisonExpr [UserIndex EgisonExpr]
   | PowerExpr EgisonExpr EgisonExpr
   | InductiveDataExpr String [EgisonExpr]
@@ -242,6 +246,9 @@
   | QuoteExpr EgisonExpr
   | QuoteFunctionExpr EgisonExpr
   
+  | WedgeExpr EgisonExpr
+  | WedgeApplyExpr EgisonExpr EgisonExpr
+
   | DoExpr [BindingExpr] EgisonExpr
   | IoExpr EgisonExpr
     
@@ -280,6 +287,7 @@
     Subscript a
   | Superscript a
   | SupSubscript a
+  | DFscript Integer Integer -- DifferentialForm
  deriving (Eq)
 
 data UserIndex a = Userscript a
@@ -469,6 +477,9 @@
 symbolScalarData :: String -> String -> EgisonValue
 symbolScalarData id name = ScalarData (Div (Plus [(Term 1 [(Symbol id name [], 1)])]) (Plus [(Term 1 [])]))
 
+getSymId :: EgisonValue -> String
+getSymId (ScalarData (Div (Plus [(Term 1 [(Symbol id name [], 1)])]) (Plus [(Term 1 [])]))) = id
+
 mathExprToEgison :: ScalarData -> EgisonValue
 mathExprToEgison (Div p1 p2) = InductiveData "Div" [(polyExprToEgison p1), (polyExprToEgison p2)]
 
@@ -840,12 +851,42 @@
                   (Just j') -> do js' <- g is js
                                   return $ j':js'
 
+appendDFscripts :: Integer -> WHNFData -> EgisonM WHNFData
+appendDFscripts id (Intermediate (ITensor (Tensor s xs is))) = do
+  let k = fromIntegral ((length s) - (length is))
+  return $ Intermediate (ITensor (Tensor s xs (is ++ (map (DFscript id) [1..k]))))
+appendDFscripts id (Value (TensorData (Tensor s xs is))) = do
+  let k = fromIntegral ((length s) - (length is))
+  return $ Value (TensorData (Tensor s xs (is ++ (map (DFscript id) [1..k]))))
+appendDFscripts _ whnf = return whnf
+
+removeDFscripts :: WHNFData -> EgisonM WHNFData
+removeDFscripts (Intermediate (ITensor (Tensor s xs is))) = do
+  let (ds, js) = partition isDF is
+  (Tensor s ys _) <- tTranspose (js ++ ds) (Tensor s xs is)
+  return (Intermediate (ITensor (Tensor s ys js)))
+ where
+  isDF (DFscript _ _) = True
+  isDF _ = False
+removeDFscripts (Value (TensorData (Tensor s xs is))) = do
+  let (ds, js) = partition isDF is
+  (Tensor s ys _) <- tTranspose (js ++ ds) (Tensor s xs is)
+  return (Value (TensorData (Tensor s ys js)))
+ where
+  isDF (DFscript _ _) = True
+  isDF _ = False
+removeDFscripts whnf = return whnf
+
 tMap :: HasTensor a => (a -> EgisonM a) -> (Tensor a) -> EgisonM (Tensor a)
-tMap f (Tensor ns xs js) = do
+tMap f (Tensor ns xs js') = do
+  let k = fromIntegral $ (length ns) - (length js')
+  let js = (js' ++ (map (DFscript 0) [1..k]))
   xs' <- mapM f (V.toList xs) >>= return . V.fromList
   t <- toTensor (V.head xs')
   case t of
-    (Tensor ns1 _ js1) ->
+    (Tensor ns1 _ js1') -> do
+      let k1 = fromIntegral $ (length ns1) - (length js1')
+      let js1 = (js1' ++ (map (DFscript 0) [1..k1]))
       tContract' $ Tensor (ns ++ ns1) (V.concat (V.toList (V.map tensorElems xs'))) (js ++ js1)
     _ -> return $ Tensor ns xs' js
 tMap f (Scalar x) = f x >>= return . Scalar
@@ -857,7 +898,11 @@
 tMapN f xs = mapM fromTensor xs >>= f >>= return . Scalar
 
 tMap2 :: HasTensor a => (a -> a -> EgisonM a) -> Tensor a -> Tensor a -> EgisonM (Tensor a)
-tMap2 f t1@(Tensor ns1 xs1 js1) t2@(Tensor ns2 xs2 js2) = do
+tMap2 f t1@(Tensor ns1 xs1 js1') t2@(Tensor ns2 xs2 js2') = do
+  let k1 = fromIntegral $ (length ns1) - (length js1')
+  let js1 = (js1' ++ (map (DFscript 0) [1..k1]))
+  let k2 = fromIntegral $ (length ns2) - (length js2')
+  let js2 = (js2' ++ (map (DFscript 0) [1..k2]))
   let (cjs, tjs1, tjs2) = h js1 js2
   t1' <- tTranspose (cjs ++ tjs1) t1
   t2' <- tTranspose (cjs ++ tjs2) t2
@@ -929,8 +974,14 @@
       | otherwise -> throwError $ InconsistentTensorSize
 
 tProduct :: HasTensor a => (a -> a -> EgisonM a) -> (Tensor a) -> (Tensor a) -> EgisonM (Tensor a)
-tProduct f t1@(Tensor ns1 xs1 js1) t2@(Tensor ns2 xs2 js2) = do
+tProduct f t1''@(Tensor ns1 xs1 js1') t2''@(Tensor ns2 xs2 js2') = do
+  let k1 = fromIntegral $ (length ns1) - (length js1')
+  let js1 = (js1' ++ (map (DFscript 0) [1..k1]))
+  let k2 = fromIntegral $ (length ns2) - (length js2')
+  let js2 = (js2' ++ (map (DFscript 0) [1..k2]))
   let (cjs1, cjs2, tjs1, tjs2) = h js1 js2
+  let t1 = (Tensor ns1 xs1 js1)
+  let t2 = (Tensor ns2 xs2 js2)
   case cjs1 of
     [] -> do
       xs' <- mapM (\is -> do let is1 = take (length ns1) is
@@ -948,7 +999,8 @@
                               rt2 <- tIntRef is t2'
                               tProduct f rt1 rt2) (enumTensorIndices cns1)
       let ret = Tensor (cns1 ++ (tSize (head rts'))) (V.concat (map tToVector rts')) ((map g cjs1) ++ tIndex (head rts'))
-      tTranspose ((map g cjs1) ++ tjs1 ++ tjs2) ret
+      ret2 <- tTranspose (uniq ((map g cjs1) ++ tjs1 ++ tjs2)) ret
+      return ret2
  where
   h :: [Index EgisonValue] -> [Index EgisonValue] -> ([Index EgisonValue], [Index EgisonValue], [Index EgisonValue], [Index EgisonValue])
   h js1 js2 = let cjs = filter (\j -> any (p j) js2) js1 in
@@ -963,6 +1015,9 @@
   g :: Index EgisonValue -> Index EgisonValue
   g (Superscript i) = (SupSubscript i)
   g (Subscript i) = (SupSubscript i)
+  uniq :: [Index EgisonValue] -> [Index EgisonValue]
+  uniq [] = []
+  uniq (x:xs) = x:(uniq (delete x xs))
 tProduct f (Scalar x) (Tensor ns xs js) = do
   xs' <- V.mapM (f x) xs
   return $ Tensor ns xs' js
@@ -997,6 +1052,7 @@
   p :: Index EgisonValue -> Index EgisonValue -> Bool
   p (Superscript i) (Superscript j) = i == j
   p (Subscript i) (Subscript j) = i == j
+  p (DFscript i1 j1) (DFscript i2 j2) = (i1 == i2) && (j1 == j2)
   p _ _ = False
 tContract' val = return val
 
@@ -1440,26 +1496,31 @@
   show (Superscript ()) = "~"
   show (Subscript ()) = "_"
   show (SupSubscript ()) = "~_"
+  show (DFscript _ _) = ""
 
 instance Show (Index String) where
   show (Superscript s) = "~" ++ s
   show (Subscript s) = "_" ++ s
   show (SupSubscript s) = "~_" ++ s
+  show (DFscript _ _) = ""
 
 instance Show (Index EgisonExpr) where
   show (Superscript i) = "~" ++ show i
   show (Subscript i) = "_" ++ show i
   show (SupSubscript i) = "~_" ++ show i
+  show (DFscript _ _) = ""
 
 instance Show (Index ScalarData) where
   show (Superscript i) = "~" ++ show i
   show (Subscript i) = "_" ++ show i
   show (SupSubscript i) = "~_" ++ show i
+  show (DFscript _ _) = ""
 
 instance Show (Index EgisonValue) where
   show (Superscript i) = "~" ++ show i
   show (Subscript i) = "_" ++ show i
   show (SupSubscript i) = "~_" ++ show i
+  show (DFscript i j) = "_d" ++ show i ++ show j
 
 instance Show (UserIndex EgisonExpr) where
   show (Userscript i) = "|" ++ show i
diff --git a/lib/core/base.egi b/lib/core/base.egi
--- a/lib/core/base.egi
+++ b/lib/core/base.egi
@@ -41,8 +41,10 @@
 (define $flip (lambda [$fn] (lambda [$x $y] (fn y x))))
 
 (define $ref
-  (lambda [$xa $i]
-    (array-ref xa i)))
+  (lambda [%xa $is]
+    (match is (list integer)
+      {[<nil> xa]
+       [<cons $i $rs> (ref (array-ref xa i) rs)]})))
 
 (define $eq?/m
   (lambda [$a $x $y]
diff --git a/lib/math/algebra/matrix.egi b/lib/math/algebra/matrix.egi
--- a/lib/math/algebra/matrix.egi
+++ b/lib/math/algebra/matrix.egi
@@ -9,7 +9,7 @@
 (define $M.b.*
   (lambda [%m1 %m2]
     (with-symbols {j}
-      (. m1~#~j m2_j))))
+      (. m1~#~j m2_j_#))))
 
 (define $M.*'
   (cambda $ms
diff --git a/lib/math/algebra/tensor.egi b/lib/math/algebra/tensor.egi
--- a/lib/math/algebra/tensor.egi
+++ b/lib/math/algebra/tensor.egi
@@ -4,6 +4,10 @@
 ;;;;;
 ;;;;;
 
+(define $tensor-order
+  (lambda [%A]
+    (length (tensor-size A))))
+
 (define $unit-tensor
   (lambda [$ns]
     (generate-tensor kronecker-delta ns)))
@@ -14,7 +18,7 @@
 
 (define $zero-tensor
   (lambda [$ns]
-    (generate-tensor 1#0 ns)))
+    (generate-tensor (cambda $xs 0) ns)))
 
 (define $b..' (lambda [%t1 %t2] (contract +' (*' t1 t2))))
 (define $b.. (lambda [%t1 %t2] (contract + (* t1 t2))))
diff --git a/lib/math/analysis/derivative.egi b/lib/math/analysis/derivative.egi
--- a/lib/math/analysis/derivative.egi
+++ b/lib/math/analysis/derivative.egi
@@ -53,6 +53,14 @@
 (define $grad ∇)
 
 (define $taylor-expansion
+  (lambda [$f $x $a]
+    (map2 *
+          (map 1#(/ (** (- x a) %1) (fact %1)) nats0)
+          (map (substitute {[x a]} $) (iterate (∂/∂ $ x) f)))))
+
+(define $maclaurin-expansion (taylor-expansion $ $ 0))
+
+(define $multivariate-taylor-expansion
   (lambda [%f %xs %as]
     (with-symbols {h}
       (let {[$hs (generate-tensor 1#h_%1 (tensor-size xs))]}
@@ -62,7 +70,7 @@
                             1#(V.substitute hs (with-symbols {i} (- xs_i as_i)) %1))
                    (iterate (compose 1#(∇ %1 xs) 1#(V.* hs %1)) f)))))))
 
-(define $maclaurin-expansion
+(define $multivariate-maclaurin-expansion
   (lambda [%f %xs]
     (multivariate-taylor-expansion f xs (tensor-map 1#0 xs))))
 
diff --git a/lib/math/common/constants.egi b/lib/math/common/constants.egi
--- a/lib/math/common/constants.egi
+++ b/lib/math/common/constants.egi
@@ -5,3 +5,5 @@
 ;;;;;
 
 (define $pi π)
+
+(define $Minkowski-metric [| [| -1 0 0 0 |] [| 0 1 0 0 |] [| 0 0 1 0 |] [| 0 0 0 1 |] |])
diff --git a/lib/math/common/functions.egi b/lib/math/common/functions.egi
--- a/lib/math/common/functions.egi
+++ b/lib/math/common/functions.egi
@@ -4,6 +4,18 @@
 ;;;;;
 ;;;;;
 
+(define $abs
+  (lambda [$x]
+    (if (rational? x)
+      (b.abs x)
+      x)))
+
+(define $neg
+  (lambda [$x]
+    (if (rational? x)
+      (b.neg x)
+      (* -1 x))))
+
 (define $exp
   (lambda [$x]
     (if (float? x)
@@ -100,3 +112,29 @@
     (* n
        (product (map (lambda [$p] (- 1 (/ 1 p)))
                      (unique (p-f n)))))))
+
+(define $ε
+  (memoized-lambda [$n]
+    (let {[[$es $os] (even-and-odd-permutations' n)]}
+      (generate-tensor
+        (cambda $is
+          (if (member? is es)
+            1
+            (if (member? is os)
+              -1
+              0)))
+        (take n (repeat1 n))))))
+
+(define $ε'
+  (memoized-lambda [$n $k]
+    (let {[[$es $os] (even-and-odd-permutations' n)]}
+      (generate-tensor
+        (cambda $is
+          (match (drop k is) (list integer)
+            {[<join _ <cons $x <join _ <cons ?1#(lt? %1 x) _>>>> 0]
+             [_ (if (member? is es)
+                  1
+                  (if (member? is os)
+                    -1
+                    0))]}))
+        (take n (repeat1 n))))))
diff --git a/lib/math/geometry/3d-euclidean-space.egi b/lib/math/geometry/3d-euclidean-space.egi
new file mode 100644
--- /dev/null
+++ b/lib/math/geometry/3d-euclidean-space.egi
@@ -0,0 +1,8 @@
+(define $coordinates {x y z})
+
+(define $metric
+  (generate-tensor
+    (match-lambda [integer integer]
+      {[[$n ,n] 1]
+       [[_ _] 0]})
+    {3 3}))
diff --git a/lib/math/geometry/4d-euclidean-space.egi b/lib/math/geometry/4d-euclidean-space.egi
new file mode 100644
--- /dev/null
+++ b/lib/math/geometry/4d-euclidean-space.egi
@@ -0,0 +1,8 @@
+(define $coordinates {x y z w})
+
+(define $metric
+  (generate-tensor
+    (match-lambda [integer integer]
+      {[[$n ,n] 1]
+       [[_ _] 0]})
+    {4 4}))
diff --git a/lib/math/geometry/differential-form.egi b/lib/math/geometry/differential-form.egi
new file mode 100644
--- /dev/null
+++ b/lib/math/geometry/differential-form.egi
@@ -0,0 +1,18 @@
+;(define $d
+;  (lambda [%X]
+;    !(∂ X)))
+
+(define $wedge
+  (lambda [%X %Y]
+    !(. X Y)))
+
+(define $Lie.wedge
+  (lambda [%X %Y]
+    (- !(. X Y) !(. Y X))))
+
+;(define $hodge
+;  (lambda [%A]
+;    (* (sqrt (abs (det g_#_#)))
+;       (foldl . ε_i.1..._i.n A_j.1..._j.k
+;              (map (lambda [$l] g~i.l~j.l) (between 1 k))))
+;    ))
diff --git a/lib/math/geometry/minkowski-space.egi b/lib/math/geometry/minkowski-space.egi
new file mode 100644
--- /dev/null
+++ b/lib/math/geometry/minkowski-space.egi
@@ -0,0 +1,12 @@
+(define $coordinates {t x y z})
+
+(define $metric
+  (generate-tensor
+    (match-lambda [integer integer]
+      {[[,1 ,1] -1]
+       [[$n ,n] 1]
+       [[_ _] 0]})
+    {4 4}))
+
+(define $.
+  (lambda [$
diff --git a/sample/bellman-ford.egi b/sample/bellman-ford.egi
new file mode 100644
--- /dev/null
+++ b/sample/bellman-ford.egi
@@ -0,0 +1,19 @@
+; initiate a distance graph
+(define $A
+  [|[| 0 19 36 66 99 65 |]
+    [| 19 0 25 59 64 31 |]
+    [| 36 25 0 84 48 28 |]
+    [| 66 59 84 0 59 29 |]
+    [| 99 64 48 59 0 9 |]
+    [| 65 31 28 29 9 0 |]|])
+
+(define $G.* (lambda [%t1 %t2] (with-symbols {i} (contract min (+ t1~#_i t2~i_#)))))
+
+(match (iterate (lambda [%P] (G.* P A)) A) (list something)
+  {[<join _ <cons $P <cons ,P _>>> P]})
+;[|[| 0 19 36 66 59 50 |]
+;  [| 19 0 25 59 40 31 |]
+;  [| 36 25 0 57 37 28 |]
+;  [| 66 59 57 0 38 29 |]
+;  [| 59 40 37 38 0 9 |]
+;  [| 50 31 28 29 9 0 |]|]
diff --git a/sample/mickey.egi b/sample/mickey.egi
new file mode 100644
--- /dev/null
+++ b/sample/mickey.egi
@@ -0,0 +1,10 @@
+(define $mickey'
+  (lambda [$cs]
+    (match cs (list char)
+      {[<snoc $z <snoc $y <snoc $x (& <snoc _ _> $hs)>>>
+        {@(mickey' hs) c#, x y z}]
+       [_ cs]})))
+
+(define $mickey (lambda [$s] (pack (mickey' (unpack s)))))
+
+(mickey "10000000000")