diff --git a/Algebra/Algorithms/Groebner.hs b/Algebra/Algorithms/Groebner.hs --- a/Algebra/Algorithms/Groebner.hs +++ b/Algebra/Algorithms/Groebner.hs @@ -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 diff --git a/Algebra/Algorithms/Groebner/Monomorphic.hs b/Algebra/Algorithms/Groebner/Monomorphic.hs --- a/Algebra/Algorithms/Groebner/Monomorphic.hs +++ b/Algebra/Algorithms/Groebner/Monomorphic.hs @@ -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 diff --git a/Algebra/Ring/Noetherian.hs b/Algebra/Ring/Noetherian.hs --- a/Algebra/Ring/Noetherian.hs +++ b/Algebra/Ring/Noetherian.hs @@ -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) diff --git a/Algebra/Ring/Polynomial.hs b/Algebra/Ring/Polynomial.hs --- a/Algebra/Ring/Polynomial.hs +++ b/Algebra/Ring/Polynomial.hs @@ -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 #-} diff --git a/Algebra/Ring/Polynomial/Monomorphic.hs b/Algebra/Ring/Polynomial/Monomorphic.hs --- a/Algebra/Ring/Polynomial/Monomorphic.hs +++ b/Algebra/Ring/Polynomial/Monomorphic.hs @@ -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 diff --git a/computational-algebra.cabal b/computational-algebra.cabal --- a/computational-algebra.cabal +++ b/computational-algebra.cabal @@ -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.*