computational-algebra 0.0.1.1 → 0.0.2.0
raw patch · 12 files changed
+775/−187 lines, 12 files
Files
- Algebra/Algorithms/Groebner.hs +35/−10
- Algebra/Algorithms/Groebner/Monomorphic.hs +151/−11
- Algebra/Internal.hs +243/−17
- Algebra/Ring/Polynomial.hs +4/−1
- Algebra/Ring/Polynomial/Monomorphic.hs +123/−62
- Algebra/Ring/Polynomial/Parser.hs +26/−24
- Example.hs +0/−59
- Monomorphic.hs +0/−1
- README.md +46/−0
- computational-algebra.cabal +4/−2
- examples/monomorphic.hs +84/−0
- examples/polymorphic.hs +59/−0
Algebra/Algorithms/Groebner.hs view
@@ -1,7 +1,17 @@ {-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts, GADTs #-} {-# LANGUAGE MultiParamTypeClasses, NoImplicitPrelude, ParallelListComp #-} {-# LANGUAGE RankNTypes, ScopedTypeVariables, TypeOperators #-}-module Algebra.Algorithms.Groebner where+module Algebra.Algorithms.Groebner (+ -- * Polynomial division+ divModPolynomial, divPolynomial, modPolynomial+ -- * Groebner basis+ , calcGroebnerBasis, calcGroebnerBasisWith+ , simpleBuchberger, reduceMinimalGroebnerBasis, minimizeGroebnerBasis+ -- * Ideal operations+ , isIdealMember, intersection, thEliminationIdeal+ , quotIdeal, quotByPrincipalIdeal+ , saturationIdeal, saturationByPrincipalIdeal+ ) where import Algebra.Internal import Algebra.Ring.Noetherian import Algebra.Ring.Polynomial@@ -9,6 +19,7 @@ import Numeric.Algebra import Prelude hiding (Num (..), recip) +-- | Calculate a polynomial quotient and remainder w.r.t. second argument. divModPolynomial :: (IsMonomialOrder order, IsPolynomial r n, Field r) => OrderedPolynomial r order n -> [OrderedPolynomial r order n] -> ([(OrderedPolynomial r order n, OrderedPolynomial r order n)], OrderedPolynomial r order n) divModPolynomial f0 fs = loop f0 zero (zip (nub fs) (repeat zero))@@ -24,12 +35,14 @@ dic' = xs ++ (g, old + q) : ys in loop (p - (q * g)) r dic' +-- | Remainder of given polynomial w.r.t. the second argument. modPolynomial :: (IsPolynomial r n, Field r, IsMonomialOrder order) => OrderedPolynomial r order n -> [OrderedPolynomial r order n] -> OrderedPolynomial r order n modPolynomial = (snd .) . divModPolynomial +-- | A Quotient of given polynomial w.r.t. the second argument. divPolynomial :: (IsPolynomial r n, Field r, IsMonomialOrder order) => OrderedPolynomial r order n -> [OrderedPolynomial r order n]@@ -40,6 +53,7 @@ infixl 7 `modPolynomial` infixl 7 `divModPolynomial` +-- | Apply Buchberger's algorithm and calculate Groebner basis for the given ideal. simpleBuchberger :: (Field k, IsPolynomial k n, IsMonomialOrder order) => Ideal (OrderedPolynomial k order n) -> [OrderedPolynomial k order n] simpleBuchberger ideal =@@ -51,40 +65,46 @@ , let q = sPolynomial f g `modPolynomial` acc, q /= zero ] +-- | Minimize the given groebner basis. minimizeGroebnerBasis :: (Field k, IsPolynomial k n, IsMonomialOrder order) => [OrderedPolynomial k order n] -> [OrderedPolynomial k order n] minimizeGroebnerBasis = foldr step [] where step x xs = injectCoeff (recip $ leadingCoeff x) * x : filter (not . (leadingMonomial x `divs`) . leadingMonomial) xs +-- | Reduce minimum Groebner basis into reduced Groebner basis. reduceMinimalGroebnerBasis :: (Field k, IsPolynomial k n, IsMonomialOrder order) => [OrderedPolynomial k order n] -> [OrderedPolynomial k order n] reduceMinimalGroebnerBasis bs = filter (/= zero) $ map red bs where red x = x `modPolynomial` delete x bs +-- | Caliculating reduced Groebner basis of the given ideal w.r.t. the specified monomial order. calcGroebnerBasisWith :: (Field k, IsPolynomial k n, IsMonomialOrder order, IsMonomialOrder order') => order -> Ideal (OrderedPolynomial k order' n) -> [OrderedPolynomial k order n]-calcGroebnerBasisWith order i = calcGroebnerBasis $ mapIdeal (changeOrder order) i+calcGroebnerBasisWith ord i = calcGroebnerBasis $ mapIdeal (changeOrder ord) i +-- | Caliculating reduced Groebner basis of the given ideal w.r.t. the graded reversed lexicographic order. calcGroebnerBasis :: (Field k, IsPolynomial k n, IsMonomialOrder order) => Ideal (OrderedPolynomial k order n) -> [OrderedPolynomial k order n] calcGroebnerBasis = reduceMinimalGroebnerBasis . minimizeGroebnerBasis . simpleBuchberger +-- | Test if the given polynomial is the member of the ideal. isIdealMember :: (IsPolynomial k n, Field k, IsMonomialOrder o) => OrderedPolynomial k o n -> Ideal (OrderedPolynomial k o n) -> Bool isIdealMember f ideal = groebnerTest f (calcGroebnerBasis ideal) +-- | Test if the given polynomial can be divided by the given polynomials. groebnerTest :: (IsPolynomial k n, Field k, IsMonomialOrder order) => OrderedPolynomial k order n -> [OrderedPolynomial k order n] -> Bool groebnerTest f fs = f `modPolynomial` fs == zero --thEliminationIdeal :: ( IsMonomialOrder ord, Field k, IsPolynomial k m, IsPolynomial k (n :+: m)- )+-- | Calculate n-th elimination ideal.+thEliminationIdeal :: ( IsMonomialOrder ord, Field k, IsPolynomial k m, IsPolynomial k (m :-: n)+ , (n :<<= m) ~ True) => SNat n- -> Ideal (OrderedPolynomial k ord (n :+: m))- -> Ideal (OrderedPolynomial k Lex m)+ -> Ideal (OrderedPolynomial k ord m)+ -> Ideal (OrderedPolynomial k Lex (m :-: n)) thEliminationIdeal n ideal = toIdeal $ [transformMonomial (dropV n) f | f <- calcGroebnerBasisWith Lex ideal , all (all (== 0) . take (toInt n) . toList . snd) $ getTerms f@@ -95,7 +115,8 @@ ( IsMonomialOrder ord, Field r, IsPolynomial r k, IsPolynomial r n , IsPolynomial r (k :+: n) )- => Vector (Ideal (OrderedPolynomial r ord n)) k -> Ideal (OrderedPolynomial r Lex n)+ => Vector (Ideal (OrderedPolynomial r ord n)) k+ -> Ideal (OrderedPolynomial r Lex n) intersection Nil = Ideal $ singletonV one intersection idsv@(_ :- _) = let sk = sLengthV idsv@@ -103,7 +124,9 @@ ts = genVars (sk %+ sn) tis = zipWith (\ideal t -> mapIdeal ((t *) . shiftR sk) ideal) (toList idsv) ts j = foldr appendIdeal (principalIdeal (one - foldr (+) zero ts)) tis- in sk `thEliminationIdeal` j+ in case plusMinusEqR sn sk of+ Eql -> case propToBoolLeq (plusLeqL sk sn) of+ LeqTrueInstance -> sk `thEliminationIdeal` j -- | Ideal quotient by a principal ideals. quotByPrincipalIdeal :: (Field k, IsPolynomial k n, IsMonomialOrder ord)@@ -114,6 +137,7 @@ case intersection (i :- (Ideal $ singletonV g) :- Nil) of Ideal gs -> Ideal $ mapV (snd . head . (`divPolynomial` [changeOrder Lex g])) gs +-- | Ideal quotient by the given ideal. quotIdeal :: forall k ord n. (IsPolynomial k n, Field k, IsMonomialOrder ord) => Ideal (OrderedPolynomial k ord n) -> Ideal (OrderedPolynomial k ord n)@@ -129,9 +153,10 @@ => Ideal (OrderedPolynomial k ord n) -> OrderedPolynomial k ord n -> Ideal (OrderedPolynomial k Lex n) saturationByPrincipalIdeal is g =- case leqSucc (sDegree g) of+ case propToClassLeq $ leqSucc (sDegree g) of LeqInstance -> sOne `thEliminationIdeal` addToIdeal (one - (castPolynomial g * var sOne)) (mapIdeal (shiftR sOne) is) +-- | Saturation ideal saturationIdeal :: forall k ord n. (IsPolynomial k n, Field k, IsMonomialOrder ord) => Ideal (OrderedPolynomial k ord n) -> Ideal (OrderedPolynomial k ord n)
Algebra/Algorithms/Groebner/Monomorphic.hs view
@@ -1,30 +1,170 @@-{-# LANGUAGE FlexibleInstances, GADTs, PolyKinds, RecordWildCards #-}-{-# LANGUAGE TypeFamilies, TypeOperators, ViewPatterns #-}-module Algebra.Algorithms.Groebner.Monomorphic where+{-# LANGUAGE ConstraintKinds, FlexibleInstances, GADTs, PolyKinds #-}+{-# LANGUAGE RecordWildCards, ScopedTypeVariables, TypeFamilies #-}+{-# LANGUAGE TypeOperators, UndecidableInstances #-}+-- | Monomorphic interface for Groenber basis.+module Algebra.Algorithms.Groebner.Monomorphic+ ( Groebnerable+ -- * Polynomial division+ , divModPolynomial, divPolynomial, modPolynomial+ , divModPolynomialWith, divPolynomialWith, modPolynomialWith+ -- * Groebner basis+ , calcGroebnerBasis, calcGroebnerBasisWith+ -- * Ideal operations+ , isIdealMember, intersection, thEliminationIdeal, eliminate+ , quotIdeal, quotByPrincipalIdeal+ , saturationIdeal, saturationByPrincipalIdeal+ -- * Re-exports+ , Lex(..), Revlex(..), Grlex(..), Grevlex(..), IsOrder, IsMonomialOrder+ ) where import qualified Algebra.Algorithms.Groebner as Gr import Algebra.Internal import Algebra.Ring.Noetherian+import Algebra.Ring.Polynomial (Grevlex (..), Grlex (..),+ IsMonomialOrder, IsOrder,+ Lex (..), Revlex (..),+ orderedBy) import qualified Algebra.Ring.Polynomial as Poly import Algebra.Ring.Polynomial.Monomorphic+import Control.Arrow import Data.List-import Monomorphic+import qualified Data.Map as M+import Numeric.Algebra+import Prelude hiding (Num (..)) -calcGroebnerBasis :: [Polyn] -> [Polyn]-calcGroebnerBasis (filter (any ((/= 0).fst)) -> []) = []-calcGroebnerBasis j =- case uniformlyPromote j :: Monomorphic (Ideal :.: Poly.Polynomial Rational) of+-- | Synonym+class (Eq r, Field r, NoetherianRing r) => Groebnerable r+instance (Eq r, Field r, NoetherianRing r) => Groebnerable r++-- | Calculate a intersection of given ideals.+intersection :: forall r. (Groebnerable r)+ => [[Polynomial r]] -> [Polynomial r]+intersection ps =+ let vars = nub $ sort $ concatMap (concatMap buildVarsList) ps+ dim = length vars+ in case promote dim of+ Monomorphic sdim ->+ case singInstance sdim of+ SingInstance ->+ case promote ps :: Monomorphic (Vector [Polynomial r]) of+ Monomorphic vec ->+ let slen = sLengthV vec+ in case singInstance slen of+ SingInstance ->+ let ids = mapV (toIdeal . map (flip orderedBy Lex . Poly.polynomial . M.mapKeys (Poly.OrderedMonomial . Poly.fromList sdim . encodeMonomList vars) . unPolynomial)) vec+ in case singInstance (slen %+ sdim) of+ SingInstance -> demoteComposed $ Gr.intersection ids++freshVar :: [Polynomial r] -> Variable+freshVar ps =+ case maximum $ concatMap buildVarsList ps of+ Variable c Nothing -> Variable c (Just 1)+ Variable c (Just n) -> Variable c (Just $ n + 1)++-- | Calculate saturation ideal by the principal ideal generated by the second argument.+saturationByPrincipalIdeal :: (Groebnerable r)+ => [Polynomial r] -> Polynomial r -> [Polynomial r]+saturationByPrincipalIdeal j g =+ let t = freshVar (g : j)+ in eliminate [t] $ (one - g * injectVar t) : j++-- | Calculate saturation ideal. The saturation of an ideal I by an ideal J is defined as follows:+-- I : J^∞ = { f ∈ k[X] | ∃ n > 0 s.t. f J^n ⊆ I }+saturationIdeal :: Groebnerable r => [Polynomial r] -> [Polynomial r] -> [Polynomial r]+saturationIdeal i g = intersection $ map (i `saturationByPrincipalIdeal`) g++-- | Calculate ideal quotient of I by principal ideal+quotByPrincipalIdeal :: Groebnerable r => [Polynomial r] -> Polynomial r -> [Polynomial r]+quotByPrincipalIdeal i g =+ map (snd . head . flip (divPolynomialWith Lex) [g]) $ intersection [i, [g]]++-- | Calculate the ideal quotient of I of J, defind as follows:+-- I : J = { f ∈ k[X] | fJ ⊆ I }+quotIdeal :: Groebnerable r => [Polynomial r] -> [Polynomial r] -> [Polynomial r]+quotIdeal i g = intersection $ map (i `quotByPrincipalIdeal`) g++divModPolynomial :: Groebnerable r+ => Polynomial r -> [Polynomial r] -> ([(Polynomial r, Polynomial r)], Polynomial r)+divModPolynomial = divModPolynomialWith Grevlex++divModPolynomialWith :: forall ord r. (IsMonomialOrder ord, Groebnerable r)+ => ord -> Polynomial r -> [Polynomial r]+ -> ([(Polynomial r, Polynomial r)], Polynomial r)+divModPolynomialWith _ f gs =+ case promoteList (f:gs) :: Monomorphic ([] :.: Poly.OrderedPolynomial r ord) of+ Monomorphic (Comp (f' : gs')) ->+ let sn = Poly.sDegree f'+ in case singInstance sn of+ SingInstance ->+ let (q, r) = Gr.divModPolynomial f' gs'+ in (map (renameVars vars . polyn . demote' *** polyn . demote') q, polyn $ demote' r)+ where+ vars = nub $ sort $ concatMap buildVarsList (f:gs)++divPolynomial :: Groebnerable r => Polynomial r -> [Polynomial r] -> [(Polynomial r, Polynomial r)]+divPolynomial = (fst .) . divModPolynomial++modPolynomial :: Groebnerable r => Polynomial r -> [Polynomial r] -> Polynomial r+modPolynomial = (snd .) . divModPolynomial++divPolynomialWith :: Groebnerable r => IsMonomialOrder ord => ord -> Polynomial r -> [Polynomial r] -> [(Polynomial r, Polynomial r)]+divPolynomialWith ord = (fst .) . divModPolynomialWith ord++modPolynomialWith :: (Groebnerable r, IsMonomialOrder ord)+ => ord -> Polynomial r -> [Polynomial r] -> Polynomial r+modPolynomialWith ord = (snd .) . divModPolynomialWith ord++calcGroebnerBasis :: Groebnerable r => [Polynomial r] -> [Polynomial r]+calcGroebnerBasis = calcGroebnerBasisWith Grevlex++calcGroebnerBasisWith :: forall ord r. (Groebnerable r, IsMonomialOrder ord)+ => ord -> [Polynomial r] -> [Polynomial r]+calcGroebnerBasisWith _ ps | any (== zero) ps = []+calcGroebnerBasisWith ord j =+ case uniformlyPromote j :: Monomorphic (Ideal :.: Poly.OrderedPolynomial r ord) of Monomorphic (Comp ideal) -> case ideal of Ideal vec -> case singInstance (Poly.sDegree (head $ toList vec)) of- SingInstance -> map (renameVars vars . polyn . demote . Monomorphic) $ Gr.calcGroebnerBasis ideal+ SingInstance -> map (renameVars vars . polyn . demote . Monomorphic) $ Gr.calcGroebnerBasisWith ord ideal where vars = nub $ sort $ concatMap buildVarsList j -isIdealMember :: Polyn -> [Polyn] -> Bool+isIdealMember :: forall r. Groebnerable r => Polynomial r -> [Polynomial r] -> Bool isIdealMember f ideal =- case promoteList (f:ideal) :: Monomorphic ([] :.: Poly.Polynomial Rational) of+ case promoteList (f:ideal) :: Monomorphic ([] :.: Poly.Polynomial r) of Monomorphic (Comp (f':ideal')) -> case singInstance (Poly.sDegree f') of SingInstance -> Gr.isIdealMember f' (toIdeal ideal') _ -> error "impossible happend!"++-- | Computes the ideal with specified variables eliminated.+eliminate :: forall r. Groebnerable r => [Variable] -> [Polynomial r] -> [Polynomial r]+eliminate elvs j =+ case promoteListWithVarOrder (els ++ rest) j :: Monomorphic ([] :.: Poly.OrderedPolynomial r Poly.Lex) of+ Monomorphic (Comp fs) ->+ case promote k of+ Monomorphic sk ->+ let sdim = Poly.sDegree $ head fs+ newDim = sMax sk sdim+ in case singInstance sdim of+ SingInstance ->+ case propToClassLeq $ maxLeqR sk sdim of+ LeqInstance ->+ case singInstance newDim of+ SingInstance ->+ let fs' = map ((flip Poly.orderedBy Poly.Lex) . Poly.scastPolynomial newDim) fs+ in case propToBoolLeq $ maxLeqL sk sdim of+ LeqTrueInstance ->+ case singInstance (newDim %- sk) of+ SingInstance ->+ map (renameVars rest) $ demoteComposed $ sk `Gr.thEliminationIdeal` toIdeal fs'+ where+ vars = nub $ sort $ concatMap buildVarsList j+ (els, rest) = partition (`elem` elvs) vars+ k = length els++-- | Computes nth elimination ideal.+thEliminationIdeal :: Groebnerable r => Int -> [Polynomial r] -> [Polynomial r]+thEliminationIdeal k j = eliminate (take k vars) j+ where+ vars = nub $ sort $ concatMap buildVarsList j
Algebra/Internal.hs view
@@ -2,7 +2,26 @@ {-# LANGUAGE DataKinds, FlexibleContexts, FlexibleInstances, GADTs #-} {-# LANGUAGE MultiParamTypeClasses, PolyKinds, StandaloneDeriving #-} {-# LANGUAGE TypeFamilies, TypeOperators #-}-module Algebra.Internal where+module Algebra.Internal ( toProxy, Nat(..), SNat(..), Vector(..), Sing(..)+ , SingInstance(..), singInstance, toInt+ , Min, Max, sMin, sMax, sZ, sS, (:+:), (%+), (:-:), (%-)+ , sZero, sOne, sTwo, sThree, Zero, One, Two, Three+ , SZero, SOne, STwo, SThree+ , lengthV, sLengthV, takeV, dropV, splitAtV, appendV+ , foldrV, foldlV, singletonV, zipWithV, toList, allV+ , mapV, headV, tailV+ , Leq(..), (:<<=), (:<=), LeqInstance(..)+ , LeqTrueInstance(..), boolToPropLeq, boolToClassLeq+ , propToClassLeq, propToBoolLeq+ , leqRefl, leqSucc, Eql(..), eqlRefl, eqlSymm+ , eqlTrans, plusZR, plusZL, eqPreservesS, plusAssociative+ , sAndPlusOne, plusCommutative, minusCongEq, minusNilpotent+ , eqSuccMinus, plusMinusEqL, plusMinusEqR, plusLeqL, plusLeqR+ , zAbsorbsMinR, zAbsorbsMinL, minLeqL, minLeqR+ , leqRhs, leqLhs, leqTrans, minComm, leqAnitsymmetric+ , maxZL, maxComm, maxZR, maxLeqL, maxLeqR+ , module Monomorphic+ ) where import Data.Proxy import Monomorphic @@ -29,6 +48,16 @@ type instance Max (S n) Z = S n type instance Max (S n) (S m) = S (Max n m) +-- | The smart constructor for @SZ@.+sZ :: SNat Z+sZ = case singInstance SZ of+ SingInstance -> SZ++-- | The smart constructor for @SS n@.+sS :: SNat n -> SNat (S n)+sS n = case singInstance n of+ SingInstance -> SS n+ type Zero = Z type One = S Z type Two = S (S Z)@@ -88,7 +117,7 @@ promote n | n < 0 = error "negative integer!" | n == 0 = Monomorphic SZ- | otherwise = withPolymorhic n $ \sn -> Monomorphic $ SS sn+ | otherwise = withPolymorhic (n - 1) $ \sn -> Monomorphic $ SS sn instance Monomorphicable (Vector a) where type MonomorphicRep (Vector a) = [a]@@ -110,6 +139,48 @@ SZ %+ n = n SS n %+ m = SS (n %+ m) +type family (n :: Nat) :-: (m :: Nat) :: Nat+type instance n :-: Z = n+type instance Z :-: m = Z+type instance S n :-: S m = n :-: m++(%-) :: (m :<<= n) ~ True => SNat n -> SNat m -> SNat (n :-: m)+n %- SZ = n+SS n %- SS m = n %- m+_ %- _ = error "impossible!"++-- | Comparison function+type family (n :: Nat) :<<= (m :: Nat) :: Bool+type instance Z :<<= n = True+type instance S n :<<= Z = False+type instance S n :<<= S m = n :<<= m++-- | Comparison witness via GADTs.+data Leq (n :: Nat) (m :: Nat) where+ ZeroLeq :: SNat m -> Leq Zero m+ SuccLeqSucc :: Leq n m -> Leq (S n) (S m)++data LeqInstance n m where+ LeqInstance :: (n :<= m) => LeqInstance n m++boolToPropLeq :: (n :<<= m) ~ True => SNat n -> SNat m -> Leq n m+boolToPropLeq SZ m = ZeroLeq m+boolToPropLeq (SS n) (SS m) = SuccLeqSucc $ boolToPropLeq n m+boolToPropLeq _ _ = error "impossible happend!"++boolToClassLeq :: (n :<<= m) ~ True => SNat n -> SNat m -> LeqInstance n m+boolToClassLeq SZ _ = LeqInstance+boolToClassLeq (SS n) (SS m) =+ case boolToClassLeq n m of+ LeqInstance -> LeqInstance+boolToClassLeq _ _ = error "impossible!"++propToClassLeq :: Leq n m -> LeqInstance n m+propToClassLeq (ZeroLeq _) = LeqInstance+propToClassLeq (SuccLeqSucc leq) =+ case propToClassLeq leq of+ LeqInstance -> LeqInstance+ appendV :: Vector a n -> Vector a m -> Vector a (n :+: m) appendV (x :- xs) ys = x :- appendV xs ys appendV Nil ys = ys@@ -138,17 +209,21 @@ (x :- xs) == (y :- ys) = x == y && xs == ys _ == _ = error "impossible!" -allV :: (a -> Bool) -> Vector a n-> Bool+allV :: (a -> Bool) -> Vector a n -> Bool allV p = foldrV ((&&) . p) False -dropV :: SNat n -> Vector a (n :+: m) -> Vector a m+dropV :: (n :<<= m) ~ True => SNat n -> Vector a m -> Vector a (m :-: n) dropV n = snd . splitAtV n +takeV :: (n :<<= m) ~ True => SNat n -> Vector a m -> Vector a n+takeV n = fst . splitAtV n+ toInt :: SNat n -> Int toInt SZ = 0 toInt (SS n) = 1 + toInt n -splitAtV :: SNat n -> Vector a (n :+: m) -> (Vector a n, Vector a m)+splitAtV :: (n :<<= m) ~ True+ => SNat n -> Vector a m -> (Vector a n, Vector a (m :-: n)) splitAtV SZ xs = (Nil, xs) splitAtV (SS n) (x :- xs) = case splitAtV n xs of@@ -178,18 +253,169 @@ case singInstance n of SingInstance -> SingInstance -data LeqInstance n m where- LeqInstance :: (n :<= m) => LeqInstance n m+leqRefl :: SNat n -> Leq n n+leqRefl SZ = ZeroLeq sZ+leqRefl (SS n) = SuccLeqSucc $ leqRefl n -leqRefl :: SNat n -> LeqInstance n n-leqRefl SZ = LeqInstance-leqRefl (SS n) =- case leqRefl n of- LeqInstance -> LeqInstance+leqSucc :: SNat n -> Leq n (S n)+leqSucc SZ = ZeroLeq sOne+leqSucc (SS n) = SuccLeqSucc $ leqSucc n -leqSucc :: SNat n -> LeqInstance n (S n)-leqSucc SZ = LeqInstance-leqSucc (SS n) =- case leqSucc n of- LeqInstance -> LeqInstance+data Eql a b where+ Eql :: Eql a a +eqlRefl :: SNat a -> Eql a a+eqlRefl _ = Eql++eqlSymm :: Eql a b -> Eql b a+eqlSymm Eql = Eql++eqlTrans :: Eql a b -> Eql b c -> Eql a c+eqlTrans Eql Eql = Eql++plusZR :: SNat n -> Eql (n :+: Z) n+plusZR SZ = Eql+plusZR (SS n) =+ case plusZR n of+ Eql -> Eql++plusZL :: SNat n -> Eql (Z :+: n) n+plusZL _ = Eql++eqPreservesS :: Eql n m -> Eql (S n) (S m)+eqPreservesS Eql = Eql++plusAssociative :: SNat n -> SNat m -> SNat l+ -> Eql (n :+: (m :+: l)) ((n :+: m) :+: l)+plusAssociative SZ _ _ = Eql+plusAssociative (SS n) m l =+ case plusAssociative n m l of+ Eql -> Eql++sAndPlusOne :: SNat n -> Eql (S n) (n :+: One)+sAndPlusOne SZ = Eql+sAndPlusOne (SS n) =+ case sAndPlusOne n of+ Eql -> Eql++plusCommutative :: SNat n -> SNat m -> Eql (n :+: m) (m :+: n)+plusCommutative SZ SZ = Eql+plusCommutative SZ (SS m) =+ case plusZR (SS m) of+ Eql -> Eql+plusCommutative (SS n) m =+ case plusCommutative n m of+ Eql -> case sAndPlusOne (m %+ n) of+ Eql -> case plusAssociative m n sOne of+ Eql -> case sAndPlusOne n of+ Eql -> Eql++minusCongEq :: Eql n m -> SNat l -> Eql (n :-: l) (m :-: l)+minusCongEq Eql _ = Eql++minusNilpotent :: SNat n -> Eql (n :-: n) Zero+minusNilpotent SZ = Eql+minusNilpotent (SS n) =+ case minusNilpotent n of+ Eql -> Eql++eqSuccMinus :: ((m :<<= n) ~ True)+ => SNat n -> SNat m -> Eql (S n :-: m) (S (n :-: m))+eqSuccMinus _ SZ = Eql+eqSuccMinus (SS n) (SS m) = case eqSuccMinus n m of Eql -> Eql+eqSuccMinus _ _ = error "impossible!"++plusMinusEqL :: SNat n -> SNat m -> Eql ((n :+: m) :-: m) n+plusMinusEqL SZ m = minusNilpotent m+plusMinusEqL (SS n) m =+ case propToBoolLeq (plusLeqR n m) of+ LeqTrueInstance -> eqlTrans (eqSuccMinus (n %+ m) m) (eqPreservesS $ plusMinusEqL n m)++plusMinusEqR :: SNat n -> SNat m -> Eql ((m :+: n) :-: m) n+plusMinusEqR n m = eqlTrans (minusCongEq (plusCommutative n m) m) (plusMinusEqL n m)++data LeqTrueInstance a b where+ LeqTrueInstance :: (a :<<= b) ~ True => LeqTrueInstance a b++propToBoolLeq :: Leq n m -> LeqTrueInstance n m+propToBoolLeq (ZeroLeq _) = LeqTrueInstance+propToBoolLeq (SuccLeqSucc leq) =+ case propToBoolLeq leq of+ LeqTrueInstance -> LeqTrueInstance+++plusLeqL :: SNat n -> SNat m -> Leq n (n :+: m)+plusLeqL SZ m = case plusZR m of Eql -> ZeroLeq m+plusLeqL (SS n) m = SuccLeqSucc $ plusLeqL n m++plusLeqR :: SNat n -> SNat m -> Leq m (n :+: m)+plusLeqR n m =+ case plusCommutative n m of+ Eql -> plusLeqL m n++zAbsorbsMinR :: SNat n -> Eql (Min n Z) Z+zAbsorbsMinR SZ = Eql+zAbsorbsMinR (SS n) =+ case zAbsorbsMinR n of+ Eql -> Eql++zAbsorbsMinL :: SNat n -> Eql (Min Z n) Z+zAbsorbsMinL SZ = Eql+zAbsorbsMinL (SS n) =+ case zAbsorbsMinL n of+ Eql -> Eql++minLeqL :: SNat n -> SNat m -> Leq (Min n m) n+minLeqL SZ m = case zAbsorbsMinL m of Eql -> ZeroLeq sZ+minLeqL n SZ = case zAbsorbsMinR n of Eql -> ZeroLeq n+minLeqL (SS n) (SS m) = SuccLeqSucc (minLeqL n m)++minLeqR :: SNat n -> SNat m -> Leq (Min n m) m+minLeqR n m = case minComm n m of Eql -> minLeqL m n++leqRhs :: Leq n m -> SNat m+leqRhs (ZeroLeq m) = m+leqRhs (SuccLeqSucc leq) = SS $ leqRhs leq++leqLhs :: Leq n m -> SNat n+leqLhs (ZeroLeq _) = SZ+leqLhs (SuccLeqSucc leq) = SS $ leqLhs leq++leqTrans :: Leq n m -> Leq m l -> Leq n l+leqTrans (ZeroLeq _) leq = ZeroLeq $ leqRhs leq+leqTrans (SuccLeqSucc nLeqm) (SuccLeqSucc mLeql) = SuccLeqSucc $ leqTrans nLeqm mLeql++minComm :: SNat n -> SNat m -> Eql (Min n m) (Min m n)+minComm SZ SZ = Eql+minComm SZ (SS _) = Eql+minComm (SS _) SZ = Eql+minComm (SS n) (SS m) = case minComm n m of Eql -> Eql++leqAnitsymmetric :: Leq n m -> Leq m n -> Eql n m+leqAnitsymmetric (ZeroLeq _) (ZeroLeq _) = Eql+leqAnitsymmetric (SuccLeqSucc leq1) (SuccLeqSucc leq2) = eqPreservesS $ leqAnitsymmetric leq1 leq2+leqAnitsymmetric _ _ = error "impossible"++maxZL :: SNat n -> Eql (Max Z n) n+maxZL SZ = Eql+maxZL (SS _) = Eql++maxComm :: SNat n -> SNat m -> Eql (Max n m) (Max m n)+maxComm SZ SZ = Eql+maxComm SZ (SS _) = Eql+maxComm (SS _) SZ = Eql+maxComm (SS n) (SS m) = case maxComm n m of Eql -> Eql++maxZR :: SNat n -> Eql (Max n Z) n+maxZR n = eqlTrans (maxComm n sZ) (maxZL n)++maxLeqL :: SNat n -> SNat m -> Leq n (Max n m)+maxLeqL SZ m = ZeroLeq (sMax sZ m)+maxLeqL n SZ = case maxZR n of+ Eql -> leqRefl n+maxLeqL (SS n) (SS m) = SuccLeqSucc $ maxLeqL n m++maxLeqR :: SNat n -> SNat m -> Leq m (Max n m)+maxLeqR n m = case maxComm n m of+ Eql -> maxLeqL m n+-- (m + S n) - m = S (m + n) - m
Algebra/Ring/Polynomial.hs view
@@ -9,7 +9,7 @@ , IsPolynomial, coeff, lcmMonomial, sPolynomial, polynomial , castMonomial, castPolynomial, toPolynomial, changeOrder , scastMonomial, scastPolynomial, OrderedPolynomial, showPolynomialWithVars- , normalize, injectCoeff, varX, var, getTerms, shiftR+ , normalize, injectCoeff, varX, var, getTerms, shiftR, orderedBy , divs, tryDiv, fromList -- , genVarsV , leadingTerm, leadingMonomial, leadingCoeff, genVars, sDegree , OrderedMonomial(..), Grevlex(..), Revlex(..), Lex(..), Grlex(..)@@ -285,6 +285,9 @@ transformMonomial :: (IsOrder o, IsPolynomial k n, IsPolynomial k m) => (Monomial n -> Monomial m) -> OrderedPolynomial k o n -> OrderedPolynomial k o m transformMonomial trans (Polynomial d) = Polynomial $ M.mapKeys (OrderedMonomial . trans . getMonomial) d++orderedBy :: IsOrder o => OrderedPolynomial k o n -> o -> OrderedPolynomial k o n+p `orderedBy` _ = p shiftR :: forall k r n ord. (Field r, IsPolynomial r n, IsPolynomial r (k :+: n), IsOrder ord) => SNat k -> OrderedPolynomial r ord n -> OrderedPolynomial r ord (k :+: n)
Algebra/Ring/Polynomial/Monomorphic.hs view
@@ -1,124 +1,185 @@-{-# LANGUAGE DataKinds, FlexibleInstances, GADTs, PolyKinds, RecordWildCards #-}-{-# LANGUAGE TypeFamilies, TypeOperators #-}+{-# LANGUAGE DataKinds, FlexibleContexts, FlexibleInstances, GADTs #-}+{-# LANGUAGE MultiParamTypeClasses, PolyKinds, RecordWildCards, TypeFamilies #-}+{-# LANGUAGE TypeOperators, ViewPatterns #-} {-# OPTIONS_GHC -fno-warn-orphans #-} module Algebra.Ring.Polynomial.Monomorphic where-import qualified Algebra.Algorithms.Groebner as Gr import Algebra.Internal import Algebra.Ring.Noetherian-import qualified Algebra.Ring.Polynomial as Poly+import qualified Algebra.Ring.Polynomial as Poly import Control.Arrow import Data.List-import qualified Data.Map as M+import qualified Data.Map as M import Data.Maybe-import Monomorphic+import qualified Numeric.Algebra as NA data Variable = Variable { varName :: Char , varIndex :: Maybe Int } deriving (Eq, Ord) +instance (Eq r, NoetherianRing r, Num r) => Num (Polynomial r) where+ fromInteger n = Polynomial $ M.singleton M.empty $ fromInteger n+ (+) = (NA.+)+ (*) = (NA.*)+ negate = NA.negate+ abs = id+ signum (normalize -> f)+ | f == NA.zero = NA.zero+ | otherwise = NA.one+ instance Show Variable where showsPrec _ v = showChar (varName v) . maybe id ((showChar '_' .) . shows) (varIndex v) -type Polyn = [(Rational, [(Variable, Integer)])]+type Monomial = M.Map Variable Integer -buildVarsList :: Polyn -> [Variable]-buildVarsList = nub . sort . concatMap (map fst . snd)+newtype Polynomial k = Polynomial { unPolynomial :: M.Map Monomial k }+ deriving (Eq, Ord) -encodeMonomList :: [Variable] -> [(Variable, Integer)] -> [Int]-encodeMonomList vars mono = map (maybe 0 fromInteger . flip lookup mono) vars+normalize :: (Eq k, NA.Monoidal k) => Polynomial k -> Polynomial k+normalize (Polynomial dic) =+ Polynomial $ M.filterWithKey (\k v -> v /= NA.zero || M.null k) $ M.mapKeysWith (NA.+) normalizeMonom dic -encodeMonomial :: [Variable] -> [(Variable, Integer)] -> Monomorphic (Vector Int)+normalizeMonom :: Monomial -> Monomial+normalizeMonom = M.filter (/= 0)++instance (Eq r, NoetherianRing r) => NoetherianRing (Polynomial r)+instance (Eq r, NoetherianRing r) => NA.Commutative (Polynomial r)+instance (Eq r, NoetherianRing r) => NA.Multiplicative (Polynomial r) where+ Polynomial (M.toList -> d1) * Polynomial (M.toList -> d2) =+ let dic = [ (M.unionWith (+) a b, r NA.* r') | (a, r) <- d1, (b, r') <- d2 ]+ in normalize $ Polynomial $ M.fromListWith (NA.+) dic++instance (Eq r, NoetherianRing r) => NA.Ring (Polynomial r)+instance (Eq r, NoetherianRing r) => NA.Group (Polynomial r) where+ negate (Polynomial dic) = Polynomial $ fmap NA.negate dic+instance (Eq r, NoetherianRing r) => NA.Rig (Polynomial r)+instance (Eq r, NoetherianRing r) => NA.Unital (Polynomial r) where+ one = Polynomial $ M.singleton M.empty NA.one+instance (Eq r, NoetherianRing r) => NA.Monoidal (Polynomial r) where+ zero = Polynomial $ M.singleton M.empty NA.zero+instance (Eq r, NoetherianRing r) => NA.LeftModule NA.Natural (Polynomial r) where+ n .* Polynomial dic = Polynomial $ fmap (n NA..*) dic +instance (Eq r, NoetherianRing r) => NA.RightModule NA.Natural (Polynomial r) where+ (*.) = flip (NA..*)+instance (Eq r, NoetherianRing r) => NA.LeftModule Integer (Polynomial r) where+ n .* Polynomial dic = Polynomial $ fmap (n NA..*) dic +instance (Eq r, NoetherianRing r) => NA.RightModule Integer (Polynomial r) where+ (*.) = flip (NA..*)+instance (Eq r, NoetherianRing r) => NA.Semiring (Polynomial r)+instance (Eq r, NoetherianRing r) => NA.Abelian (Polynomial r)+instance (Eq r, NoetherianRing r) => NA.Additive (Polynomial r) where+ (Polynomial f) + (Polynomial g) = normalize $ Polynomial $ M.unionWith (NA.+) f g++buildVarsList :: Polynomial r -> [Variable]+buildVarsList = nub . sort . concatMap M.keys . M.keys . unPolynomial++encodeMonomList :: [Variable] -> Monomial -> [Int]+encodeMonomList vars mono = map (maybe 0 fromInteger . flip M.lookup mono) vars++encodeMonomial :: [Variable] -> Monomial -> Monomorphic (Vector Int) encodeMonomial vars mono = promote $ encodeMonomList vars mono -encodePolynomial :: Polyn -> Monomorphic (Poly.Polynomial Rational)+encodePolynomial :: (Monomorphicable (Poly.Polynomial r))+ => Polynomial r -> Monomorphic (Poly.Polynomial r) encodePolynomial = promote . toPolynomialSetting -toPolynomialSetting :: Polyn -> PolynomialSetting+toPolynomialSetting :: Polynomial r -> PolynomialSetting r toPolynomialSetting p = PolySetting { polyn = p , dimension = promote $ length $ buildVarsList p } -data PolynomialSetting = PolySetting { dimension :: Monomorphic SNat- , polyn :: Polyn- } deriving (Show)+data PolynomialSetting r = PolySetting { dimension :: Monomorphic SNat+ , polyn :: Polynomial r+ } deriving (Show) +instance (Eq r, NoetherianRing r, Show r) => Show (Polynomial r) where+ show = showPolynomial -instance Poly.IsMonomialOrder ord => Monomorphicable (Poly.OrderedPolynomial Rational ord) where- type MonomorphicRep (Poly.OrderedPolynomial Rational ord) = PolynomialSetting+instance (Eq r, NoetherianRing r, Poly.IsMonomialOrder ord)+ => Monomorphicable (Poly.OrderedPolynomial r ord) where+ type MonomorphicRep (Poly.OrderedPolynomial r ord) = PolynomialSetting r promote PolySetting{..} = case dimension of Monomorphic dim ->- case singInstance dim of- SingInstance -> Monomorphic $ Poly.polynomial $ M.fromList (map ((Poly.OrderedMonomial . Poly.fromList dim . encodeMonomList vars . snd) &&& fst) polyn)+ case singInstance dim of+ SingInstance -> Monomorphic $ Poly.polynomial $ M.mapKeys (Poly.OrderedMonomial . Poly.fromList dim . encodeMonomList vars) $ unPolynomial polyn where vars = buildVarsList polyn demote (Monomorphic f) =- PolySetting { polyn = map (second $ toMonom . map toInteger . demote . Monomorphic) $ Poly.getTerms f+ PolySetting { polyn = Polynomial $ M.fromList $+ map (toMonom . map toInteger . demote . Monomorphic . snd &&& fst) $ Poly.getTerms f , dimension = Monomorphic $ Poly.sDegree f } where- toMonom = zip $ Variable 'X' Nothing : [Variable 'X' (Just i) | i <- [1..]]+ toMonom = M.fromList . zip (Variable 'X' Nothing : [Variable 'X' (Just i) | i <- [1..]]) -uniformlyPromote :: Poly.IsMonomialOrder ord- => [Polyn] -> Monomorphic (Ideal :.: Poly.OrderedPolynomial Rational ord)-uniformlyPromote ps =- case promote (length vars) of+uniformlyPromoteWithDim :: (Eq r, NoetherianRing r)+ => Poly.IsMonomialOrder ord+ => Int -> [Polynomial r] -> Monomorphic (Ideal :.: Poly.OrderedPolynomial r ord)+uniformlyPromoteWithDim d ps =+ case promote d of Monomorphic dim -> case singInstance dim of- SingInstance -> Monomorphic $ Comp $ toIdeal $ map (Poly.polynomial . M.fromList . map (Poly.OrderedMonomial . Poly.fromList dim . encodeMonomList vars . snd &&& fst)) ps+ SingInstance -> Monomorphic $ Comp $ toIdeal $ map (Poly.polynomial . M.mapKeys (Poly.OrderedMonomial . Poly.fromList dim . encodeMonomList vars) . unPolynomial) ps where vars = nub $ sort $ concatMap buildVarsList ps -instance Poly.IsMonomialOrder ord => Monomorphicable (Ideal :.: Poly.OrderedPolynomial Rational ord) where- type MonomorphicRep (Ideal :.: Poly.OrderedPolynomial Rational ord) = [Polyn]+uniformlyPromote :: (Eq r, NoetherianRing r, Poly.IsMonomialOrder ord)+ => [Polynomial r] -> Monomorphic (Ideal :.: Poly.OrderedPolynomial r ord)+uniformlyPromote ps = uniformlyPromoteWithDim (length vars) ps+ where+ vars = nub $ sort $ concatMap buildVarsList ps++instance (NoetherianRing r, Eq r, Poly.IsMonomialOrder ord)+ => Monomorphicable (Ideal :.: Poly.OrderedPolynomial r ord) where+ type MonomorphicRep (Ideal :.: Poly.OrderedPolynomial r ord) = [Polynomial r] promote = uniformlyPromote demote (Monomorphic (Comp (Ideal v))) = map (polyn . demote . Monomorphic) $ toList v -promoteList :: Poly.IsMonomialOrder ord => [Polyn] -> Monomorphic ([] :.: Poly.OrderedPolynomial Rational ord)-promoteList ps =- case promote (length vars) of- Monomorphic dim ->- case singInstance dim of- SingInstance -> Monomorphic $ Comp $ map (Poly.polynomial . M.fromList . map (Poly.OrderedMonomial . Poly.fromList dim . encodeMonomList vars . snd &&& fst)) ps+promoteList :: (Eq r, NoetherianRing r, Poly.IsMonomialOrder ord)+ => [Polynomial r] -> Monomorphic ([] :.: Poly.OrderedPolynomial r ord)+promoteList ps = promoteListWithDim (length vars) ps where vars = nub $ sort $ concatMap buildVarsList ps --{--data Equal a b where- Equal :: Equal a a--(%==) :: (a ~ b) => a -> b -> Equal a b-_ %== _ = Equal+promoteListWithVarOrder :: (Eq r, NoetherianRing r, Poly.IsMonomialOrder ord)+ => [Variable] -> [Polynomial r] -> Monomorphic ([] :.: Poly.OrderedPolynomial r ord)+promoteListWithVarOrder dic ps =+ case promote dim of+ Monomorphic sdim ->+ case singInstance sdim of+ SingInstance -> Monomorphic $ Comp $ map (Poly.polynomial . M.mapKeys (Poly.OrderedMonomial . Poly.fromList sdim . encodeMonomList vars) . unPolynomial) ps+ where+ vs0 = nub $ sort $ concatMap buildVarsList ps+ (_, rest) = partition (`elem` dic) vs0+ vars = dic ++ rest+ dim = length vars -thEliminationIdeal' :: Int -> [Polyn] -> [Polyn]-thEliminationIdeal' n [] = []-thEliminationIdeal' n ideal =- let dim = length $ nub $ sort $ concatMap buildVarsList ideal- in if n <= 0 || dim <= n- then error "Degree error!"- else case promoteList ideal of- Monomorphic (Comp is@(f:_))->- case singInstance (sDegree f) of- SingInstance ->- case promote n of- Monomorphic sn ->- case sDegree f %== (sn %+ sm) of- Equal -> demote $ Monomorphic $ Comp $ sn `thEliminationIdeal` toIdeal is--}+promoteListWithDim :: (NoetherianRing r, Eq r, Poly.IsMonomialOrder ord)+ => Int -> [Polynomial r] -> Monomorphic ([] :.: Poly.OrderedPolynomial r ord)+promoteListWithDim dim ps =+ case promote dim of+ Monomorphic sdim ->+ case singInstance sdim of+ SingInstance -> Monomorphic $ Comp $ map (Poly.polynomial . M.mapKeys (Poly.OrderedMonomial . Poly.fromList sdim . encodeMonomList vars) . unPolynomial) ps+ where+ vars = nub $ sort $ concatMap buildVarsList ps -renameVars :: [Variable] -> Polyn -> Polyn-renameVars vars = map (second $ map $ first ren)+renameVars :: [Variable] -> Polynomial r -> Polynomial r+renameVars vars = Polynomial . M.mapKeys (M.mapKeys ren) . unPolynomial where ren v = fromMaybe v $ lookup v dic dic = zip (Variable 'X' Nothing : [Variable 'X' (Just i) | i <- [1..]]) vars -showPolyn :: Polyn -> String-showPolyn f =+showPolynomial :: (Show r, Eq r, NoetherianRing r) => Polynomial r -> String+showPolynomial f = case encodePolynomial f of Monomorphic f' -> case singInstance (Poly.sDegree f') of SingInstance -> Poly.showPolynomialWithVars dic f' where dic = zip [1..] $ map show $ buildVarsList f++injectVar :: NA.Unital r => Variable -> Polynomial r+injectVar var = Polynomial $ M.singleton (M.singleton var 1) NA.one+
Algebra/Ring/Polynomial/Parser.hs view
@@ -1,12 +1,14 @@ module Algebra.Ring.Polynomial.Parser where-import Algebra.Ring.Polynomial.Monomorphic-import Control.Applicative hiding (many)-import Control.Arrow-import Data.Char-import Data.Maybe-import Data.Ratio-import Text.Parsec hiding (optional, (<|>))-import Text.Parsec.String+import Algebra.Ring.Polynomial.Monomorphic+import Control.Applicative hiding (many)+import Control.Arrow+import Data.Char+import qualified Data.Map as M+import Data.Maybe+import Data.Ratio+import qualified Numeric.Algebra as NA+import Text.Parsec hiding (optional, (<|>))+import Text.Parsec.String variable :: Parser Variable variable = Variable <$> letter <*> optional (char '_' *> index)@@ -20,18 +22,18 @@ index = digitToInt <$> digit <|> read <$ symbol '{' <*> lexeme (many1 digit) <* symbol '}' -monomial :: Parser [(Variable, Integer)]-monomial = many variableWithPower+monomial :: Parser Monomial+monomial = M.fromList <$> many variableWithPower -term :: Parser (Rational, [(Variable, Integer)])-term = signed' $ try $ (,) <$> option 1 coefficient- <*> monomial- <|> (,) <$> number <*> pure []+term :: Parser (Monomial, Rational)+term = signed' $ try $ flip (,) <$> option 1 coefficient+ <*> monomial+ <|> flip (,) <$> number <*> pure M.empty signed' p = do s <- optional sign- (c, n) <- p- return (fromMaybe 1 s * c, n)+ (n, c) <- p+ return (n, fromMaybe 1 s * c) where sign = lexeme $ char '-' *> return (negate 1) <|> char '+' *> return 1@@ -43,15 +45,15 @@ lexeme :: Parser a -> Parser a lexeme p = p <* spaces -polyOp :: Parser (Polyn -> Polyn -> Polyn)-polyOp = minusPolyn <$ symbol '-'- <|> (++) <$ symbol '+'- where- minusPolyn xs ys = xs ++ map (first negate) ys+toPolyn = normalize . Polynomial . M.fromList -expression :: Parser [(Rational, [(Variable, Integer)])]-expression = spaces *> count 1 term `chainl1` polyOp <* eof+polyOp :: Parser (Polynomial Rational -> Polynomial Rational -> Polynomial Rational)+polyOp = (NA.-) <$ symbol '-'+ <|> (NA.+) <$ symbol '+' +expression :: Parser (Polynomial Rational)+expression = (spaces *> (toPolyn <$> count 1 term) `chainl1` polyOp <* eof)+ coefficient :: Parser Rational coefficient = char '(' *> number <* char ')' <|> number@@ -83,5 +85,5 @@ float <- many1 digit return $ read $ int ++ '.':float -parsePolyn :: String -> Either ParseError Polyn+parsePolyn :: String -> Either ParseError (Polynomial Rational) parsePolyn = parse expression "polynomial"
− Example.hs
@@ -1,59 +0,0 @@-{-# OPTIONS_GHC -fno-warn-type-defaults #-}-{-# LANGUAGE ConstraintKinds, NoImplicitPrelude, TypeOperators #-}-module Example where-import Algebra.Algorithms.Groebner-import Algebra.Internal-import Algebra.Ring.Noetherian-import Algebra.Ring.Polynomial-import Data.Ratio-import Numeric.Algebra-import Prelude hiding (Fractional (..), Integral (..),- Num (..), (^), (^^))--default (Int)--(^^) :: Unital r => r -> Natural -> r-(^^) = pow--x, y, f, f1, f2 :: Polynomial (Ratio Integer) Two-x = var sOne-y = var sTwo-f = x^^2 * y + x * y^^2 + y^^2-f1 = x * y - 1-f2 = y^^2 - 1--type LexPolynomial r n = OrderedPolynomial r Lex n--heron :: Ideal (LexPolynomial (Ratio Integer) (Two :+: Two))-heron = sTwo `thEliminationIdeal` ideal- where- [x, y, a, b, c, s] = genVars (sThree %+ sThree) :: [LexPolynomial (Ratio Integer) (Three :+: Three)]- ideal = toIdeal [ 2 * s - a * y- , b^^2 - (x^^2 + y^^2)- , c^^2 - ( (a-x) ^^ 2 + y^^2)- ]--main :: IO ()-main = do- putStrLn $ unwords ["(" ++ show (x + 1) ++ ")^2", "="- , show $ (x + 1) ^^2 ]- putStrLn $ unwords ["(" ++ show (x + 1) ++ ")(" ++ show (x - 1) ++ ")", "="- , show $ (x + 1) * (x - 1) ]- putStrLn $ unwords ["(" ++ show (x - 1) ++ ")(" ++ show (y^^2 + y - 1) ++ ")", "="- , show $ (x - 1) * (y^^2 + y- 1) ]- putStrLn ""- putStrLn "*** deriving Heron's formula ***"- putStrLn "Area of triangles can be determined from following equations:"- putStrLn "\t2S = ay, b^2 = x^2 + y^2, c^2 = (a-x)^2 + y^2"- putStrLn ", where a, b, c and S stands for three lengths of the traiangle and its area, "- putStrLn "and (x, y) stands for the coordinate of one of its vertices"- putStrLn "(other two vertices are assumed to be on the origin and x-axis)."- putStrLn "Erasing x and y from the equations above, we can get Heron's formula."- putStrLn "Using elimination ideal, this can be automatically solved."- putStrLn "We calculate this with theory of Groebner basis with respect to 'lex'."- putStrLn "This might take a while. please wait..."- print heron- putStrLn "In equation above, X_1, X_2, X_3 and X_4 stands for a, b, c and S, respectively."- putStrLn "The ideal has just one polynomial `f' as its only generator."- putStrLn "Solving the equation `f = 0' assuming S > 0, we can get Heron's formula."-
Monomorphic.hs view
@@ -59,4 +59,3 @@ instance (Read (MonomorphicRep k), Monomorphicable k) => Read (Monomorphic k) where readsPrec i = map (first promote) . readsPrec i #endif-
+ README.md view
@@ -0,0 +1,46 @@+Computational Algebra Library+==============================++Installation+-------------+```{sh}+$ cabal install computational-algebra+```++If you once installed the same version of this package and want to reinstall, please run `cabal clean` first to avoid the GHC's bug.++What is this?+-------------+This library provides data-types and functions to manipulate polynomials.+This is built up with GHC's nice type features.++It contains following things:++* Compute Groebner basis using Buchberger Algorithm+* Ideal membership problem+* Elimination ideal calculation+* Ideal operations+ * Saturation Ideal, Quotient ideal,...++There are two interfaces:++* Dependently-typed I/F+ * Arity-paramaterized polynomials. It uses vector representations for monomials.+ `Algebra.Ring.Polynomial` and `Algebra.Algorithms.Groebner`.++*Monomorphic wrapper I/F+ * Not-so-dependently-typed interface to wrap dependently-typed ones. `Algebra.Ring.Polynomial.Monomorphic` and `Algebra.Algorithms.Groebner.Monomorphic`.+++For more information, please read `examples/polymorphic.hs` and `examples/monomorphic.hs`.++Known Issues+------------+Due to GHC 7.4.*'s bug, this library contains extra modules and functionalities as follows:++* `Monomorphic` data-type and his frieds+ * This is completely separeted as [`monomorphic`](http://hackage.haskell.org/package/monomorphic) package. But due to GHC 7.4.1, which is shipped with latest Haskell Platform, I include the functionality from this library for a while.+* Singleton types and functions+ * Because the [`singletons`](http://hackage.haskell.org/package/singletons) package is not available in GHC 7.4.1, I provide limited version of the functionalities of that package in `Algebra.Internal` module. After new HP released, I will entirely rewrite all source codes using `singletons`.+* Type-level natural numbers and size-parameterized vectors+ * For the similar reason, I include `SNat` and `Vector` data-type in `Algebra.Internal` module, which is separated as [`sized-vector`](http://hackage.haskell.org/package/sized-vector) package. Their proofs are so messy, so I will entirely rewrite these after new HP released with my unreleased package [`equational-reasoning`](https://github.com/konn/equational-reasoning-in-haskell), which provides the functionalities similar to Agda's EqReasoning.
computational-algebra.cabal view
@@ -2,7 +2,7 @@ -- further documentation, see http://haskell.org/cabal/users-guide/ name: computational-algebra-version: 0.0.1.1+version: 0.0.2.0 synopsis: Well-kinded computational algebra library, currently supporting Groebner basis. description: Dependently-typed computational algebra libray for Groebner basis. homepage: https://github.com/konn/computational-algebra@@ -12,6 +12,7 @@ maintainer: konn.jinro_at_gmail.com copyright: (C) Hiromi ISHII 2013 category: Math+extra-source-files: README.md, examples/*.hs build-type: Simple cabal-version: >=1.8 source-repository head@@ -25,7 +26,8 @@ , Algebra.Ring.Polynomial , Algebra.Ring.Polynomial.Monomorphic , Algebra.Ring.Polynomial.Parser- other-modules: Example, Monomorphic, Algebra.Internal+ , Algebra.Internal+ other-modules: Monomorphic build-depends: base >= 2.0 && < 5 , algebra == 3.* , tagged >= 0.4 && < 1
+ examples/monomorphic.hs view
@@ -0,0 +1,84 @@+{-# OPTIONS_GHC -fno-warn-type-defaults #-}+{-# LANGUAGE ConstraintKinds, NoImplicitPrelude, OverloadedStrings #-}+{-# LANGUAGE TypeOperators #-}+module Example where+import Algebra.Algorithms.Groebner.Monomorphic+import Algebra.Ring.Polynomial.Monomorphic+import Algebra.Ring.Polynomial.Parser+import Data.Either+import Data.List (intercalate)+import Data.Ratio+import qualified Data.Text as T+import Numeric.Algebra+import Prelude hiding+ (Fractional (..),+ Integral (..), (*),+ (+), (-), (^), (^^))+import System.IO++default (Int)++(^^) :: Unital r => r -> Natural -> r+(^^) = pow++x, y, f, f1, f2 :: Polynomial (Ratio Integer)+x = injectVar $ Variable 'x' Nothing+y = injectVar $ Variable 'y' Nothing+f = x^^2 * y + x * y^^2 + y^^2+f1 = x * y - one+f2 = y^^2 - one++heron :: [Polynomial (Ratio Integer)]+heron = eliminate [Variable 'x' Nothing, Variable 'y' Nothing] ideal+ where+ [a, b, c, s] = map (injectVar . flip Variable Nothing) "abcS"+ ideal = [ 2 * s - a * y+ , b^^2 - (x^^2 + y^^2)+ , c^^2 - ( (a-x) ^^ 2 + y^^2)+ ]++main :: IO ()+main = do+ putStrLn $ unwords ["(" ++ show (x + 1) ++ ")^2", "="+ , show $ (x + 1) ^^2 ]+ putStrLn $ unwords ["(" ++ show (x + 1) ++ ")(" ++ show (x - 1) ++ ")", "="+ , show $ (x + 1) * (x - 1) ]+ putStrLn $ unwords ["(" ++ show (x - 1) ++ ")(" ++ show (y^^2 + y - 1) ++ ")", "="+ , show $ (x - 1) * (y^^2 + y- 1) ]+ putStrLn "\n==================================================\n"+ idealMembershipDemo+ putStrLn "\n==================================================\n"+ putStrLn ""+ putStrLn "*** deriving Heron's formula ***"+ putStrLn "Area of triangles can be determined from following equations:"+ putStrLn "\t2S = ay, b^2 = x^2 + y^2, c^2 = (a-x)^2 + y^2"+ putStrLn ", where a, b, c and S stands for three lengths of the traiangle and its area, "+ putStrLn "and (x, y) stands for the coordinate of one of its vertices"+ putStrLn "(other two vertices are assumed to be on the origin and x-axis)."+ putStrLn "Erasing x and y from the equations above, we can get Heron's formula."+ putStrLn "Using elimination ideal, this can be automatically solved."+ putStrLn "We calculate this with theory of Groebner basis with respect to 'lex'."+ putStrLn "This might take a while. please wait..."+ print heron+ putStrLn "In equation above, X_1, X_2, X_3 and X_4 stands for a, b, c and S, respectively."+ putStrLn "The ideal has just one polynomial `f' as its only generator."+ putStrLn "Solving the equation `f = 0' assuming S > 0, we can get Heron's formula."++idealMembershipDemo :: IO ()+idealMembershipDemo = do+ putStrLn "======= Ideal Membership Problem ========"+ putStrLn "Enter ideal generators, separetated by comma."+ putStr "enter: "+ hFlush stdout+ src <- getLine+ let (ls, rs) = partitionEithers $ map (parsePolyn . T.unpack) $ T.splitOn "," $ T.pack src+ putStrLn "Enter the polynomial which you want to know whether it's a member of ideal above or not."+ putStr "enter: "+ hFlush stdout+ src <- getLine+ let ex = parsePolyn src+ case (ls, ex) of+ ([], Right f)+ | f `isIdealMember` rs -> putStrLn $ concat ["[YES!] ", show f, " ∈ 〈", intercalate ", " $ map show rs]+ | otherwise -> putStrLn $ concat ["[NO!] ", show f, " ∉ 〈", intercalate ", " $ map show rs]+ _ -> putStrLn "Parse error! try again." >> idealMembershipDemo
+ examples/polymorphic.hs view
@@ -0,0 +1,59 @@+{-# OPTIONS_GHC -fno-warn-type-defaults #-}+{-# LANGUAGE ConstraintKinds, NoImplicitPrelude, TypeOperators #-}+module Example where+import Algebra.Algorithms.Groebner+import Algebra.Internal+import Algebra.Ring.Noetherian+import Algebra.Ring.Polynomial+import Data.Ratio+import Numeric.Algebra+import Prelude hiding (Fractional (..), Integral (..),+ Num (..), (^), (^^))++default (Int)++(^^) :: Unital r => r -> Natural -> r+(^^) = pow++x, y, f, f1, f2 :: Polynomial (Ratio Integer) Two+x = var sOne+y = var sTwo+f = x^^2 * y + x * y^^2 + y^^2+f1 = x * y - 1+f2 = y^^2 - 1++type LexPolynomial r n = OrderedPolynomial r Lex n++heron :: Ideal (LexPolynomial (Ratio Integer) (Two :+: Two))+heron = sTwo `thEliminationIdeal` ideal+ where+ [x, y, a, b, c, s] = genVars (sThree %+ sThree) :: [LexPolynomial (Ratio Integer) (Three :+: Three)]+ ideal = toIdeal [ 2 * s - a * y+ , b^^2 - (x^^2 + y^^2)+ , c^^2 - ( (a-x) ^^ 2 + y^^2)+ ]++main :: IO ()+main = do+ putStrLn $ unwords ["(" ++ show (x + 1) ++ ")^2", "="+ , show $ (x + 1) ^^2 ]+ putStrLn $ unwords ["(" ++ show (x + 1) ++ ")(" ++ show (x - 1) ++ ")", "="+ , show $ (x + 1) * (x - 1) ]+ putStrLn $ unwords ["(" ++ show (x - 1) ++ ")(" ++ show (y^^2 + y - 1) ++ ")", "="+ , show $ (x - 1) * (y^^2 + y- 1) ]+ putStrLn ""+ putStrLn "*** deriving Heron's formula ***"+ putStrLn "Area of triangles can be determined from following equations:"+ putStrLn "\t2S = ay, b^2 = x^2 + y^2, c^2 = (a-x)^2 + y^2"+ putStrLn ", where a, b, c and S stands for three lengths of the traiangle and its area, "+ putStrLn "and (x, y) stands for the coordinate of one of its vertices"+ putStrLn "(other two vertices are assumed to be on the origin and x-axis)."+ putStrLn "Erasing x and y from the equations above, we can get Heron's formula."+ putStrLn "Using elimination ideal, this can be automatically solved."+ putStrLn "We calculate this with theory of Groebner basis with respect to 'lex'."+ putStrLn "This might take a while. please wait..."+ print heron+ putStrLn "In equation above, X_1, X_2, X_3 and X_4 stands for a, b, c and S, respectively."+ putStrLn "The ideal has just one polynomial `f' as its only generator."+ putStrLn "Solving the equation `f = 0' assuming S > 0, we can get Heron's formula."+