module Data.Polynomial.GBasis
(
Options (..)
, Strategy (..)
, defaultOptions
, basis
, basis'
, spolynomial
, reduceGBasis
) where
import qualified Data.Set as Set
import qualified Data.Heap as H
import Data.Polynomial
data Options
= Options
{ optStrategy :: Strategy
}
defaultOptions :: Options
defaultOptions =
Options
{ optStrategy = NormalStrategy
}
data Strategy
= NormalStrategy
| SugarStrategy
deriving (Eq, Ord, Show, Read, Bounded, Enum)
spolynomial
:: (Eq k, Fractional k, Ord v)
=> MonomialOrder v -> Polynomial k v -> Polynomial k v -> Polynomial k v
spolynomial cmp f g =
fromMonomial ((1,xs) `monomialDiv` (c1,xs1)) * f
fromMonomial ((1,xs) `monomialDiv` (c2,xs2)) * g
where
xs = mmLCM xs1 xs2
(c1, xs1) = leadingTerm cmp f
(c2, xs2) = leadingTerm cmp g
basis
:: forall k v. (Eq k, Fractional k, Ord k, Ord v)
=> MonomialOrder v
-> [Polynomial k v]
-> [Polynomial k v]
basis = basis' defaultOptions
basis'
:: forall k v. (Eq k, Fractional k, Ord k, Ord v)
=> Options
-> MonomialOrder v
-> [Polynomial k v]
-> [Polynomial k v]
basis' opt cmp fs =
reduceGBasis cmp $ go fs (H.fromList [item cmp fi fj | (fi,fj) <- pairs fs, checkGCD fi fj])
where
go :: [Polynomial k v] -> H.Heap (Item k v) -> [Polynomial k v]
go gs h | H.null h = gs
go gs h
| r == 0 = go gs h'
| otherwise = go (r:gs) (H.union h' (H.fromList [item cmp r g | g <- gs, checkGCD fi fj]))
where
Just (i, h') = H.viewMin h
fi = iFst i
fj = iSnd i
spoly = spolynomial cmp fi fj
r = reduce cmp spoly gs
checkGCD fi fj = mmGCD mm1 mm2 /= mmOne
where
(_, mm1) = leadingTerm cmp fi
(_, mm2) = leadingTerm cmp fj
reduceGBasis
:: forall k v. (Eq k, Ord k, Fractional k, Ord v)
=> MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]
reduceGBasis cmp ps = Set.toList $ Set.fromList $ go ps []
where
go [] qs = qs
go (p:ps) qs
| q == 0 = go ps qs
| otherwise = go ps (constant (1/c) * q : qs)
where
q = reduce cmp p (ps++qs)
(c,_) = leadingTerm cmp q
data Item k v
= Item
{ iFst :: Polynomial k v
, iSnd :: Polynomial k v
, iCmp :: MonomialOrder v
, iLCM :: MonicMonomial v
}
item :: (Eq k, Num k, Ord v) => MonomialOrder v -> Polynomial k v -> Polynomial k v -> Item k v
item cmp f g = Item f g cmp (mmLCM mm1 mm2)
where
(_, mm1) = leadingTerm cmp f
(_, mm2) = leadingTerm cmp g
instance Ord v => Ord (Item k v) where
a `compare` b = iCmp a (iLCM a) (iLCM b)
instance Ord v => Eq (Item k v) where
a == b = compare a b == EQ
pairs :: [a] -> [(a,a)]
pairs [] = []
pairs (x:xs) = [(x,y) | y <- xs] ++ pairs xs