module Algebra.Algorithms.Groebner
(
isGroebnerBasis
, calcGroebnerBasis, calcGroebnerBasisWith
, calcGroebnerBasisWithStrategy
, buchberger, syzygyBuchberger
, simpleBuchberger, primeTestBuchberger
, reduceMinimalGroebnerBasis, minimizeGroebnerBasis
, syzygyBuchbergerWithStrategy
, SelectionStrategy(..), calcWeight', GrevlexStrategy(..)
, NormalStrategy(..), SugarStrategy(..), GradedStrategy(..)
, isIdealMember, intersection, thEliminationIdeal, thEliminationIdealWith
, unsafeThEliminationIdealWith
, quotIdeal, quotByPrincipalIdeal
, saturationIdeal, saturationByPrincipalIdeal
, resultant, hasCommonFactor
, lcmPolynomial, gcdPolynomial
) where
import Algebra.Internal
import Algebra.Prelude.Core
import Algebra.Ring.Polynomial.Univariate (Unipol)
import Control.Lens ((%~), (&), _Wrapped)
import Control.Monad.Loops (whileM_)
import Control.Monad.ST (ST, runST)
import qualified Data.Foldable as H
import qualified Data.Heap as H
import qualified Data.Map as M
import Data.Singletons.Prelude (POrd (..), SEq (..))
import Data.Singletons.Prelude (Sing (SFalse, STrue), withSingI)
import Data.Singletons.Prelude.List (Length, Replicate, Sing (SCons))
import Data.Singletons.Prelude.List (sLength, sReplicate)
import Data.Sized.Builtin (toList)
import qualified Data.Sized.Builtin as V
import Data.STRef (STRef, modifySTRef, newSTRef)
import Data.STRef (readSTRef, writeSTRef)
import qualified Prelude as P
import Proof.Equational
isGroebnerBasis :: (IsOrderedPolynomial poly, Field (Coefficient poly))
=> Ideal poly -> Bool
isGroebnerBasis (nub . generators -> ideal) = all check $ combinations ideal
where
check (f, g) =
let (t, u) = (leadingMonomial f , leadingMonomial g)
in t*u == lcmMonomial t u || sPolynomial f g `modPolynomial` ideal == zero
simpleBuchberger :: (Field (Coefficient poly), IsOrderedPolynomial poly)
=> Ideal poly -> [poly]
simpleBuchberger 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
, let q = sPolynomial f g `modPolynomial` acc, q /= zero
]
primeTestBuchberger :: (Field (Coefficient poly), IsOrderedPolynomial poly)
=> Ideal poly -> [poly]
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 /= 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
buchberger :: (Field (Coefficient poly), IsOrderedPolynomial poly)
=> Ideal poly -> [poly]
buchberger = syzygyBuchberger
syzygyBuchberger :: (Field (Coefficient poly), IsOrderedPolynomial poly)
=> Ideal poly -> [poly]
syzygyBuchberger = syzygyBuchbergerWithStrategy (SugarStrategy NormalStrategy)
syzygyBuchbergerWithStrategy :: (Field (Coefficient poly), IsOrderedPolynomial poly,
SelectionStrategy (Arity poly) strategy,
Ord (Weight (Arity poly) strategy (MOrder poly)))
=> strategy -> Ideal poly -> [poly]
syzygyBuchbergerWithStrategy strategy ideal = runST $ do
let gens = zip [1..] $ filter (/= zero) $ generators ideal
gs <- newSTRef $ H.fromList [H.Entry (leadingMonomial 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
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.all ((/=k) . H.payload) rest)
[(f, h), (g, h), (h, f), (h, g)])
&& leadingMonomial h `divs` l) gs0
when (l /= f0 * g0 && not redundant) $ do
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 (leadingMonomial s) s)
len %= (*2)
map H.payload . H.toList <$> readSTRef gs
calcWeight' :: (SelectionStrategy (Arity poly) s, IsOrderedPolynomial poly)
=> s -> poly -> poly -> Weight (Arity poly) s (MOrder poly)
calcWeight' s = calcWeight (toProxy s)
class SelectionStrategy n s where
type Weight n s ord :: *
calcWeight :: (IsOrderedPolynomial poly, n ~ Arity poly)
=> Proxy s -> poly -> poly -> Weight n s (MOrder poly)
data NormalStrategy = NormalStrategy deriving (Read, Show, Eq, Ord)
instance SelectionStrategy n NormalStrategy where
type Weight n NormalStrategy ord = OrderedMonomial ord n
calcWeight _ f g = lcmMonomial (leadingMonomial f) (leadingMonomial g)
data GrevlexStrategy = GrevlexStrategy deriving (Read, Show, Eq, Ord)
instance SelectionStrategy n GrevlexStrategy where
type Weight n GrevlexStrategy ord = OrderedMonomial Grevlex n
calcWeight _ f g = changeMonomialOrderProxy Proxy $
lcmMonomial (leadingMonomial f) (leadingMonomial g)
data GradedStrategy = GradedStrategy deriving (Read, Show, Eq, Ord)
instance SelectionStrategy n GradedStrategy where
type Weight n GradedStrategy ord = OrderedMonomial (Graded ord) n
calcWeight _ f g = changeMonomialOrderProxy Proxy $
lcmMonomial (leadingMonomial f) (leadingMonomial g)
data SugarStrategy s = SugarStrategy s deriving (Read, Show, Eq, Ord)
instance SelectionStrategy n s => SelectionStrategy n (SugarStrategy s) where
type Weight n (SugarStrategy s) ord = (Int, Weight n s ord)
calcWeight (Proxy :: Proxy (SugarStrategy s)) f g = (sugar, calcWeight (Proxy :: Proxy s) f g)
where
deg' = maximum . map totalDegree . H.toList . orderedMonomials
tsgr h = deg' h totalDegree (leadingMonomial h)
sugar = max (tsgr f) (tsgr g) + totalDegree (lcmMonomial (leadingMonomial f) (leadingMonomial g))
minimizeGroebnerBasis :: (Field (Coefficient poly), IsOrderedPolynomial poly)
=> [poly] -> [poly]
minimizeGroebnerBasis bs = runST $ do
left <- newSTRef $ map monoize $ filter (/= zero) bs
right <- newSTRef []
whileM_ (not . null <$> readSTRef left) $ do
f : xs <- readSTRef left
writeSTRef left xs
ys <- readSTRef right
unless (any (\g -> leadingMonomial g `divs` leadingMonomial f) xs
|| any (\g -> leadingMonomial g `divs` leadingMonomial f) ys) $
writeSTRef right (f : ys)
readSTRef right
reduceMinimalGroebnerBasis :: (Field (Coefficient poly), IsOrderedPolynomial poly)
=> [poly] -> [poly]
reduceMinimalGroebnerBasis bs = runST $ do
left <- newSTRef bs
right <- newSTRef []
whileM_ (not . null <$> readSTRef left) $ do
f : xs <- readSTRef left
writeSTRef left xs
ys <- readSTRef right
let q = f `modPolynomial` (xs ++ ys)
if q == zero then writeSTRef right ys else writeSTRef right (q : ys)
readSTRef right
calcGroebnerBasisWith :: (IsOrderedPolynomial poly,
Field (Coefficient poly),
IsMonomialOrder (Arity poly) order)
=> order -> Ideal poly
-> [OrderedPolynomial (Coefficient poly) order (Arity poly)]
calcGroebnerBasisWith _ord = calcGroebnerBasis . mapIdeal injectVars
calcGroebnerBasisWithStrategy :: (Field (Coefficient poly), IsOrderedPolynomial poly
, SelectionStrategy (Arity poly) strategy
, Ord (Weight (Arity poly) strategy (MOrder poly)))
=> strategy -> Ideal poly -> [poly]
calcGroebnerBasisWithStrategy strategy =
reduceMinimalGroebnerBasis . minimizeGroebnerBasis . syzygyBuchbergerWithStrategy strategy
calcGroebnerBasis :: (Field (Coefficient poly), IsOrderedPolynomial poly)
=> Ideal poly -> [poly]
calcGroebnerBasis = reduceMinimalGroebnerBasis . minimizeGroebnerBasis . syzygyBuchberger
isIdealMember :: (Field (Coefficient poly), IsOrderedPolynomial poly)
=> poly -> Ideal poly -> Bool
isIdealMember f ideal = groebnerTest f (calcGroebnerBasis ideal)
groebnerTest :: (Field (Coefficient poly), IsOrderedPolynomial poly)
=> poly -> [poly] -> Bool
groebnerTest f fs = f `modPolynomial` fs == zero
newtype LengthReplicate n =
LengthReplicate { runLengthReplicate :: forall x. Sing (x :: Nat)
-> Length (Replicate n x) :~: n }
lengthReplicate :: SNat n -> SNat x -> Length (Replicate n x) :~: n
lengthReplicate = runLengthReplicate . induction base step
where
base :: LengthReplicate 0
base = LengthReplicate $ const Refl
step :: SNat n -> LengthReplicate n -> LengthReplicate (Succ n)
step n (LengthReplicate ih) = LengthReplicate $ \x ->
case (n %:+ sOne) %:== sZero of
SFalse ->
start (sLength (sReplicate (sSucc n) x))
=~= sLength (SCons x (sReplicate (sSucc n %:- sOne) x))
=~= sOne %:+ sLength (sReplicate (sSucc n %:- sOne) x)
=== sSucc (sLength (sReplicate (sSucc n %:- sOne) x))
`because` sym (succAndPlusOneL (sLength (sReplicate (sSucc n %:- sOne) x)))
=== sSucc (sLength (sReplicate (n %:+ sOne %:- sOne) x))
`because` succCong (lengthCong (replicateCong (minusCongL (succAndPlusOneR n) sOne) x))
=== sSucc (sLength (sReplicate n x))
`because` succCong (lengthCong (replicateCong (plusMinus n sOne) x))
=== sSucc n `because` succCong (ih x)
STrue -> case sCompare (n %:+ sOne) sZero of {}
lengthCong :: a :~: b -> Length a :~: Length b
lengthCong Refl = Refl
replicateCong :: a :~: b -> Sing x -> Replicate a x :~: Replicate b x
replicateCong Refl _ = Refl
thEliminationIdeal :: forall poly n.
( IsMonomialOrder (Arity poly n) (MOrder poly),
Field (Coefficient poly),
IsOrderedPolynomial poly,
(n :<= Arity poly) ~ 'True)
=> SNat n
-> Ideal poly
-> Ideal (OrderedPolynomial (Coefficient poly) (MOrder poly) (Arity poly :-. n))
thEliminationIdeal n = withSingI (sOnes n) $
withRefl (lengthReplicate n sOne) $
withKnownNat n $
withKnownNat ((sing :: SNat (Arity poly)) %:-. n) $
mapIdeal (changeOrderProxy Proxy) . thEliminationIdealWith (weightedEliminationOrder n) n
thEliminationIdealWith :: ( IsOrderedPolynomial poly,
m ~ Arity poly,
k ~ Coefficient poly, Field k,
KnownNat (m :-. n), (n :<= m) ~ 'True,
EliminationType m n ord)
=> ord
-> SNat n
-> Ideal poly
-> Ideal (OrderedPolynomial k Grevlex (m :-. n))
thEliminationIdealWith = unsafeThEliminationIdealWith
unsafeThEliminationIdealWith :: ( IsOrderedPolynomial poly,
m ~ Arity poly,
k ~ Coefficient poly,
Field k,
IsMonomialOrder m ord,
KnownNat (m :-. n), (n :<= m) ~ 'True)
=> ord
-> SNat n
-> Ideal poly
-> Ideal (OrderedPolynomial k Grevlex (m :-. n))
unsafeThEliminationIdealWith ord n ideal =
withKnownNat n $ toIdeal $ [ f & _Wrapped %~ M.mapKeys (orderMonomial Nothing . V.drop n . getMonomial)
| f <- calcGroebnerBasisWith ord ideal
, all (all (== 0) . V.takeAtMost n . getMonomial . snd) $ getTerms f
]
eliminatePadding :: (IsOrderedPolynomial poly,
IsMonomialOrder n ord,
Field (Coefficient poly),
SingI (Replicate n 1),
KnownNat n
)
=> Ideal (PadPolyL n ord poly) -> Ideal poly
eliminatePadding ideal =
toIdeal $ [ c
| f0 <- calcGroebnerBasis ideal
, let (c, m) = leadingTerm $ runPadPolyL f0
, m == one
]
intersection :: forall poly k.
( Field (Coefficient poly), IsOrderedPolynomial poly)
=> Sized k (Ideal poly)
-> Ideal poly
intersection idsv@(_ :< _) =
let sk = sizedLength idsv
in withSingI (sOnes sk) $ withKnownNat sk $
let ts = take (fromIntegral $ fromSing sk) vars
inj = padLeftPoly sk Grevlex
tis = zipWith (\ideal t -> mapIdeal ((t *) . inj) ideal) (toList idsv) ts
j = foldr appendIdeal (principalIdeal (one foldr (+) zero ts)) tis
in eliminatePadding j
intersection _ = Ideal $ singleton one
quotByPrincipalIdeal :: (Field (Coefficient poly), IsOrderedPolynomial poly)
=> Ideal poly
-> poly
-> Ideal poly
quotByPrincipalIdeal i g =
case intersection (i :< (Ideal $ singleton g) :< NilL) of
Ideal gs -> Ideal $ V.map (snd . head . (`divPolynomial` [g])) gs
quotIdeal :: forall poly l.
(IsOrderedPolynomial poly, Field (Coefficient poly))
=> Ideal poly
-> Sized l poly
-> Ideal poly
quotIdeal i g =
withKnownNat (sizedLength g) $
intersection $ V.map (i `quotByPrincipalIdeal`) g
saturationByPrincipalIdeal :: forall poly.
(IsOrderedPolynomial poly, Field (Coefficient poly))
=> Ideal poly
-> poly
-> Ideal poly
saturationByPrincipalIdeal is g =
let n = sArity' g
in withRefl (plusMinus' sOne n) $ withRefl (plusComm n sOne) $
withWitness (leqStep sOne (sOne %:+ n) n Refl) $
withWitness (lneqZero n) $
eliminatePadding $
addToIdeal (one (padLeftPoly sOne Grevlex g * var 0)) $
mapIdeal (padLeftPoly sOne Grevlex) is
saturationIdeal :: forall poly l.
(Field (Coefficient poly),
IsOrderedPolynomial poly)
=> Ideal poly
-> Sized l poly
-> Ideal poly
saturationIdeal i g =
withKnownNat (sizedLength g) $
intersection $ V.map (i `saturationByPrincipalIdeal`) g
resultant :: forall poly.
(Field (Coefficient poly),
IsOrderedPolynomial poly,
Arity poly ~ 1)
=> poly
-> poly
-> (Coefficient poly)
resultant = go one
where
go res h s
| totalDegree' s > 0 =
let r = h `modPolynomial` [s]
res' = res * negate one ^ (totalDegree' h * totalDegree' s)
* (leadingCoeff s) ^ (totalDegree' h P.- totalDegree' r)
in go res' s r
| isZero h || isZero s = zero
| totalDegree' h > 0 = (leadingCoeff s ^ totalDegree' h) * res
| otherwise = res
_ = Refl :: Arity poly :~: 1
hasCommonFactor :: (Field (Coefficient poly),
IsOrderedPolynomial poly,
Arity poly ~ 1)
=> poly
-> poly
-> Bool
hasCommonFactor f g = isZero $ resultant f g
lcmPolynomial :: forall poly.
(Field (Coefficient poly),
IsOrderedPolynomial poly)
=> poly
-> poly
-> poly
lcmPolynomial f g = head $ generators $ intersection (principalIdeal f :< principalIdeal g :< NilL)
gcdPolynomial :: (Field (Coefficient poly),
IsOrderedPolynomial poly)
=> poly
-> poly
-> poly
gcdPolynomial f g = snd $ head $ f * g `divPolynomial` [lcmPolynomial f g]