{-# LANGUAGE ConstraintKinds, DataKinds, EmptyCase, FlexibleContexts        #-}
{-# LANGUAGE FlexibleInstances, GADTs, MultiParamTypeClasses                #-}
{-# LANGUAGE NoImplicitPrelude, ParallelListComp, PolyKinds, RankNTypes     #-}
{-# LANGUAGE ScopedTypeVariables, TypeFamilies, TypeOperators, ViewPatterns #-}
{-# OPTIONS_GHC -fno-warn-type-defaults -fno-warn-orphans #-}
module Algebra.Algorithms.Groebner
       (
       -- * Groebner basis
         isGroebnerBasis
       , 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, eliminatePadding
       , quotIdeal, quotByPrincipalIdeal
       , saturationIdeal, saturationByPrincipalIdeal
       -- * Resultant
       , 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

-- | Test if the given ideal is Groebner basis, using Buchberger criteria and relatively primeness.
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

-- | The Naive buchberger's algorithm to calculate Groebner basis for the given ideal.
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
               ]

-- | Buchberger's algorithm slightly improved by discarding relatively prime pair.
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
{-# INLINE combinations #-}

-- | Calculate Groebner basis applying (modified) Buchberger's algorithm.
-- This function is same as 'syzygyBuchberger'.
buchberger :: (Field (Coefficient poly), IsOrderedPolynomial poly)
           => Ideal poly -> [poly]
buchberger = syzygyBuchberger

-- | 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 (Coefficient poly), IsOrderedPolynomial poly)
                    => Ideal poly -> [poly]
syzygyBuchberger = syzygyBuchbergerWithStrategy (SugarStrategy NormalStrategy)
{-# SPECIALISE INLINE [0]
    syzygyBuchberger :: (CoeffRing r, Field r, IsMonomialOrder n ord, KnownNat n)
                     => Ideal (OrderedPolynomial r ord n) -> [OrderedPolynomial r ord n]
 #-}
{-# SPECIALISE INLINE [0]
    syzygyBuchberger :: (CoeffRing r, Field r)
                     => Ideal (Unipol r) -> [Unipol r]
 #-}
{-# INLINE [1] syzygyBuchberger #-}

-- | apply buchberger's algorithm using given selection strategy.
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
{-# SPECIALISE INLINE [0]
 syzygyBuchbergerWithStrategy :: (Field k, CoeffRing k, KnownNat n)
                    => SugarStrategy NormalStrategy -> Ideal (OrderedPolynomial k Grevlex n) -> [OrderedPolynomial k Grevlex n]
 #-}
{-# SPECIALISE INLINE [0]
 syzygyBuchbergerWithStrategy :: (Field k, CoeffRing k)
                    => SugarStrategy NormalStrategy -> Ideal (Unipol k) -> [Unipol k]
 #-}

{-# SPECIALISE INLINE [1]
 syzygyBuchbergerWithStrategy :: (Field k, CoeffRing k, KnownNat n, IsMonomialOrder n ord)
                    => SugarStrategy NormalStrategy -> Ideal (OrderedPolynomial k ord n) -> [OrderedPolynomial k ord n]
 #-}
{-# SPECIALISE INLINE [1]
 syzygyBuchbergerWithStrategy :: (Field k, CoeffRing k, IsMonomialOrder n ord,
                                 SelectionStrategy n strategy, KnownNat n,
                                 Ord (Weight n strategy ord))
                    => strategy -> Ideal (OrderedPolynomial k ord n) -> [OrderedPolynomial k ord n]
 #-}
{-# INLINABLE [2] syzygyBuchbergerWithStrategy #-}


-- | 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 (Arity poly) s, IsOrderedPolynomial poly)
            => s -> poly -> poly -> Weight (Arity poly) s (MOrder poly)
calcWeight' s = calcWeight (toProxy s)
{-# INLINE calcWeight' #-}

-- | Type-class for selection strategies in Buchberger's algorithm.
class SelectionStrategy n s where
  type Weight n s ord :: *
  -- | Calculates the weight for the given pair of polynomial used for selection strategy.
  calcWeight :: (IsOrderedPolynomial poly, n ~ Arity poly)
             => Proxy s -> poly -> poly -> Weight n s (MOrder poly)

-- | 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 n NormalStrategy where
  type Weight n NormalStrategy ord = OrderedMonomial ord n
  calcWeight _ f g = lcmMonomial (leadingMonomial f)  (leadingMonomial g)
  {-# INLINE calcWeight #-}

-- | 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 n GrevlexStrategy where
  type Weight n GrevlexStrategy ord = OrderedMonomial Grevlex n
  calcWeight _ f g = changeMonomialOrderProxy Proxy $
                     lcmMonomial (leadingMonomial f) (leadingMonomial g)
  {-# INLINE calcWeight #-}

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 n GradedStrategy where
  type Weight n GradedStrategy ord = OrderedMonomial (Graded ord) n
  calcWeight _ f g = changeMonomialOrderProxy Proxy $
                     lcmMonomial (leadingMonomial f)  (leadingMonomial g)
  {-# INLINE calcWeight #-}


-- | 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 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))
  {-# INLINE calcWeight #-}


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

-- | Reduce minimum Groebner basis into reduced Groebner basis.
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

-- | Caliculating reduced Groebner basis of the given ideal w.r.t. the specified monomial order.
calcGroebnerBasisWith :: (IsOrderedPolynomial poly,
                          Field (Coefficient poly),
                          IsMonomialOrder (Arity poly) order)
                      => order -> Ideal poly
                      -> [OrderedPolynomial (Coefficient poly) order (Arity poly)]
calcGroebnerBasisWith _ord = calcGroebnerBasis . mapIdeal injectVars
{-# INLINE [1] calcGroebnerBasisWith #-}
{-# RULES
"calcGroebnerBasisWith/sameOrderPolyn" [~1] forall x.
  calcGroebnerBasisWith x = calcGroebnerBasis
  #-}

-- | Caliculating reduced Groebner basis of the given ideal w.r.t. the specified monomial order.
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

-- | Caliculating reduced Groebner basis of the given ideal.
calcGroebnerBasis :: (Field (Coefficient poly), IsOrderedPolynomial poly)
                  => Ideal poly -> [poly]
calcGroebnerBasis = reduceMinimalGroebnerBasis . minimizeGroebnerBasis . syzygyBuchberger
{-# SPECIALISE INLINE [0]
    calcGroebnerBasis :: (CoeffRing r, Field r, IsMonomialOrder n ord, KnownNat n)
                      => Ideal (OrderedPolynomial r ord n) -> [OrderedPolynomial r ord n]
 #-}
{-# SPECIALISE INLINE [0]
    calcGroebnerBasis :: (CoeffRing r, Field r)
                      => Ideal (Unipol r) -> [Unipol r]
 #-}
{-# INLINE [0] calcGroebnerBasis #-}


-- | Test if the given polynomial is the member of the ideal.
isIdealMember :: (Field (Coefficient poly), IsOrderedPolynomial poly)
              => poly -> Ideal poly -> Bool
isIdealMember f ideal = groebnerTest f (calcGroebnerBasis ideal)

-- | Test if the given polynomial can be divided by the given polynomials.
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

-- | Calculate n-th elimination ideal using 'WeightedEliminationOrder' ordering.
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

-- | Calculate n-th elimination ideal using the specified n-th elimination type order.
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

-- | Calculate n-th elimination ideal using the specified n-th elimination type order.
-- This function should be used carefully because it does not check whether the given ordering is
-- n-th elimintion type or not.
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
            ]

-- | An intersection ideal of given ideals (using 'WeightedEliminationOrder').
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 withRefl (plusMinus' sk sn) $
    --    withWitness (plusLeqL sk sn) $
    --    mapIdeal injectVars $
    --    coerce (cong Proxy $ minusCongL (plusComm sk sn) sk `trans` plusMinus sn sk) $
    --    thEliminationIdeal sk j
    in eliminatePadding j
intersection _ = Ideal $ singleton one

-- | Ideal quotient by a principal ideals.
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

-- | Ideal quotient by the given ideal.
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

-- | Saturation by a principal ideal.
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

-- | Saturation ideal
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

-- | Calculate resultant for given two unary polynomimals.
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
        -- to suppress "redundant" warning for univariate constraint.

-- | Determine whether two polynomials have a common factor with positive degree using resultant.
hasCommonFactor :: (Field (Coefficient poly),
                    IsOrderedPolynomial poly,
                    Arity poly ~ 1)
                => poly
                -> poly
                -> Bool
hasCommonFactor f g = isZero $ resultant f g

-- | Calculates the Least Common Multiply of the given pair of polynomials.
lcmPolynomial :: forall poly.
                 (Field (Coefficient poly),
                  IsOrderedPolynomial poly)
              => poly
              -> poly
              -> poly
lcmPolynomial f g = head $ generators $ intersection (principalIdeal f :< principalIdeal g :< NilL)

-- | Calculates the Greatest Common Divisor of the given pair of polynomials.
gcdPolynomial :: (Field (Coefficient poly),
                  IsOrderedPolynomial poly)
              => poly
              -> poly
              -> poly
gcdPolynomial f g = snd $ head $ f * g `divPolynomial` [lcmPolynomial f g]