diff --git a/Algebra/Algorithms/Groebner.hs b/Algebra/Algorithms/Groebner.hs --- a/Algebra/Algorithms/Groebner.hs +++ b/Algebra/Algorithms/Groebner.hs @@ -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 -> diff --git a/Algebra/Algorithms/Groebner/Monomorphic.hs b/Algebra/Algorithms/Groebner/Monomorphic.hs --- a/Algebra/Algorithms/Groebner/Monomorphic.hs +++ b/Algebra/Algorithms/Groebner/Monomorphic.hs @@ -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] diff --git a/Algebra/Internal.hs b/Algebra/Internal.hs --- a/Algebra/Internal.hs +++ b/Algebra/Internal.hs @@ -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 diff --git a/Algebra/Ring/Polynomial.hs b/Algebra/Ring/Polynomial.hs --- a/Algebra/Ring/Polynomial.hs +++ b/Algebra/Ring/Polynomial.hs @@ -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 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.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 diff --git a/examples/bench.hs b/examples/bench.hs --- a/examples/bench.hs +++ b/examples/bench.hs @@ -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) + ] + ] diff --git a/examples/elimination-bench.hs b/examples/elimination-bench.hs new file mode 100644 --- /dev/null +++ b/examples/elimination-bench.hs @@ -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 + ] + ] + diff --git a/examples/polymorphic.hs b/examples/polymorphic.hs --- a/examples/polymorphic.hs +++ b/examples/polymorphic.hs @@ -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) diff --git a/examples/sandpit-poly.hs b/examples/sandpit-poly.hs new file mode 100644 --- /dev/null +++ b/examples/sandpit-poly.hs @@ -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]) diff --git a/examples/sandpit.hs b/examples/sandpit.hs --- a/examples/sandpit.hs +++ b/examples/sandpit.hs @@ -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] diff --git a/examples/sugar-bench.hs b/examples/sugar-bench.hs new file mode 100644 --- /dev/null +++ b/examples/sugar-bench.hs @@ -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] + ] +