diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c) 2011-2016, Satoshi Egi
+Copyright (c) 2011-2017, Satoshi Egi
 
 Permission is hereby granted, free of charge, to any person obtaining
 a copy of this software and associated documentation files (the "Software"),
diff --git a/egison.cabal b/egison.cabal
--- a/egison.cabal
+++ b/egison.cabal
@@ -1,5 +1,5 @@
 Name:                egison
-Version:             3.6.4
+Version:             3.6.5
 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.
@@ -68,7 +68,7 @@
   location: https://github.com/egison/egison.git
   
 Library
-  Build-Depends:   base >= 4.0 && < 5, array, random, containers, unordered-containers, haskeline, transformers, mtl, parsec >= 3.0, directory, ghc, ghc-paths, text, regex-tdfa, process, vector, monad-parallel
+  Build-Depends:   base >= 4.0 && < 5, array, random, containers, unordered-containers, haskeline, transformers, mtl, parsec >= 3.0, directory, ghc, ghc-paths, text, regex-tdfa, process, vector, parallel
   Hs-Source-Dirs:  hs-src
   Exposed-Modules:
                    Language.Egison
diff --git a/hs-src/Interpreter/egison.hs b/hs-src/Interpreter/egison.hs
--- a/hs-src/Interpreter/egison.hs
+++ b/hs-src/Interpreter/egison.hs
@@ -2,7 +2,7 @@
 
 import Prelude hiding (catch)
 import Control.Exception ( AsyncException(..), catch )
-import Control.Monad.Error
+import Control.Monad.Except
 
 import qualified Data.Text as T
 import Data.Char
@@ -216,7 +216,7 @@
 
 showBanner :: IO ()
 showBanner = do
-  putStrLn $ "Egison Version " ++ showVersion version ++ " (C) 2011-2016 Satoshi Egi"
+  putStrLn $ "Egison Version " ++ showVersion version ++ " (C) 2011-2017 Satoshi Egi"
   putStrLn $ "https://www.egison.org"
   putStrLn $ "Welcome to Egison Interpreter!"
 --  putStrLn $ "** Information **"
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
@@ -43,10 +43,9 @@
 
 import Control.Arrow
 import Control.Applicative
-import Control.Monad.Error hiding (mapM)
+import Control.Monad.Except hiding (mapM)
 import Control.Monad.State hiding (mapM, state)
 import Control.Monad.Trans.Maybe
-import qualified Control.Monad.Parallel as MP
 
 import Data.Sequence (Seq, ViewL(..), ViewR(..), (><))
 import qualified Data.Sequence as Sq
@@ -127,8 +126,8 @@
   collectDefs (expr:exprs) bindings rest =
     case expr of
       Define name expr -> collectDefs exprs ((show name, expr) : bindings) rest
-      Load _ -> throwError $ strMsg "No IO support"
-      LoadFile _ -> throwError $ strMsg "No IO support"
+      Load _ -> throwError $ Default "No IO support"
+      LoadFile _ -> throwError $ Default "No IO support"
       _ -> collectDefs exprs bindings (expr : rest)
   collectDefs [] bindings rest = return (bindings, reverse rest)
 
@@ -205,7 +204,7 @@
   return . Intermediate . IArray $ Array.listArray (1, toInteger (length exprs)) refs'
 
 evalExpr env (VectorExpr exprs) = do
-  whnfs <- MP.mapM (evalExpr env) exprs
+  whnfs <- mapM (evalExpr env) exprs
   case whnfs of
     ((Intermediate (ITensor (Tensor _ _ _))):_) -> do
       ret <- mapM toTensor whnfs >>= tConcat' >>= fromTensor
@@ -217,7 +216,7 @@
   nsWhnf <- evalExpr env nsExpr
   ns <- ((fromCollection nsWhnf >>= fromMList >>= mapM evalRef >>= mapM fromWHNF) :: EgisonM [Integer])
   xsWhnf <- evalExpr env xsExpr
-  xs <- fromCollection xsWhnf >>= fromMList >>= MP.mapM evalRef
+  xs <- fromCollection xsWhnf >>= fromMList >>= mapM evalRef
   supWhnf <- evalExpr env supExpr
   sup <- fromCollection supWhnf >>= fromMList >>= mapM evalRefDeep -- >>= mapM extractScalar'
   subWhnf <- evalExpr env subExpr
@@ -358,6 +357,18 @@
   genVar :: State Int String
   genVar = modify (1+) >> gets (('#':) . show)
 
