diff --git a/Algebra/Algorithms/Groebner.hs b/Algebra/Algorithms/Groebner.hs --- a/Algebra/Algorithms/Groebner.hs +++ b/Algebra/Algorithms/Groebner.hs @@ -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 + diff --git a/computational-algebra.cabal b/computational-algebra.cabal --- a/computational-algebra.cabal +++ b/computational-algebra.cabal @@ -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.* diff --git a/examples/bench.hs b/examples/bench.hs new file mode 100644 --- /dev/null +++ b/examples/bench.hs @@ -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 + ]