egison 3.6.3 → 3.6.4

raw patch · 13 files changed

+339/−100 lines, 13 filesdep +monad-parallelPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: monad-parallel

API changes (from Hackage documentation)

Files

• egison.cabal +5/−5
• hs-src/Language/Egison/Core.hs +41/−32
• hs-src/Language/Egison/Desugar.hs +1/−1
• hs-src/Language/Egison/Parser.hs +4/−2
• hs-src/Language/Egison/Types.hs +52/−29
• lib/core/collection.egi +12/−0
• lib/math/algebra/equations.egi +2/−2
• lib/math/algebra/matrix.egi +54/−6
• lib/math/algebra/root.egi +3/−3
• lib/math/analysis/derivative.egi +5/−1
• lib/math/common/arithmetic.egi +29/−4
• lib/math/expression.egi +68/−5
• lib/math/normalize.egi +63/−10

@@ -1,5 +1,5 @@ Name:                egison-Version:             3.6.3+Version:             3.6.4 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+  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   Hs-Source-Dirs:  hs-src   Exposed-Modules:                    Language.Egison@@ -79,7 +79,7 @@                    Language.Egison.Primitives                    Language.Egison.Util   Other-modules:   Paths_egison-  ghc-prof-options:   --enable-profiling -auto-all+  ghc-options:  -O3  Test-Suite test   Type:           exitcode-stdio-1.0@@ -100,5 +100,5 @@   Build-depends:       egison, base >= 4.0 && < 5, array, containers, unordered-containers, haskeline, transformers, mtl, parsec >= 3.0, directory, ghc, ghc-paths, filepath, text, regex-tdfa, process, vector   Hs-Source-Dirs:      hs-src/Interpreter   Other-modules:       Paths_egison-  ghc-options:  -rtsopts -threaded-  ghc-prof-options:   --enable-profiling -auto-all+  ghc-options:  -O3 -threaded -eventlog -rtsopts+  ghc-prof-options:   -auto-all

