module Algebra.Ring.Polynomial.Monomial
( Monomial, OrderedMonomial(..),
IsOrder(..), IsMonomialOrder, MonomialOrder,
IsStrongMonomialOrder,
isRelativelyPrime, totalDegree, ProductOrder(..),
productOrder, productOrder', WeightProxy, WeightOrder(..),
gcdMonomial, divs, isPowerOf, tryDiv, lcmMonomial,
Lex(..), EliminationType, EliminationOrder,
WeightedEliminationOrder, eliminationOrder, weightedEliminationOrder,
lex, revlex, graded, grlex, grevlex,
weightOrder, Grevlex(..), fromList,
Revlex(..), Grlex(..), Graded(..),
castMonomial, scastMonomial, varMonom,
changeMonomialOrder, changeMonomialOrderProxy, sOnes,
withStrongMonomialOrder, cmpAnyMonomial, orderMonomial
) where
import Algebra.Internal
import AlgebraicPrelude hiding (lex)
import Control.DeepSeq (NFData (..))
import Control.Lens (Ixed (..), alaf, imap,
makeLenses, makeWrapped, (%~),
(&), (.~), _Wrapped)
import Data.Constraint ((:=>) (..), Dict (..))
import qualified Data.Constraint as C
import Data.Constraint.Forall
import qualified Data.Foldable as F
import Data.Hashable (Hashable (..))
import Data.Kind (Type)
import Data.Maybe (catMaybes)
import Data.Monoid (Dual (..))
import Data.Monoid ((<>))
import qualified Data.MonoTraversable.Unprefixed as MT
import Data.Ord (comparing)
import Data.Singletons.Prelude (POrd (..), SList, Sing ())
import Data.Singletons.Prelude (SingKind (..))
import Data.Singletons.Prelude.List (Length, Replicate, sReplicate)
import Data.Singletons.TypeLits (withKnownNat)
import qualified Data.Sized.Builtin as V
import Data.Type.Natural.Class (IsPeano (..), PeanoOrder (..))
import Data.Type.Ordinal (Ordinal (..), ordToInt)
import qualified Prelude as P
type Monomial n = Sized' n Int
newtype OrderedMonomial ordering n =
OrderedMonomial { getMonomial :: Monomial n }
deriving (NFData)
makeLenses ''OrderedMonomial
makeWrapped ''OrderedMonomial
fromList :: SNat n -> [Int] -> Monomial n
fromList len = V.fromListWithDefault len 0
type MonomialOrder n = Monomial n -> Monomial n -> Ordering
isRelativelyPrime :: OrderedMonomial ord n -> OrderedMonomial ord n -> Bool
isRelativelyPrime n m = lcmMonomial n m == n * m
totalDegree :: OrderedMonomial ord n -> Int
totalDegree = P.sum . getMonomial
lex :: MonomialOrder n
lex m n = P.foldMap (uncurry compare) $ V.zipSame m n
revlex :: MonomialOrder n
revlex xs ys = foldl (flip (<>)) EQ $ V.zipWithSame (flip compare) xs ys
graded :: MonomialOrder n -> MonomialOrder n
graded cmp xs ys = comparing F.sum xs ys <> cmp xs ys
grlex :: MonomialOrder n
grlex = graded lex
grevlex :: MonomialOrder n
grevlex = graded revlex
deriving instance Hashable (Monomial n) => Hashable (OrderedMonomial ordering n)
deriving instance (Eq (Monomial n)) => Eq (OrderedMonomial ordering n)
instance KnownNat n => Show (OrderedMonomial ord n) where
show xs =
let vs = catMaybes $ V.toList $
imap (\n i ->
if i > 0
then Just ("X_" ++ show (ordToInt n) ++ if i == 1 then "" else "^" ++ show i)
else Nothing)
$ getMonomial xs
in if null vs then "1" else unwords vs
instance Multiplicative (OrderedMonomial ord n) where
OrderedMonomial n * OrderedMonomial m = OrderedMonomial $ V.zipWithSame (+) n m
instance KnownNat n => Division (OrderedMonomial ord n) where
recip = _Wrapped %~ V.map P.negate
OrderedMonomial n / OrderedMonomial m = OrderedMonomial $ V.zipWithSame () n m
instance KnownNat n => Unital (OrderedMonomial ord n) where
one = OrderedMonomial $ fromList sing []
class IsOrder (n :: Nat) (ordering :: *) where
cmpMonomial :: Proxy ordering -> MonomialOrder n
head' :: (0 :< n) ~ 'True => Sized' n a -> a
head' = V.head
data Lex = Lex
deriving (Show, Eq, Ord)
data Revlex = Revlex
deriving (Show, Eq, Ord)
data Grevlex = Grevlex
deriving (Show, Eq, Ord)
data Grlex = Grlex
deriving (Show, Eq, Ord)
data Graded ord = Graded ord
deriving (Read, Show, Eq, Ord)
instance IsOrder n ord => IsOrder n (Graded ord) where
cmpMonomial Proxy = graded (cmpMonomial (Proxy :: Proxy ord))
instance IsMonomialOrder n ord => IsMonomialOrder n (Graded ord)
data ProductOrder (n :: Nat) (m :: Nat) (a :: *) (b :: *) where
ProductOrder :: Sing n -> Sing m -> ord -> ord' -> ProductOrder n m ord ord'
productOrder :: forall ord ord' n m. (IsOrder n ord, IsOrder m ord', KnownNat n, KnownNat m)
=> Proxy (ProductOrder n m ord ord') -> MonomialOrder (n + m)
productOrder _ mon mon' =
let n = sing :: SNat n
m = sing :: SNat m
in withWitness (plusLeqL n m) $
case (V.splitAt n mon, V.splitAt n mon') of
((xs, xs'), (ys, ys')) ->
cmpMonomial (Proxy :: Proxy ord) xs ys <>
cmpMonomial (Proxy :: Proxy ord')
(coerceLength (plusMinus' n m) xs')
(coerceLength (plusMinus' n m) ys')
productOrder' :: forall n ord ord' m.(IsOrder n ord, IsOrder m ord')
=> SNat n -> SNat m -> ord -> ord' -> MonomialOrder (n + m)
productOrder' n m _ _ =
withKnownNat n $ withKnownNat m $
productOrder (Proxy :: Proxy (ProductOrder n m ord ord'))
type WeightProxy (v :: [Nat]) = SList v
data WeightOrder (v :: [Nat]) (ord :: Type) where
WeightOrder :: SList (v :: [Nat]) -> Proxy ord -> WeightOrder v ord
calcOrderWeight :: forall vs n. (SingI vs, KnownNat n)
=> Proxy (vs :: [Nat]) -> Monomial n -> Int
calcOrderWeight Proxy = calcOrderWeight' (sing :: SList vs)
calcOrderWeight' :: forall vs n. KnownNat n => SList (vs :: [Nat]) -> Monomial n -> Int
calcOrderWeight' slst m =
let cfs = V.fromListWithDefault' (0 :: Int) $ map P.fromIntegral $ fromSing slst
in P.sum $ V.zipWithSame (*) cfs m
weightOrder :: forall n ns ord. (KnownNat n, IsOrder n ord, SingI ns)
=> Proxy (WeightOrder ns ord) -> MonomialOrder n
weightOrder Proxy m m' =
comparing (calcOrderWeight (Proxy :: Proxy ns)) m m'
<> cmpMonomial (Proxy :: Proxy ord) m m'
instance (KnownNat n, IsOrder n ord, SingI ws)
=> IsOrder n (WeightOrder ws ord) where
cmpMonomial p = weightOrder p
instance (IsOrder n ord, IsOrder m ord', KnownNat m, KnownNat n, k ~ (n + m))
=> IsOrder k (ProductOrder n m ord ord') where
cmpMonomial p = productOrder p
instance IsOrder n Grevlex where
cmpMonomial _ = grevlex
instance IsOrder n Revlex where
cmpMonomial _ = revlex
instance IsOrder n Lex where
cmpMonomial _ = lex
instance IsOrder n Grlex where
cmpMonomial _ = grlex
class IsOrder n name => IsMonomialOrder n name where
instance IsMonomialOrder n Grlex
instance IsMonomialOrder n Grevlex
instance IsMonomialOrder n Lex
instance (KnownNat n, KnownNat m, IsMonomialOrder n o, IsMonomialOrder m o', k ~ (n + m))
=> IsMonomialOrder k (ProductOrder n m o o')
instance (KnownNat k, SingI ws, IsMonomialOrder k ord)
=> IsMonomialOrder k (WeightOrder ws ord)
lcmMonomial :: OrderedMonomial ord n -> OrderedMonomial ord n -> OrderedMonomial ord n
lcmMonomial (OrderedMonomial m) (OrderedMonomial n) = OrderedMonomial $ V.zipWithSame max m n
gcdMonomial :: OrderedMonomial ord n -> OrderedMonomial ord n -> OrderedMonomial ord n
gcdMonomial (OrderedMonomial m) (OrderedMonomial n) = OrderedMonomial $ V.zipWithSame P.min m n
divs :: OrderedMonomial ord n -> OrderedMonomial ord n -> Bool
(OrderedMonomial xs) `divs` (OrderedMonomial ys) = and $ V.toList $ V.zipWith (<=) xs ys
isPowerOf :: KnownNat n => OrderedMonomial ord n -> OrderedMonomial ord n -> Bool
OrderedMonomial n `isPowerOf` OrderedMonomial m =
case V.sFindIndices (> 0) m of
[ind] -> F.sum n == V.sIndex ind n
_ -> False
tryDiv :: Field r => (r, OrderedMonomial ord n) -> (r, OrderedMonomial ord n) -> (r, OrderedMonomial ord n)
tryDiv (a, f) (b, g)
| g `divs` f = (a * recip b, OrderedMonomial $ V.zipWithSame () (getMonomial f) (getMonomial g))
| otherwise = error "cannot divide."
varMonom :: SNat n -> Ordinal n -> Monomial n
varMonom len o = V.replicate len 0 & ix o .~ 1
class (IsMonomialOrder n ord, KnownNat n) => EliminationType n m ord
instance KnownNat n => EliminationType n m Lex
instance (KnownNat n, KnownNat m, IsMonomialOrder n ord, IsMonomialOrder m ord', k ~ (n + m), KnownNat k)
=> EliminationType k n (ProductOrder n m ord ord')
instance (IsMonomialOrder k ord, ones ~ (Replicate n 1), SingI ones,
(Length ones :<= k) ~ 'True, KnownNat k)
=> EliminationType k n (WeightOrder ones ord)
type EliminationOrder n m = ProductOrder n m Grevlex Grevlex
eliminationOrder :: SNat n -> SNat m -> EliminationOrder n m
eliminationOrder n m =
withKnownNat n $ ProductOrder n m Grevlex Grevlex
sOnes :: Sing n -> Sing (Replicate n 1)
sOnes n = sReplicate n (sing :: Sing 1)
weightedEliminationOrder :: SNat n -> WeightedEliminationOrder n Grevlex
weightedEliminationOrder n =
WeightOrder (sOnes n) (Proxy :: Proxy Grevlex)
type WeightedEliminationOrder (n :: Nat) (ord :: Type) =
WeightOrder (Replicate n 1) ord
instance (Eq (Monomial n), IsOrder n name) => Ord (OrderedMonomial name n) where
OrderedMonomial m `compare` OrderedMonomial n = cmpMonomial (Proxy :: Proxy name) m n
instance Ord (Monomial n) where
compare = grevlex
castMonomial :: (KnownNat m) => OrderedMonomial o n -> OrderedMonomial o' m
castMonomial = _Wrapped %~ fromList sing . V.toList
scastMonomial :: SNat m -> OrderedMonomial o n -> OrderedMonomial o m
scastMonomial sdim = _Wrapped %~ fromList sdim . V.toList
changeMonomialOrder :: o' -> OrderedMonomial ord n -> OrderedMonomial o' n
changeMonomialOrder _ = OrderedMonomial . getMonomial
changeMonomialOrderProxy :: Proxy o' -> OrderedMonomial ord n -> OrderedMonomial o' n
changeMonomialOrderProxy _ = OrderedMonomial . getMonomial
class (IsMonomialOrder n ord) => IsMonomialOrder' ord n
instance (IsMonomialOrder n ord) => IsMonomialOrder' ord n
instance IsMonomialOrder' ord n :=> IsMonomialOrder n ord where
ins = C.Sub Dict
type IsStrongMonomialOrder ord = Forall (IsMonomialOrder' ord)
withStrongMonomialOrder :: forall ord n r proxy (proxy' :: Nat -> Type).
(IsStrongMonomialOrder ord)
=> proxy ord -> proxy' n -> (IsMonomialOrder n ord => r) -> r
withStrongMonomialOrder _ _ r = r C.\\ dict
where
ismToPrim = (ins :: IsMonomialOrder' ord n C.:- IsMonomialOrder n ord)
primeInst = inst :: Forall (IsMonomialOrder' ord) C.:- IsMonomialOrder' ord n
dict = ismToPrim `C.trans` primeInst
cmpAnyMonomial :: IsStrongMonomialOrder ord
=> Proxy ord -> Monomial n -> Monomial m -> Ordering
cmpAnyMonomial pxy t t' =
let (l, u, u') = padVecs 0 t t'
in withStrongMonomialOrder pxy l $ cmpMonomial pxy u u'
orderMonomial :: proxy ord -> Monomial n -> OrderedMonomial ord n
orderMonomial _ = OrderedMonomial