summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorHiromiIshii <>2017-01-07 06:35:00 (GMT)
committerhdiff <hdiff@hdiff.luite.com>2017-01-07 06:35:00 (GMT)
commit6e8dfcee72b4b820ac3f4c20059300ecd1200f24 (patch)
tree6094935a50fe9ad63d024376f93c1b792b097224
parentc73227e02f034e90450c65911d6385c594492ac5 (diff)
version 0.5.0.00.5.0.0
-rw-r--r--Algebra/Algorithms/Groebner.hs68
-rw-r--r--Algebra/Prelude/Core.hs4
-rw-r--r--Algebra/Ring/Polynomial.hs452
-rw-r--r--Algebra/Ring/Polynomial/Class.hs29
-rw-r--r--Algebra/Ring/Polynomial/Internal.hs549
-rw-r--r--Algebra/Ring/Polynomial/Labeled.hs59
-rw-r--r--Algebra/Ring/Polynomial/Univariate.hs20
-rw-r--r--computational-algebra.cabal66
8 files changed, 718 insertions, 529 deletions
diff --git a/Algebra/Algorithms/Groebner.hs b/Algebra/Algorithms/Groebner.hs
index 242f842..8222f9f 100644
--- a/Algebra/Algorithms/Groebner.hs
+++ b/Algebra/Algorithms/Groebner.hs
@@ -374,31 +374,42 @@ unsafeThEliminationIdealWith ord n 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.
- ( IsMonomialOrder (k + Arity poly) (MOrder poly),
- Field (Coefficient poly), IsOrderedPolynomial poly)
+ ( Field (Coefficient poly), IsOrderedPolynomial poly)
=> Sized k (Ideal poly)
-> Ideal poly
intersection idsv@(_ :< _) =
let sk = sizedLength idsv
- sn = sing :: SNat (Arity poly)
- in withSingI (sOnes sk) $ withKnownNat (sk %:+ sn) $
+ in withSingI (sOnes sk) $ withKnownNat sk $
let ts = take (fromIntegral $ fromSing sk) vars
- inj :: poly -> OrderedPolynomial (Coefficient poly) (MOrder poly) (k + Arity poly)
- inj = transformMonomial (V.append $ V.replicate sk 0) . injectVars
+ 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 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 :: (IsMonomialOrder (2 + Arity poly) (MOrder poly),
- Field (Coefficient poly), IsOrderedPolynomial poly)
+quotByPrincipalIdeal :: (Field (Coefficient poly), IsOrderedPolynomial poly)
=> Ideal poly
-> poly
-> Ideal poly
@@ -408,49 +419,38 @@ quotByPrincipalIdeal i g =
-- | Ideal quotient by the given ideal.
quotIdeal :: forall poly l.
- (IsOrderedPolynomial poly, Field (Coefficient poly),
- IsMonomialOrder (l + Arity poly) (MOrder poly),
- IsMonomialOrder (2 + Arity poly) (MOrder poly))
+ (IsOrderedPolynomial poly, Field (Coefficient poly))
=> Ideal poly
-> Sized l poly
-> Ideal poly
quotIdeal i g =
withKnownNat (sizedLength g) $
- withKnownNat (sizedLength g %:+ sArity g) $
intersection $ V.map (i `quotByPrincipalIdeal`) g
-- | Saturation by a principal ideal.
saturationByPrincipalIdeal :: forall poly.
- (IsOrderedPolynomial poly, Field (Coefficient poly),
- IsMonomialOrder (1 + Arity poly) (MOrder poly))
+ (IsOrderedPolynomial poly, Field (Coefficient poly))
=> Ideal poly
-> poly
-> Ideal poly
saturationByPrincipalIdeal is g =
let n = sArity' g
- remap :: poly -> OrderedPolynomial (Coefficient poly) (MOrder poly) (1 + Arity poly)
- remap = shiftR sOne . injectVars
- in withKnownNat (sOne %:+ n) $
- withRefl (plusMinus' sOne n) $ withRefl (plusComm n sOne) $
+ in withRefl (plusMinus' sOne n) $ withRefl (plusComm n sOne) $
withWitness (leqStep sOne (sOne %:+ n) n Refl) $
withWitness (lneqZero n) $
- mapIdeal injectVars $
- thEliminationIdeal sOne $
- addToIdeal (one - (remap g * varX)) $
- mapIdeal remap is
+ 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,
- IsMonomialOrder (l + Arity poly) (MOrder poly),
- IsMonomialOrder (1 + Arity poly) (MOrder poly))
+ IsOrderedPolynomial poly)
=> Ideal poly
-> Sized l poly
-> Ideal poly
saturationIdeal i g =
withKnownNat (sizedLength g) $
- withKnownNat (sizedLength g %:+ sArity g) $
intersection $ V.map (i `saturationByPrincipalIdeal`) g
-- | Calculate resultant for given two unary polynomimals.
@@ -488,8 +488,7 @@ 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,
- IsMonomialOrder (2 + Arity poly) (MOrder poly))
+ IsOrderedPolynomial poly)
=> poly
-> poly
-> poly
@@ -497,8 +496,7 @@ lcmPolynomial f g = head $ generators $ intersection (principalIdeal f :< princi
-- | Calculates the Greatest Common Divisor of the given pair of polynomials.
gcdPolynomial :: (Field (Coefficient poly),
- IsOrderedPolynomial poly,
- IsMonomialOrder (2 + Arity poly) (MOrder poly))
+ IsOrderedPolynomial poly)
=> poly
-> poly
-> poly
diff --git a/Algebra/Prelude/Core.hs b/Algebra/Prelude/Core.hs
index 57aa6e1..458d555 100644
--- a/Algebra/Prelude/Core.hs
+++ b/Algebra/Prelude/Core.hs
@@ -3,7 +3,7 @@ module Algebra.Prelude.Core
((%),Scalar(..),(.*.), od,Ordinal, enumOrdinal,
logBase2,ceilingLogBase2,
module AlgebraicPrelude,
- module Algebra.Ring.Polynomial,
+ module Algebra.Ring.Polynomial.Internal,
module Algebra.Ring.Ideal,
module Algebra.Normed,
module Algebra.Internal) where
@@ -11,7 +11,7 @@ module Algebra.Prelude.Core
import Algebra.Internal
import Algebra.Normed
import Algebra.Ring.Ideal
-import Algebra.Ring.Polynomial
+import Algebra.Ring.Polynomial.Internal
import Algebra.Scalar
import AlgebraicPrelude hiding (lex, (%))
diff --git a/Algebra/Ring/Polynomial.hs b/Algebra/Ring/Polynomial.hs
index f01043a..f1e4851 100644
--- a/Algebra/Ring/Polynomial.hs
+++ b/Algebra/Ring/Polynomial.hs
@@ -5,7 +5,6 @@
{-# LANGUAGE PatternGuards, PolyKinds, RankNTypes, ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving, TemplateHaskell, TypeFamilies #-}
{-# LANGUAGE TypeOperators, TypeSynonymInstances, UndecidableInstances #-}
-{-# LANGUAGE ViewPatterns #-}
{-# OPTIONS_GHC -fno-warn-orphans -fno-warn-type-defaults #-}
{-# OPTIONS_GHC -Wno-redundant-constraints #-}
module Algebra.Ring.Polynomial
@@ -20,448 +19,25 @@ module Algebra.Ring.Polynomial
mapCoeff, reversal, padeApprox,
eval, evalUnivariate,
substUnivariate, minpolRecurrent,
- IsOrder(..)
+ IsOrder(..), PadPolyL(..), padLeftPoly
) where
+import Algebra.Algorithms.Groebner
import Algebra.Internal
import Algebra.Ring.Polynomial.Class
+import Algebra.Ring.Polynomial.Internal
import Algebra.Ring.Polynomial.Monomial
-import Algebra.Scalar
-import AlgebraicPrelude
-import Control.DeepSeq (NFData)
-import Control.Lens hiding (assign)
-import qualified Data.Coerce as C
-import qualified Data.Foldable as F
-import qualified Data.HashSet as HS
-import Data.Map (Map)
-import qualified Data.Map.Strict as M
-import qualified Data.Set as Set
-import Data.Singletons.Prelude (POrd (..))
-import qualified Data.Sized.Builtin as S
-import Data.Type.Ordinal
-import qualified Numeric.Algebra as NA
-import Numeric.Algebra.Unital.UnitNormalForm (UnitNormalForm (..))
-import qualified Numeric.Algebra.Unital.UnitNormalForm as NA
-import Numeric.Domain.Integral (IntegralDomain (..))
-import qualified Numeric.Ring.Class as NA
-import Numeric.Semiring.ZeroProduct (ZeroProductSemiring)
-import qualified Prelude as P
-import Proof.Equational (symmetry)
+import AlgebraicPrelude
-instance Hashable r => Hashable (OrderedPolynomial r ord n) where
- hashWithSalt salt poly = hashWithSalt salt $ getTerms poly
-deriving instance (CoeffRing r, IsOrder n ord, Ord r) => Ord (OrderedPolynomial r ord n)
+instance {-# OVERLAPPABLE #-}
+ (KnownNat n, Eq r, DecidableUnits r, DecidableZero r, Field r,
+ IsMonomialOrder n ord, ZeroProductSemiring r)
+ => UFD (OrderedPolynomial r ord n)
--- | n-ary polynomial ring over some noetherian ring R.
-newtype OrderedPolynomial r order n = Polynomial { _terms :: Map (OrderedMonomial order n) r }
- deriving (NFData)
-type Polynomial r = OrderedPolynomial r Grevlex
-
-instance (KnownNat n, IsMonomialOrder n ord, CoeffRing r) => IsPolynomial (OrderedPolynomial r ord n) where
- type Coefficient (OrderedPolynomial r ord n) = r
- type Arity (OrderedPolynomial r ord n) = n
-
- injectCoeff r | isZero r = Polynomial M.empty
- | otherwise = Polynomial $ M.singleton one r
- {-# INLINE injectCoeff #-}
-
- sArity' = sizedLength . getMonomial . leadingMonomial
- {-# INLINE sArity' #-}
-
- mapCoeff' = mapCoeff
- {-# INLINE mapCoeff' #-}
-
- monomials = HS.fromList . map getMonomial . Set.toList . orderedMonomials
- {-# INLINE monomials #-}
-
- fromMonomial m = Polynomial $ M.singleton (OrderedMonomial m) one
- {-# INLINE fromMonomial #-}
-
- toPolynomial' (r, m) = Polynomial $ M.singleton (OrderedMonomial m) r
- {-# INLINE toPolynomial' #-}
-
- polynomial' dic = normalize $ Polynomial $ M.mapKeys OrderedMonomial dic
- {-# INLINE polynomial' #-}
-
- terms' = M.mapKeys getMonomial . terms
- {-# INLINE terms' #-}
-
- liftMap mor poly = sum $ map (uncurry (.*) . (Scalar *** extractPower)) $ getTerms poly
- where
- extractPower = runMult . ifoldMap (\ o -> Mult . pow (mor o) . fromIntegral) . getMonomial
- {-# INLINE liftMap #-}
-
-ordVec :: forall n. KnownNat n => Sized n (Ordinal n)
-ordVec = unsafeFromList' $ enumOrdinal (sing :: SNat n)
-
-instance (KnownNat n, CoeffRing r, IsMonomialOrder n ord)
- => IsOrderedPolynomial (OrderedPolynomial r ord n) where
- -- | coefficient for a degree.
- type MOrder (OrderedPolynomial r ord n) = ord
- coeff d = M.findWithDefault zero d . terms
- {-# INLINE coeff #-}
-
- terms = C.coerce
- {-# INLINE terms #-}
-
- orderedMonomials = M.keysSet . terms
- {-# INLINE orderedMonomials #-}
-
- toPolynomial (c, deg) =
- if isZero c
- then Polynomial M.empty
- else Polynomial $ M.singleton deg c
- {-# INLINE toPolynomial #-}
-
- polynomial = normalize . C.coerce
- {-# INLINE polynomial #-}
-
- leadingTerm (Polynomial d) =
- case M.maxViewWithKey d of
- Just ((deg, c), _) -> (c, deg)
- Nothing -> (zero, one)
- {-# INLINE leadingTerm #-}
-
- leadingMonomial = snd . leadingTerm
- {-# INLINE leadingMonomial #-}
-
- leadingCoeff = fst . leadingTerm
- {-# INLINE leadingCoeff #-}
-
-instance (KnownNat n, CoeffRing r, IsMonomialOrder n order)
- => Wrapped (OrderedPolynomial r order n) where
- type Unwrapped (OrderedPolynomial r order n) = Map (OrderedMonomial order n) r
- _Wrapped' = iso terms polynomial
-
-instance (KnownNat n, CoeffRing r, IsMonomialOrder n ord, t ~ OrderedPolynomial q ord' m)
- => Rewrapped (OrderedPolynomial r ord n) t
-
-castPolynomial :: (CoeffRing r, KnownNat n, KnownNat m,
- IsMonomialOrder n o, IsMonomialOrder m o')
- => OrderedPolynomial r o n
- -> OrderedPolynomial r o' m
-castPolynomial = _Wrapped %~ M.mapKeys castMonomial
-{-# INLINE castPolynomial #-}
-
-scastPolynomial :: (IsMonomialOrder n o, IsMonomialOrder m o', KnownNat m,
- CoeffRing r, KnownNat n)
- => SNat m -> OrderedPolynomial r o n -> OrderedPolynomial r o' m
-scastPolynomial _ = castPolynomial
-{-# INLINE scastPolynomial #-}
-
-mapCoeff :: (KnownNat n, CoeffRing b, IsMonomialOrder n ord)
- => (a -> b) -> OrderedPolynomial a ord n -> OrderedPolynomial b ord n
-mapCoeff f (Polynomial dic) = polynomial $ M.map f dic
-{-# INLINE mapCoeff #-}
-
-normalize :: (DecidableZero r)
- => OrderedPolynomial r order n -> OrderedPolynomial r order n
-normalize (Polynomial dic) =
- Polynomial $ M.filter (not . isZero) dic
-{-# INLINE normalize #-}
-
-
-instance (Eq r) => Eq (OrderedPolynomial r order n) where
- Polynomial f == Polynomial g = f == g
- {-# INLINE (==) #-}
-
--- -- | By Hilbert's finite basis theorem, a polynomial ring over a noetherian ring is also a noetherian ring.
--- instance (IsMonomialOrder order, CoeffRing r, KnownNat n) => Ring (OrderedPolynomial r order n) where
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Ring (OrderedPolynomial r order n) where
- fromInteger 0 = Polynomial M.empty
- fromInteger n = Polynomial $ M.singleton one (fromInteger' n)
- {-# INLINE fromInteger #-}
-
-decZero :: DecidableZero r => r -> Maybe r
-decZero n | isZero n = Nothing
- | otherwise = Just n
-{-# INLINE decZero #-}
-
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Rig (OrderedPolynomial r order n)
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Group (OrderedPolynomial r order n) where
- negate (Polynomial dic) = Polynomial $ fmap negate dic
- {-# INLINE negate #-}
-
- Polynomial f - Polynomial g = Polynomial $ M.mergeWithKey (\_ i j -> decZero (i - j)) id (fmap negate) f g
- {-# INLINE (-) #-}
-
-
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => LeftModule Integer (OrderedPolynomial r order n) where
- n .* Polynomial dic = polynomial $ fmap (n .*) dic
- {-# INLINE (.*) #-}
-
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => RightModule Integer (OrderedPolynomial r order n) where
- (*.) = flip (.*)
- {-# INLINE (*.) #-}
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Additive (OrderedPolynomial r order n) where
- (Polynomial f) + (Polynomial g) = polynomial $ M.unionWith (+) f g
- {-# INLINE (+) #-}
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Monoidal (OrderedPolynomial r order n) where
- zero = Polynomial M.empty
- {-# INLINE zero #-}
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => LeftModule Natural (OrderedPolynomial r order n) where
- n .* Polynomial dic = polynomial $ fmap (n .*) dic
- {-# INLINE (.*) #-}
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => RightModule Natural (OrderedPolynomial r order n) where
- (*.) = flip (.*)
- {-# INLINE (*.) #-}
-
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Unital (OrderedPolynomial r order n) where
- one = Polynomial $ M.singleton one one
- {-# INLINE one #-}
-
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Multiplicative (OrderedPolynomial r order n) where
- Polynomial (M.toList -> d1) * Polynomial (M.toList -> d2) =
- let dic = (one, zero) : [ (a * b, r * r') | (a, r) <- d1, (b, r') <- d2, not $ isZero (r * r')
- ]
- in polynomial $ M.fromListWith (+) dic
- {-# INLINE (*) #-}
-
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Semiring (OrderedPolynomial r order n) where
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Commutative (OrderedPolynomial r order n) where
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Abelian (OrderedPolynomial r order n) where
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => LeftModule (Scalar r) (OrderedPolynomial r order n) where
- Scalar r .* Polynomial dic = polynomial $ fmap (r*) dic
- {-# INLINE (.*) #-}
-
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => RightModule (Scalar r) (OrderedPolynomial r order n) where
- Polynomial dic *. Scalar r = polynomial $ fmap (r*) dic
- {-# INLINE (*.) #-}
-
-
-instance (IsMonomialOrder n ord, Characteristic r, KnownNat n, CoeffRing r)
- => Characteristic (OrderedPolynomial r ord n) where
- char _ = char (Proxy :: Proxy r)
- {-# INLINE char #-}
-
-instance (KnownNat n, CoeffRing r, IsMonomialOrder n order, PrettyCoeff r)
- => Show (OrderedPolynomial r order n) where
- showsPrec = showsPolynomialWith $ generate sing (\i -> "X_" ++ show (fromEnum i))
-
-showPolynomialWithVars :: (CoeffRing a, Show a, KnownNat n, IsMonomialOrder n ordering)
- => [(Int, String)] -> OrderedPolynomial a ordering n -> String
-showPolynomialWithVars dic p0@(Polynomial d)
- | isZero p0 = "0"
- | otherwise = intercalate " + " $ mapMaybe showTerm $ M.toDescList d
- where
- showTerm (getMonomial -> deg, c)
- | isZero c = Nothing
- | otherwise =
- let cstr = if (not (isZero $ c - one) || isConstantMonomial deg)
- then show c ++ " "
- else if isZero (c - one) then ""
- else if isZero (c + one)
- then if any (not . isZero) (F.toList deg) then "-" else "-1"
- else ""
- in Just $ cstr ++ unwords (mapMaybe showDeg (zip [0..] $ F.toList deg))
- showDeg (n, p) | p == 0 = Nothing
- | p == 1 = Just $ showVar n
- | otherwise = Just $ showVar n ++ "^" ++ show p
- showVar n = fromMaybe ("X_" ++ show n) $ lookup n dic
-
-isConstantMonomial :: Monomial n -> Bool
-isConstantMonomial v = all (== 0) $ F.toList v
-
--- | We provide Num instance to use trivial injection R into R[X].
--- Do not use signum or abs.
-instance (IsMonomialOrder n order, CoeffRing r, KnownNat n)
- => P.Num (OrderedPolynomial r order n) where
- (+) = (+)
- {-# INLINE (+) #-}
-
- (*) = (*)
- {-# INLINE (*) #-}
-
- fromInteger = normalize . injectCoeff . fromInteger'
- {-# INLINE fromInteger #-}
-
- signum f = if isZero f then zero else injectCoeff one
- {-# INLINE signum #-}
-
- abs = id
- {-# INLINE abs #-}
-
- negate = ((P.negate 1 :: Integer) .*)
- {-# INLINE negate #-}
-
-
-instance (CoeffRing r, KnownNat n, IsMonomialOrder n ord) => DecidableZero (OrderedPolynomial r ord n) where
- isZero (Polynomial d) = M.null d
- {-# INLINE isZero #-}
-
-instance (CoeffRing r, IsMonomialOrder 1 ord, ZeroProductSemiring r)
- => ZeroProductSemiring (OrderedPolynomial r ord 1)
-
-instance (Eq r, DecidableUnits r, DecidableZero r, Field r,
- IsMonomialOrder 1 ord, ZeroProductSemiring r)
- => DecidableAssociates (OrderedPolynomial r ord 1) where
- isAssociate = (==) `on` NA.normalize
- {-# INLINE isAssociate #-}
-
-instance (Eq r, DecidableUnits r, DecidableZero r, Field r,
- IsMonomialOrder 1 ord, ZeroProductSemiring r)
- => UnitNormalForm (OrderedPolynomial r ord 1) where
- splitUnit f
- | isZero f = (zero, f)
- | otherwise = let lc = leadingCoeff f
- in (injectCoeff lc, injectCoeff (recip lc) * f)
- {-# INLINE splitUnit #-}
-
-instance (Eq r, DecidableUnits r, DecidableZero r, Field r,
- IsMonomialOrder 1 ord, ZeroProductSemiring r)
- => GCDDomain (OrderedPolynomial r ord 1)
-instance (Eq r, DecidableUnits r, DecidableZero r, Field r,
- IsMonomialOrder 1 ord, ZeroProductSemiring r)
- => UFD (OrderedPolynomial r ord 1)
-instance (Eq r, DecidableUnits r, DecidableZero r, Field r,
- IsMonomialOrder 1 ord, ZeroProductSemiring r)
- => PID (OrderedPolynomial r ord 1)
-instance (Eq r, DecidableUnits r, DecidableZero r, Field r, IsMonomialOrder 1 ord, ZeroProductSemiring r) => Euclidean (OrderedPolynomial r ord 1) where
- f0 `divide` g = step f0 zero
- where
- lm = leadingMonomial g
- step p quo
- | isZero p = (quo, p)
- | lm `divs` leadingMonomial p =
- let q = toPolynomial $ leadingTerm p `tryDiv` leadingTerm g
- in step (p - (q * g)) (quo + q)
- | otherwise = (quo, p)
- degree f | isZero f = Nothing
- | otherwise = Just $ P.fromIntegral $ totalDegree' f
-
-
-instance (Eq r, DecidableUnits r, DecidableZero r, KnownNat n,
- Field r, IsMonomialOrder n ord, ZeroProductSemiring r)
- => ZeroProductSemiring (OrderedPolynomial r ord n)
-
-instance (Eq r, DecidableUnits r, DecidableZero r, KnownNat n,
- Field r, IsMonomialOrder n ord, ZeroProductSemiring r)
- => IntegralDomain (OrderedPolynomial r ord n) where
- p `divides` q = isZero $ p `modPolynomial` [q]
- p `maybeQuot` q =
- if isZero q
- then Nothing
- else let (r, s) = p `divModPolynomial` [q]
- in if isZero s
- then Just $ snd $ head r
- else Nothing
-
-instance (CoeffRing r, IsMonomialOrder n ord, DecidableUnits r, KnownNat n) => DecidableUnits (OrderedPolynomial r ord n) where
- isUnit f =
- let (lc, lm) = leadingTerm f
- in lm == one && isUnit lc
- recipUnit f | isUnit f = injectCoeff <$> recipUnit (leadingCoeff f)
- | otherwise = Nothing
-
-varX :: forall r n order. (CoeffRing r, KnownNat n, IsMonomialOrder n order, (0 :< n) ~ 'True)
- => OrderedPolynomial r order n
-varX = var OZ
-
--- | Substitute univariate polynomial using Horner's rule
-substUnivariate :: (Module (Scalar r) b, Unital b, CoeffRing r, IsMonomialOrder 1 order)
- => b -> OrderedPolynomial r order 1 -> b
-substUnivariate u f =
- let n = totalDegree' f
- in foldr (\a b -> Scalar a .* one + b * u)
- (Scalar (coeff (OrderedMonomial $ singleton $ fromIntegral n) f) .* one)
- [ coeff (OrderedMonomial $ singleton $ fromIntegral i) f | i <- [0 .. n P.- 1] ]
-
-evalUnivariate :: (CoeffRing b, IsMonomialOrder 1 order) => b -> OrderedPolynomial b order 1 -> b
-evalUnivariate u f =
- let n = totalDegree' f
- in if n == 0
- then coeff one f
- else foldr1 (\a b -> a + b * u) [ coeff (OrderedMonomial $ singleton $ fromIntegral i) f | i <- [0 .. n] ]
-
--- | Evaluate polynomial at some point.
-eval :: (CoeffRing r, IsMonomialOrder n order, KnownNat n)
- => Sized n r -> OrderedPolynomial r order n -> r
-eval = substWith (*)
-
--- evalOn :: forall k a order . (SingI k, CoeffRing a, IsMonomialOrder order)
--- => OrderedPolynomial a order k -> RepArgs k a a
--- evalOn p = fromNAry $ (fromVecFun (flip eval p) :: NAry k a a)
-
--- | @substVar n f@ substitutes @n@-th variable with polynomial @f@,
--- without changing arity.
-substVar :: (CoeffRing r, KnownNat n, IsMonomialOrder n ord, (1 :<= n) ~ 'True)
- => Ordinal n
- -> OrderedPolynomial r ord n
- -> OrderedPolynomial r ord n
- -> OrderedPolynomial r ord n
-substVar p val =
- liftMap (\o -> if o == p then val else var o)
-
-allVars :: forall k ord n . (IsMonomialOrder n ord, CoeffRing k, KnownNat n)
- => Sized n (OrderedPolynomial k ord n)
-allVars = unsafeFromList' vars
-
-changeOrder :: (CoeffRing k, Eq (Monomial n), IsMonomialOrder n o, IsMonomialOrder n o', KnownNat n)
- => o' -> OrderedPolynomial k o n -> OrderedPolynomial k o' n
-changeOrder _ = _Wrapped %~ M.mapKeys (OrderedMonomial . getMonomial)
-
-changeOrderProxy :: (CoeffRing k, Eq (Monomial n), IsMonomialOrder n o,
- IsMonomialOrder n o', KnownNat n)
- => Proxy o' -> OrderedPolynomial k o n -> OrderedPolynomial k o' n
-changeOrderProxy _ = _Wrapped %~ M.mapKeys (OrderedMonomial . getMonomial)
-
-getTerms :: OrderedPolynomial k order n -> [(k, OrderedMonomial order n)]
-getTerms = map (snd &&& fst) . M.toDescList . _terms
-
-transformMonomial :: (IsMonomialOrder m o, CoeffRing k, KnownNat m)
- => (Monomial n -> Monomial m) -> OrderedPolynomial k o n -> OrderedPolynomial k o m
-transformMonomial tr (Polynomial d) =
- polynomial $ M.mapKeys (OrderedMonomial . tr . getMonomial) d
-
-orderedBy :: OrderedPolynomial k o n -> o -> OrderedPolynomial k o n
-p `orderedBy` _ = p
-
-shiftR :: forall k r n ord. (CoeffRing r, KnownNat n, IsMonomialOrder n ord,
- IsMonomialOrder (k + n) ord)
- => SNat k -> OrderedPolynomial r ord n -> OrderedPolynomial r ord (k :+ n)
-shiftR k = withKnownNat (k %:+ (sing :: SNat n)) $
- withKnownNat k $ transformMonomial (S.append (fromList k []))
-
--- | Calculate the homogenized polynomial of given one, with additional variable is the last variable.
-homogenize :: forall k ord n.
- (CoeffRing k, KnownNat n, IsMonomialOrder (n+1) ord, IsMonomialOrder n ord)
- => OrderedPolynomial k ord n -> OrderedPolynomial k ord (n + 1)
-homogenize f =
- withKnownNat (sSucc (sing :: SNat n)) $
- let g = substWith (.*.) (S.init allVars) f
- d = fromIntegral (totalDegree' g)
- in mapMonomialMonotonic (\m -> m & _Wrapped.ix maxBound .~ d - P.sum (m^._Wrapped)) g
-
-unhomogenize :: forall k ord n.
- (CoeffRing k, KnownNat n, IsMonomialOrder n ord,
- IsMonomialOrder (n+1) ord)
- => OrderedPolynomial k ord (Succ n) -> OrderedPolynomial k ord n
-unhomogenize f =
- withKnownNat (sSucc (sing :: SNat n)) $
- substWith (.*.)
- (coerceLength (symmetry $ succAndPlusOneR (sing :: SNat n)) $
- allVars `S.append` S.singleton one)
- f
-
-reversal :: (CoeffRing k, IsMonomialOrder 1 o)
- => Int -> OrderedPolynomial k o 1 -> OrderedPolynomial k o 1
-reversal k = transformMonomial (S.map (k - ))
-
-padeApprox :: (Field r, DecidableUnits r, CoeffRing r, ZeroProductSemiring r,
- IsMonomialOrder 1 order)
- => Natural -> Natural -> OrderedPolynomial r order 1
- -> (OrderedPolynomial r order 1, OrderedPolynomial r order 1)
-padeApprox k nmk g =
- let (r, _, t) = last $ filter ((< P.fromIntegral k) . totalDegree' . view _1) $ euclid (pow varX (k+nmk)) g
- in (r, t)
-
-
-minpolRecurrent :: forall k. (Eq k, ZeroProductSemiring k, DecidableUnits k, DecidableZero k, Field k)
- => Natural -> [k] -> Polynomial k 1
-minpolRecurrent n xs =
- let h = sum $ zipWith (\a b -> injectCoeff a * b) xs [pow varX i | i <- [0.. pred (2 * n)]]
- :: Polynomial k 1
- (s, t) = padeApprox n n h
- d = fromIntegral $ max (1 + totalDegree' s) (totalDegree' t)
- in reversal d (recip (coeff one t) .*. t)
+instance {-# OVERLAPPABLE #-}
+ (KnownNat n, Eq r, DecidableUnits r, DecidableZero r, Field r,
+ IsMonomialOrder n ord, ZeroProductSemiring r)
+ => GCDDomain (OrderedPolynomial r ord n) where
+ gcd = gcdPolynomial
+ lcm = lcmPolynomial
diff --git a/Algebra/Ring/Polynomial/Class.hs b/Algebra/Ring/Polynomial/Class.hs
index 205e577..d7e8073 100644
--- a/Algebra/Ring/Polynomial/Class.hs
+++ b/Algebra/Ring/Polynomial/Class.hs
@@ -16,6 +16,9 @@ module Algebra.Ring.Polynomial.Class
showsPolynomialWith, showPolynomialWith,
-- * Polynomial division
divModPolynomial, divPolynomial, modPolynomial
+ -- * Default instances
+ , isUnitDefault, recipUnitDefault, isAssociateDefault
+ , splitUnitDefault
) where
import Algebra.Internal
import Algebra.Normed
@@ -32,7 +35,7 @@ import qualified Data.HashSet as HS
import Data.Int
import qualified Data.List as L
import qualified Data.Map.Strict as M
-import Data.Maybe (catMaybes, fromMaybe)
+import Data.Maybe (catMaybes, fromJust, fromMaybe)
import qualified Data.Ratio as R
import qualified Data.Set as S
import Data.Singletons.Prelude (SingKind (..))
@@ -552,3 +555,27 @@ divPolynomial = (fst .) . divModPolynomial
infixl 7 `divPolynomial`
infixl 7 `modPolynomial`
infixl 7 `divModPolynomial`
+
+isUnitDefault :: (DecidableUnits r, Coefficient poly ~ r, IsPolynomial poly)
+ => poly -> Bool
+isUnitDefault p = totalDegree' p == 0 && isUnit (constantTerm p)
+
+recipUnitDefault :: (DecidableUnits r, Coefficient poly ~ r, IsPolynomial poly)
+ => poly -> Maybe poly
+recipUnitDefault p
+ | isUnitDefault p = fmap injectCoeff $ recipUnit $ constantTerm p
+ | otherwise = Nothing
+
+isAssociateDefault :: (UnitNormalForm r, Coefficient poly ~ r, IsOrderedPolynomial poly)
+ => poly -> poly -> Bool
+isAssociateDefault p q =
+ let up = fromJust $ recipUnit $ leadingUnit $ fst $ leadingTerm p
+ uq = fromJust $ recipUnit $ leadingUnit $ fst $ leadingTerm q
+ in (up !* q) == (uq !* p)
+
+splitUnitDefault :: (UnitNormalForm r, Coefficient poly ~ r, IsOrderedPolynomial poly)
+ => poly -> (poly, poly)
+splitUnitDefault f =
+ let u = leadingUnit $ leadingCoeff f
+ u' = fromJust $ recipUnit u
+ in (injectCoeff u, u' !* f)
diff --git a/Algebra/Ring/Polynomial/Internal.hs b/Algebra/Ring/Polynomial/Internal.hs
new file mode 100644
index 0000000..87288b7
--- /dev/null
+++ b/Algebra/Ring/Polynomial/Internal.hs
@@ -0,0 +1,549 @@
+{-# LANGUAGE ConstraintKinds, DataKinds, ExplicitNamespaces #-}
+{-# LANGUAGE FlexibleContexts, FlexibleInstances, GADTs #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving, LiberalTypeSynonyms #-}
+{-# LANGUAGE MultiParamTypeClasses, NoMonomorphismRestriction #-}
+{-# LANGUAGE PatternGuards, PolyKinds, RankNTypes, ScopedTypeVariables #-}
+{-# LANGUAGE StandaloneDeriving, TemplateHaskell, TypeFamilies #-}
+{-# LANGUAGE TypeOperators, TypeSynonymInstances, UndecidableInstances #-}
+{-# LANGUAGE ViewPatterns #-}
+{-# OPTIONS_GHC -fno-warn-orphans -fno-warn-type-defaults #-}
+{-# OPTIONS_GHC -Wno-redundant-constraints #-}
+module Algebra.Ring.Polynomial.Internal
+ ( module Algebra.Ring.Polynomial.Monomial,
+ module Algebra.Ring.Polynomial.Class,
+ Polynomial,
+ transformMonomial,
+ castPolynomial, changeOrder, changeOrderProxy,
+ scastPolynomial, OrderedPolynomial(..),
+ allVars, substVar, homogenize, unhomogenize,
+ normalize, varX, getTerms, shiftR, orderedBy,
+ mapCoeff, reversal, padeApprox,
+ eval, evalUnivariate,
+ substUnivariate, minpolRecurrent,
+ IsOrder(..),
+ PadPolyL(..),
+ padLeftPoly
+ ) where
+import Algebra.Internal
+import Algebra.Ring.Polynomial.Class
+import Algebra.Ring.Polynomial.Monomial
+import Algebra.Scalar
+
+import AlgebraicPrelude
+import Control.DeepSeq (NFData)
+import Control.Lens hiding (assign)
+import qualified Data.Coerce as C
+import qualified Data.Foldable as F
+import qualified Data.HashSet as HS
+import Data.Map (Map)
+import qualified Data.Map.Strict as M
+import qualified Data.Set as Set
+import Data.Singletons.Prelude (POrd (..))
+import Data.Singletons.Prelude.List (Replicate)
+import qualified Data.Sized.Builtin as S
+import Data.Type.Ordinal
+import qualified Numeric.Algebra as NA
+import Numeric.Algebra.Unital.UnitNormalForm (UnitNormalForm (..))
+import Numeric.Domain.Integral (IntegralDomain (..))
+import Numeric.Semiring.ZeroProduct (ZeroProductSemiring)
+import qualified Prelude as P
+import Proof.Equational (symmetry)
+
+instance Hashable r => Hashable (OrderedPolynomial r ord n) where
+ hashWithSalt salt poly = hashWithSalt salt $ getTerms poly
+
+deriving instance (CoeffRing r, IsOrder n ord, Ord r) => Ord (OrderedPolynomial r ord n)
+
+-- | n-ary polynomial ring over some noetherian ring R.
+newtype OrderedPolynomial r order n = Polynomial { _terms :: Map (OrderedMonomial order n) r }
+ deriving (NFData)
+type Polynomial r = OrderedPolynomial r Grevlex
+
+instance (KnownNat n, IsMonomialOrder n ord, CoeffRing r) => IsPolynomial (OrderedPolynomial r ord n) where
+ type Coefficient (OrderedPolynomial r ord n) = r
+ type Arity (OrderedPolynomial r ord n) = n
+
+ injectCoeff r | isZero r = Polynomial M.empty
+ | otherwise = Polynomial $ M.singleton one r
+ {-# INLINE injectCoeff #-}
+
+ sArity' = sizedLength . getMonomial . leadingMonomial
+ {-# INLINE sArity' #-}
+
+ mapCoeff' = mapCoeff
+ {-# INLINE mapCoeff' #-}
+
+ monomials = HS.fromList . map getMonomial . Set.toList . orderedMonomials
+ {-# INLINE monomials #-}
+
+ fromMonomial m = Polynomial $ M.singleton (OrderedMonomial m) one
+ {-# INLINE fromMonomial #-}
+
+ toPolynomial' (r, m) = Polynomial $ M.singleton (OrderedMonomial m) r
+ {-# INLINE toPolynomial' #-}
+
+ polynomial' dic = normalize $ Polynomial $ M.mapKeys OrderedMonomial dic
+ {-# INLINE polynomial' #-}
+
+ terms' = M.mapKeys getMonomial . terms
+ {-# INLINE terms' #-}
+
+ liftMap mor poly = sum $ map (uncurry (.*) . (Scalar *** extractPower)) $ getTerms poly
+ where
+ extractPower = runMult . ifoldMap (\ o -> Mult . pow (mor o) . fromIntegral) . getMonomial
+ {-# INLINE liftMap #-}
+
+ordVec :: forall n. KnownNat n => Sized n (Ordinal n)
+ordVec = unsafeFromList' $ enumOrdinal (sing :: SNat n)
+
+instance (KnownNat n, CoeffRing r, IsMonomialOrder n ord)
+ => IsOrderedPolynomial (OrderedPolynomial r ord n) where
+ -- | coefficient for a degree.
+ type MOrder (OrderedPolynomial r ord n) = ord
+ coeff d = M.findWithDefault zero d . terms
+ {-# INLINE coeff #-}
+
+ terms = C.coerce
+ {-# INLINE terms #-}
+
+ orderedMonomials = M.keysSet . terms
+ {-# INLINE orderedMonomials #-}
+
+ toPolynomial (c, deg) =
+ if isZero c
+ then Polynomial M.empty
+ else Polynomial $ M.singleton deg c
+ {-# INLINE toPolynomial #-}
+
+ polynomial = normalize . C.coerce
+ {-# INLINE polynomial #-}
+
+ leadingTerm (Polynomial d) =
+ case M.maxViewWithKey d of
+ Just ((deg, c), _) -> (c, deg)
+ Nothing -> (zero, one)
+ {-# INLINE leadingTerm #-}
+
+ leadingMonomial = snd . leadingTerm
+ {-# INLINE leadingMonomial #-}
+
+ leadingCoeff = fst . leadingTerm
+ {-# INLINE leadingCoeff #-}
+
+instance (KnownNat n, CoeffRing r, IsMonomialOrder n order)
+ => Wrapped (OrderedPolynomial r order n) where
+ type Unwrapped (OrderedPolynomial r order n) = Map (OrderedMonomial order n) r
+ _Wrapped' = iso terms polynomial
+
+instance (KnownNat n, CoeffRing r, IsMonomialOrder n ord, t ~ OrderedPolynomial q ord' m)
+ => Rewrapped (OrderedPolynomial r ord n) t
+
+castPolynomial :: (CoeffRing r, KnownNat n, KnownNat m,
+ IsMonomialOrder n o, IsMonomialOrder m o')
+ => OrderedPolynomial r o n
+ -> OrderedPolynomial r o' m
+castPolynomial = _Wrapped %~ M.mapKeys castMonomial
+{-# INLINE castPolynomial #-}
+
+scastPolynomial :: (IsMonomialOrder n o, IsMonomialOrder m o', KnownNat m,
+ CoeffRing r, KnownNat n)
+ => SNat m -> OrderedPolynomial r o n -> OrderedPolynomial r o' m
+scastPolynomial _ = castPolynomial
+{-# INLINE scastPolynomial #-}
+
+mapCoeff :: (KnownNat n, CoeffRing b, IsMonomialOrder n ord)
+ => (a -> b) -> OrderedPolynomial a ord n -> OrderedPolynomial b ord n
+mapCoeff f (Polynomial dic) = polynomial $ M.map f dic
+{-# INLINE mapCoeff #-}
+
+normalize :: (DecidableZero r)
+ => OrderedPolynomial r order n -> OrderedPolynomial r order n
+normalize (Polynomial dic) =
+ Polynomial $ M.filter (not . isZero) dic
+{-# INLINE normalize #-}
+
+
+instance (Eq r) => Eq (OrderedPolynomial r order n) where
+ Polynomial f == Polynomial g = f == g
+ {-# INLINE (==) #-}
+
+-- -- | By Hilbert's finite basis theorem, a polynomial ring over a noetherian ring is also a noetherian ring.
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Ring (OrderedPolynomial r order n) where
+ fromInteger 0 = Polynomial M.empty
+ fromInteger n = Polynomial $ M.singleton one (fromInteger' n)
+ {-# INLINE fromInteger #-}
+
+decZero :: DecidableZero r => r -> Maybe r
+decZero n | isZero n = Nothing
+ | otherwise = Just n
+{-# INLINE decZero #-}
+
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Rig (OrderedPolynomial r order n)
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Group (OrderedPolynomial r order n) where
+ negate (Polynomial dic) = Polynomial $ fmap negate dic
+ {-# INLINE negate #-}
+
+ Polynomial f - Polynomial g = Polynomial $ M.mergeWithKey (\_ i j -> decZero (i - j)) id (fmap negate) f g
+ {-# INLINE (-) #-}
+
+
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => LeftModule Integer (OrderedPolynomial r order n) where
+ n .* Polynomial dic = polynomial $ fmap (n .*) dic
+ {-# INLINE (.*) #-}
+
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => RightModule Integer (OrderedPolynomial r order n) where
+ (*.) = flip (.*)
+ {-# INLINE (*.) #-}
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Additive (OrderedPolynomial r order n) where
+ (Polynomial f) + (Polynomial g) = polynomial $ M.unionWith (+) f g
+ {-# INLINE (+) #-}
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Monoidal (OrderedPolynomial r order n) where
+ zero = Polynomial M.empty
+ {-# INLINE zero #-}
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => LeftModule Natural (OrderedPolynomial r order n) where
+ n .* Polynomial dic = polynomial $ fmap (n .*) dic
+ {-# INLINE (.*) #-}
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => RightModule Natural (OrderedPolynomial r order n) where
+ (*.) = flip (.*)
+ {-# INLINE (*.) #-}
+
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Unital (OrderedPolynomial r order n) where
+ one = Polynomial $ M.singleton one one
+ {-# INLINE one #-}
+
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Multiplicative (OrderedPolynomial r order n) where
+ Polynomial (M.toList -> d1) * Polynomial (M.toList -> d2) =
+ let dic = (one, zero) : [ (a * b, r * r') | (a, r) <- d1, (b, r') <- d2, not $ isZero (r * r')
+ ]
+ in polynomial $ M.fromListWith (+) dic
+ {-# INLINE (*) #-}
+
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Semiring (OrderedPolynomial r order n) where
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Commutative (OrderedPolynomial r order n) where
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => Abelian (OrderedPolynomial r order n) where
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => LeftModule (Scalar r) (OrderedPolynomial r order n) where
+ Scalar r .* Polynomial dic = polynomial $ fmap (r*) dic
+ {-# INLINE (.*) #-}
+
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n) => RightModule (Scalar r) (OrderedPolynomial r order n) where
+ Polynomial dic *. Scalar r = polynomial $ fmap (r*) dic
+ {-# INLINE (*.) #-}
+
+
+instance (IsMonomialOrder n ord, Characteristic r, KnownNat n, CoeffRing r)
+ => Characteristic (OrderedPolynomial r ord n) where
+ char _ = char (Proxy :: Proxy r)
+ {-# INLINE char #-}
+
+instance (KnownNat n, CoeffRing r, IsMonomialOrder n order, PrettyCoeff r)
+ => Show (OrderedPolynomial r order n) where
+ showsPrec = showsPolynomialWith $ generate sing (\i -> "X_" ++ show (fromEnum i))
+
+showPolynomialWithVars :: (CoeffRing a, Show a, KnownNat n, IsMonomialOrder n ordering)
+ => [(Int, String)] -> OrderedPolynomial a ordering n -> String
+showPolynomialWithVars dic p0@(Polynomial d)
+ | isZero p0 = "0"
+ | otherwise = intercalate " + " $ mapMaybe showTerm $ M.toDescList d
+ where
+ showTerm (getMonomial -> deg, c)
+ | isZero c = Nothing
+ | otherwise =
+ let cstr = if (not (isZero $ c - one) || isConstantMonomial deg)
+ then show c ++ " "
+ else if isZero (c - one) then ""
+ else if isZero (c + one)
+ then if any (not . isZero) (F.toList deg) then "-" else "-1"
+ else ""
+ in Just $ cstr ++ unwords (mapMaybe showDeg (zip [0..] $ F.toList deg))
+ showDeg (n, p) | p == 0 = Nothing
+ | p == 1 = Just $ showVar n
+ | otherwise = Just $ showVar n ++ "^" ++ show p
+ showVar n = fromMaybe ("X_" ++ show n) $ lookup n dic
+
+isConstantMonomial :: Monomial n -> Bool
+isConstantMonomial v = all (== 0) $ F.toList v
+
+-- | We provide Num instance to use trivial injection R into R[X].
+-- Do not use signum or abs.
+instance (IsMonomialOrder n order, CoeffRing r, KnownNat n)
+ => P.Num (OrderedPolynomial r order n) where
+ (+) = (+)
+ {-# INLINE (+) #-}
+
+ (*) = (*)
+ {-# INLINE (*) #-}
+
+ fromInteger = normalize . injectCoeff . fromInteger'
+ {-# INLINE fromInteger #-}
+
+ signum f = if isZero f then zero else injectCoeff one
+ {-# INLINE signum #-}
+
+ abs = id
+ {-# INLINE abs #-}
+
+ negate = ((P.negate 1 :: Integer) .*)
+ {-# INLINE negate #-}
+
+instance (CoeffRing r, KnownNat n, IsMonomialOrder n ord) => DecidableZero (OrderedPolynomial r ord n) where
+ isZero (Polynomial d) = M.null d
+ {-# INLINE isZero #-}
+
+instance (Eq r, KnownNat n, Euclidean r, IsMonomialOrder n ord)
+ => DecidableAssociates (OrderedPolynomial r ord n) where
+ isAssociate = isAssociateDefault
+ {-# INLINE isAssociate #-}
+
+instance (Eq r, Euclidean r, KnownNat n,
+ IsMonomialOrder n ord)
+ => UnitNormalForm (OrderedPolynomial r ord n) where
+ splitUnit = splitUnitDefault
+ {-# INLINE splitUnit #-}
+
+instance {-# OVERLAPPING #-}
+ (Eq r, DecidableUnits r, DecidableZero r, Field r,
+ IsMonomialOrder 1 ord, ZeroProductSemiring r)
+ => GCDDomain (OrderedPolynomial r ord 1)
+instance {-# OVERLAPPING #-}
+ (Eq r, DecidableUnits r, DecidableZero r, Field r,
+ IsMonomialOrder 1 ord, ZeroProductSemiring r)
+ => UFD (OrderedPolynomial r ord 1)
+instance (Eq r, DecidableUnits r, DecidableZero r, Field r,
+ IsMonomialOrder 1 ord, ZeroProductSemiring r)
+ => PID (OrderedPolynomial r ord 1)
+instance (Eq r, DecidableUnits r, DecidableZero r, Field r, IsMonomialOrder 1 ord, ZeroProductSemiring r) => Euclidean (OrderedPolynomial r ord 1) where
+ f0 `divide` g = step f0 zero
+ where
+ lm = leadingMonomial g
+ step p quo
+ | isZero p = (quo, p)
+ | lm `divs` leadingMonomial p =
+ let q = toPolynomial $ leadingTerm p `tryDiv` leadingTerm g
+ in step (p - (q * g)) (quo + q)
+ | otherwise = (quo, p)
+ degree f | isZero f = Nothing
+ | otherwise = Just $ P.fromIntegral $ totalDegree' f
+
+
+instance (Eq r, DecidableUnits r, DecidableZero r, KnownNat n,
+ Field r, IsMonomialOrder n ord, ZeroProductSemiring r)
+ => ZeroProductSemiring (OrderedPolynomial r ord n)
+
+instance {-# OVERLAPPING #-}
+ (Eq r, DecidableUnits r, DecidableZero r,
+ Field r, IsMonomialOrder 1 ord, ZeroProductSemiring r)
+ => IntegralDomain (OrderedPolynomial r ord 1)
+
+instance (Eq r, DecidableUnits r, DecidableZero r, KnownNat n,
+ Field r, IsMonomialOrder n ord, ZeroProductSemiring r)
+ => IntegralDomain (OrderedPolynomial r ord n) where
+ p `divides` q = isZero $ p `modPolynomial` [q]
+ p `maybeQuot` q =
+ if isZero q
+ then Nothing
+ else let (r, s) = p `divModPolynomial` [q]
+ in if isZero s
+ then Just $ snd $ head r
+ else Nothing
+
+instance (CoeffRing r, IsMonomialOrder n ord, DecidableUnits r, KnownNat n)
+ => DecidableUnits (OrderedPolynomial r ord n) where
+ isUnit = isUnitDefault
+ recipUnit = recipUnitDefault
+
+varX :: forall r n order. (CoeffRing r, KnownNat n, IsMonomialOrder n order, (0 :< n) ~ 'True)
+ => OrderedPolynomial r order n
+varX = var OZ
+
+-- | Substitute univariate polynomial using Horner's rule
+substUnivariate :: (Module (Scalar r) b, Unital b, CoeffRing r, IsMonomialOrder 1 order)
+ => b -> OrderedPolynomial r order 1 -> b
+substUnivariate u f =
+ let n = totalDegree' f
+ in foldr (\a b -> Scalar a .* one + b * u)
+ (Scalar (coeff (OrderedMonomial $ singleton $ fromIntegral n) f) .* one)
+ [ coeff (OrderedMonomial $ singleton $ fromIntegral i) f | i <- [0 .. n P.- 1] ]
+
+evalUnivariate :: (CoeffRing b, IsMonomialOrder 1 order) => b -> OrderedPolynomial b order 1 -> b
+evalUnivariate u f =
+ let n = totalDegree' f
+ in if n == 0
+ then coeff one f
+ else foldr1 (\a b -> a + b * u) [ coeff (OrderedMonomial $ singleton $ fromIntegral i) f | i <- [0 .. n] ]
+
+-- | Evaluate polynomial at some point.
+eval :: (CoeffRing r, IsMonomialOrder n order, KnownNat n)
+ => Sized n r -> OrderedPolynomial r order n -> r
+eval = substWith (*)
+
+-- evalOn :: forall k a order . (SingI k, CoeffRing a, IsMonomialOrder order)
+-- => OrderedPolynomial a order k -> RepArgs k a a
+-- evalOn p = fromNAry $ (fromVecFun (flip eval p) :: NAry k a a)
+
+-- | @substVar n f@ substitutes @n@-th variable with polynomial @f@,
+-- without changing arity.
+substVar :: (CoeffRing r, KnownNat n, IsMonomialOrder n ord, (1 :<= n) ~ 'True)
+ => Ordinal n
+ -> OrderedPolynomial r ord n
+ -> OrderedPolynomial r ord n
+ -> OrderedPolynomial r ord n
+substVar p val =
+ liftMap (\o -> if o == p then val else var o)
+
+allVars :: forall k ord n . (IsMonomialOrder n ord, CoeffRing k, KnownNat n)
+ => Sized n (OrderedPolynomial k ord n)
+allVars = unsafeFromList' vars
+
+changeOrder :: (CoeffRing k, Eq (Monomial n), IsMonomialOrder n o, IsMonomialOrder n o', KnownNat n)
+ => o' -> OrderedPolynomial k o n -> OrderedPolynomial k o' n
+changeOrder _ = _Wrapped %~ M.mapKeys (OrderedMonomial . getMonomial)
+
+changeOrderProxy :: (CoeffRing k, Eq (Monomial n), IsMonomialOrder n o,
+ IsMonomialOrder n o', KnownNat n)
+ => Proxy o' -> OrderedPolynomial k o n -> OrderedPolynomial k o' n
+changeOrderProxy _ = _Wrapped %~ M.mapKeys (OrderedMonomial . getMonomial)
+
+getTerms :: OrderedPolynomial k order n -> [(k, OrderedMonomial order n)]
+getTerms = map (snd &&& fst) . M.toDescList . _terms
+
+transformMonomial :: (IsMonomialOrder m o, CoeffRing k, KnownNat m)
+ => (Monomial n -> Monomial m) -> OrderedPolynomial k o n -> OrderedPolynomial k o m
+transformMonomial tr (Polynomial d) =
+ polynomial $ M.mapKeys (OrderedMonomial . tr . getMonomial) d
+
+orderedBy :: OrderedPolynomial k o n -> o -> OrderedPolynomial k o n
+p `orderedBy` _ = p
+
+shiftR :: forall k r n ord. (CoeffRing r, KnownNat n, IsMonomialOrder n ord,
+ IsMonomialOrder (k + n) ord)
+ => SNat k -> OrderedPolynomial r ord n -> OrderedPolynomial r ord (k :+ n)
+shiftR k = withKnownNat (k %:+ (sing :: SNat n)) $
+ withKnownNat k $ transformMonomial (S.append (fromList k []))
+
+-- | Calculate the homogenized polynomial of given one, with additional variable is the last variable.
+homogenize :: forall k ord n.
+ (CoeffRing k, KnownNat n, IsMonomialOrder (n+1) ord, IsMonomialOrder n ord)
+ => OrderedPolynomial k ord n -> OrderedPolynomial k ord (n + 1)
+homogenize f =
+ withKnownNat (sSucc (sing :: SNat n)) $
+ let g = substWith (.*.) (S.init allVars) f
+ d = fromIntegral (totalDegree' g)
+ in mapMonomialMonotonic (\m -> m & _Wrapped.ix maxBound .~ d - P.sum (m^._Wrapped)) g
+
+unhomogenize :: forall k ord n.
+ (CoeffRing k, KnownNat n, IsMonomialOrder n ord,
+ IsMonomialOrder (n+1) ord)
+ => OrderedPolynomial k ord (Succ n) -> OrderedPolynomial k ord n
+unhomogenize f =
+ withKnownNat (sSucc (sing :: SNat n)) $
+ substWith (.*.)
+ (coerceLength (symmetry $ succAndPlusOneR (sing :: SNat n)) $
+ allVars `S.append` S.singleton one)
+ f
+
+reversal :: (CoeffRing k, IsMonomialOrder 1 o)
+ => Int -> OrderedPolynomial k o 1 -> OrderedPolynomial k o 1
+reversal k = transformMonomial (S.map (k - ))
+
+padeApprox :: (Field r, DecidableUnits r, CoeffRing r, ZeroProductSemiring r,
+ IsMonomialOrder 1 order)
+ => Natural -> Natural -> OrderedPolynomial r order 1
+ -> (OrderedPolynomial r order 1, OrderedPolynomial r order 1)
+padeApprox k nmk g =
+ let (r, _, t) = last $ filter ((< P.fromIntegral k) . totalDegree' . view _1) $ euclid (pow varX (k+nmk)) g
+ in (r, t)
+
+
+minpolRecurrent :: forall k. (Eq k, ZeroProductSemiring k, DecidableUnits k, DecidableZero k, Field k)
+ => Natural -> [k] -> Polynomial k 1
+minpolRecurrent n xs =
+ let h = sum $ zipWith (\a b -> injectCoeff a * b) xs [pow varX i | i <- [0.. pred (2 * n)]]
+ :: Polynomial k 1
+ (s, t) = padeApprox n n h
+ d = fromIntegral $ max (1 + totalDegree' s) (totalDegree' t)
+ in reversal d (recip (coeff one t) .*. t)
+
+newtype PadPolyL n ord poly = PadPolyL { runPadPolyL :: OrderedPolynomial poly (Graded ord) n }
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => Additive (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => LeftModule Natural (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => RightModule Natural (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => Monoidal (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => LeftModule Integer (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => RightModule Integer (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => Group (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => Abelian (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => Multiplicative (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => Unital (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => Commutative (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => Eq (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => Semiring (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => Rig (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => Ring (PadPolyL n ord poly)
+deriving instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => DecidableZero (PadPolyL n ord poly)
+instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly,
+ LeftModule (Scalar r) poly)
+ => LeftModule (Scalar r) (PadPolyL n ord poly) where
+ r .* PadPolyL f = PadPolyL $ mapCoeff' (r .*) f
+ {-# INLINE (.*) #-}
+instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly,
+ RightModule (Scalar r) poly)
+ => RightModule (Scalar r) (PadPolyL n ord poly) where
+ PadPolyL f *. r = PadPolyL $ mapCoeff' (*. r) f
+ {-# INLINE (*.) #-}
+
+instance (KnownNat n, IsMonomialOrder n ord, IsPolynomial poly)
+ => IsPolynomial (PadPolyL n ord poly) where
+ type Coefficient (PadPolyL n ord poly) = Coefficient poly
+ type Arity (PadPolyL n ord poly) = n + Arity poly
+ sArity _ = sing
+ liftMap f = subst $ S.generate sing f
+ subst vec (PadPolyL f) =
+ let sn = sing :: Sing n
+ sm = sing :: Sing (Arity poly)
+ in withWitness (plusLeqL sn sm) $
+ withRefl (plusMinus' sn sm) $
+ case S.splitAt sn vec of
+ (ls, rs) -> substWith (\ g a -> a * subst rs g) ls f
+ injectCoeff = PadPolyL . injectCoeff . injectCoeff
+ fromMonomial m =
+ let sn = sing :: Sing n
+ sm = sing :: Sing (Arity poly)
+ in withWitness (plusLeqL sn sm) $
+ withRefl (plusMinus' sn sm) $
+ case S.splitAt sn m of
+ (ls, rs) -> PadPolyL $ fromMonomial ls * injectCoeff (fromMonomial rs)
+ terms' (PadPolyL m) =
+ M.fromList
+ [ (ls S.++ rs, k)
+ | (ls, pol) <- M.toList $ terms' m
+ , (rs, k) <- M.toList $ terms' pol
+ ]
+
+instance (SingI (Replicate n 1), KnownNat n, IsMonomialOrder n ord, IsOrderedPolynomial poly)
+ => IsOrderedPolynomial (PadPolyL n ord poly) where
+ type MOrder (PadPolyL n ord poly) =
+ ProductOrder n (Arity poly) (Graded ord) (MOrder poly)
+ leadingTerm (PadPolyL f) =
+ let (p, OrderedMonomial ls) = leadingTerm f
+ (k, OrderedMonomial rs) = leadingTerm p
+ in (k, OrderedMonomial $ ls S.++ rs)
+
+padLeftPoly :: (IsMonomialOrder n ord, IsPolynomial poly)
+ => Sing n -> ord -> poly -> PadPolyL n ord poly
+padLeftPoly n _ = withKnownNat n $ PadPolyL . injectCoeff
diff --git a/Algebra/Ring/Polynomial/Labeled.hs b/Algebra/Ring/Polynomial/Labeled.hs
index 3ebeebf..649e108 100644
--- a/Algebra/Ring/Polynomial/Labeled.hs
+++ b/Algebra/Ring/Polynomial/Labeled.hs
@@ -1,9 +1,9 @@
{-# LANGUAGE CPP, ConstraintKinds, DataKinds, EmptyCase, FlexibleContexts #-}
-{-# LANGUAGE FlexibleInstances, GADTs, KindSignatures, IncoherentInstances #-}
-{-# LANGUAGE MultiParamTypeClasses, PolyKinds, RankNTypes #-}
-{-# LANGUAGE ScopedTypeVariables, StandaloneDeriving, TemplateHaskell #-}
-{-# LANGUAGE TypeFamilies, TypeInType, TypeOperators, UndecidableInstances #-}
-{-# LANGUAGE UndecidableSuperClasses, OverloadedLabels #-}
+{-# LANGUAGE FlexibleInstances, GADTs, GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE IncoherentInstances, KindSignatures, MultiParamTypeClasses #-}
+{-# LANGUAGE OverloadedLabels, PolyKinds, RankNTypes, ScopedTypeVariables #-}
+{-# LANGUAGE StandaloneDeriving, TemplateHaskell, TypeFamilies, TypeInType #-}
+{-# LANGUAGE TypeOperators, UndecidableInstances, UndecidableSuperClasses #-}
module Algebra.Ring.Polynomial.Labeled
(IsUniqueList, LabPolynomial(..),
LabPolynomial', LabUnipol,
@@ -11,13 +11,14 @@ module Algebra.Ring.Polynomial.Labeled
canonicalMap',
IsSubsetOf) where
import Algebra.Internal
-import Algebra.Ring.Polynomial.Class
import Algebra.Ring.Polynomial
import Algebra.Ring.Polynomial.Univariate
import Algebra.Scalar
-import qualified Prelude as P
+import AlgebraicPrelude
+import Control.Lens (each, (%~), (&))
import Data.Function (on)
+import qualified Data.List as L
import Data.Singletons.Prelude
import Data.Singletons.Prelude.Enum (SEnum (..))
import Data.Singletons.Prelude.List hiding (Group)
@@ -25,12 +26,9 @@ import qualified Data.Sized.Builtin as S
import Data.Type.Natural.Class (IsPeano (..), sOne)
import Data.Type.Ordinal
import GHC.Exts (Constraint)
-import qualified Data.List as L
-import Numeric.Algebra hiding (Order (..))
-import Numeric.Decidable.Zero
-import Prelude hiding (Integral (..), Num (..),
- product, sum)
-import GHC.OverloadedLabels (IsLabel(..))
+import GHC.OverloadedLabels (IsLabel (..))
+import qualified Numeric.Algebra as NA
+import qualified Prelude as P
type family UniqueList' (x :: Symbol) (xs :: [Symbol]) :: Constraint where
UniqueList' x '[] = ()
@@ -132,7 +130,7 @@ instance (Wraps vars poly, Monoidal poly) => Monoidal (LabPolynomial poly vars)
instance (Wraps vars poly, Semiring poly) => Semiring (LabPolynomial poly vars)
instance (Wraps vars poly, Rig poly) => Rig (LabPolynomial poly vars)
instance (Wraps vars poly, Ring poly) => Ring (LabPolynomial poly vars) where
- fromInteger n = LabelPolynomial (fromInteger n :: poly)
+ fromInteger n = LabelPolynomial (NA.fromInteger n :: poly)
{-# INLINE fromInteger #-}
instance (Wraps vars poly, LeftModule (Scalar r) poly) => LeftModule (Scalar r) (LabPolynomial poly vars) where
@@ -223,6 +221,39 @@ class (All (FlipSym0 @@ ElemSym0 @@ ys) xs ~ 'True) => IsSubsetOf (xs :: [a])
_suppress _ _ = id
instance (All (FlipSym0 @@ ElemSym0 @@ ys) xs ~ 'True) => IsSubsetOf (xs :: [a]) (ys :: [a])
+instance (ZeroProductSemiring poly , Wraps vars poly) => ZeroProductSemiring (LabPolynomial poly vars)
+instance (IntegralDomain poly , Wraps vars poly) => IntegralDomain (LabPolynomial poly vars) where
+ divides = divides `on` unLabelPolynomial
+ maybeQuot f g = LabelPolynomial <$> maybeQuot (unLabelPolynomial f) (unLabelPolynomial g)
+instance (UFD poly , Wraps vars poly) => UFD (LabPolynomial poly vars)
+instance (PID poly , Wraps vars poly) => PID (LabPolynomial poly vars) where
+ egcd (LabelPolynomial f) (LabelPolynomial g) =
+ egcd f g & each %~ LabelPolynomial
+instance (GCDDomain poly , Wraps vars poly) => GCDDomain (LabPolynomial poly vars) where
+ gcd f g = LabelPolynomial $ gcd (unLabelPolynomial f) (unLabelPolynomial g)
+ reduceFraction f g =
+ reduceFraction (unLabelPolynomial f) (unLabelPolynomial g)
+ & each %~ LabelPolynomial
+ lcm f g = LabelPolynomial $ lcm (unLabelPolynomial f) (unLabelPolynomial g)
+instance (UnitNormalForm poly , Wraps vars poly) => UnitNormalForm (LabPolynomial poly vars) where
+ splitUnit = (each %~ LabelPolynomial) . splitUnit . unLabelPolynomial
+instance (DecidableUnits poly , Wraps vars poly) => DecidableUnits (LabPolynomial poly vars) where
+ isUnit = isUnit . unLabelPolynomial
+ recipUnit = fmap LabelPolynomial . recipUnit . unLabelPolynomial
+ LabelPolynomial f ^? n = LabelPolynomial <$> (f ^? n)
+
+instance (DecidableAssociates poly , Wraps vars poly)
+ => DecidableAssociates (LabPolynomial poly vars) where
+ isAssociate = isAssociate `on` unLabelPolynomial
+
+instance (Euclidean poly , Wraps vars poly)
+ => Euclidean (LabPolynomial poly vars) where
+ degree = degree . unLabelPolynomial
+ divide (LabelPolynomial f) (LabelPolynomial g) =
+ divide f g & each %~ LabelPolynomial
+ quot f g = LabelPolynomial $ quot (unLabelPolynomial f) (unLabelPolynomial g)
+ rem f g = LabelPolynomial $ rem (unLabelPolynomial f) (unLabelPolynomial g)
+
-- | So unsafe! Don't expose it!
permute0 :: (SEq k) => SList (xs :: [k]) -> SList (ys :: [k]) -> Sized (Length xs) Integer
permute0 SNil _ = S.NilL
diff --git a/Algebra/Ring/Polynomial/Univariate.hs b/Algebra/Ring/Polynomial/Univariate.hs
index 39c802d..e3e981a 100644
--- a/Algebra/Ring/Polynomial/Univariate.hs
+++ b/Algebra/Ring/Polynomial/Univariate.hs
@@ -7,7 +7,7 @@
module Algebra.Ring.Polynomial.Univariate
(Unipol(), naiveMult, karatsuba,
divModUnipolByMult, divModUnipol,
- mapCoeffUnipol,
+ mapCoeffUnipol, liftMapUnipol,
module Algebra.Ring.Polynomial.Class,
module Algebra.Ring.Polynomial.Monomial) where
import Algebra.Prelude.Core
@@ -326,12 +326,7 @@ instance CoeffRing r => IsPolynomial (Unipol r) where
arity _ = 1
constantTerm = IM.findWithDefault zero 0 . runUnipol
{-# INLINE constantTerm #-}
- liftMap g f@(Unipol dic) =
- let u = g 0
- n = maybe 0 (fst . fst) $ IM.maxViewWithKey $ runUnipol f
- in foldr (\a b -> a .*. one + b * u)
- (IM.findWithDefault zero n dic .*. one)
- [IM.findWithDefault zero k dic | k <- [0..n-1]]
+ liftMap = liftMapUnipol
{-# INLINABLE liftMap #-}
fromMonomial = Unipol . flip IM.singleton one . SV.head
{-# INLINE fromMonomial #-}
@@ -371,3 +366,14 @@ instance (CoeffRing r, PrettyCoeff r) => Show (Unipol r) where
mapCoeffUnipol :: DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol f (Unipol a) =
Unipol $ IM.mapMaybe (decZero . f) a
+{-# INLINE mapCoeffUnipol #-}
+
+liftMapUnipol :: (Module (Scalar k) r, Monoidal k, Unital r)
+ => (Ordinal 1 -> r) -> Unipol k -> r
+liftMapUnipol g f@(Unipol dic) =
+ let u = g 0
+ n = maybe 0 (fst . fst) $ IM.maxViewWithKey $ runUnipol f
+ in foldr (\a b -> a .*. one + b * u)
+ (IM.findWithDefault zero n dic .*. one)
+ [IM.findWithDefault zero k dic | k <- [0..n-1]]
+{-# INLINE liftMapUnipol #-}
diff --git a/computational-algebra.cabal b/computational-algebra.cabal
index bd45f71..7f67eef 100644
--- a/computational-algebra.cabal
+++ b/computational-algebra.cabal
@@ -1,10 +1,10 @@
name: computational-algebra
-version: 0.4.0.0
+version: 0.5.0.0
cabal-version: >=1.10
build-type: Simple
license: BSD3
license-file: LICENSE
-copyright: (C) Hiromi ISHII 2013
+copyright: (C) Hiromi ISHII 2017
maintainer: konn.jinro_at_gmail.com
homepage: https://github.com/konn/computational-algebra
synopsis: Well-kinded computational algebra library, currently supporting Groebner basis.
@@ -89,7 +89,8 @@ library
control-monad-loop ==0.1.*,
primes >=0.2.1 && <0.3,
singletons ==2.2.*,
- arithmoi >=0.4.3.0 && <0.5
+ arithmoi >=0.4.3.0 && <0.5,
+ ghc-typelits-knownnat >=0.2.2 && <0.3
default-language: Haskell2010
default-extensions: CPP DataKinds PolyKinds GADTs
MultiParamTypeClasses TypeFamilies FlexibleContexts
@@ -98,14 +99,15 @@ library
Algebra.Algorithms.FGLM
Algebra.Field.Galois.Conway
Algebra.Field.Galois.Internal
- ghc-options: -O2 -Wall -Wno-unused-top-binds
+ Algebra.Ring.Polynomial.Internal
+ ghc-options: -O2 -Wall -Wno-unused-top-binds -fplugin GHC.TypeLits.KnownNat.Solver
executable groebner-prof
main-is: groebner-prof.hs
buildable: False
build-depends:
base >=4.9.0.0 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
deepseq >=1.4.2.0 && <1.5
default-language: Haskell2010
extensions: NoImplicitPrelude
@@ -121,7 +123,7 @@ executable solve
algebra >=4.1 && <4.4,
base >=4 && <4.10,
type-natural >=0.7.1.2 && <0.8,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
random >=1.0 && <1.2,
hmatrix >=0.17.0.2 && <0.18,
matrix ==0.3.*,
@@ -139,7 +141,7 @@ executable algebraic
build-depends:
base >=4.9.0.0 && <4.10,
algebraic-prelude >=0.1.0.1 && <0.2,
- computational-algebra >=0.4.0.0 && <0.5
+ computational-algebra >=0.5.0.0 && <0.6
default-language: Haskell2010
hs-source-dirs: examples
ghc-options: -Wall -O2 -threaded
@@ -156,7 +158,7 @@ executable ipsolve
reflection >=1.4 && <2.2,
base >=4 && <4.10,
type-natural >=0.7.1.2 && <0.8,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
random >=1.0 && <1.2,
hmatrix >=0.17.0.2 && <0.18,
matrix ==0.3.*,
@@ -179,7 +181,7 @@ executable faugere-prof
algebra ==4.3.*,
base >=4 && <4.10,
type-natural >=0.7.1.2 && <0.8,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
random >=1.0 && <1.2,
hmatrix >=0.17.0.2 && <0.18,
matrix ==0.3.*,
@@ -202,7 +204,7 @@ executable hensel-prof
algebra ==4.3.*,
base >=4 && <4.10,
type-natural >=0.7.1.2 && <0.8,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
random >=1.0 && <1.2,
hmatrix >=0.17.0.2 && <0.18,
matrix ==0.3.*,
@@ -221,7 +223,7 @@ executable sandpit-poly
build-depends:
semigroups >=0.15.2 && <0.19,
constraints >=0.3 && <0.9,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
base >=4 && <4.10,
type-natural >=0.7.1.2 && <0.8,
algebra ==4.3.*,
@@ -234,7 +236,7 @@ executable quotient
build-depends:
semigroups >=0.15.2 && <0.19,
constraints >=0.3 && <0.9,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
base >=4 && <4.10,
type-natural >=0.7.1.2 && <0.8,
algebra ==4.3.*,
@@ -252,7 +254,7 @@ test-suite test-multi-table
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
hspec >=1.9.5 && <2.3,
lazysmallcheck ==0.6.*,
@@ -283,7 +285,7 @@ test-suite singular-test
MonadRandom >=0.1 && <0.5,
QuickCheck >=2.6 && <2.9,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
deepseq >=1.3 && <1.5,
hspec >=2.2.4 && <2.3,
monomorphic >=0.0.3 && <0.1,
@@ -316,7 +318,7 @@ test-suite monomial-order-test
MonadRandom >=0.1 && <0.5,
QuickCheck >=2.6 && <2.9,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
deepseq >=1.3 && <1.5,
hspec >=2.2.4 && <2.3,
@@ -341,7 +343,7 @@ test-suite linear-test
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
deepseq >=1.3 && <1.5,
hspec >=2.2.4 && <2.3,
@@ -373,7 +375,7 @@ test-suite matrix-test
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
deepseq >=1.3 && <1.5,
hspec >=2.2.4 && <2.3,
@@ -404,7 +406,7 @@ test-suite specs
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
deepseq >=1.3 && <1.5,
hspec >=2.2.4 && <2.3,
@@ -446,7 +448,7 @@ test-suite new-div-test
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
deepseq >=1.3 && <1.5,
hspec >=2.2.4 && <2.3,
@@ -474,7 +476,7 @@ benchmark unipol-bench
constraints >=0.3 && <0.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,
@@ -494,7 +496,7 @@ benchmark normal-bench
constraints >=0.3 && <0.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,
@@ -517,7 +519,7 @@ benchmark elimination-bench
constraints >=0.3 && <0.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,
@@ -540,7 +542,7 @@ benchmark quotient-bench-randomized
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,
@@ -573,7 +575,7 @@ benchmark monomial-order-bench
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,
@@ -609,7 +611,7 @@ benchmark linear-bench
algebra ==4.3.*,
sized >=0.2.1.0 && <0.3,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,
@@ -645,7 +647,7 @@ benchmark division-bench
algebra ==4.3.*,
sized >=0.2.1.0 && <0.3,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,
@@ -680,7 +682,7 @@ benchmark sugar-paper-bench
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,
@@ -715,7 +717,7 @@ benchmark solve-bench
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,
@@ -751,7 +753,7 @@ benchmark coercion-bench
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,
@@ -785,7 +787,7 @@ benchmark faugere4-bench
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,
@@ -819,7 +821,7 @@ benchmark unipol-mult-bench
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,
@@ -853,7 +855,7 @@ benchmark unipol-div-bench
QuickCheck >=2.6 && <2.9,
algebra ==4.3.*,
base >=4 && <4.10,
- computational-algebra >=0.4.0.0 && <0.5,
+ computational-algebra >=0.5.0.0 && <0.6,
containers ==0.5.*,
criterion >=0.8.1.0 && <1.2,
deepseq >=1.3 && <1.5,