@@ -46,6 +46,7 @@ import Control.Monad.Error 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@@ -204,11 +205,11 @@   return . Intermediate . IArray $ Array.listArray (1, toInteger (length exprs)) refs'  evalExpr env (VectorExpr exprs) = do-  whnfs <- mapM (evalExpr env) exprs+  whnfs <- MP.mapM (evalExpr env) exprs   case whnfs of-    [whnf] -> return $ whnf     ((Intermediate (ITensor (Tensor _ _ _))):_) -> do-      mapM toTensor whnfs >>= tConcat' >>= fromTensor+      ret <- mapM toTensor whnfs >>= tConcat' >>= fromTensor+      return ret     _ -> do       fromTensor (Tensor [fromIntegral (length whnfs)] (V.fromList whnfs) []) @@ -216,7 +217,7 @@   nsWhnf <- evalExpr env nsExpr   ns <- ((fromCollection nsWhnf >>= fromMList >>= mapM evalRef >>= mapM fromWHNF) :: EgisonM [Integer])   xsWhnf <- evalExpr env xsExpr-  xs <- fromCollection xsWhnf >>= fromMList >>= mapM evalRef+  xs <- fromCollection xsWhnf >>= fromMList >>= MP.mapM evalRef   supWhnf <- evalExpr env supExpr   sup <- fromCollection supWhnf >>= fromMList >>= mapM evalRefDeep -- >>= mapM extractScalar'   subWhnf <- evalExpr env subExpr@@ -276,31 +277,32 @@                       Subscript n -> evalExprDeep env n >>= return . Subscript                       SupSubscript n -> evalExprDeep env n >>= return . SupSubscript               ) indices-  case tensor of-    (Value (ScalarData (Div (Plus [(Term 1 [(Symbol id name [], 1)])]) (Plus [(Term 1 [])])))) -> do-      js2 <- mapM (\i -> case i of-                           Superscript n -> evalExprDeep env n >>= extractScalar >>= return . Superscript-                           Subscript n -> evalExprDeep env n >>= extractScalar >>= return . Subscript-                           SupSubscript n -> evalExprDeep env n >>= extractScalar >>= return . SupSubscript-                  ) indices-      return $ Value (ScalarData (Div (Plus [(Term 1 [(Symbol id name js2, 1)])]) (Plus [(Term 1 [])])))-    (Value (ScalarData _)) -> do-      return $ tensor-    (Value (TensorData (Tensor ns xs _))) -> do-      tref js (Tensor ns xs js) >>= toTensor >>= tContract' >>= fromTensor >>= return . Value-    (Intermediate (ITensor (Tensor ns xs _))) -> do-      tref js (Tensor ns xs js) >>= toTensor >>= tContract' >>= fromTensor-    _ -> do-      js2 <- mapM (\i -> case i of-                           Superscript n -> evalExprDeep env n >>= extractScalar >>= return . Superscript-                           Subscript n -> evalExprDeep env n >>= extractScalar >>= return . Subscript-                           SupSubscript n -> evalExprDeep env n >>= extractScalar >>= return . SupSubscript-                  ) indices-      refArray tensor (map (\j -> case j of-                                    Superscript k -> ScalarData k-                                    Subscript k -> ScalarData k-                                    SupSubscript k -> ScalarData k-                            ) js2)+  ret <- case tensor of+      (Value (ScalarData (Div (Plus [(Term 1 [(Symbol id name [], 1)])]) (Plus [(Term 1 [])])))) -> do+        js2 <- mapM (\i -> case i of+                             Superscript n -> evalExprDeep env n >>= extractScalar >>= return . Superscript+                             Subscript n -> evalExprDeep env n >>= extractScalar >>= return . Subscript+                             SupSubscript n -> evalExprDeep env n >>= extractScalar >>= return . SupSubscript+                    ) indices+        return $ Value (ScalarData (Div (Plus [(Term 1 [(Symbol id name js2, 1)])]) (Plus [(Term 1 [])])))+      (Value (ScalarData _)) -> do+        return $ tensor+      (Value (TensorData (Tensor ns xs _))) -> do+        tref js (Tensor ns xs js) >>= toTensor >>= tContract' >>= fromTensor >>= return . Value+      (Intermediate (ITensor (Tensor ns xs _))) -> do+        tref js (Tensor ns xs js) >>= toTensor >>= tContract' >>= fromTensor+      _ -> do+        js2 <- mapM (\i -> case i of+                             Superscript n -> evalExprDeep env n >>= extractScalar >>= return . Superscript+                             Subscript n -> evalExprDeep env n >>= extractScalar >>= return . Subscript+                             SupSubscript n -> evalExprDeep env n >>= extractScalar >>= return . SupSubscript+                    ) indices+        refArray tensor (map (\j -> case j of+                                      Superscript k -> ScalarData k+                                      Subscript k -> ScalarData k+                                      SupSubscript k -> ScalarData k+                              ) js2)+  return ret  where   f :: Index a -> Index ()   f (Superscript _) = Superscript ()@@ -546,7 +548,7 @@   size'' <- collectionToList size'   ns <- (mapM fromEgison size'') :: EgisonM [Integer]   fn <- evalExpr env fnExpr-  xs <-  mapM (\ms -> applyFunc env fn (Value (makeTuple ms))) (map (\ms -> map toEgison ms) (enumTensorIndices ns))+  xs <-  MP.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) [])@@ -669,9 +671,16 @@ evalWHNF (Intermediate (ITuple [ref])) = evalRefDeep ref evalWHNF (Intermediate (ITuple refs)) = Tuple <$> mapM evalRefDeep refs evalWHNF (Intermediate (ITensor (Tensor ns whnfs js))) = do-  vals <- mapM evalWHNF whnfs-  return $ TensorData $ Tensor ns vals js+  vals <- MP.mapM evalWHNF (V.toList whnfs)+  return $ TensorData $ Tensor ns (V.fromList vals) js+--  vals <- mapM evalWHNF whnfs+--  return $ TensorData $ Tensor ns vals js evalWHNF coll = Collection <$> (fromCollection coll >>= fromMList >>= mapM evalRefDeep . Sq.fromList)++fib :: Integer -> Integer+fib 0 = 1+fib 1 = 1+fib n = (fib (n - 1)) + (fib (n - 2))  applyFunc :: Env -> WHNFData -> WHNFData -> EgisonM WHNFData applyFunc _ (Value (PartialFunc env n body)) arg = do

@@ -143,7 +143,7 @@ desugar (PowerExpr expr1 expr2) = do   expr1' <- desugar expr1   expr2' <- desugar expr2-  return $ ApplyExpr (VarExpr "**'") (TupleExpr [expr1', expr2'])+  return $ ApplyExpr (VarExpr "**") (TupleExpr [expr1', expr2'])  desugar (ArrayBoundsExpr expr) = do   expr' <- desugar expr

@@ -155,13 +155,15 @@ exprs = endBy expr whiteSpace  expr :: Parser EgisonExpr-expr = P.lexeme lexer (do expr0 <- expr'+expr = P.lexeme lexer (do expr0 <- expr' <|> quoteExpr'                           expr1 <- option expr0 $ IndexedExpr 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)                           option expr2 $ PowerExpr expr1 <$> (try $ char '^' >> expr'))                           +quoteExpr' :: Parser EgisonExpr+quoteExpr' = char '\'' >> QuoteExpr <$> expr'  expr' :: Parser EgisonExpr expr' = (try partialExpr@@ -175,7 +177,7 @@              <|> try tupleExpr              <|> try hashExpr              <|> collectionExpr-             <|> quoteExpr+--             <|> quoteExpr              <|> quoteFunctionExpr              <|> parens (ifExpr                          <|> lambdaExpr

