computational-algebra 0.0.3.2 → 0.1.0.0
raw patch · 3 files changed
+122/−17 lines, 3 filesdep +heapsdep +monad-loops
Dependencies added: heaps, monad-loops
Files
- Algebra/Algorithms/Groebner.hs +84/−16
- computational-algebra.cabal +3/−1
- examples/bench.hs +35/−0
Algebra/Algorithms/Groebner.hs view
@@ -1,23 +1,32 @@-{-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts, GADTs #-}-{-# LANGUAGE MultiParamTypeClasses, NoImplicitPrelude, ParallelListComp #-}-{-# LANGUAGE RankNTypes, ScopedTypeVariables, TypeOperators #-}+{-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts, GADTs #-}+{-# LANGUAGE MultiParamTypeClasses, NoImplicitPrelude, ParallelListComp #-}+{-# LANGUAGE RankNTypes, ScopedTypeVariables, TemplateHaskell, TypeOperators #-} module Algebra.Algorithms.Groebner ( -- * Polynomial division divModPolynomial, divPolynomial, modPolynomial -- * Groebner basis , calcGroebnerBasis, calcGroebnerBasisWith- , simpleBuchberger, reduceMinimalGroebnerBasis, minimizeGroebnerBasis+ , buchberger, syzygyBuchberger, simpleBuchberger, primeTestBuchberger+ , reduceMinimalGroebnerBasis, minimizeGroebnerBasis -- * Ideal operations , isIdealMember, intersection, thEliminationIdeal , quotIdeal, quotByPrincipalIdeal , saturationIdeal, saturationByPrincipalIdeal ) where-import Algebra.Internal-import Algebra.Ring.Noetherian-import Algebra.Ring.Polynomial-import Data.List-import Numeric.Algebra-import Prelude hiding (Num (..), recip)+import Algebra.Internal+import Algebra.Ring.Noetherian+import Algebra.Ring.Polynomial+import Control.Applicative+import Control.Monad+import Control.Monad.Loops+import Control.Monad.ST+import qualified Data.Foldable as H+import Data.Function+import qualified Data.Heap as H+import Data.List+import Data.STRef+import Numeric.Algebra+import Prelude hiding (Num (..), recip) -- | Calculate a polynomial quotient and remainder w.r.t. second argument. divModPolynomial :: (IsMonomialOrder order, IsPolynomial r n, Field r)@@ -53,7 +62,7 @@ infixl 7 `modPolynomial` infixl 7 `divModPolynomial` --- | Apply Buchberger's algorithm and calculate Groebner basis for the given ideal.+-- | The Naive buchberger's algorithm to calculate Groebner basis for the given ideal. simpleBuchberger :: (Field k, IsPolynomial k n, IsMonomialOrder order) => Ideal (OrderedPolynomial k order n) -> [OrderedPolynomial k order n] simpleBuchberger ideal =@@ -65,6 +74,63 @@ , let q = sPolynomial f g `modPolynomial` acc, q /= zero ] +-- | Buchberger's algorithm slightly improved by discarding relatively prime pair.+primeTestBuchberger :: (Field r, IsPolynomial r n, IsMonomialOrder order)+ => Ideal (OrderedPolynomial r order n) -> [OrderedPolynomial r order n]+primeTestBuchberger ideal =+ let gs = nub $ generators ideal+ in fst $ until (null . snd) (\(ggs, acc) -> let cur = nub $ ggs ++ acc in+ (cur, calc cur)) (gs, calc gs)+ where+ calc acc = [ q | f <- acc, g <- acc, f /= g+ , let f0 = leadingMonomial f, let g0 = leadingMonomial g+ , lcmMonomial f0 g0 /= zipWithV (+) f0 g0+ , let q = sPolynomial f g `modPolynomial` acc, q /= zero+ ]++(.=) :: STRef s a -> a -> ST s ()+x .= v = writeSTRef x v++(%=) :: STRef s a -> (a -> a) -> ST s ()+x %= f = modifySTRef x f++combinations :: [a] -> [(a, a)]+combinations xs = concat $ zipWith (map . (,)) xs $ drop 1 $ tails xs++-- | Calculate Groebner basis applying (modified) Buchberger's algorithm.+-- This function is same as 'syzygyBuchberger'.+buchberger :: (Field r, IsPolynomial r n, IsMonomialOrder order)+ => Ideal (OrderedPolynomial r order n) -> [OrderedPolynomial r order n]+buchberger = syzygyBuchberger++-- | Buchberger's algorithm greately improved using the syzygy theory.+-- Utilizing priority queues, this function reduces division complexity and comparison time.+-- If you don't have strong reason to avoid this function, this function is recommended to use.+syzygyBuchberger :: (Field r, IsPolynomial r n, IsMonomialOrder order)+ => Ideal (OrderedPolynomial r order n) -> [OrderedPolynomial r order n]+syzygyBuchberger ideal = runST $ do+ let gens = generators ideal+ lcmm = lcmMonomial `on` leadingMonomial+ gs <- newSTRef $ H.fromList [H.Entry (leadingMonomial g) g | g <- gens]+ b <- newSTRef $ H.fromList $ [H.Entry (lcmm f g) (f, g) | (f, g) <- combinations gens]+ whileM_ (not . H.null <$> readSTRef b) $ do+ Just (H.Entry _ (f, g), rest) <- H.viewMin <$> readSTRef b+ gs0 <- readSTRef gs+ b .= rest+ let f0 = leadingMonomial f+ g0 = leadingMonomial g+ l = lcmMonomial f0 g0+ redundant = H.any (\(H.Entry _ h) -> h `notElem` [f, g]+ && all (\k -> H.any ((==k) . H.payload) rest) [(f, h), (g, h), (h, f), (h, g)]+ && leadingMonomial h `divs` l) gs0+ when (l /= zipWithV (+) f0 g0 && not redundant) $ do+ let qs = (H.toList gs0)+ s = sPolynomial f g `modPolynomial` map H.payload qs+ when (s /= zero) $ do+ b %= H.union (H.fromList [H.Entry (lcmm q s) (q, s) | H.Entry _ q <- qs])+ gs %= H.insert (H.Entry (leadingMonomial s) s)+ map H.payload . H.toList <$> readSTRef gs+ -- | Minimize the given groebner basis. minimizeGroebnerBasis :: (Field k, IsPolynomial k n, IsMonomialOrder order) => [OrderedPolynomial k order n] -> [OrderedPolynomial k order n]@@ -84,10 +150,10 @@ => order -> Ideal (OrderedPolynomial k order' n) -> [OrderedPolynomial k order n] calcGroebnerBasisWith ord i = calcGroebnerBasis $ mapIdeal (changeOrder ord) i --- | Caliculating reduced Groebner basis of the given ideal w.r.t. the graded reversed lexicographic order.+-- | Caliculating reduced Groebner basis of the given ideal. calcGroebnerBasis :: (Field k, IsPolynomial k n, IsMonomialOrder order) => Ideal (OrderedPolynomial k order n) -> [OrderedPolynomial k order n]-calcGroebnerBasis = reduceMinimalGroebnerBasis . minimizeGroebnerBasis . simpleBuchberger+calcGroebnerBasis = reduceMinimalGroebnerBasis . minimizeGroebnerBasis . syzygyBuchberger -- | Test if the given polynomial is the member of the ideal. isIdealMember :: (IsPolynomial k n, Field k, IsMonomialOrder o)@@ -106,9 +172,10 @@ -> Ideal (OrderedPolynomial k ord m) -> Ideal (OrderedPolynomial k Lex (m :-: n)) thEliminationIdeal n ideal =- toIdeal $ [transformMonomial (dropV n) f | f <- calcGroebnerBasisWith Lex ideal- , all (all (== 0) . take (toInt n) . toList . snd) $ getTerms f- ]+ toIdeal $ [ transformMonomial (dropV n) f+ | f <- calcGroebnerBasisWith Lex ideal+ , all (all (== 0) . take (toInt n) . toList . snd) $ getTerms f+ ] -- | An intersection ideal of given ideals. intersection :: forall r k n ord.@@ -166,3 +233,4 @@ SingInstance -> case singInstance (sLengthV g %+ (sing :: SNat n)) of SingInstance -> intersection $ mapV (i `saturationByPrincipalIdeal`) g+
computational-algebra.cabal view
@@ -2,7 +2,7 @@ -- further documentation, see http://haskell.org/cabal/users-guide/ name: computational-algebra-version: 0.0.3.2+version: 0.1.0.0 synopsis: Well-kinded computational algebra library, currently supporting Groebner basis. description: Dependently-typed computational algebra libray for Groebner basis. homepage: https://github.com/konn/computational-algebra@@ -34,5 +34,7 @@ , lens == 3.* , containers >= 0.4 && < 0.6 , peggy == 0.3.*+ , monad-loops == 0.3.*+ , heaps == 0.2.* if impl(ghc >= 7.6.1) build-depends: monomorphic == 0.0.*
+ examples/bench.hs view
@@ -0,0 +1,35 @@+{-# LANGUAGE DataKinds, FlexibleContexts, FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses, OverloadedStrings, PolyKinds #-}+{-# LANGUAGE TemplateHaskell, UndecidableInstances #-}+import Algebra.Algorithms.Groebner.Monomorphic+import Algebra.Ring.Polynomial.Monomorphic+import Control.DeepSeq+import Criterion.Types+import Numeric.Algebra (LeftModule (..))+import qualified Numeric.Algebra as NA+import Progression.Main++x, y, z, w, s, a, b, c :: Polynomial Rational+[x, y, z, w, s, a, b, c] =+ map (injectVar . flip Variable Nothing) "xyzwSabc"++instance NFData Variable where+ rnf (Variable x y) = rnf x `seq` rnf y `seq` ()++instance NFData (Polynomial Rational) where+ rnf (Polynomial dic) = rnf dic++ideal1, ideal2, ideal3, ideal4 :: [Polynomial Rational]+ideal1 = [x^2 + y^2 + z^2 - 1, x^2 + y^2 + z^2 - 2*x, 2*x -3*y - z]+ideal2 = [x^2 * y - 2*x*y - 4*z - 1, z-y^2, x^3 - 4*z*y]+ideal3 = [ 2 * s - a * y, b^2 - (x^2 + y^2), c^2 - ( (a-x) ^ 2 + y^2)+ ]+ideal4 = [ x^5 + y^4 + z^3 - 1, x^3 + y^3 + z^2 - 1]++main :: IO ()+main =+ defaultMain $ bgroup "groebner"+ [ bench "grevlex" $ nf calcGroebnerBasis ideal4+ , bench "grlex" $ nf (calcGroebnerBasisWith Grlex) ideal4+ , bench "lex" $ nf (calcGroebnerBasisWith Lex) ideal4+ ]