computational-algebra 0.1.2.0 → 0.1.3.0
raw patch · 11 files changed
+515/−126 lines, 11 files
Files
- Algebra/Algorithms/Groebner.hs +130/−30
- Algebra/Algorithms/Groebner/Monomorphic.hs +29/−12
- Algebra/Internal.hs +4/−1
- Algebra/Ring/Polynomial.hs +88/−41
- computational-algebra.cabal +1/−1
- examples/bench.hs +86/−21
- examples/elimination-bench.hs +68/−0
- examples/polymorphic.hs +0/−3
- examples/sandpit-poly.hs +36/−0
- examples/sandpit.hs +19/−17
- examples/sugar-bench.hs +54/−0
Algebra/Algorithms/Groebner.hs view
@@ -1,13 +1,21 @@-{-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts, GADTs #-}-{-# LANGUAGE MultiParamTypeClasses, NoImplicitPrelude, ParallelListComp #-}-{-# LANGUAGE RankNTypes, ScopedTypeVariables, TemplateHaskell, TypeOperators #-}+{-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts, FlexibleInstances #-}+{-# LANGUAGE GADTs, MultiParamTypeClasses, NoImplicitPrelude #-}+{-# LANGUAGE ParallelListComp, RankNTypes, ScopedTypeVariables #-}+{-# LANGUAGE StandaloneDeriving, TemplateHaskell, TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# OPTIONS_GHC -fno-warn-type-defaults -fno-warn-orphans #-} module Algebra.Algorithms.Groebner ( -- * Polynomial division divModPolynomial, divPolynomial, modPolynomial -- * Groebner basis- , calcGroebnerBasis, calcGroebnerBasisWith- , buchberger, syzygyBuchberger, simpleBuchberger, primeTestBuchberger+ , calcGroebnerBasis, calcGroebnerBasisWith, calcGroebnerBasisWithStrategy+ , buchberger, syzygyBuchberger+ , simpleBuchberger, primeTestBuchberger , reduceMinimalGroebnerBasis, minimizeGroebnerBasis+ -- ** Selection Strategies+ , syzygyBuchbergerWithStrategy+ , SelectionStrategy(..), calcWeight', GrevlexStrategy(..)+ , NormalStrategy(..), SugarStrategy(..), GradedStrategy(..) -- * Ideal operations , isIdealMember, intersection, thEliminationIdeal, thEliminationIdealWith , unsafeThEliminationIdealWith@@ -96,6 +104,20 @@ (%=) :: STRef s a -> (a -> a) -> ST s () x %= f = modifySTRef x f +padVec :: a -> Vector a n -> Vector a m -> (Vector a (Max n m), Vector a (Max n m))+padVec _ Nil Nil = (Nil, Nil)+padVec def (x :- xs) (y :- ys) =+ case padVec def xs ys of+ (xs', ys') -> (x :- xs', y :- ys')+padVec def (x :- xs) Nil =+ case padVec def xs Nil of+ (xs', ys') -> case maxZR (sLengthV xs) of+ Eql -> (x :- xs', def :- ys')+padVec def Nil (y :- ys) =+ case padVec def Nil ys of+ (xs', ys') -> case maxZL (sLengthV ys) of+ Eql -> (def :- xs', y :- ys')+ combinations :: [a] -> [(a, a)] combinations xs = concat $ zipWith (map . (,)) xs $ drop 1 $ tails xs @@ -105,34 +127,107 @@ => Ideal (OrderedPolynomial r order n) -> [OrderedPolynomial r order n] buchberger = syzygyBuchberger --- | Buchberger's algorithm greately improved using the syzygy theory.+-- | Buchberger's algorithm greately improved using the syzygy theory with the sugar strategy. -- 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]+syzygyBuchberger = syzygyBuchbergerWithStrategy (SugarStrategy NormalStrategy)++(=@=) :: Vector Int n -> Vector Int m -> Bool+Nil =@= Nil = True+(x :- xs) =@= (y :- ys) = x == y && xs =@= ys+Nil =@= (0 :- ys) = Nil =@= ys+(0 :- xs) =@= Nil = xs =@= Nil+_ =@= _ = False++instance Eq (Monomorphic (OrderedMonomial ord)) where+ Monomorphic xs == Monomorphic ys = getMonomial xs =@= getMonomial ys++instance IsMonomialOrder ord => Ord (Monomorphic (OrderedMonomial ord)) where+ compare (Monomorphic (OrderedMonomial m)) (Monomorphic (OrderedMonomial m')) =+ let (mm, mm') = padVec 0 m m'+ in cmpMonomial (Proxy :: Proxy ord) mm mm'++-- | apply buchberger's algorithm using given selection strategy.+syzygyBuchbergerWithStrategy :: ( Field r, IsPolynomial r n, IsMonomialOrder order, SelectionStrategy strategy+ , Ord (Weight strategy order))+ => strategy -> Ideal (OrderedPolynomial r order n) -> [OrderedPolynomial r order n]+syzygyBuchbergerWithStrategy strategy ideal = runST $ do+ let gens = zip [1..] $ generators ideal+ gs <- newSTRef $ H.fromList [H.Entry (leadingOrderedMonomial g) g | (_, g) <- gens]+ b <- newSTRef $ H.fromList $ [H.Entry (calcWeight' strategy f g, j) (f, g) | ((_, f), (j, g)) <- combinations gens]+ len <- newSTRef (genericLength gens :: Integer) whileM_ (not . H.null <$> readSTRef b) $ do- Just (H.Entry _ (f, g), rest) <- H.viewMin <$> readSTRef b+ 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)]+ 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)+ len0 <- readSTRef len+ 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 (calcWeight' strategy q s, j) (q, s) | H.Entry _ q <- qs | j <- [len0+1..]])+ gs %= H.insert (H.Entry (leadingOrderedMonomial s) s)+ len %= (*2) map H.payload . H.toList <$> readSTRef gs +-- | Calculate the weight of given polynomials w.r.t. the given strategy.+-- Buchberger's algorithm proccesses the pair with the most least weight first.+-- This function requires the @Ord@ instance for the weight; this constraint is not required+-- in the 'calcWeight' because of the ease of implementation. So use this function.+calcWeight' :: (SelectionStrategy s, IsPolynomial r n, IsMonomialOrder ord, Ord (Weight s ord))+ => s -> OrderedPolynomial r ord n -> OrderedPolynomial r ord n -> Weight s ord+calcWeight' s = calcWeight (toProxy s)++-- | Type-class for selection strategies in Buchberger's algorithm.+class SelectionStrategy s where+ type Weight s ord :: *+ calcWeight :: (IsPolynomial r n, IsMonomialOrder ord)+ => Proxy s -> OrderedPolynomial r ord n -> OrderedPolynomial r ord n -> Weight s ord++-- | Buchberger's normal selection strategy. This selects the pair with+-- the least LCM(LT(f), LT(g)) w.r.t. current monomial ordering.+data NormalStrategy = NormalStrategy deriving (Read, Show, Eq, Ord)++instance SelectionStrategy NormalStrategy where+ type Weight NormalStrategy ord = Monomorphic (OrderedMonomial ord)+ calcWeight _ f g = Monomorphic $+ OrderedMonomial (lcmMonomial (leadingMonomial f) (leadingMonomial g))+ `asTypeOf` leadingOrderedMonomial f++-- | Choose the pair with the least LCM(LT(f), LT(g)) w.r.t. 'Grevlex' order.+data GrevlexStrategy = GrevlexStrategy deriving (Read, Show, Eq, Ord)++instance SelectionStrategy GrevlexStrategy where+ type Weight GrevlexStrategy ord = Monomorphic (OrderedMonomial Grevlex)+ calcWeight _ f g = Monomorphic $ OrderedMonomial $ lcmMonomial (leadingMonomial f) (leadingMonomial g)++data GradedStrategy = GradedStrategy deriving (Read, Show, Eq, Ord)++-- | Choose the pair with the least LCM(LT(f), LT(g)) w.r.t. graded current ordering.+instance SelectionStrategy GradedStrategy where+ type Weight GradedStrategy ord = Monomorphic (OrderedMonomial (Graded ord))+ calcWeight _ f g = Monomorphic $ OrderedMonomial (lcmMonomial (leadingMonomial f) (leadingMonomial g))++-- | Sugar strategy. This chooses the pair with the least phantom homogenized degree and then break the tie with the given strategy (say @s@).+data SugarStrategy s = SugarStrategy s deriving (Read, Show, Eq, Ord)++instance SelectionStrategy s => SelectionStrategy (SugarStrategy s) where+ type Weight (SugarStrategy s) ord = (Int, Weight s ord)+ calcWeight (Proxy :: Proxy (SugarStrategy s)) f g = (sugar, calcWeight (Proxy :: Proxy s) f g)+ where+ deg' = maximum . map (totalDegree . snd) . getTerms+ tsgr h = deg' h - totalDegree (leadingMonomial h)+ sugar = max (tsgr f) (tsgr g) + totalDegree (lcmMonomial (leadingMonomial f) (leadingMonomial g))+ -- | Minimize the given groebner basis. minimizeGroebnerBasis :: (Field k, IsPolynomial k n, IsMonomialOrder order) => [OrderedPolynomial k order n] -> [OrderedPolynomial k order n]@@ -142,7 +237,7 @@ -- | Reduce minimum Groebner basis into reduced Groebner basis. reduceMinimalGroebnerBasis :: (Field k, IsPolynomial k n, IsMonomialOrder order)- => [OrderedPolynomial k order n] -> [OrderedPolynomial k order n]+ => [OrderedPolynomial k order n] -> [OrderedPolynomial k order n] reduceMinimalGroebnerBasis bs = filter (/= zero) $ map red bs where red x = x `modPolynomial` delete x bs@@ -150,8 +245,15 @@ -- | Caliculating reduced Groebner basis of the given ideal w.r.t. the specified monomial order. calcGroebnerBasisWith :: (Field k, IsPolynomial k n, IsMonomialOrder order, IsMonomialOrder order') => order -> Ideal (OrderedPolynomial k order' n) -> [OrderedPolynomial k order n]-calcGroebnerBasisWith ord i = calcGroebnerBasis $ mapIdeal (changeOrder ord) i+calcGroebnerBasisWith ord i = calcGroebnerBasis $ mapIdeal (changeOrder ord) i +-- | Caliculating reduced Groebner basis of the given ideal w.r.t. the specified monomial order.+calcGroebnerBasisWithStrategy :: ( Field k, IsPolynomial k n, IsMonomialOrder order+ , SelectionStrategy strategy, Ord (Weight strategy order))+ => strategy -> Ideal (OrderedPolynomial k order n) -> [OrderedPolynomial k order n]+calcGroebnerBasisWithStrategy strategy =+ reduceMinimalGroebnerBasis . minimizeGroebnerBasis . syzygyBuchbergerWithStrategy strategy+ -- | 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]@@ -167,7 +269,7 @@ => OrderedPolynomial k order n -> [OrderedPolynomial k order n] -> Bool groebnerTest f fs = f `modPolynomial` fs == zero --- | Calculate n-th elimination ideal using 'Lex' ordering.+-- | Calculate n-th elimination ideal using 'WeightedEliminationOrder' ordering. thEliminationIdeal :: ( IsMonomialOrder ord, Field k, IsPolynomial k m, IsPolynomial k (m :-: n) , (n :<<= m) ~ True) => SNat n@@ -176,8 +278,7 @@ thEliminationIdeal n = case singInstance n of SingInstance ->- case weightVInstance n of- ToWeightVectorInstance -> mapIdeal (changeOrderProxy Proxy) . thEliminationIdealWith (weightedEliminationOrder n) n+ mapIdeal (changeOrderProxy Proxy) . thEliminationIdealWith (weightedEliminationOrder n) n -- | Calculate n-th elimination ideal using the specified n-th elimination type order. thEliminationIdealWith :: ( IsMonomialOrder ord, Field k, IsPolynomial k m, IsPolynomial k (m :-: n)@@ -209,8 +310,7 @@ , all (all (== 0) . take (toInt n) . toList . snd) $ getTerms f ] ---- | An intersection ideal of given ideals.+-- | An intersection ideal of given ideals (using 'WeightedEliminationOrder'). intersection :: forall r k n ord. ( IsMonomialOrder ord, Field r, IsPolynomial r k, IsPolynomial r n , IsPolynomial r (k :+: n)@@ -251,16 +351,16 @@ -- | Saturation by a principal ideal. saturationByPrincipalIdeal :: (Field k, IsPolynomial k n, IsMonomialOrder ord) => Ideal (OrderedPolynomial k ord n)- -> OrderedPolynomial k ord n -> Ideal (OrderedPolynomial k Lex n)+ -> OrderedPolynomial k ord n -> Ideal (OrderedPolynomial k ord n) saturationByPrincipalIdeal is g = case propToClassLeq $ leqSucc (sDegree g) of- LeqInstance -> thEliminationIdealWith Lex sOne $ addToIdeal (one - (castPolynomial g * var sOne)) (mapIdeal (shiftR sOne) is)+ LeqInstance -> thEliminationIdeal sOne $ addToIdeal (one - (castPolynomial g * var sOne)) (mapIdeal (shiftR sOne) is) -- | Saturation ideal saturationIdeal :: forall k ord n. (IsPolynomial k n, Field k, IsMonomialOrder ord) => Ideal (OrderedPolynomial k ord n) -> Ideal (OrderedPolynomial k ord n)- -> Ideal (OrderedPolynomial k Lex n)+ -> Ideal (OrderedPolynomial k ord n) saturationIdeal i (Ideal g) = case singInstance (sLengthV g) of SingInstance ->
Algebra/Algorithms/Groebner/Monomorphic.hs view
@@ -1,6 +1,6 @@-{-# LANGUAGE ConstraintKinds, FlexibleInstances, GADTs, PolyKinds #-}-{-# LANGUAGE RecordWildCards, ScopedTypeVariables, TypeFamilies #-}-{-# LANGUAGE TypeOperators, UndecidableInstances #-}+{-# LANGUAGE ConstraintKinds, FlexibleContexts, FlexibleInstances, GADTs #-}+{-# LANGUAGE PolyKinds, RecordWildCards, ScopedTypeVariables, TypeFamilies #-}+{-# LANGUAGE TypeOperators, UndecidableInstances #-} -- | Monomorphic interface for Groenber basis. module Algebra.Algorithms.Groebner.Monomorphic ( Groebnerable@@ -9,7 +9,7 @@ , divModPolynomialWith, divPolynomialWith, modPolynomialWith -- * Groebner basis , calcGroebnerBasis, calcGroebnerBasisWith- , syzygyBuchberger, syzygyBuchbergerWith+ , syzygyBuchberger, syzygyBuchbergerWith, syzygyBuchbergerWithStrategy , primeTestBuchberger, primeTestBuchbergerWith , simpleBuchberger, simpleBuchbergerWith -- * Ideal operations@@ -17,8 +17,15 @@ , quotIdeal, quotByPrincipalIdeal , saturationIdeal, saturationByPrincipalIdeal -- * Re-exports- , Lex(..), Revlex(..), Grlex(..), Grevlex(..), IsOrder, IsMonomialOrder+ , Lex(..), Revlex(..), Grlex(..), Grevlex(..), IsOrder(..), IsMonomialOrder+ , SelectionStrategy(..), NormalStrategy(..), SugarStrategy(..), Gr.GrevlexStrategy(..)+ , Gr.GradedStrategy(..)+ , calcWeight' ) where+import Algebra.Algorithms.Groebner (NormalStrategy (..),+ SelectionStrategy (..),+ SugarStrategy (..),+ calcWeight') import qualified Algebra.Algorithms.Groebner as Gr import Algebra.Internal import Algebra.Ring.Noetherian@@ -97,7 +104,7 @@ in case singInstance sn of SingInstance -> let (q, r) = Gr.divModPolynomial f' gs'- in (map (renameVars vars . polyn . demote' *** polyn . demote') q, polyn $ demote' r)+ in (map (renameVars vars . polyn . demote' *** renameVars vars . polyn . demote') q, renameVars vars $ polyn $ demote' r) where vars = nub $ sort $ concatMap buildVarsList (f:gs) @@ -165,19 +172,27 @@ syzygyBuchberger :: (Groebnerable r) => [Polynomial r] -> [Polynomial r] syzygyBuchberger = syzygyBuchbergerWith Grevlex -syzygyBuchbergerWith :: forall ord r. (Groebnerable r, IsMonomialOrder ord)- => ord -> [Polynomial r] -> [Polynomial r]-syzygyBuchbergerWith _ ps | any (== zero) ps = []-syzygyBuchbergerWith ord j =+syzygyBuchbergerWithStrategy :: forall strategy ord r.+ ( Groebnerable r, IsMonomialOrder ord+ , Gr.SelectionStrategy strategy, Ord (Gr.Weight strategy ord))+ => strategy -> ord -> [Polynomial r] -> [Polynomial r]+syzygyBuchbergerWithStrategy _ _ ps | any (== zero) ps = []+syzygyBuchbergerWithStrategy strategy _ j = case uniformlyPromote j :: Monomorphic (Ideal :.: Poly.OrderedPolynomial r ord) of Monomorphic (Comp ideal) -> case ideal of Ideal vec -> case singInstance (Poly.sDegree (head $ toList vec)) of- SingInstance -> map (renameVars vars . polyn . demote . Monomorphic) $ Gr.syzygyBuchberger ideal+ SingInstance -> map (renameVars vars . polyn . demote . Monomorphic) $ Gr.syzygyBuchbergerWithStrategy strategy ideal where vars = nub $ sort $ concatMap buildVarsList j ++syzygyBuchbergerWith :: forall ord r. (Groebnerable r, IsMonomialOrder ord)+ => ord -> [Polynomial r] -> [Polynomial r]+syzygyBuchbergerWith _ ps | any (== zero) ps = []+syzygyBuchbergerWith ord j = syzygyBuchbergerWithStrategy (SugarStrategy NormalStrategy) ord j+ isIdealMember :: forall r. Groebnerable r => Polynomial r -> [Polynomial r] -> Bool isIdealMember f ideal = case promoteList (f:ideal) :: Monomorphic ([] :.: Poly.Polynomial r) of@@ -214,7 +229,9 @@ k = length els eliminate :: forall r. Groebnerable r => [Variable] -> [Polynomial r] -> [Polynomial r]-eliminate = eliminateWith Lex+eliminate vs j = eliminateWith Lex vs j+ where+ vars = nub $ sort $ concatMap buildVarsList j -- | Computes nth elimination ideal. thEliminationIdeal :: Groebnerable r => Int -> [Polynomial r] -> [Polynomial r]
Algebra/Internal.hs view
@@ -10,7 +10,7 @@ , SZero, SOne, STwo, SThree , lengthV, sLengthV, takeV, dropV, splitAtV, appendV , foldrV, foldlV, singletonV, zipWithV, toList, allV- , mapV, headV, tailV, splitAtLess+ , mapV, headV, tailV, splitAtLess, takeVAtMost , Leq(..), (:<<=), (:<=), LeqInstance(..) , LeqTrueInstance(..), boolToPropLeq, boolToClassLeq , propToClassLeq, propToBoolLeq@@ -218,6 +218,9 @@ takeV :: (n :<<= m) ~ True => SNat n -> Vector a m -> Vector a n takeV n = fst . splitAtV n++takeVAtMost :: SNat n -> Vector a m -> Vector a (Min n m)+takeVAtMost = (fst .) . splitAtLess splitAtLess :: SNat n -> Vector a m -> (Vector a (Min n m), Vector a (m :-: n)) splitAtLess SZ v = case zAbsorbsMinL (sLengthV v) of
Algebra/Ring/Polynomial.hs view
@@ -1,21 +1,22 @@ {-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts, FlexibleInstances #-}-{-# LANGUAGE GADTs, GeneralizedNewtypeDeriving, MultiParamTypeClasses #-}-{-# LANGUAGE OverlappingInstances, PolyKinds, ScopedTypeVariables #-}-{-# LANGUAGE StandaloneDeriving, TypeFamilies, TypeOperators #-}-{-# LANGUAGE UndecidableInstances, ViewPatterns #-}+{-# LANGUAGE GADTs, GeneralizedNewtypeDeriving, LiberalTypeSynonyms #-}+{-# LANGUAGE MultiParamTypeClasses, OverlappingInstances, PolyKinds #-}+{-# LANGUAGE RankNTypes, ScopedTypeVariables, StandaloneDeriving #-}+{-# LANGUAGE TypeFamilies, TypeOperators, UndecidableInstances, ViewPatterns #-} {-# OPTIONS_GHC -fno-warn-orphans -fno-warn-type-defaults #-} module Algebra.Ring.Polynomial ( Polynomial, Monomial, MonomialOrder, EliminationType, EliminationOrder , WeightedEliminationOrder, eliminationOrder, weightedEliminationOrder , lex, revlex, graded, grlex, grevlex, productOrder, productOrder'- , transformMonomial, WeightProxy(..), weightOrder+ , transformMonomial, WeightProxy(..), weightOrder, totalDegree , IsPolynomial, coeff, lcmMonomial, sPolynomial, polynomial , castMonomial, castPolynomial, toPolynomial, changeOrder, changeOrderProxy- , scastMonomial, scastPolynomial, OrderedPolynomial, showPolynomialWithVars, showPolynomialWith, showRational, ToWeightVectorInstance(..), weightVInstance+ , scastMonomial, scastPolynomial, OrderedPolynomial, showPolynomialWithVars, showPolynomialWith, showRational , normalize, injectCoeff, varX, var, getTerms, shiftR, orderedBy- , divs, tryDiv, fromList -- , genVarsV- , leadingTerm, leadingMonomial, leadingCoeff, genVars, sDegree- , OrderedMonomial(..), Grevlex(..), Revlex(..), Lex(..), Grlex(..)+ , divs, tryDiv, fromList, Coefficient(..),ToWeightVector(..)+ , leadingTerm, leadingMonomial, leadingOrderedMonomial, leadingCoeff, genVars, sDegree+ , OrderedMonomial(..), OrderedMonomial'(..), Grevlex(..)+ , Revlex(..), Lex(..), Grlex(..), Graded(..) , ProductOrder (..), WeightOrder(..) , IsOrder(..), IsMonomialOrder) where import Algebra.Internal@@ -30,30 +31,60 @@ import Data.Ord import Data.Proxy import Data.Ratio-import Numeric.Algebra hiding (Order(..))+import Numeric.Algebra hiding (Order(..), sum) import Prelude hiding (lex, (*), (+), (-), (^), (^^), recip, negate) import qualified Prelude as P+import Data.Function -- | N-ary Monomial. IntMap contains degrees for each x_i. type Monomial (n :: Nat) = Vector Int n +-- | Monomorphic representation for monomial.+newtype OrderedMonomial' ord = OM' { getMonomial' :: [Int] }+ deriving (Read, Show, Eq)++instance Monomorphicable (OrderedMonomial ord) where+ type MonomorphicRep (OrderedMonomial ord) = OrderedMonomial' ord+ promote (OM' xs) =+ case promote xs of+ Monomorphic m -> Monomorphic $ OrderedMonomial m+ demote (Monomorphic (OrderedMonomial m)) = OM' (demote (Monomorphic m))++instance IsMonomialOrder ord => Ord (OrderedMonomial' ord) where+ compare = cmpMonomial' (Proxy :: Proxy ord) `on` getMonomial'++instance Monomorphicable (Vector Int) where+ type MonomorphicRep (Vector Int) = [Int]+ promote [] = Monomorphic Nil+ promote (n : ns) =+ case promote ns of+ Monomorphic ns' -> Monomorphic (n :- ns')+ demote (Monomorphic Nil) = []+ demote (Monomorphic (n :- ns)) = n : demote (Monomorphic ns)+ -- | convert NAry list into Monomial. fromList :: SNat n -> [Int] -> Monomial n fromList SZ _ = Nil fromList (SS n) [] = 0 :- fromList n [] fromList (SS n) (x : xs) = x :- fromList n xs +-- | apply monomial ordering to monomorphic monomials.+cmpMonomial' :: IsMonomialOrder ord => Proxy ord -> [Int] -> [Int] -> Ordering+cmpMonomial' pxy xs ys =+ withPolymorhic (max (length xs) (length ys)) $ \n ->+ cmpMonomial pxy (fromList n xs) (fromList n ys)+ -- | Monomial order (of degree n). This should satisfy following laws: -- (1) Totality: forall a, b (a < b || a == b || b < a) -- (2) Additivity: a <= b ==> a + c <= b + c -- (3) Non-negative: forall a, 0 <= a-type MonomialOrder n = Monomial n -> Monomial n -> Ordering+type MonomialOrder = forall n. Monomial n -> Monomial n -> Ordering totalDegree :: Monomial n -> Int-totalDegree = foldrV (+) 0+totalDegree = foldlV (+) 0 -- | Lexicographical order. This *is* a monomial order.-lex :: MonomialOrder n+lex :: MonomialOrder lex Nil Nil = EQ lex (x :- xs) (y :- ys) = x `compare` y <> xs `lex` ys lex _ _ = error "cannot happen"@@ -69,11 +100,11 @@ graded cmp xs ys = comparing totalDegree xs ys <> cmp xs ys -- | Graded lexicographical order. This *is* a monomial order.-grlex :: MonomialOrder n+grlex :: MonomialOrder grlex = graded lex -- | Graded reversed lexicographical order. This *is* a monomial order.-grevlex :: MonomialOrder n+grevlex :: MonomialOrder grevlex = graded revlex -- | A wrapper for monomials with a certain (monomial) order.@@ -85,20 +116,36 @@ -- | Class to lookup ordering from its (type-level) name. class IsOrder (ordering :: *) where- cmpMonomial :: Proxy ordering -> MonomialOrder n+ cmpMonomial :: Proxy ordering -> MonomialOrder -- * Names for orderings. -- We didn't choose to define one single type for ordering names for the extensibility.-data Grevlex = Grevlex -- Graded reversed lexicographical order- deriving (Show, Eq, Ord)-data Lex = Lex -- Lexicographical order+-- | Lexicographical order+data Lex = Lex deriving (Show, Eq, Ord)-data Grlex = Grlex -- Graded lexicographical order- deriving (Show, Eq, Ord)-data Revlex = Revlex -- Reversed lexicographical order++-- | Reversed lexicographical order+data Revlex = Revlex deriving (Show, Eq, Ord) -data ProductOrder n a b where+-- | Graded reversed lexicographical order. Same as @Graded Revlex@.+data Grevlex = Grevlex+ deriving (Show, Eq, Ord)++-- | Graded lexicographical order. Same as @Graded Lex@.+data Grlex = Grlex+ deriving (Show, Eq, Ord)++-- | Graded order from another monomial order.+data Graded ord = Graded ord+ deriving (Read, Show, Eq, Ord)++instance IsOrder ord => IsOrder (Graded ord) where+ cmpMonomial Proxy = graded (cmpMonomial (Proxy :: Proxy ord))++instance IsMonomialOrder ord => IsMonomialOrder (Graded ord)++data ProductOrder (n :: Nat) (a :: *) (b :: *) where ProductOrder :: SNat n -> ord -> ord' -> ProductOrder n ord ord' productOrder :: forall ord ord' n m. (IsOrder ord, IsOrder ord', Sing n)@@ -180,33 +227,29 @@ type EliminationOrder n = ProductOrder n Grevlex Grevlex -data ToWeightVectorInstance n where- ToWeightVectorInstance :: (EliminationType n (WeightedEliminationOrder n), ToWeightVector (EWeight n)) => ToWeightVectorInstance n--weightVInstance :: SNat n -> ToWeightVectorInstance n-weightVInstance SZ = ToWeightVectorInstance-weightVInstance (SS n) =- case weightVInstance n of- ToWeightVectorInstance -> ToWeightVectorInstance- eliminationOrder :: SNat n -> EliminationOrder n eliminationOrder n = case singInstance n of SingInstance -> ProductOrder n Grevlex Grevlex -weightedEliminationOrder :: SNat n -> WeightedEliminationOrder n-weightedEliminationOrder n = WeightOrder (mkWeight n) Grevlex--mkWeight :: SNat n -> WeightProxy (EWeight n)-mkWeight SZ = NilWeight-mkWeight (SS n) = ConsWeight sOne $ mkWeight n-+weightedEliminationOrder :: SNat n -> WeightedEliminationOrder n Grevlex+weightedEliminationOrder n = WEOrder n (Proxy :: Proxy Grevlex) type family EWeight n :: [Nat] type instance EWeight Z = '[] type instance EWeight (S n) = One ': EWeight n -type WeightedEliminationOrder n = WeightOrder (EWeight n) Grevlex+data WeightedEliminationOrder (n :: Nat) (ord :: *) where+ WEOrder :: SNat n -> Proxy ord -> WeightedEliminationOrder n ord +instance (Sing n, IsMonomialOrder ord) => IsOrder (WeightedEliminationOrder n ord) where+ cmpMonomial Proxy m m' = comparing calc m m' <> cmpMonomial (Proxy :: Proxy ord) m m'+ where+ calc = foldlV (+) 0 . takeVAtMost (sing :: SNat n)++instance (Sing n, IsMonomialOrder ord) => IsMonomialOrder (WeightedEliminationOrder n ord)++instance (Sing n, IsMonomialOrder ord) => EliminationType n (WeightedEliminationOrder n ord) where+ -- | Special ordering for ordered-monomials. instance (Eq (Monomial n), IsOrder name) => Ord (OrderedMonomial name n) where OrderedMonomial m `compare` OrderedMonomial n = cmpMonomial (Proxy :: Proxy name) m n@@ -250,7 +293,7 @@ normalize = unwrapped %~ M.insertWith (+) (OrderedMonomial $ fromList sing []) zero . M.filter (/= zero) instance (Eq r, IsOrder order, IsPolynomial r n) => Eq (OrderedPolynomial r order n) where- (normalize -> Polynomial f) == (normalize -> Polynomial g) = f == g+ Polynomial f == Polynomial g = f == g injectCoeff :: (IsPolynomial r n) => r -> OrderedPolynomial r order n injectCoeff r = Polynomial $ M.singleton (OrderedMonomial $ fromList sing []) r@@ -385,6 +428,10 @@ leadingMonomial :: (IsOrder order, IsPolynomial r n) => OrderedPolynomial r order n -> Monomial n leadingMonomial = snd . leadingTerm++leadingOrderedMonomial :: (IsOrder order, IsPolynomial r n)+ => OrderedPolynomial r order n -> OrderedMonomial order n+leadingOrderedMonomial = OrderedMonomial . leadingMonomial leadingCoeff :: (IsOrder order, IsPolynomial r n) => OrderedPolynomial r order n -> r leadingCoeff = fst . leadingTerm
computational-algebra.cabal view
@@ -2,7 +2,7 @@ -- further documentation, see http://haskell.org/cabal/users-guide/ name: computational-algebra-version: 0.1.2.0+version: 0.1.3.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
examples/bench.hs view
@@ -1,13 +1,13 @@-{-# 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+{-# LANGUAGE DataKinds, FlexibleContexts, FlexibleInstances, GADTs #-}+{-# LANGUAGE MultiParamTypeClasses, OverloadedStrings, PolyKinds #-}+{-# LANGUAGE QuasiQuotes, TemplateHaskell, UndecidableInstances #-}+{-# OPTIONS_GHC -fno-warn-type-defaults -fno-warn-orphans #-}+import Algebra.Algorithms.Groebner.Monomorphic+import Algebra.Ring.Polynomial.Monomorphic+import Control.DeepSeq+import Control.Parallel.Strategies+import Criterion.Main+import SingularBench x, y, z, w, s, a, b, c :: Polynomial Rational [x, y, z, w, s, a, b, c] = map (injectVar . flip Variable Nothing) "xyzwSabc"@@ -18,17 +18,82 @@ 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)+i1, i2, i3, i4 :: [Polynomial Rational]+i1 = [x^2 + y^2 + z^2 - 1, x^2 + y^2 + z^2 - 2*x, 2*x -3*y - z]+i2 = [x^2 * y - 2*x*y - 4*z - 1, z-y^2, x^3 - 4*z*y]+i3 = [ 2 * s - a * y, b^2 - (x^2 + y^2), c^2 - ( (a-x) ^ 2 + y^2) ]-ideal4 = [ z^5 + y^4 + x^3 - 1, z^3 + y^3 + x^2 - 1]+i4 = [ z^5 + y^4 + x^3 - 1, z^3 + y^3 + x^2 - 1] main :: IO ()-main =- defaultMain $ bgroup "groebner"- [ bench "simple" $ nf (simpleBuchbergerWith Lex) ideal3- , bench "relprime" $ nf (primeTestBuchbergerWith Lex) ideal3- , bench "syzygy" $ nf (syzygyBuchbergerWith Lex) ideal3- ]+main = do+ ideal1 <- return $! (i1 `using` rdeepseq)+ ideal2 <- return $! (i2 `using` rdeepseq)+ ideal3 <- return $! (i3 `using` rdeepseq)+ ideal4 <- return $! (i4 `using` rdeepseq)+ defaultMain $ [bgroup "lex01"+ [ bench "simple" $ nf (simpleBuchbergerWith Lex) ideal1+ , bench "relprime" $ nf (primeTestBuchbergerWith Lex) ideal1+ , bench "syzygy" $ nf (syzygyBuchbergerWithStrategy NormalStrategy Lex) ideal1+ , bench "syz+sugar" $ nf (syzygyBuchbergerWith Lex) ideal1+ ]+ ,bgroup "grlex01"+ [ bench "simple" $ nf (simpleBuchbergerWith Grlex) ideal1+ , bench "relprime" $ nf (primeTestBuchbergerWith Grlex) ideal1+ , bench "syzygy" $ nf (syzygyBuchbergerWithStrategy NormalStrategy Grlex) ideal1+ , bench "syz+sugar" $ nf (syzygyBuchbergerWith Grlex) ideal1+ ]+ ,bgroup "grevlex01"+ [ bench "simple" $ nf (simpleBuchbergerWith Grevlex) ideal1+ , bench "relprime" $ nf (primeTestBuchbergerWith Grevlex) ideal1+ , bench "syzygy" $ nf (syzygyBuchbergerWithStrategy NormalStrategy Grevlex) ideal1+ , bench "syz+sugar" $ nf (syzygyBuchbergerWith Grevlex) ideal1+ ]+ ,bgroup "grlex02"+ [ bench "simple" $ nf (simpleBuchbergerWith Grlex) ideal2+ , bench "relprime" $ nf (primeTestBuchbergerWith Grlex) ideal2+ , bench "syzygy" $ nf (syzygyBuchbergerWithStrategy NormalStrategy Grlex) ideal2+ , bench "syz+sugar" $ nf (syzygyBuchbergerWith Grlex) ideal2+ ]+ ,bgroup "grevlex02"+ [ bench "simple" $ nf (simpleBuchbergerWith Grevlex) ideal2+ , bench "relprime" $ nf (primeTestBuchbergerWith Grevlex) ideal2+ , bench "syzygy" $ nf (syzygyBuchbergerWithStrategy NormalStrategy Grevlex) ideal2+ , bench "syz+sugar" $ nf (syzygyBuchbergerWith Grevlex) ideal2+ -- , bench "singular" $ nfIO (singularWith Grevlex ideal2)+ ]+ ,bgroup "lex03"+ [ bench "simple" $ nf (simpleBuchbergerWith Lex) ideal3+ , bench "relprime" $ nf (primeTestBuchbergerWith Lex) ideal4+ , bench "syzygy" $ nf (syzygyBuchbergerWithStrategy NormalStrategy Lex) ideal4+ , bench "syz+sugar" $ nf (syzygyBuchbergerWith Lex) ideal4+ -- , bench "singular" $ nfIO (singularWith Lex ideal3)+ ]+ ,bgroup "grlex03"+ [ bench "simple" $ nf (simpleBuchbergerWith Grlex) ideal3+ , bench "relprime" $ nf (primeTestBuchbergerWith Grlex) ideal3+ , bench "syzygy" $ nf (syzygyBuchbergerWithStrategy NormalStrategy Grlex) ideal3+ , bench "syz+sugar" $ nf (syzygyBuchbergerWith Grlex) ideal3+ -- , bench "singular" $ nfIO (singularWith Grlex ideal3)+ ]+ ,bgroup "grevlex03"+ [ bench "simple" $ nf (simpleBuchbergerWith Grevlex) ideal3+ , bench "relprime" $ nf (primeTestBuchbergerWith Grevlex) ideal3+ , bench "syzygy" $ nf (syzygyBuchbergerWithStrategy NormalStrategy Grevlex) ideal3+ , bench "syz+sugar" $ nf (syzygyBuchbergerWith Grevlex) ideal3+ ]+ ,bgroup "grlex04"+ [ bench "simple" $ nf (simpleBuchbergerWith Grlex) ideal4+ , bench "relprime" $ nf (primeTestBuchbergerWith Grlex) ideal4+ , bench "syzygy" $ nf (syzygyBuchbergerWithStrategy NormalStrategy Grlex) ideal+ , bench "syz+sugar" $ nf (syzygyBuchbergerWith Grlex) ideal4+ -- , bench "singular" $ nfIO (singularWith Grlex ideal4)+ ]+ ,bgroup "grevlex04"+ [ bench "simple" $ nf (simpleBuchbergerWith Grevlex) ideal4+ , bench "relprime" $ nf (primeTestBuchbergerWith Grevlex) ideal4+ , bench "syzygy" $ nf (syzygyBuchbergerWithStrategy NormalStrategy Grevlex) ideal4+ , bench "syz+sugar" $ nf (syzygyBuchbergerWith Grevlex) ideal4+ -- , bench "singular" $ nfIO (singularWith Grevlex ideal4)+ ]+ ]
+ examples/elimination-bench.hs view
@@ -0,0 +1,68 @@+{-# LANGUAGE DataKinds, FlexibleContexts, FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses, OverloadedStrings, PolyKinds #-}+{-# LANGUAGE TemplateHaskell, UndecidableInstances #-}+import Algebra.Algorithms.Groebner.Monomorphic+import Algebra.Internal+import Algebra.Ring.Polynomial (eliminationOrder,+ weightedEliminationOrder)+import Algebra.Ring.Polynomial.Monomorphic+import Control.DeepSeq+import Control.Parallel.Strategies+import Criterion.Main+import Criterion.Types+import Numeric.Algebra (LeftModule (..))+import qualified Numeric.Algebra as NA++x, y, z, w, s, a, b, c, t, u, v :: Polynomial Rational+[x, y, z, w, s, a, b, c, t, u, v] = map (injectVar . flip Variable Nothing) "xyzwSabctuv"++instance NFData Variable where+ rnf (Variable x y) = rnf x `seq` rnf y `seq` ()++instance NFData (Polynomial Rational) where+ rnf (Polynomial dic) = rnf dic++i1, i2, i3, i4, i5 :: [Polynomial Rational]+i1 = [x - (t + u), y - (t^2 + 2*t*u), z - (t^3 + 3*t^2*u)]+i2 = [t^2 + x^2+y^2+z^2, t^2 + 2*x^2 - x*y -z^2, t+ y^3-z^3]+i3 = [ 2 * s - a * y, b^2 - (x^2 + y^2), c^2 - ( (a-x) ^ 2 + y^2)+ ]+i4 = [ x - (3*u + 3*u*v^2 - u^3), y - (3*v + 3*u^2*v - v^3), z - (3*u^2 - 3*v^2)]+i5 = [t^2+x^2+y^2+z^2, t^2+2*x^2-x*y-z^2, t+y^3-z^3]+++main :: IO ()+main = do+ ideal1 <- return $! (i1 `using` rdeepseq)+ ideal2 <- return $! (i2 `using` rdeepseq)+ ideal3 <- return $! (i3 `using` rdeepseq)+ ideal4 <- return $! (i4 `using` rdeepseq)+ ideal5 <- return $! (i5 `using` rdeepseq)+ [var_x, var_y, var_t] <- return $! (map (flip Variable Nothing) "xyt" `using` rdeepseq)+ defaultMain $ [ bgroup "heron"+ [ bench "lex" $ nf (eliminateWith Lex [var_x, var_y]) ideal3+ , bench "product" $ nf (eliminateWith (eliminationOrder sTwo) [var_x, var_y]) ideal3+ , bench "weight" $ nf (eliminateWith (weightedEliminationOrder sTwo) [var_x, var_y]) ideal3+ ]+ , bgroup "implicit01"+ [ bench "lex" $ nf (eliminateWith Lex [var_t]) ideal5+ , bench "product" $ nf (eliminateWith (eliminationOrder sOne) [var_t]) ideal5+ , bench "weight" $ nf (eliminateWith (weightedEliminationOrder sOne) [var_t]) ideal5+ ]+ , bgroup "implicit02"+ [ bench "lex" $ nf (eliminateWith Lex [var_t]) ideal2+ , bench "product" $ nf (eliminateWith (eliminationOrder sOne) [var_t]) ideal2+ , bench "weight" $ nf (eliminateWith (weightedEliminationOrder sOne) [var_t]) ideal2+ ]+ , bgroup "implicit03"+ [ bench "lex" $ nf (eliminateWith Lex [var_t]) ideal1+ , bench "product" $ nf (eliminateWith (eliminationOrder sOne) [var_t]) ideal1+ , bench "weight" $ nf (eliminateWith (weightedEliminationOrder sOne) [var_t]) ideal1+ ]+ , bgroup "implicit04"+ [ bench "lex" $ nf (eliminateWith Lex [var_t]) ideal4+ , bench "product" $ nf (eliminateWith (eliminationOrder sOne) [var_t]) ideal4+ , bench "weight" $ nf (eliminateWith (weightedEliminationOrder sOne) [var_t]) ideal4+ ]+ ]+
examples/polymorphic.hs view
@@ -24,9 +24,6 @@ type LexPolynomial r n = OrderedPolynomial r Lex n -heron :: Ideal (LexPolynomial (Ratio Integer) (Two :+: Two))-heron = sTwo `thEliminationIdeal` heronIdeal- heronIdeal :: Ideal (Polynomial (Ratio Integer) (Three :+: Three)) heronIdeal = toIdeal [ 2 * s - a * y , b^^2 - (x^^2 + y^^2)
+ examples/sandpit-poly.hs view
@@ -0,0 +1,36 @@+{-# LANGUAGE DataKinds, OverlappingInstances, PolyKinds #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}+module Main (module Algebra.Algorithms.Groebner, module Algebra.Ring.Polynomial+ , module Data.Ratio, module Main, module Algebra.Internal+ ) where+import Algebra.Algorithms.Groebner+import Algebra.Internal+import Algebra.Ring.Noetherian+import Algebra.Ring.Polynomial+import Data.Ratio+import qualified Numeric.Algebra as NA++c, s, x, y :: Polynomial Rational (S Three)+[c, s, x, y] = genVars (sS sThree)++(.+), (.*), (.-) :: Polynomial Rational (S Three) -> Polynomial Rational (S Three) -> Polynomial Rational (S Three)+(.+) = (NA.+)+(.*) = (NA.*)+(.-) = (NA.-)++infixl 6 .+, .-+infixl 7 .*++(^^^) :: Polynomial Rational (S Three) -> NA.Natural -> Polynomial Rational (S Three)+(^^^) = NA.pow++fromRight :: Either t t1 -> t1+fromRight (Right a) = a++{-+parse :: String -> Polynomial Rational+parse = fromRight . parsePolyn+-}++main :: IO ()+main = print $ thEliminationIdealWith (weightedEliminationOrder sTwo) sTwo (toIdeal [x-c^3, y-s^3, c^2+s^2-1])
examples/sandpit.hs view
@@ -1,25 +1,25 @@ module Main (module Algebra.Algorithms.Groebner.Monomorphic, module Algebra.Ring.Polynomial- , module Algebra.Ring.Polynomial.Parser, module Algebra.Ring.Polynomial.Monomorphic- , module Data.Ratio, module Main, module Algebra.Internal) where-import Algebra.Algorithms.Groebner.Monomorphic-import Algebra.Ring.Polynomial (Grevlex (..),- Grlex (..), Lex (..),- ProductOrder (..),- WeightOrder (..),- WeightProxy (..),- eliminationOrder,- weightedEliminationOrder)-import Algebra.Ring.Polynomial.Monomorphic-import Algebra.Ring.Polynomial.Parser-import Data.Ratio+ {- , module Algebra.Ring.Polynomial.Parser-} , module Algebra.Ring.Polynomial.Monomorphic+ , module Data.Ratio, module Main, module Algebra.Internal+ ) where+import Algebra.Algorithms.Groebner.Monomorphic import Algebra.Internal-import qualified Numeric.Algebra as NA+import Algebra.Ring.Polynomial (Grevlex (..), Grlex (..),+ Lex (..), ProductOrder (..),+ WeightOrder (..),+ WeightProxy (..),+ eliminationOrder,+ weightedEliminationOrder)+import Algebra.Ring.Polynomial.Monomorphic+-- import Algebra.Ring.Polynomial.Parser+import Data.Ratio+import qualified Numeric.Algebra as NA var_x, var_y, var_z, var_t, var_u :: Variable-[var_x, var_y, var_z, var_t, var_u] = map (flip Variable Nothing) "xyztu"+[var_c, var_s, var_x, var_y, var_z, var_t, var_u] = map (flip Variable Nothing) "csxyztu" x, y, z, t, u :: Polynomial Rational-[x, y, z, t, u] = map injectVar [var_x, var_y, var_z, var_t, var_u]+[c, s, x, y, z, t, u] = map injectVar [var_c, var_s, var_x, var_y, var_z, var_t, var_u] (.+), (.*), (.-) :: Polynomial Rational -> Polynomial Rational -> Polynomial Rational (.+) = (NA.+)@@ -35,8 +35,10 @@ fromRight :: Either t t1 -> t1 fromRight (Right a) = a +{- parse :: String -> Polynomial Rational parse = fromRight . parsePolyn+-} main :: IO ()-main = return ()+main = print $ eliminate [var_s, var_c] [x-c^3, y-s^3, c^2+s^2-1]
+ examples/sugar-bench.hs view
@@ -0,0 +1,54 @@+{-# LANGUAGE DataKinds, FlexibleContexts, FlexibleInstances, GADTs #-}+{-# LANGUAGE MultiParamTypeClasses, OverloadedStrings, PolyKinds #-}+{-# LANGUAGE QuasiQuotes, TemplateHaskell, UndecidableInstances #-}+{-# OPTIONS_GHC -fno-warn-type-defaults -fno-warn-orphans #-}+import Algebra.Algorithms.Groebner.Monomorphic+import Algebra.Ring.Polynomial.Monomorphic+import Control.DeepSeq+import Control.Parallel.Strategies+import Criterion.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++i1, i2, i3, i4 :: [Polynomial Rational]+i1 = [x^2 + y^2 + z^2 - 1, x^2 + y^2 + z^2 - 2*x, 2*x -3*y - z]+i2 = [x^2 * y - 2*x*y - 4*z - 1, z-y^2, x^3 - 4*z*y]+i3 = [ 2 * s - a * y, b^2 - (x^2 + y^2), c^2 - ( (a-x) ^ 2 + y^2)+ ]+i4 = [ z^5 + y^4 + x^3 - 1, z^3 + y^3 + x^2 - 1]++mkTestCases :: (Eq r, Show a, r ~ Rational) => a -> [Polynomial r] -> [Benchmark]+mkTestCases num ideal = [ mkTC ("lex0" ++ show num) ideal Lex+ , mkTC ("grlex0" ++ show num) ideal Grlex+ , mkTC ("grevlex0" ++ show num) ideal Grevlex+ ]++mkTC :: (r ~ Rational, IsMonomialOrder ord) => String -> [Polynomial r] -> ord -> Benchmark+mkTC name ideal ord =+ bgroup name [ bench "syzygy" $ nf (syzygyBuchbergerWithStrategy NormalStrategy ord) ideal+ , bench "syz+grevsel" $ nf (syzygyBuchbergerWithStrategy GrevlexStrategy ord) ideal+ , bench "syz+grad" $ nf (syzygyBuchbergerWithStrategy GradedStrategy ord) ideal+ , bench "syz+sugar" $ nf (syzygyBuchbergerWithStrategy (SugarStrategy NormalStrategy) ord) ideal+ , bench "syz+grsugar" $ nf (syzygyBuchbergerWithStrategy (SugarStrategy GradedStrategy) ord) ideal+ ]++main :: IO ()+main = do+ ideal1 <- return $! (i1 `using` rdeepseq)+ ideal2 <- return $! (i2 `using` rdeepseq)+ ideal3 <- return $! (i3 `using` rdeepseq)+ ideal4 <- return $! (i4 `using` rdeepseq)+ defaultMain $ concat $+ [ mkTestCases 1 ideal1+ , [mkTC "grlex02" ideal2 Grlex, mkTC "grevlex02" ideal2 Grevlex]+ , mkTestCases 3 ideal3+ , [mkTC "grlex04" ideal4 Grlex, mkTC "grevlex04" ideal4 Grevlex]+ ]+