@@ -153,6 +153,7 @@ 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@@ -426,10 +427,6 @@   tensorElems (TensorData (Tensor _ xs _)) = xs   tensorSize (TensorData (Tensor ns _ _)) = ns   tensorIndices (TensorData (Tensor _ _ js)) = js-  fromTensor (Tensor [] xs []) =-    if V.length xs == 1-      then return $ V.head xs-      else throwError $ InconsistentTensorIndex   fromTensor t@(Tensor _ _ _) = return $ TensorData t   fromTensor (Scalar x) = return x   toTensor (TensorData t) = return t@@ -440,10 +437,6 @@   tensorElems (Intermediate (ITensor (Tensor _ xs _))) = xs   tensorSize (Intermediate (ITensor (Tensor ns _ _))) = ns   tensorIndices (Intermediate (ITensor (Tensor _ _ js))) = js-  fromTensor (Tensor [] xs []) =-    if V.length xs == 1-      then return $ V.head xs-      else throwError $ InconsistentTensorIndex   fromTensor t@(Tensor _ _ _) = return $ Intermediate $ ITensor t   fromTensor (Scalar x) = return x   toTensor (Intermediate (ITensor t)) = return t@@ -714,18 +707,36 @@ tref ((Subscript (ScalarData (Div (Plus [(Term m [])]) (Plus [(Term 1 [])])))):ms) t = tIntRef' m t >>= toTensor >>= tref ms tref ((Superscript (ScalarData (Div (Plus [(Term m [])]) (Plus [(Term 1 [])])))):ms) t = tIntRef' m t >>= toTensor >>= tref ms tref ((SupSubscript (ScalarData (Div (Plus [(Term m [])]) (Plus [(Term 1 [])])))):ms) t = tIntRef' m t >>= toTensor >>= tref ms-tref ((Subscript (Tuple [ScalarData (Div (Plus [(Term m [])]) (Plus [(Term 1 [])])), ScalarData (Div (Plus [(Term n [])]) (Plus [(Term 1 [])]))])):ms) t = do-  ts <- mapM (\i -> tIntRef' i t >>= toTensor >>= tref ms >>= toTensor) [m..n]-  symId <- fresh-  tConcat (Subscript (symbolScalarData "" (":::" ++ symId))) ts >>= fromTensor-tref ((Superscript (Tuple [ScalarData (Div (Plus [(Term m [])]) (Plus [(Term 1 [])])), ScalarData (Div (Plus [(Term n [])]) (Plus [(Term 1 [])]))])):ms) t = do-  ts <- mapM (\i -> tIntRef' i t >>= toTensor >>= tref ms >>= toTensor) [m..n]-  symId <- fresh-  tConcat (Superscript (symbolScalarData "" (":::" ++ symId))) ts >>= fromTensor-tref ((SupSubscript (Tuple [ScalarData (Div (Plus [(Term m [])]) (Plus [(Term 1 [])])), ScalarData (Div (Plus [(Term n [])]) (Plus [(Term 1 [])]))])):ms) t = do-  ts <- mapM (\i -> tIntRef' i t >>= toTensor >>= tref ms >>= toTensor) [m..n]-  symId <- fresh-  tConcat (SupSubscript (symbolScalarData "" (":::" ++ symId))) ts >>= fromTensor+tref ((Subscript (Tuple [mVal, nVal])):ms) t@(Tensor is _ _) = do+  m <- fromEgison mVal+  n <- fromEgison nVal+  if m > n+    then do+      fromTensor (Tensor (take (length is) (repeat 0)) V.empty [])+    else do+      ts <- mapM (\i -> tIntRef' i t >>= toTensor >>= tref ms >>= toTensor) [m..n]+      symId <- fresh+      tConcat (Subscript (symbolScalarData "" (":::" ++ symId))) ts >>= fromTensor+tref ((Superscript (Tuple [mVal, nVal])):ms) t@(Tensor is _ _) = do+  m <- fromEgison mVal+  n <- fromEgison nVal+  if m > n+    then do+      fromTensor (Tensor (take (length is) (repeat 0)) V.empty [])+    else do+      ts <- mapM (\i -> tIntRef' i t >>= toTensor >>= tref ms >>= toTensor) [m..n]+      symId <- fresh+      tConcat (Superscript (symbolScalarData "" (":::" ++ symId))) ts >>= fromTensor+tref ((SupSubscript (Tuple [mVal, nVal])):ms) t@(Tensor is _ _) = do+  m <- fromEgison mVal+  n <- fromEgison nVal+  if m > n+    then do+      fromTensor (Tensor (take (length is) (repeat 0)) V.empty [])+    else do+      ts <- mapM (\i -> tIntRef' i t >>= toTensor >>= tref ms >>= toTensor) [m..n]+      symId <- fresh+      tConcat (SupSubscript (symbolScalarData "" (":::" ++ symId))) ts >>= fromTensor tref (s:ms) (Tensor (n:ns) xs js) = do   let yss = split (product ns) xs   ts <- mapM (\ys -> tref ms (Tensor ns ys (cdr js))) yss@@ -755,7 +766,7 @@  tMap :: HasTensor a => (a -> EgisonM a) -> (Tensor a) -> EgisonM (Tensor a) tMap f (Tensor ns xs js) = do-  xs' <- V.mapM f xs+  xs' <- MP.mapM f (V.toList xs) >>= return . V.fromList   t <- toTensor (V.head xs')   case t of     (Tensor ns1 _ js1) ->@@ -765,9 +776,9 @@  tMapN :: HasTensor a => ([a] -> EgisonM a) -> [Tensor a] -> EgisonM (Tensor a) tMapN f ts@((Tensor ns xs js):_) = do-  xs' <- mapM (\is -> mapM (tIntRef is) ts >>= mapM fromTensor >>= f) (enumTensorIndices ns)+  xs' <- MP.mapM (\is -> mapM (tIntRef is) ts >>= mapM fromTensor >>= f) (enumTensorIndices ns)   return $ Tensor ns (V.fromList xs') js-tMapN f xs = mapM fromTensor xs >>= f >>= return . Scalar+tMapN f xs = MP.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@@ -777,7 +788,7 @@   let cns = take (length cjs) (tSize t1')   rts1 <- mapM (flip tIntRef t1') (enumTensorIndices cns)   rts2 <- mapM (flip tIntRef t2') (enumTensorIndices cns)-  rts' <- mapM (\(t1, t2) -> tProduct f t1 t2) (zip rts1 rts2)+  rts' <- MP.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@@ -929,13 +940,18 @@                  hs:(split w ts)  tConcat :: HasTensor a => Index EgisonValue -> [Tensor a] -> EgisonM (Tensor a)+tConcat s ((Tensor ns@(0:_) _ js):_) = do+  return $ Tensor (0:ns) V.empty (s:js) tConcat s ts@((Tensor ns _ js):_) = return $ Tensor ((fromIntegral (length ts)):ns) (V.concat (map tToVector ts)) (s:js) tConcat s ts = do   ts' <- mapM getScalar ts   return $ Tensor [fromIntegral (length ts)] (V.fromList ts') [s]  tConcat' :: HasTensor a => [Tensor a] -> EgisonM (Tensor a)-tConcat' ts@((Tensor ns _ _):_) = return $ Tensor ((fromIntegral (length ts)):ns) (V.concat (map tToVector ts)) []+tConcat' ((Tensor ns@(0:_) _ _):_) = do+  return $ Tensor (0:ns) V.empty []+tConcat' ts@((Tensor ns v _):_) = do+  return $ Tensor ((fromIntegral (length ts)):ns) (V.concat (map tToVector ts)) [] tConcat' ts = do   ts' <- mapM getScalar ts   return $ Tensor [fromIntegral (length ts)] (V.fromList ts') []@@ -1010,6 +1026,7 @@   show (ScalarData mExpr) = show mExpr --  show (TensorData (Scalar x)) = "invalid scalar:" ++ show x   show (TensorData (Tensor [_] xs js)) = "[| " ++ unwords (map show (V.toList xs)) ++ " |]" ++ concat (map show js)+  show (TensorData (Tensor [0, 0] _ js)) = "[| [|  |] |]" ++ concat (map show js)   show (TensorData (Tensor [i, j] xs js)) = "[| " ++ f (fromIntegral j) (V.toList xs) ++ "|]" ++ concat (map show js)    where     f j [] = ""@@ -1258,7 +1275,8 @@   show (Intermediate (IIntHash _)) = "{|...|}"    show (Intermediate (ICharHash _)) = "{|...|}"    show (Intermediate (IStrHash _)) = "{|...|}" -  show (Intermediate (ITensor _)) = "[|...|]" +--  show (Intermediate (ITensor _)) = "[|...|]" +  show (Intermediate (ITensor (Tensor ns xs _))) = "[|" ++ show (length ns) ++ show (V.length xs) ++ "|]"   instance Show Object where   show (Thunk _) = "#<thunk>"@@ -1457,8 +1475,9 @@ --  newtype EgisonM a = EgisonM {-    unEgisonM :: ErrorT EgisonError (FreshT IO) a-  } deriving (Functor, Applicative, Monad, MonadIO, MonadError EgisonError, MonadFresh)+    unEgisonM :: (ErrorT EgisonError (FreshT IO) a)+  } deriving (Functor, Applicative, Monad, MonadIO, MonadError EgisonError, MonadFresh, MP.MonadParallel)+--  } deriving (Functor, Applicative, Monad, MonadIO, MonadError EgisonError, MonadFresh)  runEgisonM :: EgisonM a -> FreshT IO (Either EgisonError a) runEgisonM = runErrorT . unEgisonM@@ -1490,7 +1509,8 @@   return result    newtype FreshT m a = FreshT { unFreshT :: StateT Int m a }-  deriving (Functor, Applicative, Monad, MonadState Int, MonadTrans)+  deriving (Functor, Applicative, Monad, MonadState Int, MonadTrans, MP.MonadParallel)+--  deriving (Functor, Applicative, Monad, MonadState Int, MonadTrans)  type Fresh = FreshT Identity @@ -1523,6 +1543,9 @@  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 seed = flip (runStateT . unFreshT) seed

@@ -359,6 +359,18 @@         {s @(between (+ s 1) e)}         {})))) +(define $L./+  (lambda [$xs $ys]+    (if (lt? (length xs) (length ys))+      [{} xs]+      (match [ys xs] [(list math-expr) (list math-expr)]+        {+         [[<cons $y $yrs> <cons $x $xrs>]+          (let {[[$zs $rs] (L./ {@(map2 - (take (length yrs) xrs) (map (* $ (/ x y)) yrs))+                                 @(drop (length yrs) xrs)} ys)]}+            [{(/ x y) @zs} rs])]+         }))))+ ;;; ;;; Multiset ;;;

