diff --git a/egison.cabal b/egison.cabal --- a/egison.cabal +++ b/egison.cabal @@ -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 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 @@ -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 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 @@ -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 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 @@ -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 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 @@ -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 _) = "#" @@ -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 diff --git a/lib/core/collection.egi b/lib/core/collection.egi --- a/lib/core/collection.egi +++ b/lib/core/collection.egi @@ -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)] + { + [[ ] + (let {[[$zs $rs] (L./ {@(map2 - (take (length yrs) xrs) (map (* $ (/ x y)) yrs)) + @(drop (length yrs) xrs)} ys)]} + [{(/ x y) @zs} rs])] + })))) + ;;; ;;; Multiset ;;; diff --git a/lib/math/algebra/equations.egi b/lib/math/algebra/equations.egi --- a/lib/math/algebra/equations.egi +++ b/lib/math/algebra/equations.egi @@ -26,7 +26,7 @@ (define $q-f (lambda [$f $x] - (match (coefficients x f) (list math-expr) + (match (coefficients f x) (list math-expr) {[>>> (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) {[>>>> (c-f' a_3 a_2 a_1 a_0)]}))) 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 @@ -31,10 +31,14 @@ (define $M.inverse (lambda [%m] - (match (tensor-size m) (list integer) - {[>> - (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 + {[ + (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 @@ {[> {[tgt_1_1 tgt_1_[2 n] tgt_[2 m]_1 tgt_[2 m]_[2 n]]}] [_ {}]})]}] + [ [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) - {[>> + {[>> 1] + [>> (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)}))) 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 @@ -17,9 +17,9 @@ [
> (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)))) 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 @@ -21,7 +21,11 @@ (sum (map 2#(* (capply `g|%1 args) (∂/∂ %2 x)) (zip nats args)))] ; quote - [ (∂/∂ g x)] + [ + (let {[$g' (∂/∂ g x)]} + (if (monomial? g') + g' + 'g'))] ; term (constant) [,0 0] [(* _ ,1) 0] 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 @@ -24,7 +24,7 @@ (match [x1 x2] [math-expr math-expr] {[[
] (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 {[
> - (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 +; [ _>> 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))]))) + diff --git a/lib/math/expression.egi b/lib/math/expression.egi --- a/lib/math/expression.egi +++ b/lib/math/expression.egi @@ -85,6 +85,13 @@ [ [x (product' (map 2#(**' %1 %2) rs))]])] [_ {}]})]}] + [ [math-expr mult-expr] + {[$tgt (match tgt math-expr + {[ + (match-all xs (list [math-expr integer]) + [> + [x (product' (map 2#(**' %1 %2) {@hs [x (- n k)] @ts}))]])] + [_ {}]})]}] [ [math-expr integer mult-expr] {[$tgt (match tgt math-expr {[ @@ -242,6 +249,10 @@ xs))]})) {@ts1 @ts2})]})))) +(define $find-symbols-from-poly + (lambda [$poly] + (match-all poly math-expr + [ $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 + [ + (capply g (map expand-all args))] + ; quote + [ g] + ; term (multiplication) + [ + (* a (product (map 2#(** (expand-all %1) (expand-all %2)) ps)))] + ; polynomial + [ (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 + [ + (capply g (map expand-all' args))] + ; quote + [ g] + ; term (multiplication) + [ + (*' a (product' (map 2#(**' (expand-all' %1) (expand-all' %2)) ps)))] + ; polynomial + [ (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 [
> _>> _> 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 [
$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 [
> _>> _> (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 + [
>> _>> _> + (foldl *' a (map 2#(**' %1 %2) ts)) + ])) (denominator f)))) diff --git a/lib/math/normalize.egi b/lib/math/normalize.egi --- a/lib/math/normalize.egi +++ b/lib/math/normalize.egi @@ -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] + {[[
+ (|
+
)] + (rewrite-rule-for-cos-and-sin-expr (/' (+' (*' a (-' (cos (/ θ 2))^2 (sin (/ θ 2))^2) mr) pr1) pr2))] + [[
+ (|
+
)] + (rewrite-rule-for-cos-and-sin-expr (/' (+' (*' (*' a 2) (*' (cos (/ θ 2)) (sin (/ θ 2))) mr) pr1) pr2))] + [[
+ (|
+
)] + (rewrite-rule-for-cos-and-sin-expr (/' pr2 (+' (*' a (-' (cos (/ θ 2))^2 (sin (/ θ 2))^2) mr) pr1)))] + [[
+ (|
+
)] + (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]}))) ;;