computational-algebra 0.2.0.0 → 0.3.0.0
raw patch · 6 files changed
+87/−49 lines, 6 filesdep ~type-natural
Dependency ranges changed: type-natural
Files
- Algebra/Algorithms/Groebner.hs +36/−17
- Algebra/Algorithms/Groebner/Monomorphic.hs +27/−7
- Algebra/Ring/Noetherian.hs +4/−5
- Algebra/Ring/Polynomial.hs +15/−13
- Algebra/Ring/Polynomial/Monomorphic.hs +3/−5
- computational-algebra.cabal +2/−2
Algebra/Algorithms/Groebner.hs view
@@ -1,7 +1,7 @@-{-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts, FlexibleInstances #-}-{-# LANGUAGE GADTs, MultiParamTypeClasses, NoImplicitPrelude #-}-{-# LANGUAGE ParallelListComp, RankNTypes, ScopedTypeVariables #-}-{-# LANGUAGE StandaloneDeriving, TemplateHaskell, TypeFamilies, TypeOperators #-}+{-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts, FlexibleInstances #-}+{-# LANGUAGE GADTs, MultiParamTypeClasses, NoImplicitPrelude #-}+{-# LANGUAGE ParallelListComp, RankNTypes, ScopedTypeVariables #-}+{-# LANGUAGE TemplateHaskell, TypeFamilies, TypeOperators #-} {-# OPTIONS_GHC -fno-warn-type-defaults -fno-warn-orphans #-} module Algebra.Algorithms.Groebner ( -- * Polynomial division@@ -20,6 +20,8 @@ , unsafeThEliminationIdealWith , quotIdeal, quotByPrincipalIdeal , saturationIdeal, saturationByPrincipalIdeal+ -- * Resultant+ , resultant, hasCommonFactor ) where import Algebra.Internal import Algebra.Ring.Noetherian@@ -33,15 +35,14 @@ import qualified Data.Heap as H import Data.List import Data.Maybe-import Data.Proxy import Data.STRef import Data.Type.Monomorphic import Data.Type.Natural hiding (max, one, zero) import Data.Vector.Sized hiding (all, drop, foldr, head, map, take, zipWith) import qualified Data.Vector.Sized as V-import Numeric.Algebra-import Prelude hiding (Num (..), recip)+import Numeric.Algebra hiding ((>))+import Prelude hiding (Num (..), recip, (^)) import Proof.Equational -- | Calculate a polynomial quotient and remainder w.r.t. second argument.@@ -162,7 +163,7 @@ syzygyBuchbergerWithStrategy strategy ideal = runST $ do let gens = zip [1..] $ generators ideal gs <- newSTRef $ H.fromList [H.Entry (leadingOrderedMonomial g) g | (_, g) <- gens]- b <- newSTRef $ H.fromList $ [H.Entry (calcWeight' strategy f g, j) (f, g) | ((_, f), (j, g)) <- combinations gens]+ b <- newSTRef $ H.fromList [H.Entry (calcWeight' strategy f g, j) (f, g) | ((_, f), (j, g)) <- combinations gens] len <- newSTRef (genericLength gens :: Integer) whileM_ (not . H.null <$> readSTRef b) $ do Just (H.Entry _ (f, g), rest) <- H.viewMin <$> readSTRef b@@ -237,16 +238,15 @@ minimizeGroebnerBasis :: (Field k, IsPolynomial k n, IsMonomialOrder order) => [OrderedPolynomial k order n] -> [OrderedPolynomial k order n] minimizeGroebnerBasis bs = runST $ do- left <- newSTRef bs+ left <- newSTRef $ map monoize $ filter (/= zero) bs right <- newSTRef [] whileM_ (not . null <$> readSTRef left) $ do f : xs <- readSTRef left writeSTRef left xs ys <- readSTRef right- if any (\g -> leadingMonomial g `divs` leadingMonomial f) xs ||- any (\g -> leadingMonomial g `divs` leadingMonomial f) ys- then writeSTRef right ys- else writeSTRef right (monoize f : ys)+ unless (any (\g -> leadingMonomial g `divs` leadingMonomial f) xs+ || any (\g -> leadingMonomial g `divs` leadingMonomial f) ys) $+ writeSTRef right (f : ys) readSTRef right -- | Reduce minimum Groebner basis into reduced Groebner basis.@@ -263,9 +263,6 @@ if q == zero then writeSTRef right ys else writeSTRef right (q : ys) readSTRef right --- foldr step [] [f, g, h]--- f `step` (g `step` (h `step` []))- monoize :: (Field k, IsPolynomial k n, IsMonomialOrder order) => OrderedPolynomial k order n -> OrderedPolynomial k order n monoize f = injectCoeff (recip $ leadingCoeff f) * f@@ -381,7 +378,7 @@ => Ideal (OrderedPolynomial k ord n) -> OrderedPolynomial k ord n -> Ideal (OrderedPolynomial k ord n) saturationByPrincipalIdeal is g =- case propToClassLeq $ leqSucc (sDegree g) of+ case propToClassLeq $ leqSucc (sArity g) of LeqInstance -> thEliminationIdeal sOne $ addToIdeal (one - (castPolynomial g * var sOne)) (mapIdeal (shiftR sOne) is) -- | Saturation ideal@@ -394,3 +391,25 @@ SingInstance -> case singInstance (sLength g %+ (sing :: SNat n)) of SingInstance -> intersection $ V.map (i `saturationByPrincipalIdeal`) g++-- | Calculate resultant for given two unary polynomimals.+resultant :: forall k ord . (Eq k, NoetherianRing k, Field k, IsMonomialOrder ord)+ => OrderedPolynomial k ord One+ -> OrderedPolynomial k ord One+ -> k+resultant = go one+ where+ go res h s+ | totalDegree' s > 0 = let r = h `modPolynomial` [s]+ res' = res * negate one ^ (totalDegree' h * totalDegree' s)+ * (leadingCoeff s) ^ (totalDegree' h - totalDegree' r)+ in go res' s r+ | h == zero || s == zero = zero+ | totalDegree' h > 0 = (leadingCoeff s ^ totalDegree' h) * res+ | otherwise = res++hasCommonFactor :: forall k ord . (NoetherianRing k, Eq k, Field k, IsMonomialOrder ord)+ => OrderedPolynomial k ord One+ -> OrderedPolynomial k ord One+ -> Bool+hasCommonFactor f g = resultant f g == zero
Algebra/Algorithms/Groebner/Monomorphic.hs view
@@ -17,6 +17,8 @@ , isIdealMember, intersection, thEliminationIdeal, eliminate, thEliminationIdealWith, eliminateWith , quotIdeal, quotByPrincipalIdeal , saturationIdeal, saturationByPrincipalIdeal+ -- * Resultant+ , resultant, hasCommonFactor -- * Re-exports , Lex(..), Revlex(..), Grlex(..), Grevlex(..), IsOrder(..), IsMonomialOrder , SelectionStrategy(..), NormalStrategy(..), SugarStrategy(..), Gr.GrevlexStrategy(..)@@ -108,7 +110,7 @@ divModPolynomialWith _ f gs = case promoteList (f:gs) :: Monomorphic ([] :.: Poly.OrderedPolynomial r ord) of Monomorphic (Comp (f' : gs')) ->- let sn = Poly.sDegree f'+ let sn = Poly.sArity f' in case singInstance sn of SingInstance -> let (q, r) = Gr.divModPolynomial f' gs'@@ -140,7 +142,7 @@ Monomorphic (Comp ideal) -> case ideal of Ideal vec ->- case singInstance (Poly.sDegree (head $ toList vec)) of+ case singInstance (Poly.sArity (head $ toList vec)) of SingInstance -> map (renameVars vars . polyn . demote . Monomorphic) $ Gr.calcGroebnerBasisWith ord ideal where vars = nub $ sort $ concatMap buildVarsList j@@ -156,7 +158,7 @@ Monomorphic (Comp ideal) -> case ideal of Ideal vec ->- case singInstance (Poly.sDegree (head $ toList vec)) of+ case singInstance (Poly.sArity (head $ toList vec)) of SingInstance -> map (renameVars vars . polyn . demote . Monomorphic) $ Gr.simpleBuchberger ideal where vars = nub $ sort $ concatMap buildVarsList j@@ -172,7 +174,7 @@ Monomorphic (Comp ideal) -> case ideal of Ideal vec ->- case singInstance (Poly.sDegree (head $ toList vec)) of+ case singInstance (Poly.sArity (head $ toList vec)) of SingInstance -> map (renameVars vars . polyn . demote . Monomorphic) $ Gr.primeTestBuchberger ideal where vars = nub $ sort $ concatMap buildVarsList j@@ -190,7 +192,7 @@ Monomorphic (Comp ideal) -> case ideal of Ideal vec ->- case singInstance (Poly.sDegree (head $ toList vec)) of+ case singInstance (Poly.sArity (head $ toList vec)) of SingInstance -> map (renameVars vars . polyn . demote . Monomorphic) $ Gr.syzygyBuchbergerWithStrategy strategy ideal where vars = nub $ sort $ concatMap buildVarsList j@@ -205,7 +207,7 @@ isIdealMember f ideal = case promoteList (f:ideal) :: Monomorphic ([] :.: Poly.Polynomial r) of Monomorphic (Comp (f':ideal')) ->- case singInstance (Poly.sDegree f') of+ case singInstance (Poly.sArity f') of SingInstance -> Gr.isIdealMember f' (toIdeal ideal') _ -> error "impossible happend!" @@ -217,7 +219,7 @@ Monomorphic (Comp fs) -> case promote k of Monomorphic sk ->- let sdim = Poly.sDegree $ head fs+ let sdim = Poly.sArity $ head fs newDim = sMax sk sdim in case singInstance sdim of SingInstance ->@@ -247,3 +249,21 @@ thEliminationIdealWith ord k j = eliminateWith ord (take k vars) j where vars = nub $ sort $ concatMap buildVarsList j++-- | Calculates resultants for given two unary-polynomials.+resultant :: forall r. Groebnerable r+ => Polynomial r -> Polynomial r -> r+resultant f g =+ let vars = nub $ buildVarsList f ++ buildVarsList g+ in case vars of+ [_] ->+ let f' = Poly.polynomial $ M.mapKeys (Poly.OrderedMonomial . Poly.fromList sOne . encodeMonomList vars) $+ unPolynomial f+ g' = Poly.polynomial $ M.mapKeys (Poly.OrderedMonomial . Poly.fromList sOne . encodeMonomList vars) $+ unPolynomial g+ in Gr.resultant (f' `orderedBy` Grevlex) g'+ _ -> error "currently supports only unary polynomial."++-- | Determin if given two unary polynomials have common factor.+hasCommonFactor :: (Eq r, Division r, NoetherianRing r) => Polynomial r -> Polynomial r -> Bool+hasCommonFactor f g = resultant f g == zero
Algebra/Ring/Noetherian.hs view
@@ -4,20 +4,17 @@ {-# OPTIONS_GHC -fno-warn-orphans #-} module Algebra.Ring.Noetherian ( NoetherianRing, Ideal(..), addToIdeal, toIdeal, appendIdeal , generators, filterIdeal, mapIdeal, principalIdeal) where-import Algebra.Internal import qualified Data.Complex as C import Data.Function import Data.Ord import Data.Ratio+import Data.Vector.Sized (Vector (..))+import qualified Data.Vector.Sized as V import Numeric.Algebra-import qualified Numeric.Algebra as NA import qualified Numeric.Algebra.Complex as NA import Prelude hiding (negate, subtract, (*), (+), (-)) import qualified Prelude as P-import qualified Data.Vector.Sized as V-import Data.Vector.Sized (Vector(..))-import Data.Type.Natural class (Commutative r, Ring r) => NoetherianRing r where @@ -36,6 +33,8 @@ instance Integral n => TriviallyInvolutive (Ratio n) instance (P.Num n) => P.Num (NA.Complex n) where+ abs = error "unimplemented"+ signum = error "unimplemented" fromInteger n = NA.Complex (P.fromInteger n) 0 negate (NA.Complex x y) = NA.Complex (P.negate x) (P.negate y) NA.Complex x y + NA.Complex z w = NA.Complex (x P.+ y) (z P.+ w)
Algebra/Ring/Polynomial.hs view
@@ -9,13 +9,13 @@ ( Polynomial, Monomial, MonomialOrder, EliminationType, EliminationOrder , WeightedEliminationOrder, eliminationOrder, weightedEliminationOrder , lex, revlex, graded, grlex, grevlex, productOrder, productOrder'- , transformMonomial, WeightProxy(..), weightOrder, totalDegree+ , transformMonomial, WeightProxy(..), weightOrder, totalDegree, totalDegree' , IsPolynomial, coeff, lcmMonomial, sPolynomial, polynomial , castMonomial, castPolynomial, toPolynomial, changeOrder, changeOrderProxy , scastMonomial, scastPolynomial, OrderedPolynomial, showPolynomialWithVars, showPolynomialWith, showRational , normalize, injectCoeff, varX, var, getTerms, shiftR, orderedBy , divs, tryDiv, fromList, Coefficient(..),ToWeightVector(..)- , leadingTerm, leadingMonomial, leadingOrderedMonomial, leadingCoeff, genVars, sDegree+ , leadingTerm, leadingMonomial, leadingOrderedMonomial, leadingCoeff, genVars, sArity , OrderedMonomial(..), OrderedMonomial'(..), Grevlex(..) , Revlex(..), Lex(..), Grlex(..), Graded(..) , ProductOrder (..), WeightOrder(..)@@ -31,7 +31,6 @@ import Data.Maybe import Data.Monoid import Data.Ord-import Data.Proxy import Data.Ratio import Data.Type.Monomorphic import Data.Type.Natural hiding (max, one, promote, zero)@@ -90,6 +89,9 @@ totalDegree = V.foldl (+) 0 {-# INLINE totalDegree #-} +totalDegree' :: OrderedPolynomial k ord n -> Int+totalDegree' = maximum . (0:) . map (totalDegree . snd) . getTerms+ -- | Lexicographical order. This *is* a monomial order. lex :: MonomialOrder lex Nil Nil = EQ@@ -295,7 +297,7 @@ castMonomial = unwrapped %~ fromList sing . V.toList scastMonomial :: (n :<= m) => SNat m -> OrderedMonomial o n -> OrderedMonomial o m-scastMonomial snat = unwrapped %~ fromList snat . V.toList+scastMonomial sdim = unwrapped %~ fromList sdim . V.toList castPolynomial :: (IsPolynomial r n, IsPolynomial r m, SingRep m, IsOrder o, IsOrder o', n :<= m) => OrderedPolynomial r o n@@ -503,14 +505,14 @@ seed = cycle $ 1 : replicate (n - 1) 0 in map (\m -> Polynomial $ M.singleton (OrderedMonomial $ fromList sn $ take n (drop (n-m) seed)) one) [0..n-1] -sDegree :: OrderedPolynomial k ord n -> SNat n-sDegree (Polynomial dic) = V.sLength $ getMonomial $ fst $ M.findMin dic+sArity :: OrderedPolynomial k ord n -> SNat n+sArity (Polynomial dic) = V.sLength $ getMonomial $ fst $ M.findMin dic {-# RULES-"sDegree/zero" forall (v :: OrderedPolynomial k ord Z). sDegree v = SZ-"sDegree/one" forall (v :: OrderedPolynomial k ord (S Z)). sDegree v = SS SZ-"sDegree/two" forall (v :: OrderedPolynomial k ord (S (S Z))). sDegree v = SS (SS SZ)-"sDegree/three" forall (v :: OrderedPolynomial k ord (S (S (S Z)))). sDegree v = SS (SS (sS SZ))-"sDegree/four" forall (v :: OrderedPolynomial k ord (S (S (S (S Z))))). sDegree v = SS (SS (SS (SS SZ)))-"sDegree/five" forall (v :: OrderedPolynomial k ord (S (S (S (S (S Z)))))). sDegree v = SS (SS (SS (SS (SS SZ))))-"sDegree/sing" forall (v :: SingRep n => OrderedPolynomial k ord n). sDegree (v :: OrderedPolynomial k ord n) = sing :: SNat n+"sArity/zero" forall (v :: OrderedPolynomial k ord Z). sArity v = SZ+"sArity/one" forall (v :: OrderedPolynomial k ord (S Z)). sArity v = SS SZ+"sArity/two" forall (v :: OrderedPolynomial k ord (S (S Z))). sArity v = SS (SS SZ)+"sArity/three" forall (v :: OrderedPolynomial k ord (S (S (S Z)))). sArity v = SS (SS (sS SZ))+"sArity/four" forall (v :: OrderedPolynomial k ord (S (S (S (S Z))))). sArity v = SS (SS (SS (SS SZ)))+"sArity/five" forall (v :: OrderedPolynomial k ord (S (S (S (S (S Z)))))). sArity v = SS (SS (SS (SS (SS SZ))))+"sArity/sing" forall (v :: SingRep n => OrderedPolynomial k ord n). sArity (v :: OrderedPolynomial k ord n) = sing :: SNat n #-}
Algebra/Ring/Polynomial/Monomorphic.hs view
@@ -3,7 +3,6 @@ {-# LANGUAGE TypeOperators, ViewPatterns, OverlappingInstances #-} {-# OPTIONS_GHC -fno-warn-orphans #-} module Algebra.Ring.Polynomial.Monomorphic where-import Algebra.Internal import Algebra.Ring.Noetherian import qualified Algebra.Ring.Polynomial as Poly import Control.Arrow@@ -14,7 +13,6 @@ import Data.Type.Monomorphic import qualified Numeric.Algebra as NA import Data.Ratio-import Data.Singletons hiding (promote) import qualified Data.Vector.Sized as V data Variable = Variable { varName :: Char@@ -116,7 +114,7 @@ demote (Monomorphic f) = PolySetting { polyn = Polynomial $ M.fromList $ map (toMonom . map toInteger . demote . Monomorphic . snd &&& fst) $ Poly.getTerms f- , dimension = Monomorphic $ Poly.sDegree f+ , dimension = Monomorphic $ Poly.sArity f } where toMonom = M.fromList . zip (Variable 'X' Nothing : [Variable 'X' (Just i) | i <- [1..]])@@ -183,7 +181,7 @@ showPolynomial f = case encodePolynomial f of Monomorphic f' ->- case singInstance (Poly.sDegree f') of+ case singInstance (Poly.sArity f') of SingInstance -> Poly.showPolynomialWithVars dic f' where dic = zip [1 :: Int ..] $ map show $ buildVarsList f@@ -192,7 +190,7 @@ showRatPolynomial f = case encodePolynomial f of Monomorphic f' ->- case singInstance (Poly.sDegree f') of+ case singInstance (Poly.sArity f') of SingInstance -> Poly.showPolynomialWith dic Poly.showRational f' where dic = zip [1 :: Int ..] $ map show $ buildVarsList f
computational-algebra.cabal view
@@ -2,7 +2,7 @@ -- further documentation, see http://haskell.org/cabal/users-guide/ name: computational-algebra-version: 0.2.0.0+version: 0.3.0.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@@ -36,7 +36,7 @@ , peggy == 0.3.* , monad-loops >= 0.3 && <0.5 , heaps == 0.2.*- , type-natural >= 0.0.2.1+ , type-natural == 0.0.3.* , sized-vector == 0.0.* , singletons >= 0.8 , equational-reasoning == 0.0.*