@@ -26,7 +26,7 @@  (define $q-f   (lambda [$f $x]-    (match (coefficients x f) (list math-expr)+    (match (coefficients f x) (list math-expr)       {[<cons $a_0 <cons $a_1 <cons $a_2 <nil>>>>         (q-f' a_2 a_1 a_0)]}))) @@ -42,7 +42,7 @@  (define $c-f   (lambda [$f $x]-    (match (coefficients x f) (list math-expr)+    (match (coefficients f x) (list math-expr)       {[<cons $a_0 <cons $a_1 <cons $a_2 <cons $a_3 <nil>>>>>         (c-f' a_3 a_2 a_1 a_0)]}))) 

@@ -31,10 +31,14 @@  (define $M.inverse   (lambda [%m]-    (match (tensor-size m) (list integer)-      {[<cons ,2 <cons ,2 <nil>>>-        (T.map (/ $ (M.det m)) (tensor {2 2} {m_2_2 (* -1 m_1_2) (* -1 m_2_1) m_1_1}))]-       [_ undefined]})))+    (let {[$d (M.det m)]}+      (generate-tensor+        2#(match m matrix+            {[<cons ,%2 ,%1 _ $A $B $C $D>+              (if (even? (+ %1 %2))+                (/ (M.det (M.join A B C D)) d)+                (* -1 (/ (M.det (M.join A B C D)) d)))]})+        (tensor-size m)))))  (define $trace (lambda [%t] (with-symbols {i} (contract + t~i_i)))) @@ -45,12 +49,44 @@                {[<cons $m <cons $n _>>                  {[tgt_1_1 tgt_1_[2 n] tgt_[2 m]_1 tgt_[2 m]_[2 n]]}]                 [_ {}]})]}]+     [<cons ,$i ,$j $ $ $ $ $> [math-expr matrix matrix matrix matrix]+      {[$tgt+        (let* {[$ns (tensor-size tgt)]+               [$m (nth 1 ns)]+               [$n (nth 2 ns)]}+          {[tgt_i_j+            tgt_[1 (- i 1)]_[1 (- j 1)]+            tgt_[1 (- i 1)]_[(+ j 1) n]+            tgt_[(+ i 1) m]_[1 (- j 1)]+            tgt_[(+ i 1) m]_[(+ j 1) n]+            ]})]}]      [,$val []       {[$tgt (if (eq? val tgt) {[]} {})]}]      [$ [something]       {[$tgt {tgt}]}]      })) +(define $M.join+  (lambda [%A %B %C %D]+    (let* {[$as (tensor-size A)]+           [$a1 (nth 1 as)] [$a2 (nth 2 as)]+           [$bs (tensor-size B)]+           [$b1 (nth 1 bs)] [$b2 (nth 2 bs)]+           [$cs (tensor-size C)]+           [$c1 (nth 1 cs)] [$c2 (nth 2 cs)]+           [$ds (tensor-size D)]+           [$d1 (nth 1 ds)] [$d2 (nth 2 ds)]+           [$m1 (max a1 b1)] [$m2 (max a2 c2)]+           [$n1 (max c1 d1)] [$n2 (max b2 d2)]+           }+      (generate-tensor+        2#(match [%1 %2] [integer integer]+            {[[?(lte? $ a1) ?(lte? $ a2)] A_%1_%2]+             [[?(lte? $ m1) _] B_%1_(- %2 a2)]+             [[_ ?(lte? $ m2)] C_(- %1 a1)_%2]+             [[_ _] D_(- %1 m1)_(- %2 m2)]})+        {(+ m1 n1) (+ m2 n2)}))))+ ;; ;; Determinant ;;@@ -58,7 +94,8 @@ (define $even-and-odd-permutations   (lambda [$n]     (match n integer-      {[,2 [{{1 2}} {{2 1}}]]+      {[,1 [{{1}} {}]]+       [,2 [{{1 2}} {{2 1}}]]        [_ (let* {[[$es $os] (even-and-odd-permutations (- n 1))]                  [$es' (map 1#{@%1 n} es)]                  [$os' (map 1#{@%1 n} os)]}@@ -82,7 +119,8 @@ (define $M.determinant   (lambda [%m]     (match (tensor-size m) (list integer)-      {[<cons $n <cons ,n <nil>>>+      {[<cons ,0 <cons ,0 <nil>>> 1]+       [<cons $n <cons ,n <nil>>>         (let {[[$es $os] (even-and-odd-permutations n)]}           (- (sum (map (lambda [$e]                          (product (map2 (lambda [$i $j] m_i_j)@@ -143,3 +181,13 @@                              [m_3_1 x_3_1 b_3_1] [m_3_2 x_3_2 b_3_2] [m_3_3 x_3_3 c_3_3]})]}           [(substitute ret L) (substitute ret U)])]        [_ undefined]})))++;;+;; Utility+;;++(define $generate-matrix-from-quadratic-expr+  (lambda [$f $xs]+    (generate-tensor+      2#(coefficient2 f (nth %1 xs) (nth %2 xs))+      {(length xs) (length xs)})))

@@ -17,9 +17,9 @@          [<div <poly $xs> <poly $ys>>           (let {[$xd (reduce gcd xs)]                 [$yd (reduce gcd ys)]}-            (let {[[$a $r] (from-monomial (rt-monomial n (/ xd yd)))]}-              (*' a-                 (rt'' n (*' (** r n) (/' (sum' (map (/' $ xd) xs)) (sum' (map (/' $ yd) ys)))))+            (let {[$d (rt-monomial n (/ xd yd))]}+              (*' d+                 (rt'' n (*' (/' (sum' (map (/' $ xd) xs)) (sum' (map (/' $ yd) ys)))))                  )))]})       (rt'' n x)))) 