+evalExpr env (TransposeExpr vars expr) = do
+  syms <- evalExpr env vars >>= collectionToList
+  whnf <- evalExpr env expr
+  case whnf of
+    (Intermediate (ITensor t)) -> do
+      t' <- tTranspose' syms t
+      return (Intermediate (ITensor t'))
+    (Value (TensorData t)) -> do
+      t' <- tTranspose' syms t
+      return (Value (TensorData t'))
+    _ -> return whnf
+
 evalExpr env (WithSymbolsExpr vars expr) = do
   symId <- fresh
   syms <- mapM (\var -> (newEvaluatedObjectRef (Value (symbolScalarData symId var)))) vars
@@ -467,7 +478,7 @@
             case result of
               MCons bindings _ -> evalExpr (extendEnv env bindings) expr
               MNil -> cont
-      foldr tryMatchClause (throwError $ strMsg "failed pattern match") clauses
+      foldr tryMatchClause (throwError $ Default "failed pattern match") clauses
 
 evalExpr env (SeqExpr expr1 expr2) = do
   evalExprDeep env expr1
@@ -548,7 +559,7 @@
   size'' <- collectionToList size'
   ns <- (mapM fromEgison size'') :: EgisonM [Integer]
   fn <- evalExpr env fnExpr
-  xs <-  MP.mapM (\ms -> applyFunc env fn (Value (makeTuple ms))) (map (\ms -> map toEgison ms) (enumTensorIndices ns))
+  xs <-  mapM (\ms -> applyFunc env fn (Value (makeTuple ms))) (map (\ms -> map toEgison ms) (enumTensorIndices ns))
   case (ns, xs) of
     ([1], x:[]) -> return $ x
     _ -> fromTensor (Tensor ns (V.fromList xs) [])
@@ -621,6 +632,19 @@
     yRef <- newEvaluatedObjectRef y
     applyFunc env fn (Intermediate (ITuple [xRef, yRef]))
 
+evalExpr env (ParExpr expr1 expr2) = undefined
+evalExpr env (PseqExpr expr1 expr2) = undefined
+
+evalExpr env (PmapExpr fnExpr cExpr) = do
+  fn <- evalExpr env fnExpr
+  xs <- evalExpr env cExpr >>= collectionToList
+  ys <- parallelMapM (applyFunc' env fn) xs
+  return $ Value $ Collection (Sq.fromList ys)
+ where
+  applyFunc' :: Env -> WHNFData -> EgisonValue -> EgisonM EgisonValue
+  applyFunc' env fn x = applyFunc env fn (Value x) >>= evalWHNF
+  
+
 evalExpr _ SomethingExpr = return $ Value Something
 evalExpr _ UndefinedExpr = return $ Value Undefined
 evalExpr _ expr = throwError $ NotImplemented ("evalExpr for " ++ show expr)
@@ -671,7 +695,7 @@
 evalWHNF (Intermediate (ITuple [ref])) = evalRefDeep ref
 evalWHNF (Intermediate (ITuple refs)) = Tuple <$> mapM evalRefDeep refs
 evalWHNF (Intermediate (ITensor (Tensor ns whnfs js))) = do
-  vals <- MP.mapM evalWHNF (V.toList whnfs)
+  vals <- mapM evalWHNF (V.toList whnfs)
   return $ TensorData $ Tensor ns (V.fromList vals) js
 --  vals <- mapM evalWHNF whnfs
 --  return $ TensorData $ Tensor ns vals js
@@ -977,7 +1001,7 @@
   let env' = extendEnvForNonLinearPatterns env bindings loops
   case pattern of
     NotPat _ -> throwError $ EgisonBug "should not reach here (not pattern)"
-    VarPat _ -> throwError $ strMsg $ "cannot use variable except in pattern function:" ++ show pattern
+    VarPat _ -> throwError $ Default $ "cannot use variable except in pattern function:" ++ show pattern
 
     LetPat bindings' pattern' ->
       let extractBindings ([name], expr) =
@@ -1020,7 +1044,7 @@
           return $ msingleton $ MState env ((LoopPatContext (name, startNumRef) endsRef endPat pat pat'):loops) bindings ((MAtom ContPat target matcher):trees)
     ContPat ->
       case loops of
-        [] -> throwError $ strMsg "cannot use cont pattern except in loop pattern"
+        [] -> throwError $ Default "cannot use cont pattern except in loop pattern"
         LoopPatContext (name, startNumRef) endsRef endPat pat pat' : loops' -> do
           startNumWhnf <- evalRef startNumRef
           startNum <- fromWHNF startNumWhnf :: (EgisonM Integer)
@@ -1066,7 +1090,7 @@
               if not (length patterns == length matchers) then throwError $ ArgumentsNum (length patterns) (length matchers) else return ()
               let trees' = zipWith3 MAtom patterns targets matchers ++ trees
               return $ msingleton $ MState env loops bindings trees'
-            _ ->  throwError $ strMsg $ "should not reach here. matcher: " ++ show matcher ++ ", pattern:  " ++ show pattern
+            _ ->  throwError $ Default $ "should not reach here. matcher: " ++ show matcher ++ ", pattern:  " ++ show pattern
 
         Something ->
           case pattern of
@@ -1102,18 +1126,18 @@
                   keys <- return $ HL.keys hash
                   vals <- mapM (newEvaluatedObjectRef . Value) $ HL.elems hash
                   updateHash indices (Intermediate $ IIntHash $ HL.fromList $ zip keys vals)
-                updateHash _ v = throwError $ strMsg $ "expected hash value: " ++ show v
+                updateHash _ v = throwError $ Default $ "expected hash value: " ++ show v
                 subst :: (Eq a) => a -> b -> [(a, b)] -> [(a, b)]
                 subst k nv ((k', v'):xs) | k == k'   = (k', nv):(subst k nv xs)
                                          | otherwise = (k', v'):(subst k nv xs)
                 subst _ _ [] = []
-            IndexedPat pattern indices -> throwError $ strMsg ("invalid indexed-pattern: " ++ show pattern) 
+            IndexedPat pattern indices -> throwError $ Default ("invalid indexed-pattern: " ++ show pattern) 
             TuplePat patterns -> do
               targets <- fromTupleWHNF target
               if not (length patterns == length targets) then throwError $ ArgumentsNum (length patterns) (length targets) else return ()
               let trees' = zipWith3 MAtom patterns targets (take (length patterns) (repeat Something)) ++ trees
               return $ msingleton $ MState env loops bindings trees'
-            _ -> throwError $ strMsg "something can only match with a pattern variable"
+            _ -> throwError $ Default "something can only match with a pattern variable"
         _ ->  throwError $ EgisonBug $ "should not reach here. matcher: " ++ show matcher ++ ", pattern:  " ++ show pattern
 
 inductiveMatch :: Env -> EgisonPattern -> WHNFData -> Matcher ->
@@ -1137,8 +1161,8 @@
               evalExpr env expr >>= fromCollection
             _ -> cont
       _ -> cont
-  failPPPatternMatch = throwError $ strMsg "failed primitive pattern pattern match"
-  failPDPatternMatch = throwError $ strMsg "failed primitive data pattern match"
+  failPPPatternMatch = throwError $ Default "failed primitive pattern pattern match"
+  failPDPatternMatch = throwError $ Default "failed primitive data pattern match"
 
 primitivePatPatternMatch :: Env -> PrimitivePatPattern -> EgisonPattern ->
                             MatchM ([EgisonPattern], [Binding])
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
@@ -25,17 +25,17 @@
 import Data.List (span)
 import Data.Set (Set)
 import Data.Char (toUpper)
-import Control.Monad.Error
+import Control.Monad.Except
 import Control.Monad.Reader
 import Language.Egison.Types
 
 type Subst = [(String, EgisonExpr)]
 
-newtype DesugarM a = DesugarM { unDesugarM :: ReaderT Subst (ErrorT EgisonError Fresh) a }
+newtype DesugarM a = DesugarM { unDesugarM :: ReaderT Subst (ExceptT EgisonError Fresh) a }
   deriving (Functor, Applicative, Monad, MonadError EgisonError, MonadFresh, MonadReader Subst)
 
 runDesugarM :: DesugarM a -> Fresh (Either EgisonError a)
-runDesugarM = runErrorT . flip runReaderT [] . unDesugarM
+runDesugarM = runExceptT . flip runReaderT [] . unDesugarM
 
 desugarTopExpr :: EgisonTopExpr -> EgisonM EgisonTopExpr
 desugarTopExpr (Define name expr) = do
@@ -302,6 +302,26 @@
   t1Expr' <- desugar t1Expr
   t2Expr' <- desugar t2Expr
   return $ TensorMap2Expr fnExpr' t1Expr' t2Expr'
+
+desugar (TransposeExpr vars expr) = do
+  vars' <- desugar vars
+  expr' <- desugar expr
+  return $ TransposeExpr vars' expr'
+
+desugar (ParExpr expr1 expr2) = do
+  expr1' <- desugar expr1
+  expr2' <- desugar expr2
+  return $ ParExpr expr1' expr2'
+
+desugar (PseqExpr expr1 expr2) = do
+  expr1' <- desugar expr1
+  expr2' <- desugar expr2
+  return $ PseqExpr expr1' expr2'
+
+desugar (PmapExpr expr1 expr2) = do
+  expr1' <- desugar expr1
+  expr2' <- desugar expr2
+  return $ PmapExpr expr1' expr2'
 
 desugar (ApplyExpr expr0 expr1) = do
   expr0' <- desugar expr0
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
@@ -26,7 +26,7 @@
 
 import Prelude hiding (mapM)
 import Control.Monad.Identity hiding (mapM)
-import Control.Monad.Error hiding (mapM)
+import Control.Monad.Except hiding (mapM)
 import Control.Monad.State hiding (mapM)
 import Control.Applicative ((<$>), (<*>), (*>), (<*), pure)
 
@@ -98,7 +98,7 @@
 loadFile :: FilePath -> EgisonM [EgisonTopExpr]
 loadFile file = do
   doesExist <- liftIO $ doesFileExist file
-  unless doesExist $ throwError $ strMsg ("file does not exist: " ++ file)
+  unless doesExist $ throwError $ Default ("file does not exist: " ++ file)
   input <- liftIO $ readFile file
   exprs <- readTopExprs $ shebang input
   concat <$> mapM  recursiveLoad exprs
@@ -125,8 +125,8 @@
 --
 topExpr :: Parser EgisonTopExpr
 topExpr = try (Test <$> expr)
-      <|> try (parens (defineExpr
-                   <|> redefineExpr
+      <|> try defineExpr
+      <|> try (parens (redefineExpr
                    <|> testExpr
                    <|> executeExpr
                    <|> loadFileExpr
@@ -134,7 +134,28 @@
       <?> "top-level expression"
 
 defineExpr :: Parser EgisonTopExpr
-defineExpr = keywordDefine >> Define <$> varNameWithIndexType <*> expr
+defineExpr = try (parens (keywordDefine >> Define <$> varNameWithIndexType <*> expr))
+         <|> try (parens (do keywordDefine
+                             (VarWithIndices name is) <- varNameWithIndices
+                             body <- expr
+                             return $ Define (Var name (map f is)) (WithSymbolsExpr (map g is) (TransposeExpr (CollectionExpr (map (ElementExpr . h) is)) body))))
+--defineExpr = try 
+--                 (do keywordDefine
+--                     (VarWithIndices name is) <- varNameWithIndices
+--                     body <- expr
+--                     return $ Define (Var name (map f is)) (WithSymbolsExpr (map g is) (TransposeExpr (CollectionExpr (map (ElementExpr . h) is)) body)))
+--            <|> (keywordDefine >> Define <$> varNameWithIndexType <*> expr)
+--defineExpr = (keywordDefine >> Define <$> varNameWithIndexType <*> expr)
+ where
+  f (Superscript _) = Superscript ()
+  f (Subscript _) = Subscript ()
+  f (SupSubscript _) = SupSubscript ()
+  g (Superscript i) = i
+  g (Subscript i) = i
+  g (SupSubscript i) = i
+  h (Superscript i) = (VarExpr i)
+  h (Subscript i) = (VarExpr i)
+  h (SupSubscript i) = (VarExpr i)
 
 redefineExpr :: Parser EgisonTopExpr
 redefineExpr = (keywordRedefine <|> keywordSet) >> Redefine <$> varNameWithIndexType <*> expr
@@ -216,6 +237,10 @@
                          <|> tensorContractExpr
                          <|> tensorMapExpr
                          <|> tensorMap2Expr
+                         <|> transposeExpr
+                         <|> parExpr
+                         <|> pseqExpr
+                         <|> pmapExpr
                          )
              <?> "expression")
 
@@ -425,6 +450,17 @@
 indexType = try (char '~' >> return (Superscript ()))
         <|> try (char '_' >> return (Subscript ()))
 
+varNameWithIndices :: Parser VarWithIndices
+varNameWithIndices = P.lexeme lexer (do
+  char '$'
+  name <- ident
+  is <- many indexForVar
+  return $ VarWithIndices name is)
+
+indexForVar :: Parser (Index String)
+indexForVar = try (char '~' >> Superscript <$> ident)
+        <|> try (char '_' >> Subscript <$> ident)
+
 argNames :: Parser [Arg]
 argNames = return <$> argName
             <|> brackets (sepEndBy argName whiteSpace) 
@@ -454,16 +490,19 @@
   case vars of
     [] -> return . ApplyExpr func . TupleExpr $ rights args
     _ | all null vars ->
-        let genVar = modify (1+) >> gets (VarExpr . (':':) . show)
-            args' = evalState (mapM (either (const genVar) return) args) 0
-        in return . LambdaExpr (map ScalarArg (annonVars $ length vars)) . ApplyExpr func $ TupleExpr args'
+        let args' = rights args
+            args'' = map f (zip args (annonVars 1 (length args)))
+            args''' = map (VarExpr . (either id id)) args''
+        in return $ ApplyExpr (LambdaExpr (map ScalarArg (rights args'')) (LambdaExpr (map ScalarArg (lefts args'')) $ ApplyExpr func $ TupleExpr args''')) $ TupleExpr args'
       | all (not . null) vars ->
         let ns = Set.fromList $ map read vars
             n = Set.size ns
         in if Set.findMin ns == 1 && Set.findMax ns == n
              then
-               let args' = map (either (VarExpr . (':':)) id) args
-               in return . LambdaExpr (map ScalarArg (annonVars n)) . ApplyExpr func $ TupleExpr args'
+               let args' = rights args
+                   args'' = map g (zip args (annonVars (n + 1) (length args)))
+                   args''' = map (VarExpr . (either id id)) args''
+               in return $ ApplyExpr (LambdaExpr (map ScalarArg (rights args'')) (LambdaExpr (map ScalarArg (annonVars 1 n)) $ ApplyExpr func $ TupleExpr args''')) $ TupleExpr args'
              else fail "invalid partial application"
       | otherwise -> fail "invalid partial application"
  where
@@ -471,7 +510,11 @@
   arg = try (Right <$> expr)
          <|> char '$' *> (Left <$> option "" index)
   index = (:) <$> satisfy (\c -> '1' <= c && c <= '9') <*> many digit
-  annonVars n = take n $ map ((':':) . show) [1..]
+  annonVars m n = take n $ map ((':':) . show) [m..]
+  f ((Left _), var) = Left var
+  f ((Right _), var) = Right var
+  g ((Left arg), _) = Left (':':arg)
+  g ((Right _), var) = Right var
 
 partialExpr :: Parser EgisonExpr
 partialExpr = PartialExpr <$> read <$> index <*> (char '#' >> expr)
@@ -517,6 +560,18 @@
 tensorMap2Expr :: Parser EgisonExpr
 tensorMap2Expr = keywordTensorMap2 >> TensorMap2Expr <$> expr <*> expr <*> expr
 
+transposeExpr :: Parser EgisonExpr
+transposeExpr = keywordTranspose >> TransposeExpr <$> expr <*> expr
+
+parExpr :: Parser EgisonExpr
+parExpr = keywordPar >> ParExpr <$> expr <*> expr
+
+pseqExpr :: Parser EgisonExpr
+pseqExpr = keywordPseq >> PseqExpr <$> expr <*> expr
+
+pmapExpr :: Parser EgisonExpr
+pmapExpr = keywordPmap >> PmapExpr <$> expr <*> expr
+
 -- Patterns
 
 pattern :: Parser EgisonPattern
@@ -749,6 +804,10 @@
   , "contract"
   , "tensor-map"
   , "tensor-map2"
+  , "transpose"
+  , "par"
+  , "pseq"
+  , "pmap"
   , "something"
   , "undefined"]
   
@@ -824,6 +883,10 @@
 keywordTensorContract       = reserved "contract"
 keywordTensorMap            = reserved "tensor-map"
 keywordTensorMap2           = reserved "tensor-map2"
+keywordTranspose            = reserved "transpose"
+keywordPar                  = reserved "par"
+keywordPseq                 = reserved "pseq"
+keywordPmap                 = reserved "pmap"
 
 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
@@ -11,7 +11,7 @@
 module Language.Egison.Primitives (primitiveEnv, primitiveEnvNoIO) where
 
 import Control.Arrow
-import Control.Monad.Error
+import Control.Monad.Except
 import Control.Monad.Trans.Maybe
 import Control.Applicative ((<$>), (<*>), (*>), (<*), pure)
 
@@ -145,6 +145,8 @@
              , ("from-math-expr", fromScalarData)
              , ("to-math-expr", toScalarData)
              , ("to-math-expr'", toScalarData)
+             , ("decons-user-scripts", deconsUserScripts)
+             , ("append-user-scripts", appendUserScripts)
 
              , ("modulo",    integerBinaryOp mod)
              , ("quotient",   integerBinaryOp quot)
@@ -351,6 +353,18 @@
  where
   toScalarData' val = egisonToScalarData val >>= return . ScalarData . mathNormalize'
 
+appendUserScripts :: PrimitiveFunc
+appendUserScripts = twoArgs $ appendUserScripts'
+ where
+  appendUserScripts' v (Collection is) = do
+    let is' = map Userscript (toList is)
+    return $ UserIndexedData v is'
+
+deconsUserScripts :: PrimitiveFunc
+deconsUserScripts = oneArg $ deconsUserScripts'
+ where
+  deconsUserScripts' (UserIndexedData v is) = return $ Tuple [v, Collection (Sq.fromList (map (\(Userscript i) -> i) is))]
+  deconsUserScripts' v = return $ Tuple [v, Collection (Sq.fromList [])]
 --
 -- Pred
 --
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
@@ -44,6 +44,7 @@
     , tIndex
     , tref
     , enumTensorIndices
+    , tTranspose'
     , tMap
     , tMap2
     , tMapN
@@ -82,6 +83,7 @@
     -- * Environment
     , Env (..)
     , Var (..)
+    , VarWithIndices (..)
     , Binding (..)
     , nullEnv
     , extendEnv
@@ -99,6 +101,7 @@
     , liftError
     -- * Monads
     , EgisonM (..)
+    , parallelMapM
     , runEgisonM
     , liftEgisonM
     , fromEgisonM
@@ -144,16 +147,16 @@
 import Prelude hiding (foldr, mappend, mconcat)
 
 import Control.Exception
+import Control.Parallel
 import Data.Typeable
 
 import Control.Applicative
-import Control.Monad.Error
+import Control.Monad.Except
 import Control.Monad.State
 import Control.Monad.Reader (ReaderT)
 import Control.Monad.Writer (WriterT)
 import Control.Monad.Identity
 import Control.Monad.Trans.Maybe
-import qualified Control.Monad.Parallel as MP
 
 import Data.Monoid (Monoid)
 import qualified Data.HashMap.Lazy as HL
@@ -166,7 +169,7 @@
 import Data.HashMap.Strict (HashMap)
 import qualified Data.HashMap.Strict as HashMap
 
-import Data.List (intercalate, sort, sortBy, findIndex, splitAt, (\\), elem, delete, deleteBy, any)
+import Data.List (intercalate, sort, sortBy, find, findIndex, splitAt, (\\), elem, delete, deleteBy, any)
 import Data.Text (Text)
 import qualified Data.Text as T
 
@@ -253,11 +256,16 @@
   | ArrayBoundsExpr EgisonExpr
   | ArrayRefExpr EgisonExpr EgisonExpr
 
+  | ParExpr EgisonExpr EgisonExpr
+  | PseqExpr EgisonExpr EgisonExpr
+  | PmapExpr EgisonExpr EgisonExpr
+
   | GenerateTensorExpr EgisonExpr EgisonExpr
   | TensorExpr EgisonExpr EgisonExpr EgisonExpr EgisonExpr
   | TensorContractExpr EgisonExpr EgisonExpr
   | TensorMapExpr EgisonExpr EgisonExpr
   | TensorMap2Expr EgisonExpr EgisonExpr EgisonExpr
+  | TransposeExpr EgisonExpr EgisonExpr
 
   | SomethingExpr
   | UndefinedExpr
@@ -401,8 +409,19 @@
   (Term a []) == (Term b [])
     | a /= b =  False
     | otherwise = True
+  (Term a (((Quote x),n):xs)) == (Term b ys)
+    | (a /= b) && (a /= (negate b)) =  False
+    | otherwise = case findIndex ((==) ((Quote x),n)) ys of
+                    Just i -> let (hs, _:ts) = splitAt i ys in
+                                (Term a xs) == (Term b (hs ++ ts))
+                    Nothing -> case findIndex ((==) ((Quote (mathNegate x)),n)) ys of
+                                 Just i -> let (hs, _:ts) = splitAt i ys in
+                                             if even n
+                                               then (Term a xs) == (Term b (hs ++ ts))
+                                               else (Term (negate a) xs) == (Term b (hs ++ ts))
+                                 Nothing -> False
   (Term a (x:xs)) == (Term b ys)
-    | a /= b =  False
+    | (a /= b) && (a /= (negate b)) =  False
     | otherwise = case findIndex ((==) x) ys of
                     Just i -> let (hs, _:ts) = splitAt i ys in
                                 (Term a xs) == (Term b (hs ++ ts))
@@ -517,39 +536,63 @@
 mathNormalize' :: ScalarData -> ScalarData
 mathNormalize' mExpr = mathDivide (mathRemoveZero (mathFold (mathRemoveZeroSymbol mExpr)))
 
+termsGcd :: [TermExpr] -> TermExpr
+termsGcd (t:ts) = f t ts
+ where
+  f :: TermExpr -> [TermExpr] -> TermExpr
+  f ret [] =  ret
+  f (Term a xs) ((Term b ys):ts) =
+    f (Term (gcd a b) (g xs ys)) ts
+  g :: [(SymbolExpr, Integer)] -> [(SymbolExpr, Integer)] -> [(SymbolExpr, Integer)]
+  g [] ys = []
+  g ((x, n):xs) ys = let (z, m) = h (x, n) ys in
+    if m == 0 then g xs ys else ((z, m):(g xs ys))
+  h :: (SymbolExpr, Integer) -> [(SymbolExpr, Integer)] -> (SymbolExpr, Integer)
+  h (x, n) [] = (x, 0)
+  h ((Quote x), n) (((Quote y), m):ys) = if x == y
+                                        then ((Quote x), (min n m))
+                                        else if x == (mathNegate y)
+                                             then ((Quote x), (min n m))
+                                             else h ((Quote x), n) ys
+  h (x, n) ((y, m):ys) = if x == y
+                         then (x, (min n m))
+                         else h (x, n) ys
+
 mathDivide :: ScalarData -> ScalarData
 mathDivide (Div (Plus ts1) (Plus [])) = (Div (Plus ts1) (Plus []))
 mathDivide (Div (Plus []) (Plus ts2)) = (Div (Plus []) (Plus ts2))
-mathDivide (Div (Plus ts1) (Plus [(Term a xs)])) =
-  case f (Term a xs) ts1 of
-    Nothing -> (Div (Plus ts1) (Plus [(Term a xs)]))
-    Just ts1' -> (Div (Plus ts1') (Plus [(Term 1 [])]))
- where
-  f :: TermExpr -> [TermExpr] -> Maybe [TermExpr]
-  f _ [] = Just []
-  f t1 (t:ts) = do t' <- f' t1 t
-                   ts' <- f t1 ts
-                   return (t':ts')
-  f' :: TermExpr -> TermExpr -> Maybe TermExpr
-  f' (Term a xs) (Term b ys) =
-    if b `mod` a == 0
-      then do ys' <- g xs ys
-              return (Term (b `quot` a) ys')
-      else Nothing
-  g :: [(SymbolExpr, Integer)] -> [(SymbolExpr, Integer)] -> Maybe [(SymbolExpr, Integer)]
-  g [] ys = Just ys
-  g ((x,n):xs) ys = do
-    ys' <- g' (x,n) ys
-    g xs ys'
-  g' :: (SymbolExpr, Integer) -> [(SymbolExpr, Integer)] -> Maybe [(SymbolExpr, Integer)]
-  g' (x, n) [] = Nothing
-  g' (x, n) ((y, m):ys) = do
-    if (x == y && n <= m)
-      then Just ((y, (m - n)):ys)
-      else do ys' <- g' (x,n) ys
-              return ((y,m):ys')
-mathDivide (Div (Plus ts1) (Plus ts2)) = (Div (Plus ts1) (Plus ts2))
+mathDivide (Div (Plus ts1) (Plus ts2)) =
+  let z = termsGcd (ts1 ++ ts2) in
+  case z of
+--    (Term 1 []) -> (Div (Plus ts1) (Plus ts2))
+    (Term c zs) -> case ts2 of
+      [(Term a _)] -> if a < 0
+                      then (Div (Plus (map (\t -> mathDivideTerm t (Term (-1 * c) zs)) ts1)) (Plus (map (\t -> mathDivideTerm t (Term (-1 * c) zs)) ts2)))
+                      else (Div (Plus (map (\t -> mathDivideTerm t z) ts1)) (Plus (map (\t -> mathDivideTerm t z) ts2)))
+      _ -> (Div (Plus (map (\t -> mathDivideTerm t z) ts1)) (Plus (map (\t -> mathDivideTerm t z) ts2)))
 
+mathDivideTerm :: TermExpr -> TermExpr -> TermExpr
+mathDivideTerm (Term a xs) (Term b ys) =
+  let (sgn, zs) = f 1 xs ys in
+  (Term (sgn * (div a b)) zs)
+ where
+  f :: Integer -> [(SymbolExpr, Integer)] -> [(SymbolExpr, Integer)] -> (Integer, [(SymbolExpr, Integer)])
+  f sgn xs [] = (sgn, xs)
+  f sgn xs ((y, n):ys) =
+    let (sgns, zs) = unzip (map (\(x, m) -> (g (x, m) (y, n))) xs) in
+    f (sgn * (product sgns)) zs ys
+  g :: (SymbolExpr, Integer) -> (SymbolExpr, Integer) -> (Integer, (SymbolExpr, Integer))
+  g ((Quote x), n) ((Quote y), m) =
+    if x == y
+    then (1, ((Quote x), (n - m)))
+    else if x == (mathNegate y)
+         then if even m then (1, ((Quote x), (n - m))) else (-1, ((Quote x), (n - m)))
+         else (1, ((Quote x), n))
+  g (x, n) (y, m) =
+    if x == y
+    then (1, (x, (n - m)))
+    else (1, (x, n))
+              
 mathRemoveZeroSymbol :: ScalarData -> ScalarData
 mathRemoveZeroSymbol (Div (Plus ts1) (Plus ts2)) =
   let p x = case x of
@@ -574,19 +617,26 @@
 mathSymbolFold (Div (Plus ts1) (Plus ts2)) = Div (Plus (map f ts1)) (Plus (map f ts2))
  where
   f :: TermExpr -> TermExpr
-  f (Term a xs) = Term a (g [] xs)
-  g :: [(SymbolExpr, Integer)] -> [(SymbolExpr, Integer)] -> [(SymbolExpr, Integer)]
+  f (Term a xs) = let (ys, sgns) = unzip $ g [] xs
+                    in Term ((product sgns) * a) ys
+  g :: [((SymbolExpr, Integer),Integer)] -> [(SymbolExpr, Integer)] -> [((SymbolExpr, Integer),Integer)]
   g ret [] = ret
   g ret ((x, n):xs) =
     if (any (p (x, n)) ret)
       then g (map (h (x, n)) ret) xs
-      else g (ret ++ [(x, n)]) xs
-  p :: (SymbolExpr, Integer) -> (SymbolExpr, Integer) -> Bool
-  p (x, _) (y, _) = x == y
-  h :: (SymbolExpr, Integer) -> (SymbolExpr, Integer) -> (SymbolExpr, Integer)
-  h (x, n) (y, m) = if x == y
-                     then (y, m + n)
-                     else (y, m)
+      else g (ret ++ [((x, n), 1)]) xs
+  p :: (SymbolExpr, Integer) -> ((SymbolExpr, Integer), Integer) -> Bool
+  p ((Quote x), _) (((Quote y), _),_) = (x == y) || ((mathNegate x) == y)
+  p (x, _) ((y, _),_) = x == y
+  h :: (SymbolExpr, Integer) -> ((SymbolExpr, Integer), Integer) -> ((SymbolExpr, Integer), Integer)
+  h ((Quote x), n) (((Quote y), m), sgn) = if x == y
+                                      then (((Quote y), m + n), sgn)
+                                      else if x == (mathNegate y)
+                                            then if even n then (((Quote y), m + n), sgn) else (((Quote y), m + n), -1 * sgn)
+                                            else (((Quote y), m), sgn)
+  h (x, n) ((y, m), sgn) = if x == y
+                             then ((y, m + n), sgn)
+                             else ((y, m), sgn)
 
 mathTermFold :: ScalarData -> ScalarData
 mathTermFold (Div (Plus ts1) (Plus ts2)) = Div (Plus (f ts1)) (Plus (f ts2))
@@ -596,25 +646,32 @@
   f' :: [TermExpr] -> [TermExpr] -> [TermExpr]
   f' ret [] = ret
   f' ret ((Term a xs):ts) =
-    if any (\(Term _ ys) -> (p xs ys)) ret
+    if any (\(Term _ ys) -> (fst (p 1 xs ys))) ret
       then f' (map (g (Term a xs)) ret) ts
       else f' (ret ++ [(Term a xs)]) ts
   g :: TermExpr -> TermExpr -> TermExpr
-  g (Term a xs) (Term b ys) = if p xs ys
-                                then (Term (a + b) ys)
-                                else Term b ys
-  p :: [(SymbolExpr, Integer)] -> [(SymbolExpr, Integer)] -> Bool
-  p [] [] = True
-  p [] _ = False
-  p ((x, n):xs) ys =
-    let (b, ys') = q (x, n) [] ys in
+  g (Term a xs) (Term b ys) = let (c, sgn) = p 1 xs ys in
+                                if c
+                                  then (Term ((sgn * a) + b) ys)
+                                  else Term b ys
+  p :: Integer -> [(SymbolExpr, Integer)] -> [(SymbolExpr, Integer)] -> (Bool, Integer)
+  p sgn [] [] = (True, sgn)
+  p sgn [] _ = (False, 0)
+  p sgn ((x, n):xs) ys =
+    let (b, ys', sgn2) = q (x, n) [] ys in
       if b 
-        then p xs ys'
-        else False
-  q :: (SymbolExpr, Integer) -> [(SymbolExpr, Integer)] -> [(SymbolExpr, Integer)] -> (Bool, [(SymbolExpr, Integer)])
-  q _ _ [] = (False, [])
+        then p (sgn * sgn2) xs ys'
+        else (False, 0)
+  q :: (SymbolExpr, Integer) -> [(SymbolExpr, Integer)] -> [(SymbolExpr, Integer)] -> (Bool, [(SymbolExpr, Integer)], Integer)
+  q _ _ [] = (False, [], 1)
+  q ((Quote x), n) ret (((Quote y), m):ys) = if (x == y) && (n == m)
+                                               then (True, (ret ++ ys), 1)
+                                               else if ((mathNegate x) == y) && (n == m)
+                                                      then if even n then (True, (ret ++ ys), 1) else (True, (ret ++ ys), -1)
+                                                      else q ((Quote x), n) (ret ++ [((Quote y), m)]) ys
+  q ((Quote x), n) ret ((y,m):ys) = q ((Quote x), n) (ret ++ [(y, m)]) ys
   q (x, n) ret ((y, m):ys) = if (x == y) && (n == m)
-                               then (True, (ret ++ ys))
+                               then (True, (ret ++ ys), 1)
                                else q (x, n) (ret ++ [(y, m)]) ys
 
 --
@@ -630,6 +687,9 @@
 mathMult :: ScalarData -> ScalarData -> ScalarData
 mathMult (Div m1 n1) (Div m2 n2) = mathNormalize' $ Div (mathMultPoly m1 m2) (mathMultPoly n1 n2)
 
+mathMult' :: ScalarData -> ScalarData -> ScalarData
+mathMult' (Div m1 n1) (Div m2 n2) = Div (mathMultPoly m1 m2) (mathMultPoly n1 n2)
+
 mathMultPoly :: PolyExpr -> PolyExpr -> PolyExpr
 mathMultPoly (Plus []) (Plus _) = Plus []
 mathMultPoly (Plus _) (Plus []) = Plus []
@@ -690,7 +750,7 @@
         let ys = V.take w (V.drop (w * (fromIntegral (i - 1))) xs) in
           fromTensor $ Tensor ns ys (cdr js)
    else throwError $ TensorIndexOutOfBounds i n
-tIntRef' i _ = throwError $ strMsg "More indices than the order of the tensor"
+tIntRef' i _ = throwError $ Default "More indices than the order of the tensor"
  
 tIntRef :: HasTensor a => [Integer] -> (Tensor a) -> EgisonM (Tensor a)
 tIntRef [] (Tensor [] xs _)
@@ -741,7 +801,7 @@
   let yss = split (product ns) xs
   ts <- mapM (\ys -> tref ms (Tensor ns ys (cdr js))) yss
   mapM toTensor ts >>= tConcat s >>= fromTensor
-tref _ t = throwError $ strMsg "More indices than the order of the tensor"
+tref _ t = throwError $ Default "More indices than the order of the tensor"
 
 enumTensorIndices :: [Integer] -> [[Integer]]
 enumTensorIndices [] = [[]]
@@ -764,9 +824,25 @@
   xs' <- mapM (transIndex js is) (enumTensorIndices ns') >>= mapM (flip tIntRef t) >>= mapM fromTensor >>= return . V.fromList
   return $ Tensor ns' xs' is
 
+tTranspose' :: HasTensor a => [EgisonValue] -> (Tensor a) -> EgisonM (Tensor a)
+tTranspose' is t@(Tensor ns xs js) = do
+  is' <- g is js
+  tTranspose is' t
+ where
+  f :: (Index EgisonValue) -> EgisonValue
+  f (Subscript i) = i
+  f (Superscript i) = i
+  f (SupSubscript i) = i
+  g :: [EgisonValue] -> [Index EgisonValue] -> EgisonM [Index EgisonValue]
+  g [] js = return []
+  g (i:is) js = case find (\j -> i == (f j)) js of
+                  Nothing ->  throwError $ InconsistentTensorIndex
+                  (Just j') -> do js' <- g is js
+                                  return $ j':js'
+
 tMap :: HasTensor a => (a -> EgisonM a) -> (Tensor a) -> EgisonM (Tensor a)
 tMap f (Tensor ns xs js) = do
-  xs' <- MP.mapM f (V.toList xs) >>= return . V.fromList
+  xs' <- mapM f (V.toList xs) >>= return . V.fromList
   t <- toTensor (V.head xs')
   case t of
     (Tensor ns1 _ js1) ->
@@ -776,9 +852,9 @@
 
 tMapN :: HasTensor a => ([a] -> EgisonM a) -> [Tensor a] -> EgisonM (Tensor a)
 tMapN f ts@((Tensor ns xs js):_) = do
-  xs' <- MP.mapM (\is -> mapM (tIntRef is) ts >>= mapM fromTensor >>= f) (enumTensorIndices ns)
+  xs' <- mapM (\is -> mapM (tIntRef is) ts >>= mapM fromTensor >>= f) (enumTensorIndices ns)
   return $ Tensor ns (V.fromList xs') js
-tMapN f xs = MP.mapM fromTensor xs >>= f >>= return . Scalar
+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
@@ -788,7 +864,7 @@
   let cns = take (length cjs) (tSize t1')
   rts1 <- mapM (flip tIntRef t1') (enumTensorIndices cns)
   rts2 <- mapM (flip tIntRef t2') (enumTensorIndices cns)
-  rts' <- MP.mapM (\(t1, t2) -> tProduct f t1 t2) (zip rts1 rts2)
+  rts' <- mapM (\(t1, t2) -> tProduct f t1 t2) (zip rts1 rts2)
   let ret = Tensor (cns ++ (tSize (head rts'))) (V.concat (map tToVector rts')) (cjs ++ tIndex (head rts'))
   tTranspose (uniq (tDiagIndex (js1 ++ js2))) ret
  where
@@ -972,7 +1048,7 @@
 
 getScalar :: (Tensor a) -> EgisonM a
 getScalar (Scalar x) = return x
-getScalar _ = throwError $ strMsg "Inconsitent Tensor order"
+getScalar _ = throwError $ Default "Inconsitent Tensor order"
 
 findPairs :: (a -> a -> Bool) -> [a] -> [(Int, Int)]
 findPairs p xs = reverse $ findPairs' 0 p xs
@@ -1109,8 +1185,8 @@
                                                     else showFFloat Nothing y "i"
 
 showTSV :: EgisonValue -> String
-showTSV (Tuple (val:vals)) = foldl (\r x -> r ++ "\t" ++ x) (show val) (map showTSV vals)
-showTSV (Collection vals) = intercalate "\t" (map showTSV (toList vals))
+showTSV (Tuple (val:vals)) = foldl (\r x -> r ++ "\t" ++ x) (show val) (map show vals)
+showTSV (Collection vals) = intercalate "\t" (map show (toList vals))
 showTSV val = show val
 
 instance Eq EgisonValue where
@@ -1351,14 +1427,25 @@
  deriving (Eq)
 type Binding = (String, ObjectRef)
 
+data VarWithIndices = VarWithIndices String [Index String]
+ deriving (Eq)
+
 instance Show Var where
   show (Var x is) = x ++ concat (map show is)
 
+instance Show VarWithIndices where
+  show (VarWithIndices x is) = x ++ concat (map show is)
+
 instance Show (Index ()) where
   show (Superscript ()) = "~"
   show (Subscript ()) = "_"
   show (SupSubscript ()) = "~_"
 
+instance Show (Index String) where
+  show (Superscript s) = "~" ++ s
+  show (Subscript s) = "_" ++ s
+  show (SupSubscript s) = "~_" ++ s
+
 instance Show (Index EgisonExpr) where
   show (Superscript i) = "~" ++ show i
   show (Subscript i) = "_" ++ show i
@@ -1463,10 +1550,6 @@
 
 instance Exception EgisonError
 
-instance Error EgisonError where
-  noMsg = Default "An error has occurred"
-  strMsg = Default
-
 liftError :: (MonadError e m) => Either e a -> m a
 liftError = either throwError return
 
@@ -1475,15 +1558,30 @@
 --
 
 newtype EgisonM a = EgisonM {
-    unEgisonM :: (ErrorT EgisonError (FreshT IO) a)
-  } deriving (Functor, Applicative, Monad, MonadIO, MonadError EgisonError, MonadFresh, MP.MonadParallel)
---  } deriving (Functor, Applicative, Monad, MonadIO, MonadError EgisonError, MonadFresh)
+    unEgisonM :: (ExceptT EgisonError (FreshT IO) a)
+  } deriving (Functor, Applicative, Monad, MonadIO, MonadError EgisonError, MonadFresh)
 
+parallelMapM :: (a -> EgisonM b) -> [a] -> EgisonM [b]
+parallelMapM f [] = return []
+parallelMapM f (x:xs) = do
+    let y = unsafePerformEgison (0,1) $ f x
+    let ys = unsafePerformEgison (0,1) $ parallelMapM f xs
+    y `par` (ys `pseq` return (y:ys))
+
+unsafePerformEgison :: (Int, Int) -> EgisonM a -> a
+unsafePerformEgison (x, y) ma =
+  let ((Right ret), _) = unsafePerformIO $ runFreshT (x, y + 1) $ runEgisonM ma in
+  ret
+--    f' :: (Either EgisonError a) -> (Either EgisonError b) -> EgisonM c
+--    f' (Right x) (Right y) = f x y
+--    f' (Left e) _ = liftError (Left e)
+--    f' _ (Left e) = liftError (Left e)
+
 runEgisonM :: EgisonM a -> FreshT IO (Either EgisonError a)
-runEgisonM = runErrorT . unEgisonM
+runEgisonM = runExceptT . unEgisonM
 
 liftEgisonM :: Fresh (Either EgisonError a) -> EgisonM a
-liftEgisonM m = EgisonM $ ErrorT $ FreshT $ do
+liftEgisonM m = EgisonM $ ExceptT $ FreshT $ do
   s <- get
   (a, s') <- return $ runFresh s m
   put s'
@@ -1492,13 +1590,13 @@
 fromEgisonM :: EgisonM a -> IO (Either EgisonError a)
 fromEgisonM = modifyCounter . runEgisonM
 
-counter :: IORef Int
-counter = unsafePerformIO (newIORef 0)
+counter :: IORef (Int, Int)
+counter = unsafePerformIO (newIORef (0, 0))
 
-readCounter :: IO Int
+readCounter :: IO (Int, Int)
 readCounter = readIORef counter
 
-updateCounter :: Int -> IO ()
+updateCounter :: (Int, Int) -> IO ()
 updateCounter = writeIORef counter
 
 modifyCounter :: FreshT IO a -> IO a
@@ -1508,9 +1606,8 @@
   updateCounter seed
   return result  
 
-newtype FreshT m a = FreshT { unFreshT :: StateT Int m a }
-  deriving (Functor, Applicative, Monad, MonadState Int, MonadTrans, MP.MonadParallel)
---  deriving (Functor, Applicative, Monad, MonadState Int, MonadTrans)
+newtype FreshT m a = FreshT { unFreshT :: StateT (Int, Int) m a }
+  deriving (Functor, Applicative, Monad, MonadState (Int, Int), MonadTrans)
 
 type Fresh = FreshT Identity
 
@@ -1518,8 +1615,8 @@
   fresh :: m String
 
 instance (Applicative m, Monad m) => MonadFresh (FreshT m) where
-  fresh = FreshT $ do counter <- get; modify (+ 1)
-                      return $ "$_" ++ show counter
+  fresh = FreshT $ do (x, y) <- get; modify (\(x,y) -> (x + 1, y))
+                      return $ "$_" ++ (show x) ++ (show y)
 
 instance (MonadError e m) => MonadError e (FreshT m) where
   throwError = lift . throwError
@@ -1532,7 +1629,7 @@
 instance (MonadFresh m) => MonadFresh (StateT s m) where
   fresh = lift $ fresh
 
-instance (MonadFresh m, Error e) => MonadFresh (ErrorT e m) where
+instance (MonadFresh m) => MonadFresh (ExceptT e m) where
   fresh = lift $ fresh
 
 instance (MonadFresh m, Monoid e) => MonadFresh (ReaderT e m) where
@@ -1544,13 +1641,10 @@
 instance MonadIO (FreshT IO) where
   liftIO = lift
 
-instance (MP.MonadParallel m, Error e) => MP.MonadParallel (ErrorT e m)
-instance (MP.MonadParallel m) => MP.MonadParallel (StateT s m)
-
-runFreshT :: Monad m => Int -> FreshT m a -> m (a, Int)
+runFreshT :: Monad m => (Int, Int) -> FreshT m a -> m (a, (Int, Int))
 runFreshT seed = flip (runStateT . unFreshT) seed
 
-runFresh :: Int -> Fresh a -> (a, Int)
+runFresh :: (Int, Int) -> Fresh a -> (a, (Int, Int))
 runFresh seed m = runIdentity $ flip runStateT seed $ unFreshT m
 
 
diff --git a/hs-src/Language/Egison/Util.hs b/hs-src/Language/Egison/Util.hs
--- a/hs-src/Language/Egison/Util.hs
+++ b/hs-src/Language/Egison/Util.hs
@@ -11,7 +11,7 @@
 import Data.List
 import Text.Regex.TDFA
 import System.Console.Haskeline hiding (handle, catch, throwTo)
-import Control.Monad.Error (liftIO)
+import Control.Monad.Except (liftIO)
 
 import Language.Egison.Types
 import Language.Egison.Parser
diff --git a/lib/core/base.egi b/lib/core/base.egi
--- a/lib/core/base.egi
+++ b/lib/core/base.egi
@@ -38,6 +38,8 @@
     (lambda $x
       (foldl 2#(%2 %1) x fs))))
 
+(define $flip (lambda [$fn] (lambda [$x $y] (fn y x))))
+
 (define $ref
   (lambda [$xa $i]
     (array-ref xa i)))
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
@@ -93,10 +93,16 @@
 
 (define $even-and-odd-permutations
   (lambda [$n]
+    (let {[[$es $os] (even-and-odd-permutations' n)]}
+      [(map 1#(lambda [$i] (nth i %1)) es)
+       (map 1#(lambda [$i] (nth i %1)) os)])))
+
+(define $even-and-odd-permutations'
+  (lambda [$n]
     (match n integer
       {[,1 [{{1}} {}]]
        [,2 [{{1 2}} {{2 1}}]]
-       [_ (let* {[[$es $os] (even-and-odd-permutations (- n 1))]
+       [_ (let* {[[$es $os] (even-and-odd-permutations' (- n 1))]
                  [$es' (map 1#{@%1 n} es)]
                  [$os' (map 1#{@%1 n} os)]}
             [{@es'
@@ -121,7 +127,7 @@
     (match (tensor-size m) (list integer)
       {[<cons ,0 <cons ,0 <nil>>> 1]
        [<cons $n <cons ,n <nil>>>
-        (let {[[$es $os] (even-and-odd-permutations n)]}
+        (let {[[$es $os] (even-and-odd-permutations' n)]}
           (- (sum (map (lambda [$e]
                          (product (map2 (lambda [$i $j] m_i_j)
                                         (between 1 n)
diff --git a/lib/math/algebra/root.egi b/lib/math/algebra/root.egi
--- a/lib/math/algebra/root.egi
+++ b/lib/math/algebra/root.egi
@@ -14,13 +14,14 @@
       (match x math-expr
         {[,0 0]
          [?monomial? (rt-monomial n x)]
-         [<div <poly $xs> <poly $ys>>
-          (let {[$xd (reduce gcd xs)]
-                [$yd (reduce gcd ys)]}
-            (let {[$d (rt-monomial n (/ xd yd))]}
-              (*' d
-                 (rt'' n (*' (/' (sum' (map (/' $ xd) xs)) (sum' (map (/' $ yd) ys)))))
-                 )))]})
+;         [<div <poly $xs> <poly $ys>>
+;          (let {[$xd (reduce gcd xs)]
+;                [$yd (reduce gcd ys)]}
+;            (let {[$d (rt-monomial n (/ xd yd))]}
+;              (*' d
+;                 (rt'' n (*' (/' (sum' (map (/' $ xd) xs)) (sum' (map (/' $ yd) ys)))))
+;                 )))]
+         [_ (rt'' n x)]})
       (rt'' n x))))
 
 (define $rt-monomial
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
@@ -19,6 +19,8 @@
 (define $b..' (lambda [%t1 %t2] (contract +' (*' t1 t2))))
 (define $b.. (lambda [%t1 %t2] (contract + (* t1 t2))))
 
+(define $.' (cambda $xs (foldl b..' (car xs) (cdr xs))))
+
 (define $.
   (cambda $xs
     (match xs (list something)
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
@@ -18,14 +18,15 @@
        [(,sqrt $g) (* (/ 1 (* 2 (sqrt g))) (∂/∂ g x))]
        [(,** $g $h) (* f (∂/∂ (* (log g) h) x))]
        [<apply $g $args>
-        (sum (map 2#(* (capply `g|%1 args) (∂/∂ %2 x))
+        (sum (map 2#(* (capply `(add-user-script g %1) args) (∂/∂ %2 x))
                   (zip nats args)))]
        ; quote
        [<quote $g>
         (let {[$g' (∂/∂ g x)]}
           (if (monomial? g')
             g'
-            'g'))]
+            (let {[$d (capply gcd (from-poly g'))]}
+              (*' d '(map-poly (/' $ d) g')))))]
        ; term (constant)
        [,0 0]
        [(* _ ,1) 0]
@@ -57,10 +58,15 @@
       (let {[$hs (generate-tensor 1#h_%1 (tensor-size xs))]}
         (map2 *
               (map 1#(/ 1 (fact %1)) nats0)
-              (map (compose (V.substitute xs as $)
-                            (V.substitute hs (with-symbols {i} (- xs_i as_i)) $))
+              (map (compose 1#(V.substitute xs as %1)
+                            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
   (lambda [%f %xs]
     (multivariate-taylor-expansion f xs (tensor-map 1#0 xs))))
+
+(define $add-user-script
+  (lambda [$f $i]
+    (let {[[$g $is] (decons-user-scripts f)]}
+      (append-user-scripts g (sort {@is i})))))
diff --git a/lib/math/common/arithmetic.egi b/lib/math/common/arithmetic.egi
--- a/lib/math/common/arithmetic.egi
+++ b/lib/math/common/arithmetic.egi
@@ -4,69 +4,30 @@
 ;;;;;
 ;;;;;
 
-(define $to-math-expr (macro [$arg] (math-normalize b.+ 0 (apply to-math-expr' arg))))
+(define $to-math-expr (macro [$arg] (math-normalize1 (apply to-math-expr' arg))))
 
 (define $+' (cambda $xs (foldl b.+' (car xs) (cdr xs))))
 (define $-' (cambda $xs (foldl b.-' (car xs) (cdr xs))))
 (define $*' (cambda $xs (foldl b.*' (car xs) (cdr xs))))
 (define $/' b./')
-(define $.' (cambda $xs (foldl b..' (car xs) (cdr xs))))
 
-(define $b.+
-  (lambda [$x1 $x2]
-    (match [x1 x2] [math-expr math-expr]
-      {[[<div $p1 (& !,1 $p2)> <div $q1 ,p2>]
-        (b./ (b.+ p1 q1) p2)]
-       [[_ _] (reduce-fraction (math-normalize b.+' x1 x2))]})))
-
-(define $b.-
-  (lambda [$x1 $x2]
-    (match [x1 x2] [math-expr math-expr]
-      {[[<div $p1 (& !,1 $p2)> <div $q1 ,p2>]
-        (b./ (b.- p1 q1) p2)]
-       [[_ _] (reduce-fraction (math-normalize b.+' x1 (b.*' -1 x2)))]})))
-
-(define $b.*
-  (lambda [$x1 $x2]
-    (reduce-fraction (math-normalize b.*' x1 x2))))
-
-(define $b./
-  (lambda [$x1 $x2]
-    (reduce-fraction (math-normalize b.*' x1 (b./' 1 x2)))))
-
-(define $+ (cambda $xs (foldl b.+ (car xs) (cdr xs))))
-(define $- (cambda $xs (foldl b.- (car xs) (cdr xs))))
-(define $* (cambda $xs (foldl b.* (car xs) (cdr xs))))
+(define $+ (cambda $xs (math-normalize1 (capply +' xs))))
+(define $- (cambda $xs (math-normalize1 (capply -' xs))))
+(define $* (cambda $xs (math-normalize1 (capply *' xs))))
 (define $/ b./)
 
-
-(define $reduce-fraction
-  (lambda [$mexpr]
-    (match mexpr math-expr
-      {[<div <poly $ts1>
-             <poly $ts2>>
-        (let* {[$d1 (capply gcd ts1)]
-               [$d2 (capply gcd ts2)]
-               [$d (capply gcd {@ts1 @ts2})]}
-          (if (eq? (-' (sum' (map (/' $ d1) ts1)) (sum' (map (/' $ d2) ts2))) 0)
-            (/' (/' d1 d) (/' d2 d))
-            (/' (sum' (map (/' $ d) ts1)) (sum' (map (/' $ d) ts2)))))]})))
+;(define $+ +')
+;(define $- -')
+;(define $* *')
+;(define $/ b./')
 
-;            (let* {[$e1 (expand-all' (sum' (map (/' $ d) ts1)))]
-;                   [$qexprs (match-all (/' d2 d) math-expr
-;                              [<term _ <cons <quote $qexpr> _>> qexpr])]
-;                   [$dps (filter 1#(eq? 0 (2#%2 (P./ e1 %1 (car (find-symbols-from-poly %1)))))
-;                                 qexprs)]}
-;              (if (eq? dps {})
-;                (/' (sum' (map (/' $ d) ts1)) (sum' (map (/' $ d) ts2)))
-;                (let {[$dp (car dps)]}
-;                  (reduce-fraction (/' (2#%1 (P./ e1 dp (car (find-symbols-from-poly dp))))
-;                                       (sum' (map (/' $ (*' d 'dp)) ts2))))
-;                  )))))]})))
+(define $reduce-fraction id)
 
 (define $sum
   (lambda [$xs]
-    (foldl + 0 xs)))
+    (if (empty? xs)
+      0
+      (capply + xs))))
 
 (define $sum'
   (lambda [$xs]
@@ -74,7 +35,9 @@
 
 (define $product
   (lambda [$xs]
-    (foldl * 1 xs)))
+    (if (empty? xs)
+      1
+      (capply * xs))))
 
 (define $product'
   (lambda [$xs]
@@ -82,7 +45,7 @@
 
 (define $power
   (lambda [$x $n]
-    (foldl * 1 (take n (repeat1 x)))))
+    (math-normalize1 (power' x n))))
 
 (define $power'
   (lambda [$x $n]
@@ -143,4 +106,3 @@
            [[$zs $rs] (L./ as bs)]}
       [(sum' (map2 2#(*' %1 (**' x %2)) (reverse zs) nats0))
        (sum' (map2 2#(*' %1 (**' x %2)) (reverse rs) nats0))])))
-        
diff --git a/lib/math/expression.egi b/lib/math/expression.egi
--- a/lib/math/expression.egi
+++ b/lib/math/expression.egi
@@ -177,6 +177,20 @@
       {[<div $p1 $p2>
         (/' (fn p1) (fn p2))]})))
 
+(define $from-poly
+  (lambda [$mexpr]
+    (match mexpr math-expr
+      {[<div <poly $ts1> $q>
+        (map (lambda [$t1] (/' t1 q))
+             ts1)]})))
+
+(define $map-poly
+  (lambda [$fn $mexpr]
+    (match mexpr math-expr
+      {[<div <poly $ts1> $q>
+        (foldl +' 0 (map (lambda [$t1] (fn (/' t1 q)))
+                         ts1))]})))
+
 (define $map-terms
   (lambda [$fn $mexpr]
     (match mexpr math-expr
@@ -227,6 +241,22 @@
                                               (if (eq? f g)
                                                 #t
                                                 (any id (map (contain-function? f $) args)))]
+                                             [_ #f]})
+                                        xs))]}))
+                     {@ts1 @ts2})]}))))
+
+(define $contain-function-with-order?
+  (lambda [$f $n $mexpr]
+    (any id (match mexpr math-expr
+              {[<div <poly $ts1> <poly $ts2>>
+                (map (lambda [$term]
+                       (match term term-expr
+                         {[<term _ $xs>
+                           (any id (map 2#(match %1 symbol-expr
+                                            {[<apply $g $args>
+                                              (if (and (eq? f g) (gte? %2 n))
+                                                #t
+                                                (any id (map (contain-function-with-order? f n $) args)))]
                                              [_ #f]})
                                         xs))]}))
                      {@ts1 @ts2})]}))))
diff --git a/lib/math/normalize.egi b/lib/math/normalize.egi
--- a/lib/math/normalize.egi
+++ b/lib/math/normalize.egi
@@ -4,35 +4,25 @@
 ;;;;;
 ;;;;;
 
-(define $math-normalize
-  (lambda [$fn $x1 $x2]
-    (if (or (integer? x1) (integer? x2))
-      (fn x1 x2)
-
-;      (fn x1 x2))))
-      
-;      (let* {[$tmp (fn x1 x2)]
-;             [$ret ((capply compose (map 2#%1 (filter 2#(%2 x1 x2) rewrite-rules))) tmp)]}
-;        (if (eq? tmp ret)
-;          ret
-;          (debug ((capply compose (map 2#%1 (filter 2#(%2 x1 x2) rewrite-rules))) (debug ((debug fn) x1 x2)))))))))
-
-      (let* {[$tmp (fn  x1 x2)]
-             [$ret ((capply compose (map 2#%1 (filter 2#(%2 x1 x2) rewrite-rules))) tmp)]}
-        ret))))
+(define $math-normalize1
+  (lambda [$x]
+    (if (integer? x)
+      x
+      (let {[$ret ((capply compose (map 2#%1 (filter 2#(%2 x) rewrite-rules1))) x)]} ret))))
+;      (let {[$ret ((capply compose (map 2#%1 (filter 2#(%2 fn x) rewrite-rules1))) (debug x))]} (debug ret)))))
 
-(define $rewrite-rules
+(define $rewrite-rules1
   {
-   [id 2##t]
-   [rewrite-rule-for-i 2#(and (contain-symbol? i %1) (contain-symbol? i %2))]
-   [rewrite-rule-for-w-term 2#(and (contain-symbol? w %1) (contain-symbol? w %2))]
-   [rewrite-rule-for-rtu-term 2#(and (contain-function? rtu %1) (contain-function? rtu %2))]
-   [rewrite-rule-for-w-poly 2#(and (contain-symbol? w %1) (contain-symbol? w %2))]
-   [rewrite-rule-for-rtu-poly 2#(and (contain-function? rtu %1) (contain-function? rtu %2))]
-   [rewrite-rule-for-sqrt 2#(and (contain-function? sqrt %1) (contain-function? sqrt %2))]
-   [rewrite-rule-for-rt 2#(and (contain-function? rt %1) (contain-function? rt %2))]
-   [rewrite-rule-for-cos-and-sin 2#(or (or (contain-function? cos %1) (contain-function? sin %1)) (or (contain-function? cos %2) (contain-function? sin %2)))]
-   [rewrite-rule-for-∂/∂ 2#(and (contain-function-with-index? %1) (contain-function-with-index? %2))]
+   [id 1##t]
+   [rewrite-rule-for-i 1#(contain-symbol? i %1)]
+   [rewrite-rule-for-w-term 1#(contain-symbol? w %1)]
+   [rewrite-rule-for-rtu-term 1#(contain-function? rtu %1)]
+   [rewrite-rule-for-w-poly 1#(contain-symbol? w %1)]
+   [rewrite-rule-for-rtu-poly 1#(contain-function? rtu %1)]
+   [rewrite-rule-for-sqrt 1#(contain-function? sqrt %1)]
+   [rewrite-rule-for-rt 1#(contain-function? rt %1)]
+;   [rewrite-rule-for-cos-and-sin 1#(or (contain-function-with-order? cos 2 %1) (contain-function-with-order? sin 2 %1))]
+   [rewrite-rule-for-cos-to-sin 1#(contain-function-with-order? cos 2 %1)]
    })
 
 ;;
@@ -105,9 +95,10 @@
 (define $rewrite-rule-for-rtu-poly'
   (lambda [$poly]
     (match poly math-expr
-      {[(+ (* $a (,rtu $n)^,(- n 1) $mr)
+      {
+       [(+ (* $a (,rtu $n)^,1 $mr)
            (loop $i [2 ,(- n 1)]
-             (+ (* ,a ,(rtu n)^,(- n i) ,mr) ...)
+             (+ (* ,a ,(rtu n)^,i ,mr) ...)
              $pr))
         (rewrite-rule-for-rtu-poly' (+' pr (*' -1 a mr)))]
        [_ poly]})))
@@ -210,6 +201,15 @@
 ;           $pr)
 ;        (rewrite-rule-for-cos-and-sin-poly (+' pr (*' (* 2 a) (sin (/ θ 2))^2 mr)))]
        [_ poly]})))
+
+(define $rewrite-rule-for-cos-to-sin 1#(map-terms rewrite-rule-for-cos-to-sin-term' %1))
+
+(define $rewrite-rule-for-cos-to-sin-term'
+  (lambda [$term]
+    (match term math-expr
+      {[(* $a (,cos $θ)^,2 $mr)
+        (*' a (-' 1 (sin θ)^2) (rewrite-rule-for-cos-to-sin-term' mr))]
+       [_ term]})))
 
 ;;
 ;; d