@@ -21,7 +21,11 @@         (sum (map 2#(* (capply `g|%1 args) (∂/∂ %2 x))                   (zip nats args)))]        ; quote-       [<quote $g> (∂/∂ g x)]+       [<quote $g>+        (let {[$g' (∂/∂ g x)]}+          (if (monomial? g')+            g'+            'g'))]        ; term (constant)        [,0 0]        [(* _ ,1) 0]

@@ -24,7 +24,7 @@     (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))]})))+       [[_ _] (reduce-fraction (math-normalize b.+' x1 (b.*' -1 x2)))]})))  (define $b.*   (lambda [$x1 $x2]@@ -32,7 +32,7 @@  (define $b./   (lambda [$x1 $x2]-    (reduce-fraction (math-normalize b./' 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))))@@ -45,9 +45,25 @@     (match mexpr math-expr       {[<div <poly $ts1>              <poly $ts2>>-        (let {[$d (capply gcd {@ts2 @ts1})]}-          (/' (sum' (map (/' $ d) ts1)) (sum' (map (/' $ d) 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)))))]}))) +;            (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 $sum   (lambda [$xs]     (foldl + 0 xs)))@@ -119,3 +135,12 @@        [[,0 _] y]        [[_ ?(gte? $ x)] (b.gcd' (modulo y x) x)]        [[_ _] (b.gcd' y x)]})))++(define $P./+  (lambda [$fx $gx $x]+    (let* {[$as (reverse (coefficients fx x))]+           [$bs (reverse (coefficients gx x))]+           [[$zs $rs] (L./ as bs)]}+      [(sum' (map2 2#(*' %1 (**' x %2)) (reverse zs) nats0))+       (sum' (map2 2#(*' %1 (**' x %2)) (reverse rs) nats0))])))+        

@@ -85,6 +85,13 @@                    [<cons $x $rs>                     [x (product' (map 2#(**' %1 %2) rs))]])]                 [_ {}]})]}]+     [<ncons $ ,$k $> [math-expr mult-expr]+      {[$tgt (match tgt math-expr+               {[<term _ $xs>+                 (match-all xs (list [math-expr integer])+                   [<join $hs <cons [$x (& ?(gte? $ k) $n)] $ts>>+                    [x (product' (map 2#(**' %1 %2) {@hs [x (- n k)] @ts}))]])]+                [_ {}]})]}]      [<ncons $ $ $> [math-expr integer mult-expr]       {[$tgt (match tgt math-expr                {[<term _ $xs>@@ -242,6 +249,10 @@                                         xs))]}))                      {@ts1 @ts2})]}))))   +(define $find-symbols-from-poly+  (lambda [$poly]+    (match-all poly math-expr+      [<poly <cons <term _ <cons (& <symbol _> $s) _>> _>> s])))  ;;; ;;; Substitute@@ -267,29 +278,81 @@   (lambda [%xs %ys $mexpr]     (substitute (zip (tensor-to-list xs) (tensor-to-list ys)) mexpr))) +(define $expand-all+  (lambda [$mexpr]+    (match mexpr math-expr+      {+       [?symbol? mexpr]+       ; function application+       [<apply $g $args>+        (capply g (map expand-all args))]+       ; quote+       [<quote $g> g]+       ; term (multiplication)+       [<term $a $ps>+        (* a (product (map 2#(** (expand-all %1) (expand-all %2)) ps)))]+       ; polynomial+       [<poly $ts> (sum (map (expand-all $) ts))]+       ; quotient+       [(/ $p1 $p2)+        (let {[$p1' (expand-all p1)]+              [$p2' (expand-all p2)]}+          (/ p1' p2'))]+       })))++(define $expand-all'+  (lambda [$mexpr]+    (match mexpr math-expr+      {+       [?symbol? mexpr]+       ; function application+       [<apply $g $args>+        (capply g (map expand-all' args))]+       ; quote+       [<quote $g> g]+       ; term (multiplication)+       [<term $a $ps>+        (*' a (product' (map 2#(**' (expand-all' %1) (expand-all' %2)) ps)))]+       ; polynomial+       [<poly $ts> (sum' (map (expand-all' $) ts))]+       ; quotient+       [(/ $p1 $p2)+        (let {[$p1' (expand-all' p1)]+              [$p2' (expand-all' p2)]}+          (/' p1' p2'))]+       })))+ ;;; ;;; Coefficient ;;; (define $coefficients-  (lambda [$x $f]+  (lambda [$f $x]     (let {[$m (capply max {0 @(match-all f math-expr                                 [<div <poly <cons <term $a <ncons ,x $k $ts>> _>> _> k])})]}-      (map (coefficient x $ f) (between 0 m)))))+      (map (coefficient f x $) (between 0 m)))))  (define $coefficient-  (lambda [$x $m $f]+  (lambda [$f $x $m]     (if (eq? m 0)       (/ (sum (match-all f math-expr                 [<div <poly <cons <term $a (& !<cons ,x _> $ts)> _>> _>                  (foldl *' a (map 2#(**' %1 %2) ts))]))          (denominator f))-      (coefficient' x m f))))+      (coefficient' f x m))))  (define $coefficient'-  (lambda [$x $m $f]+  (lambda [$f $x $m]     (/ (sum (match-all f math-expr               [<div <poly <cons <term $a <ncons ,x $k $ts>> _>> _>                (if (eq? m k)                  (foldl *' a (map 2#(**' %1 %2) ts))                  0)]))+       (denominator f))))++(define $coefficient2+  (lambda [$f $x $y]+    (/ (sum (match-all f math-expr+              [<div <poly <cons <term $a <cons ,x <cons ,y $ts>>> _>> _>+               (foldl *' a (map 2#(**' %1 %2) ts))+               ]))        (denominator f))))

@@ -6,10 +6,21 @@  (define $math-normalize   (lambda [$fn $x1 $x2]-    (if (and (integer? x1) (integer? x2))+    (if (or (integer? x1) (integer? x2))       (fn x1 x2)-      ((capply compose (map 2#%1 (filter 2#(%2 x1 x2) rewrite-rules))) (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 $rewrite-rules   {    [id 2##t]@@ -20,8 +31,8 @@    [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 (contain-function? cos %1) (contain-function? sin %1) (contain-function? cos %2) (contain-function? sin %2))]-;   [rewrite-rule-for-∂/∂ 2#(and (contain-function-with-index? %1) (contain-function-with-index? %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))]    })  ;;@@ -111,7 +122,7 @@   (lambda [$term]     (match term math-expr       {[(* $a (,sqrt $x) (,sqrt ,x) $r)-        (rewrite-rule-for-sqrt-term (*' a x r))]+        (rewrite-rule-for-sqrt (*' a x r))]        [(* $a (,sqrt (& ?term? $x)) (,sqrt (& ?term? $y)) $r)         (let* {[$d (gcd x y)]                [[$a1 $x1] (from-monomial (/ x d))]@@ -139,8 +150,38 @@ ;; cos, sin ;; -(define $rewrite-rule-for-cos-and-sin (map-polys rewrite-rule-for-cos-and-sin-poly $))+;(define $rewrite-rule-for-cos-and-sin 1#(rewrite-rule-for-cos-and-sin-expr (map-polys rewrite-rule-for-cos-and-sin-poly %1)))+(define $rewrite-rule-for-cos-and-sin 1#(map-polys rewrite-rule-for-cos-and-sin-poly %1)) +(define $rewrite-rule-for-cos-and-sin-expr+  (lambda [$expr]+    (match [expr expr] [math-expr math-expr]+      {[[<div (+ (* $a (,cos $θ) $mr)+                 $pr1)+              $pr2>+         (| <div (+ (* _ (| (,cos ,(/ θ 2)) (,sin ,(/ θ 2))) _) _) _>+            <div _ (+ (* _ (| (,cos ,(/ θ 2)) (,sin ,(/ θ 2))) _) _)>)]+        (rewrite-rule-for-cos-and-sin-expr (/' (+' (*' a (-' (cos (/ θ 2))^2 (sin (/ θ 2))^2) mr) pr1) pr2))]+       [[<div (+ (* $a (,sin $θ) $mr)+                 $pr1)+              $pr2>+         (| <div (+ (* _ (| (,cos ,(/ θ 2)) (,sin ,(/ θ 2))) _) _) _>+            <div _ (+ (* _ (| (,cos ,(/ θ 2)) (,sin ,(/ θ 2))) _) _)>)]+        (rewrite-rule-for-cos-and-sin-expr (/' (+' (*' (*' a 2) (*' (cos (/ θ 2)) (sin (/ θ 2))) mr) pr1) pr2))]+       [[<div $pr2+              (+ (* $a (,cos $θ) $mr)+                 $pr1)>+         (| <div (+ (* _ (| (,cos ,(/ θ 2)) (,sin ,(/ θ 2))) _) _) _>+            <div _ (+ (* _ (| (,cos ,(/ θ 2)) (,sin ,(/ θ 2))) _) _)>)]+        (rewrite-rule-for-cos-and-sin-expr (/' pr2 (+' (*' a (-' (cos (/ θ 2))^2 (sin (/ θ 2))^2) mr) pr1)))]+       [[<div $pr2+              (+ (* $a (,sin $θ) $mr)+                 $pr1)>+         (| <div (+ (* _ (| (,cos ,(/ θ 2)) (,sin ,(/ θ 2))) _) _) _>+            <div _ (+ (* _ (| (,cos ,(/ θ 2)) (,sin ,(/ θ 2))) _) _)>)]+        (rewrite-rule-for-cos-and-sin-expr (/' pr2 (+' (*' (*' a 2) (*' (cos (/ θ 2)) (sin (/ θ 2))) mr) pr1)))]+       [_ expr]})))+ (define $rewrite-rule-for-cos-and-sin-poly   (lambda [$poly]     (match poly math-expr@@ -148,14 +189,26 @@            (* ,a (,sin ,θ)^,2 ,mr)            $pr)         (rewrite-rule-for-cos-and-sin-poly (+' pr (*' a mr)))]+;       [(+ (* $a (,cos $θ)^,2 $mr)+;           (* ,(* -1 a) (,sin ,θ)^,2 ,mr)+;           $pr)+;        (rewrite-rule-for-cos-and-sin-poly (+' pr (*' a (cos (* 2 θ)) mr)))]+;       [(+ (* $a (,cos $θ) (,sin ,θ) $mr)+;           $pr)+;        (rewrite-rule-for-cos-and-sin-poly (+' pr (*' (/ a 2) (sin (* 2 θ)) mr)))]        [(+ (* $a $mr)-           (* ,(*' -1 a) (,sin $θ)^,2 ,mr)+           (* ,(* -1 a) (,sin $θ)^,2 ,mr)            $pr)-        (rewrite-rule-for-cos-and-sin-poly (+' (*' a (cos θ)^2 mr)))]+        (rewrite-rule-for-cos-and-sin-poly (+' pr (*' a (cos θ)^2 mr)))]        [(+ (* $a $mr)-           (* ,(*' -1 a) (,cos $θ)^,2 ,mr)+           (* ,(* -1 a) (,cos $θ)^,2 ,mr)            $pr)-        (rewrite-rule-for-cos-and-sin-poly (+' (*' a (sin θ)^2)))]+        (rewrite-rule-for-cos-and-sin-poly (+' pr (*' a (sin θ)^2 mr)))]+       ; (- 1 (cos (* 2 θ))) = (* 2 (sin θ)^2)+;       [(+ (* $a $mr)+;           (* ,(* -1 a) (,cos $θ) ,mr)+;           $pr)+;        (rewrite-rule-for-cos-and-sin-poly (+' pr (*' (* 2 a) (sin (/ θ 2))^2 mr)))]        [_ poly]})))  ;;