-- Copyright (c) David Amos, 2011. All rights reserved.

{-# LANGUAGE TupleSections, NoMonomorphismRestriction #-}

-- |A module providing an efficient implementation of the Buchberger algorithm for calculating the (reduced) Groebner basis for an ideal,
-- together with some straightforward applications.
module Math.CommutativeAlgebra.GroebnerBasis where

import Data.List as L
import qualified Data.IntMap as IM
import qualified Data.Set as S

import Math.Core.Utils
import Math.Core.Field
import Math.Algebras.VectorSpace
import Math.Algebras.Structures
import Math.CommutativeAlgebra.Polynomial

-- Sources:
-- Cox, Little, O'Shea: Ideals, Varieties and Algorithms
-- Giovini, Mora, Niesi, Robbiano, Traverso, "One sugar cube please, or Selection strategies in the Buchberger algorithm"


sPoly :: Vect b3 b -> Vect b3 b -> Vect b3 b
sPoly f :: Vect b3 b
f g :: Vect b3 b
g = let d :: (b, b1)
d = (b, b3) -> (b, b3) -> (b, b1)
forall a b1 b2 b3.
(Monomial a, Num b1) =>
(a, b2) -> (a, b3) -> (a, b1)
tgcd (Vect b3 b -> (b, b3)
forall k b. Vect k b -> (b, k)
lt Vect b3 b
f) (Vect b3 b -> (b, b3)
forall k b. Vect k b -> (b, k)
lt Vect b3 b
g)
            in (Vect b3 b -> (b, b3)
forall k b. Vect k b -> (b, k)
lt Vect b3 b
g (b, b3) -> (b, b3) -> (b, b3)
forall a b.
(Monomial a, Fractional b) =>
(a, b) -> (a, b) -> (a, b)
`tdiv` (b, b3)
forall b1. Num b1 => (b, b1)
d) (b, b3) -> Vect b3 b -> Vect b3 b
forall b k. (Mon b, Num k) => (b, k) -> Vect k b -> Vect k b
*-> Vect b3 b
f Vect b3 b -> Vect b3 b -> Vect b3 b
forall a. Num a => a -> a -> a
- (Vect b3 b -> (b, b3)
forall k b. Vect k b -> (b, k)
lt Vect b3 b
f (b, b3) -> (b, b3) -> (b, b3)
forall a b.
(Monomial a, Fractional b) =>
(a, b) -> (a, b) -> (a, b)
`tdiv` (b, b3)
forall b1. Num b1 => (b, b1)
d) (b, b3) -> Vect b3 b -> Vect b3 b
forall b k. (Mon b, Num k) => (b, k) -> Vect k b -> Vect k b
*-> Vect b3 b
g
-- The point about the s-poly is that it cancels out the leading terms of the two polys, exposing their second terms


isGB :: [Vect b3 b] -> Bool
isGB fs :: [Vect b3 b]
fs = (Vect b3 b -> Bool) -> [Vect b3 b] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\h :: Vect b3 b
h -> Vect b3 b
h Vect b3 b -> [Vect b3 b] -> Vect b3 b
forall k m.
(Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
Vect k m -> [Vect k m] -> Vect k m
%% [Vect b3 b]
fs Vect b3 b -> Vect b3 b -> Bool
forall a. Eq a => a -> a -> Bool
== 0) ((Vect b3 b -> Vect b3 b -> Vect b3 b) -> [Vect b3 b] -> [Vect b3 b]
forall a a. (a -> a -> a) -> [a] -> [a]
pairWith Vect b3 b -> Vect b3 b -> Vect b3 b
forall b3 b.
(Eq b3, Ord b, Algebra b3 b, Fractional b3, Monomial b) =>
Vect b3 b -> Vect b3 b -> Vect b3 b
sPoly [Vect b3 b]
fs)


-- Initial, naive version
-- Cox p87
gb1 :: [Vect b3 b] -> [Vect b3 b]
gb1 fs :: [Vect b3 b]
fs = [Vect b3 b] -> [Vect b3 b] -> [Vect b3 b]
forall b b3.
(Fractional b3, Monomial b, Ord b, Algebra b3 b, Eq b3) =>
[Vect b3 b] -> [Vect b3 b] -> [Vect b3 b]
gb' [Vect b3 b]
fs ((Vect b3 b -> Vect b3 b -> Vect b3 b) -> [Vect b3 b] -> [Vect b3 b]
forall a a. (a -> a -> a) -> [a] -> [a]
pairWith Vect b3 b -> Vect b3 b -> Vect b3 b
forall b3 b.
(Eq b3, Ord b, Algebra b3 b, Fractional b3, Monomial b) =>
Vect b3 b -> Vect b3 b -> Vect b3 b
sPoly [Vect b3 b]
fs) where
    gb' :: [Vect b3 b] -> [Vect b3 b] -> [Vect b3 b]
gb' gs :: [Vect b3 b]
gs (h :: Vect b3 b
h:hs :: [Vect b3 b]
hs) = let h' :: Vect b3 b
h' = Vect b3 b
h Vect b3 b -> [Vect b3 b] -> Vect b3 b
forall k m.
(Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
Vect k m -> [Vect k m] -> Vect k m
%% [Vect b3 b]
gs in
                    if Vect b3 b
h' Vect b3 b -> Vect b3 b -> Bool
forall a. Eq a => a -> a -> Bool
== 0 then [Vect b3 b] -> [Vect b3 b] -> [Vect b3 b]
gb' [Vect b3 b]
gs [Vect b3 b]
hs else [Vect b3 b] -> [Vect b3 b] -> [Vect b3 b]
gb' (Vect b3 b
h'Vect b3 b -> [Vect b3 b] -> [Vect b3 b]
forall a. a -> [a] -> [a]
:[Vect b3 b]
gs) ([Vect b3 b]
hs [Vect b3 b] -> [Vect b3 b] -> [Vect b3 b]
forall a. [a] -> [a] -> [a]
++ (Vect b3 b -> Vect b3 b) -> [Vect b3 b] -> [Vect b3 b]
forall a b. (a -> b) -> [a] -> [b]
map (Vect b3 b -> Vect b3 b -> Vect b3 b
forall b3 b.
(Eq b3, Ord b, Algebra b3 b, Fractional b3, Monomial b) =>
Vect b3 b -> Vect b3 b -> Vect b3 b
sPoly Vect b3 b
h') [Vect b3 b]
gs)
    gb' gs :: [Vect b3 b]
gs [] = [Vect b3 b]
gs

-- [f xi xj | xi <- xs, xj <- xs, i < j]
pairWith :: (a -> a -> a) -> [a] -> [a]
pairWith f :: a -> a -> a
f (x :: a
x:xs :: [a]
xs) = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
f a
x) [a]
xs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ (a -> a -> a) -> [a] -> [a]
pairWith a -> a -> a
f [a]
xs
pairWith _ [] = []


-- Cox p89-90
reduce :: [Vect k b] -> [Vect k b]
reduce gs :: [Vect k b]
gs = [Vect k b] -> [Vect k b] -> [Vect k b]
forall b k.
(Ord k, Fractional k, Monomial b, Ord b, Algebra k b) =>
[Vect k b] -> [Vect k b] -> [Vect k b]
reduce' [Vect k b]
gs [] where
    reduce' :: [Vect k b] -> [Vect k b] -> [Vect k b]
reduce' (l :: Vect k b
l:ls :: [Vect k b]
ls) rs :: [Vect k b]
rs = let l' :: Vect k b
l' = Vect k b
l Vect k b -> [Vect k b] -> Vect k b
forall k m.
(Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
Vect k m -> [Vect k m] -> Vect k m
%% ([Vect k b]
rs [Vect k b] -> [Vect k b] -> [Vect k b]
forall a. [a] -> [a] -> [a]
++ [Vect k b]
ls) in
                        if Vect k b
l' Vect k b -> Vect k b -> Bool
forall a. Eq a => a -> a -> Bool
== 0 then [Vect k b] -> [Vect k b] -> [Vect k b]
reduce' [Vect k b]
ls [Vect k b]
rs else [Vect k b] -> [Vect k b] -> [Vect k b]
reduce' [Vect k b]
ls (Vect k b -> Vect k b
forall b k.
(Eq k, Ord b, Show b, Algebra k b, Fractional k) =>
Vect k b -> Vect k b
toMonic Vect k b
l'Vect k b -> [Vect k b] -> [Vect k b]
forall a. a -> [a] -> [a]
:[Vect k b]
rs)
    reduce' [] rs :: [Vect k b]
rs = [Vect k b] -> [Vect k b]
forall a. Ord a => [a] -> [a]
L.sort [Vect k b]
rs
-- when using an elimination order, the elimination ideal will be at the end

-- Cox et al p106-7
-- No need to calculate an spoly fi fj if
-- 1. the lm fi and lm fj are coprime, or
-- 2. there exists some fk, with (i,k) (j,k) already considered, and lm fk divides lcm (lm fi) (lm fj) 
-- some slight inefficiencies from looking up fi, fj repeatedly
gb2 :: [Vect k b] -> [Vect k b]
gb2 fs :: [Vect k b]
fs = [Vect k b] -> [Vect k b]
forall k b.
(Fractional k, Monomial b, Ord k, Ord b, Algebra k b) =>
[Vect k b] -> [Vect k b]
reduce ([Vect k b] -> [Vect k b]) -> [Vect k b] -> [Vect k b]
forall a b. (a -> b) -> a -> b
$ [Vect k b] -> [(Int, Int)] -> Int -> [Vect k b]
forall m b.
(Fractional b, Monomial m, Ord m, Algebra b m, Eq b) =>
[Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
gb' [Vect k b]
fs ([Int] -> [(Int, Int)]
forall t. [t] -> [(t, t)]
pairs [1..Int
s]) Int
s where
    s :: Int
s = [Vect k b] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Vect k b]
fs
    gb' :: [Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
gb' gs :: [Vect b m]
gs ((i :: Int
i,j :: Int
j):ps :: [(Int, Int)]
ps) t :: Int
t =
        let fi :: Vect b m
fi = [Vect b m]
gs[Vect b m] -> Int -> Vect b m
forall a. [a] -> Int -> a
!Int
i; fj :: Vect b m
fj = [Vect b m]
gs[Vect b m] -> Int -> Vect b m
forall a. [a] -> Int -> a
!Int
j in
        if m -> m -> Bool
forall m. Monomial m => m -> m -> Bool
mcoprime (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
fi) (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
fj) Bool -> Bool -> Bool
|| Vect b m -> Vect b m -> Bool
forall b b. Vect b m -> Vect b m -> Bool
criterion Vect b m
fi Vect b m
fj
        then [Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
gb' [Vect b m]
gs [(Int, Int)]
ps Int
t
        else let h :: Vect b m
h = Vect b m -> Vect b m -> Vect b m
forall b3 b.
(Eq b3, Ord b, Algebra b3 b, Fractional b3, Monomial b) =>
Vect b3 b -> Vect b3 b -> Vect b3 b
sPoly Vect b m
fi Vect b m
fj Vect b m -> [Vect b m] -> Vect b m
forall k m.
(Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
Vect k m -> [Vect k m] -> Vect k m
%% [Vect b m]
gs in
             if Vect b m
h Vect b m -> Vect b m -> Bool
forall a. Eq a => a -> a -> Bool
== 0 then [Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
gb' [Vect b m]
gs [(Int, Int)]
ps Int
t else [Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
gb' ([Vect b m]
gs[Vect b m] -> [Vect b m] -> [Vect b m]
forall a. [a] -> [a] -> [a]
++[Vect b m
h]) ([(Int, Int)]
ps [(Int, Int)] -> [(Int, Int)] -> [(Int, Int)]
forall a. [a] -> [a] -> [a]
++ [ (Int
i,Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+1) | Int
i <- [1..Int
t] ]) (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+1)
        where
            criterion :: Vect b m -> Vect b m -> Bool
criterion fi :: Vect b m
fi fj :: Vect b m
fj = let l :: m
l = m -> m -> m
forall m. Monomial m => m -> m -> m
mlcm (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
fi) (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
fj) in (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (m -> Int -> Bool
test m
l) [1..Int
t]
            test :: m -> Int -> Bool
test l :: m
l k :: Int
k = Int
k Int -> [Int] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Int
i,Int
j]
                    Bool -> Bool -> Bool
&& Int -> Int -> (Int, Int)
forall b. Ord b => b -> b -> (b, b)
ordpair Int
i Int
k (Int, Int) -> [(Int, Int)] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [(Int, Int)]
ps
                    Bool -> Bool -> Bool
&& Int -> Int -> (Int, Int)
forall b. Ord b => b -> b -> (b, b)
ordpair Int
j Int
k (Int, Int) -> [(Int, Int)] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [(Int, Int)]
ps
                    Bool -> Bool -> Bool
&& Vect b m -> m
forall b c. Vect b c -> c
lm ([Vect b m]
gs[Vect b m] -> Int -> Vect b m
forall a. [a] -> Int -> a
!Int
k) m -> m -> Bool
forall m. Monomial m => m -> m -> Bool
`mdivides` m
l
    gb' gs :: [Vect b m]
gs [] _ = [Vect b m]
gs

xs :: [a]
xs ! :: [a] -> Int -> a
! i :: Int
i = [a]
xs [a] -> Int -> a
forall a. [a] -> Int -> a
!! (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-1) -- in other words, index the list from 1 not 0


-- Doesn't result in any speedup
gb2a :: [Vect k b] -> [Vect k b]
gb2a fs :: [Vect k b]
fs = [Vect k b] -> [Vect k b]
forall k b.
(Fractional k, Monomial b, Ord k, Ord b, Algebra k b) =>
[Vect k b] -> [Vect k b]
reduce ([Vect k b] -> [Vect k b]) -> [Vect k b] -> [Vect k b]
forall a b. (a -> b) -> a -> b
$ IntMap (Vect k b) -> [(Int, Int)] -> Int -> [Vect k b]
forall m b.
(Fractional b, Monomial m, Ord m, Algebra b m, Eq b) =>
IntMap (Vect b m) -> [(Int, Int)] -> Int -> [Vect b m]
gb' IntMap (Vect k b)
fs' ([Int] -> [(Int, Int)]
forall t. [t] -> [(t, t)]
pairs [1..Int
s]) Int
s
    where fs' :: IntMap (Vect k b)
fs' = [(Int, Vect k b)] -> IntMap (Vect k b)
forall a. [(Int, a)] -> IntMap a
IM.fromList ([(Int, Vect k b)] -> IntMap (Vect k b))
-> [(Int, Vect k b)] -> IntMap (Vect k b)
forall a b. (a -> b) -> a -> b
$ [Int] -> [Vect k b] -> [(Int, Vect k b)]
forall a b. [a] -> [b] -> [(a, b)]
zip [1..] ([Vect k b] -> [(Int, Vect k b)])
-> [Vect k b] -> [(Int, Vect k b)]
forall a b. (a -> b) -> a -> b
$ (Vect k b -> Bool) -> [Vect k b] -> [Vect k b]
forall a. (a -> Bool) -> [a] -> [a]
filter (Vect k b -> Vect k b -> Bool
forall a. Eq a => a -> a -> Bool
/= 0) [Vect k b]
fs
          s :: Int
s = IntMap (Vect k b) -> Int
forall a. IntMap a -> Int
IM.size IntMap (Vect k b)
fs'
          gb' :: IntMap (Vect b m) -> [(Int, Int)] -> Int -> [Vect b m]
gb' gs :: IntMap (Vect b m)
gs ((i :: Int
i,j :: Int
j):ps :: [(Int, Int)]
ps) t :: Int
t =
              let fi :: Vect b m
fi = IntMap (Vect b m)
gs IntMap (Vect b m) -> Int -> Vect b m
forall a. IntMap a -> Int -> a
IM.! Int
i; fj :: Vect b m
fj = IntMap (Vect b m)
gs IntMap (Vect b m) -> Int -> Vect b m
forall a. IntMap a -> Int -> a
IM.! Int
j in
              if m -> m -> Bool
forall m. Monomial m => m -> m -> Bool
mcoprime (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
fi) (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
fj) Bool -> Bool -> Bool
|| Vect b m -> Vect b m -> Bool
forall b b. Vect b m -> Vect b m -> Bool
criterion Vect b m
fi Vect b m
fj
              then IntMap (Vect b m) -> [(Int, Int)] -> Int -> [Vect b m]
gb' IntMap (Vect b m)
gs [(Int, Int)]
ps Int
t
              else let h :: Vect b m
h = Vect b m -> Vect b m -> Vect b m
forall b3 b.
(Eq b3, Ord b, Algebra b3 b, Fractional b3, Monomial b) =>
Vect b3 b -> Vect b3 b -> Vect b3 b
sPoly Vect b m
fi Vect b m
fj Vect b m -> [Vect b m] -> Vect b m
forall k m.
(Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
Vect k m -> [Vect k m] -> Vect k m
%% IntMap (Vect b m) -> [Vect b m]
forall a. IntMap a -> [a]
IM.elems IntMap (Vect b m)
gs in
                   if Vect b m
h Vect b m -> Vect b m -> Bool
forall a. Eq a => a -> a -> Bool
== 0
                   then IntMap (Vect b m) -> [(Int, Int)] -> Int -> [Vect b m]
gb' IntMap (Vect b m)
gs [(Int, Int)]
ps Int
t
                   else let t' :: Int
t' = Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+1
                            gs' :: IntMap (Vect b m)
gs' = Int -> Vect b m -> IntMap (Vect b m) -> IntMap (Vect b m)
forall a. Int -> a -> IntMap a -> IntMap a
IM.insert Int
t' Vect b m
h IntMap (Vect b m)
gs
                            ps' :: [(Int, Int)]
ps' = [(Int, Int)]
ps [(Int, Int)] -> [(Int, Int)] -> [(Int, Int)]
forall a. [a] -> [a] -> [a]
++ (Int -> (Int, Int)) -> [Int] -> [(Int, Int)]
forall a b. (a -> b) -> [a] -> [b]
map (,Int
t') [1..Int
t]
                        in IntMap (Vect b m) -> [(Int, Int)] -> Int -> [Vect b m]
gb' IntMap (Vect b m)
gs' [(Int, Int)]
ps' Int
t'
              where criterion :: Vect b m -> Vect b m -> Bool
criterion fi :: Vect b m
fi fj :: Vect b m
fj = let l :: m
l = m -> m -> m
forall m. Monomial m => m -> m -> m
mlcm (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
fi) (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
fj) in (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (m -> Int -> Bool
test m
l) [1..Int
t]
                    test :: m -> Int -> Bool
test l :: m
l k :: Int
k = Int
k Int -> [Int] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Int
i,Int
j]
                            Bool -> Bool -> Bool
&& Int -> Int -> (Int, Int)
forall b. Ord b => b -> b -> (b, b)
ordpair Int
i Int
k (Int, Int) -> [(Int, Int)] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [(Int, Int)]
ps
                            Bool -> Bool -> Bool
&& Int -> Int -> (Int, Int)
forall b. Ord b => b -> b -> (b, b)
ordpair Int
j Int
k (Int, Int) -> [(Int, Int)] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [(Int, Int)]
ps
                            Bool -> Bool -> Bool
&& Vect b m -> m
forall b c. Vect b c -> c
lm (IntMap (Vect b m)
gs IntMap (Vect b m) -> Int -> Vect b m
forall a. IntMap a -> Int -> a
IM.! Int
k) m -> m -> Bool
forall m. Monomial m => m -> m -> Bool
`mdivides` m
l
          gb' gs :: IntMap (Vect b m)
gs [] _ = IntMap (Vect b m) -> [Vect b m]
forall a. IntMap a -> [a]
IM.elems IntMap (Vect b m)
gs


-- version of gb2 where we eliminate pairs as they're created, rather than as they're processed
gb3 :: [Vect k b] -> [Vect k b]
gb3 fs :: [Vect k b]
fs = [Vect k b] -> [Vect k b]
forall k b.
(Fractional k, Monomial b, Ord k, Ord b, Algebra k b) =>
[Vect k b] -> [Vect k b]
reduce ([Vect k b] -> [Vect k b]) -> [Vect k b] -> [Vect k b]
forall a b. (a -> b) -> a -> b
$ [Vect k b] -> [Vect k b] -> [(Int, Int)] -> Int -> [Vect k b]
forall b m.
(Fractional b, Ord m, Algebra b m, Eq b, Monomial m) =>
[Vect b m] -> [Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
gb1' [] [Vect k b]
fs [] 0
    where
    gb1' :: [Vect b m] -> [Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
gb1' gs :: [Vect b m]
gs (f :: Vect b m
f:fs :: [Vect b m]
fs) ps :: [(Int, Int)]
ps t :: Int
t = let ps' :: [(Int, Int)]
ps' = [Vect b m] -> [(Int, Int)] -> Vect b m -> Int -> [(Int, Int)]
forall m b b.
Monomial m =>
[Vect b m] -> [(Int, Int)] -> Vect b m -> Int -> [(Int, Int)]
updatePairs [Vect b m]
gs [(Int, Int)]
ps Vect b m
f Int
t
                          in [Vect b m] -> [Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
gb1' ([Vect b m]
gs [Vect b m] -> [Vect b m] -> [Vect b m]
forall a. [a] -> [a] -> [a]
++ [Vect b m
f]) [Vect b m]
fs [(Int, Int)]
ps' (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+1)
    gb1' ls :: [Vect b m]
ls [] ps :: [(Int, Int)]
ps t :: Int
t = [Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
forall m b.
(Fractional b, Ord m, Algebra b m, Eq b, Monomial m) =>
[Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
gb2' [Vect b m]
ls [(Int, Int)]
ps Int
t
    gb2' :: [Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
gb2' gs :: [Vect b m]
gs ((i :: Int
i,j :: Int
j):ps :: [(Int, Int)]
ps) t :: Int
t =
        let h :: Vect b m
h = Vect b m -> Vect b m -> Vect b m
forall b3 b.
(Eq b3, Ord b, Algebra b3 b, Fractional b3, Monomial b) =>
Vect b3 b -> Vect b3 b -> Vect b3 b
sPoly ([Vect b m]
gs[Vect b m] -> Int -> Vect b m
forall a. [a] -> Int -> a
!Int
i) ([Vect b m]
gs[Vect b m] -> Int -> Vect b m
forall a. [a] -> Int -> a
!Int
j) Vect b m -> [Vect b m] -> Vect b m
forall k m.
(Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
Vect k m -> [Vect k m] -> Vect k m
%% [Vect b m]
gs in
        if Vect b m
h Vect b m -> Vect b m -> Bool
forall a. Eq a => a -> a -> Bool
== 0
        then [Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
gb2' [Vect b m]
gs [(Int, Int)]
ps Int
t
        else let ps' :: [(Int, Int)]
ps' = [Vect b m] -> [(Int, Int)] -> Vect b m -> Int -> [(Int, Int)]
forall m b b.
Monomial m =>
[Vect b m] -> [(Int, Int)] -> Vect b m -> Int -> [(Int, Int)]
updatePairs [Vect b m]
gs ((Int
i,Int
j)(Int, Int) -> [(Int, Int)] -> [(Int, Int)]
forall a. a -> [a] -> [a]
:[(Int, Int)]
ps) Vect b m
h Int
t in [Vect b m] -> [(Int, Int)] -> Int -> [Vect b m]
gb2' ([Vect b m]
gs[Vect b m] -> [Vect b m] -> [Vect b m]
forall a. [a] -> [a] -> [a]
++[Vect b m
h]) [(Int, Int)]
ps' (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+1)
    gb2' gs :: [Vect b m]
gs [] _ = [Vect b m]
gs
    updatePairs :: [Vect b m] -> [(Int, Int)] -> Vect b m -> Int -> [(Int, Int)]
updatePairs gs :: [Vect b m]
gs ps :: [(Int, Int)]
ps f :: Vect b m
f t :: Int
t =
        [(Int, Int)
p | p :: (Int, Int)
p@(i :: Int
i,j :: Int
j) <- [(Int, Int)]
ps,
             Bool -> Bool
not (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
f m -> m -> Bool
forall m. Monomial m => m -> m -> Bool
`mdivides` m -> m -> m
forall m. Monomial m => m -> m -> m
mlcm (Vect b m -> m
forall b c. Vect b c -> c
lm ([Vect b m]
gs[Vect b m] -> Int -> Vect b m
forall a. [a] -> Int -> a
!Int
i)) (Vect b m -> m
forall b c. Vect b c -> c
lm ([Vect b m]
gs[Vect b m] -> Int -> Vect b m
forall a. [a] -> Int -> a
!Int
j)))] 
     [(Int, Int)] -> [(Int, Int)] -> [(Int, Int)]
forall a. [a] -> [a] -> [a]
++ [ (Int
i,Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+1) | (gi :: Vect b m
gi,i :: Int
i) <- [Vect b m] -> [Int] -> [(Vect b m, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Vect b m]
gs [1..Int
t],
                    Bool -> Bool
not (m -> m -> Bool
forall m. Monomial m => m -> m -> Bool
mcoprime (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
gi) (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
f)),
                    Bool -> Bool
not (m -> Int -> Bool
criterion (m -> m -> m
forall m. Monomial m => m -> m -> m
mlcm (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
gi) (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
f)) Int
i) ]
        where criterion :: m -> Int -> Bool
criterion l :: m
l i :: Int
i = (m -> Bool) -> [m] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (m -> m -> Bool
forall m. Monomial m => m -> m -> Bool
`mdivides` m
l) [Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
gk | (gk :: Vect b m
gk,k :: Int
k) <- [Vect b m] -> [Int] -> [(Vect b m, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Vect b m]
gs [1..Int
t], Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i, Int -> Int -> (Int, Int)
forall b. Ord b => b -> b -> (b, b)
ordpair Int
i Int
k (Int, Int) -> [(Int, Int)] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [(Int, Int)]
ps]


-- Cox et al 108
-- 1. list smallest fs first, as more likely to reduce
-- 2. order the pairs with smallest lcm fi fj first ("normal selection strategy")
gb4 :: [Vect k b] -> [Vect k b]
gb4 fs :: [Vect k b]
fs = [Vect k b] -> [Vect k b]
forall k b.
(Fractional k, Monomial b, Ord k, Ord b, Algebra k b) =>
[Vect k b] -> [Vect k b]
reduce ([Vect k b] -> [Vect k b]) -> [Vect k b] -> [Vect k b]
forall a b. (a -> b) -> a -> b
$ [Vect k b] -> [Vect k b] -> [(b, (Int, Int))] -> Int -> [Vect k b]
forall b a.
(Fractional b, Algebra b a, Eq b, Ord a, Monomial a) =>
[Vect b a] -> [Vect b a] -> [(a, (Int, Int))] -> Int -> [Vect b a]
gb1' [] [Vect k b]
fs' [] 0
    where fs' :: [Vect k b]
fs' = [Vect k b] -> [Vect k b]
forall a. [a] -> [a]
reverse ([Vect k b] -> [Vect k b]) -> [Vect k b] -> [Vect k b]
forall a b. (a -> b) -> a -> b
$ [Vect k b] -> [Vect k b]
forall a. Ord a => [a] -> [a]
L.sort ([Vect k b] -> [Vect k b]) -> [Vect k b] -> [Vect k b]
forall a b. (a -> b) -> a -> b
$ (Vect k b -> Bool) -> [Vect k b] -> [Vect k b]
forall a. (a -> Bool) -> [a] -> [a]
filter (Vect k b -> Vect k b -> Bool
forall a. Eq a => a -> a -> Bool
/=0) [Vect k b]
fs
          gb1' :: [Vect b a] -> [Vect b a] -> [(a, (Int, Int))] -> Int -> [Vect b a]
gb1' gs :: [Vect b a]
gs (f :: Vect b a
f:fs :: [Vect b a]
fs) ps :: [(a, (Int, Int))]
ps t :: Int
t = [Vect b a] -> [Vect b a] -> [(a, (Int, Int))] -> Int -> [Vect b a]
gb1' ([Vect b a]
gs [Vect b a] -> [Vect b a] -> [Vect b a]
forall a. [a] -> [a] -> [a]
++ [Vect b a
f]) [Vect b a]
fs [(a, (Int, Int))]
ps' (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+1)
              where ps' :: [(a, (Int, Int))]
ps' = [Vect b a]
-> [(a, (Int, Int))] -> Vect b a -> Int -> [(a, (Int, Int))]
forall a b b b.
(Ord b, Ord a, Enum b, Num b, Monomial a) =>
[Vect b a] -> [(a, (b, b))] -> Vect b a -> b -> [(a, (b, b))]
updatePairs [Vect b a]
gs [(a, (Int, Int))]
ps Vect b a
f Int
t
          gb1' ls :: [Vect b a]
ls [] ps :: [(a, (Int, Int))]
ps t :: Int
t = [Vect b a] -> [(a, (Int, Int))] -> Int -> [Vect b a]
forall a b.
(Fractional b, Algebra b a, Eq b, Ord a, Monomial a) =>
[Vect b a] -> [(a, (Int, Int))] -> Int -> [Vect b a]
gb2' [Vect b a]
ls [(a, (Int, Int))]
ps Int
t
          gb2' :: [Vect b a] -> [(a, (Int, Int))] -> Int -> [Vect b a]
gb2' gs :: [Vect b a]
gs ((l :: a
l,(i :: Int
i,j :: Int
j)):ps :: [(a, (Int, Int))]
ps) t :: Int
t =
              let h :: Vect b a
h = Vect b a -> Vect b a -> Vect b a
forall b3 b.
(Eq b3, Ord b, Algebra b3 b, Fractional b3, Monomial b) =>
Vect b3 b -> Vect b3 b -> Vect b3 b
sPoly ([Vect b a]
gs[Vect b a] -> Int -> Vect b a
forall a. [a] -> Int -> a
!Int
i) ([Vect b a]
gs[Vect b a] -> Int -> Vect b a
forall a. [a] -> Int -> a
!Int
j) Vect b a -> [Vect b a] -> Vect b a
forall k m.
(Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
Vect k m -> [Vect k m] -> Vect k m
%% [Vect b a]
gs in
              if Vect b a
h Vect b a -> Vect b a -> Bool
forall a. Eq a => a -> a -> Bool
== 0
              then [Vect b a] -> [(a, (Int, Int))] -> Int -> [Vect b a]
gb2' [Vect b a]
gs [(a, (Int, Int))]
ps Int
t
              else let ps' :: [(a, (Int, Int))]
ps' = [Vect b a]
-> [(a, (Int, Int))] -> Vect b a -> Int -> [(a, (Int, Int))]
forall a b b b.
(Ord b, Ord a, Enum b, Num b, Monomial a) =>
[Vect b a] -> [(a, (b, b))] -> Vect b a -> b -> [(a, (b, b))]
updatePairs [Vect b a]
gs ((a
l,(Int
i,Int
j))(a, (Int, Int)) -> [(a, (Int, Int))] -> [(a, (Int, Int))]
forall a. a -> [a] -> [a]
:[(a, (Int, Int))]
ps) Vect b a
h Int
t in [Vect b a] -> [(a, (Int, Int))] -> Int -> [Vect b a]
gb2' ([Vect b a]
gs[Vect b a] -> [Vect b a] -> [Vect b a]
forall a. [a] -> [a] -> [a]
++[Vect b a
h]) [(a, (Int, Int))]
ps' (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+1)
          gb2' gs :: [Vect b a]
gs [] _ = [Vect b a]
gs
          updatePairs :: [Vect b a] -> [(a, (b, b))] -> Vect b a -> b -> [(a, (b, b))]
updatePairs gs :: [Vect b a]
gs ps :: [(a, (b, b))]
ps f :: Vect b a
f t :: b
t =
              let oldps :: [(a, (b, b))]
oldps = [(a, (b, b))
p | p :: (a, (b, b))
p@(l :: a
l,(i :: b
i,j :: b
j)) <- [(a, (b, b))]
ps, Bool -> Bool
not (Vect b a -> a
forall b c. Vect b c -> c
lm Vect b a
f a -> a -> Bool
forall m. Monomial m => m -> m -> Bool
`mdivides` a
l)]
                  newps :: [(a, (b, b))]
newps = ((a, (b, b)) -> (a, (b, b)) -> Ordering)
-> [(a, (b, b))] -> [(a, (b, b))]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (((a, (b, b)) -> (a, (b, b)) -> Ordering)
-> (a, (b, b)) -> (a, (b, b)) -> Ordering
forall a b c. (a -> b -> c) -> b -> a -> c
flip (a, (b, b)) -> (a, (b, b)) -> Ordering
forall a b1 b2. Ord a => (a, b1) -> (a, b2) -> Ordering
cmpfst)
                          [ (a
l,(b
i,b
tb -> b -> b
forall a. Num a => a -> a -> a
+1)) | (gi :: Vect b a
gi,i :: b
i) <- [Vect b a] -> [b] -> [(Vect b a, b)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Vect b a]
gs [1..b
t], let l :: a
l = a -> a -> a
forall m. Monomial m => m -> m -> m
mlcm (Vect b a -> a
forall b c. Vect b c -> c
lm Vect b a
gi) (Vect b a -> a
forall b c. Vect b c -> c
lm Vect b a
f),
                                          Bool -> Bool
not (a -> a -> Bool
forall m. Monomial m => m -> m -> Bool
mcoprime (Vect b a -> a
forall b c. Vect b c -> c
lm Vect b a
gi) (Vect b a -> a
forall b c. Vect b c -> c
lm Vect b a
f)),
                                          Bool -> Bool
not (a -> b -> Bool
criterion a
l b
i) ]
              in ((a, (b, b)) -> (a, (b, b)) -> Ordering)
-> [(a, (b, b))] -> [(a, (b, b))] -> [(a, (b, b))]
forall a. (a -> a -> Ordering) -> [a] -> [a] -> [a]
mergeBy (((a, (b, b)) -> (a, (b, b)) -> Ordering)
-> (a, (b, b)) -> (a, (b, b)) -> Ordering
forall a b c. (a -> b -> c) -> b -> a -> c
flip (a, (b, b)) -> (a, (b, b)) -> Ordering
forall a b1 b2. Ord a => (a, b1) -> (a, b2) -> Ordering
cmpfst) [(a, (b, b))]
oldps [(a, (b, b))]
newps
              where criterion :: a -> b -> Bool
criterion l :: a
l i :: b
i = (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (a -> a -> Bool
forall m. Monomial m => m -> m -> Bool
`mdivides` a
l) [Vect b a -> a
forall b c. Vect b c -> c
lm Vect b a
gk | (gk :: Vect b a
gk,k :: b
k) <- [Vect b a] -> [b] -> [(Vect b a, b)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Vect b a]
gs [1..b
t], b
k b -> b -> Bool
forall a. Eq a => a -> a -> Bool
/= b
i, b -> b -> (b, b)
forall b. Ord b => b -> b -> (b, b)
ordpair b
i b
k (b, b) -> [(b, b)] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` ((a, (b, b)) -> (b, b)) -> [(a, (b, b))] -> [(b, b)]
forall a b. (a -> b) -> [a] -> [b]
map (a, (b, b)) -> (b, b)
forall a b. (a, b) -> b
snd [(a, (b, b))]
ps]


mergeBy :: (a -> a -> Ordering) -> [a] -> [a] -> [a]
mergeBy cmp :: a -> a -> Ordering
cmp (t :: a
t:ts :: [a]
ts) (u :: a
u:us :: [a]
us) =
    case a -> a -> Ordering
cmp a
t a
u of
    LT -> a
t a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a -> Ordering) -> [a] -> [a] -> [a]
mergeBy a -> a -> Ordering
cmp [a]
ts (a
ua -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
us)
    EQ -> a
t a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a -> Ordering) -> [a] -> [a] -> [a]
mergeBy a -> a -> Ordering
cmp [a]
ts (a
ua -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
us)
    GT -> a
u a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a -> Ordering) -> [a] -> [a] -> [a]
mergeBy a -> a -> Ordering
cmp (a
ta -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
ts) [a]
us
mergeBy _ ts :: [a]
ts us :: [a]
us = [a]
ts [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a]
us -- one of them is null


-- Giovini et al
-- The point of sugar is, given fi, fj, to give an upper bound on the degree of sPoly fi fj without having to calculate it
-- We can then select by preference pairs with lower sugar, expecting therefore that the s-polys will have lower degree

-- It is only for Lex ordering that sugar seems to give a substantial improvement
-- gb4 is fine for Grevlex

-- !! can probably get rid of the Ord k requirement in the following

-- |Given a list of polynomials over a field, return a Groebner basis for the ideal generated by the polynomials.
gb :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
    [Vect k m] -> [Vect k m]
gb :: [Vect k m] -> [Vect k m]
gb fs :: [Vect k m]
fs =
    let fs' :: [Vect k m]
fs' = [Vect k m] -> [Vect k m]
forall a. [a] -> [a]
reverse ([Vect k m] -> [Vect k m]) -> [Vect k m] -> [Vect k m]
forall a b. (a -> b) -> a -> b
$ [Vect k m] -> [Vect k m]
forall a. Ord a => [a] -> [a]
sort ([Vect k m] -> [Vect k m]) -> [Vect k m] -> [Vect k m]
forall a b. (a -> b) -> a -> b
$ (Vect k m -> Vect k m) -> [Vect k m] -> [Vect k m]
forall a b. (a -> b) -> [a] -> [b]
map Vect k m -> Vect k m
forall b k.
(Eq k, Ord b, Show b, Algebra k b, Fractional k) =>
Vect k b -> Vect k b
toMonic ([Vect k m] -> [Vect k m]) -> [Vect k m] -> [Vect k m]
forall a b. (a -> b) -> a -> b
$ (Vect k m -> Bool) -> [Vect k m] -> [Vect k m]
forall a. (a -> Bool) -> [a] -> [a]
filter (Vect k m -> Vect k m -> Bool
forall a. Eq a => a -> a -> Bool
/=0) [Vect k m]
fs
    in [Vect k m] -> [Vect k m]
forall k b.
(Fractional k, Monomial b, Ord k, Ord b, Algebra k b) =>
[Vect k b] -> [Vect k b]
reduce ([Vect k m] -> [Vect k m]) -> [Vect k m] -> [Vect k m]
forall a b. (a -> b) -> a -> b
$ [Vect k m]
-> [Vect k m] -> [((Int, m), (Int, Int))] -> Int -> [Vect k m]
forall b b.
(Algebra b b, Fractional b, Eq b, Ord b, Monomial b) =>
[Vect b b]
-> [Vect b b] -> [((Int, b), (Int, Int))] -> Int -> [Vect b b]
gb1' [] [Vect k m]
fs' [] 0 where
    gb1' :: [Vect b b]
-> [Vect b b] -> [((Int, b), (Int, Int))] -> Int -> [Vect b b]
gb1' gs :: [Vect b b]
gs (f :: Vect b b
f:fs :: [Vect b b]
fs) ps :: [((Int, b), (Int, Int))]
ps t :: Int
t = [Vect b b]
-> [Vect b b] -> [((Int, b), (Int, Int))] -> Int -> [Vect b b]
gb1' ([Vect b b]
gs [Vect b b] -> [Vect b b] -> [Vect b b]
forall a. [a] -> [a] -> [a]
++ [Vect b b
f]) [Vect b b]
fs [((Int, b), (Int, Int))]
ps' (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+1)
        where ps' :: [((Int, b), (Int, Int))]
ps' = [Vect b b]
-> [((Int, b), (Int, Int))]
-> Vect b b
-> Int
-> [((Int, b), (Int, Int))]
forall b b b.
(Ord b, Monomial b) =>
[Vect b b]
-> [((Int, b), (Int, Int))]
-> Vect b b
-> Int
-> [((Int, b), (Int, Int))]
updatePairs [Vect b b]
gs [((Int, b), (Int, Int))]
ps Vect b b
f (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+1)
    gb1' ls :: [Vect b b]
ls [] ps :: [((Int, b), (Int, Int))]
ps t :: Int
t = [Vect b b] -> [((Int, b), (Int, Int))] -> Int -> [Vect b b]
forall b b.
(Algebra b b, Fractional b, Eq b, Ord b, Monomial b) =>
[Vect b b] -> [((Int, b), (Int, Int))] -> Int -> [Vect b b]
gb2' [Vect b b]
ls [((Int, b), (Int, Int))]
ps Int
t
    gb2' :: [Vect b b] -> [((Int, b), (Int, Int))] -> Int -> [Vect b b]
gb2' gs :: [Vect b b]
gs (p :: ((Int, b), (Int, Int))
p@(_,(i :: Int
i,j :: Int
j)):ps :: [((Int, b), (Int, Int))]
ps) t :: Int
t =
        if Vect b b
h Vect b b -> Vect b b -> Bool
forall a. Eq a => a -> a -> Bool
== 0
        then [Vect b b] -> [((Int, b), (Int, Int))] -> Int -> [Vect b b]
gb2' [Vect b b]
gs [((Int, b), (Int, Int))]
ps Int
t
        else [Vect b b] -> [((Int, b), (Int, Int))] -> Int -> [Vect b b]
gb2' ([Vect b b]
gs[Vect b b] -> [Vect b b] -> [Vect b b]
forall a. [a] -> [a] -> [a]
++[Vect b b
h]) [((Int, b), (Int, Int))]
ps' (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+1)
        where h :: Vect b b
h = Vect b b -> Vect b b
forall b k.
(Eq k, Ord b, Show b, Algebra k b, Fractional k) =>
Vect k b -> Vect k b
toMonic (Vect b b -> Vect b b) -> Vect b b -> Vect b b
forall a b. (a -> b) -> a -> b
$ Vect b b -> Vect b b -> Vect b b
forall b3 b.
(Eq b3, Ord b, Algebra b3 b, Fractional b3, Monomial b) =>
Vect b3 b -> Vect b3 b -> Vect b3 b
sPoly ([Vect b b]
gs[Vect b b] -> Int -> Vect b b
forall a. [a] -> Int -> a
!Int
i) ([Vect b b]
gs[Vect b b] -> Int -> Vect b b
forall a. [a] -> Int -> a
!Int
j) Vect b b -> [Vect b b] -> Vect b b
forall k m.
(Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
Vect k m -> [Vect k m] -> Vect k m
%% [Vect b b]
gs
              ps' :: [((Int, b), (Int, Int))]
ps' = [Vect b b]
-> [((Int, b), (Int, Int))]
-> Vect b b
-> Int
-> [((Int, b), (Int, Int))]
forall b b b.
(Ord b, Monomial b) =>
[Vect b b]
-> [((Int, b), (Int, Int))]
-> Vect b b
-> Int
-> [((Int, b), (Int, Int))]
updatePairs [Vect b b]
gs (((Int, b), (Int, Int))
p((Int, b), (Int, Int))
-> [((Int, b), (Int, Int))] -> [((Int, b), (Int, Int))]
forall a. a -> [a] -> [a]
:[((Int, b), (Int, Int))]
ps) Vect b b
h (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+1)
    gb2' gs :: [Vect b b]
gs [] _ = [Vect b b]
gs
    updatePairs :: [Vect b b]
-> [((Int, b), (Int, Int))]
-> Vect b b
-> Int
-> [((Int, b), (Int, Int))]
updatePairs gs :: [Vect b b]
gs ps :: [((Int, b), (Int, Int))]
ps gk :: Vect b b
gk k :: Int
k =
        let newps :: [((Int, b), (Int, Int))]
newps = [let l :: b
l = b -> b -> b
forall m. Monomial m => m -> m -> m
mlcm (Vect b b -> b
forall b c. Vect b c -> c
lm Vect b b
gi) (Vect b b -> b
forall b c. Vect b c -> c
lm Vect b b
gk) in ((Vect b b -> Vect b b -> b -> Int
forall m m m b b.
(Monomial m, Monomial m, Monomial m) =>
Vect b m -> Vect b m -> m -> Int
sugar Vect b b
gi Vect b b
gk b
l, b
l), (Int
i,Int
k)) | (gi :: Vect b b
gi,i :: Int
i) <- [Vect b b] -> [Int] -> [(Vect b b, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Vect b b]
gs [1..Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-1] ]
            ps' :: [((Int, b), (Int, Int))]
ps' = [((Int, b), (Int, Int))
p | p :: ((Int, b), (Int, Int))
p@((sij :: Int
sij,tij :: b
tij),(i :: Int
i,j :: Int
j)) <- [((Int, b), (Int, Int))]
ps,
                       let ((sik :: Int
sik,tik :: b
tik),_) = [((Int, b), (Int, Int))]
newps [((Int, b), (Int, Int))] -> Int -> ((Int, b), (Int, Int))
forall a. [a] -> Int -> a
! Int
i, let ((sjk :: Int
sjk,tjk :: b
tjk),_) = [((Int, b), (Int, Int))]
newps [((Int, b), (Int, Int))] -> Int -> ((Int, b), (Int, Int))
forall a. [a] -> Int -> a
! Int
j,
                       Bool -> Bool
not ( (b
tik b -> b -> Bool
forall m. Monomial m => m -> m -> Bool
`mproperlydivides` b
tij) Bool -> Bool -> Bool
&& (b
tjk b -> b -> Bool
forall m. Monomial m => m -> m -> Bool
`mproperlydivides` b
tij) ) ] -- sloppy variant
            newps' :: [((Int, b), (Int, Int))]
newps' = [((Int, b), (Int, Int))]
-> [((Int, b), (Int, Int))] -> [((Int, b), (Int, Int))]
forall a a b.
Eq a =>
[((a, a), (Int, b))]
-> [((a, a), (Int, b))] -> [((a, a), (Int, b))]
discard1 [] [((Int, b), (Int, Int))]
newps
            newps'' :: [((Int, b), (Int, Int))]
newps'' = (((Int, b), (Int, Int)) -> ((Int, b), (Int, Int)) -> Ordering)
-> [((Int, b), (Int, Int))] -> [((Int, b), (Int, Int))]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy ((((Int, b), (Int, Int)) -> ((Int, b), (Int, Int)) -> Ordering)
-> ((Int, b), (Int, Int)) -> ((Int, b), (Int, Int)) -> Ordering
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Int, b), (Int, Int)) -> ((Int, b), (Int, Int)) -> Ordering
forall a b c a a.
(Ord a, Ord b, Ord c, Num a) =>
((a, b), (a, c)) -> ((a, b), (a, c)) -> Ordering
cmpSug) ([((Int, b), (Int, Int))] -> [((Int, b), (Int, Int))])
-> [((Int, b), (Int, Int))] -> [((Int, b), (Int, Int))]
forall a b. (a -> b) -> a -> b
$ [((Int, b), (Int, Int))]
-> [((Int, b), (Int, Int))] -> [((Int, b), (Int, Int))]
forall m a a a.
(Monomial m, Eq a) =>
[((a, m), (a, a))] -> [((a, m), (a, a))] -> [((a, m), (a, a))]
discard2 [] ([((Int, b), (Int, Int))] -> [((Int, b), (Int, Int))])
-> [((Int, b), (Int, Int))] -> [((Int, b), (Int, Int))]
forall a b. (a -> b) -> a -> b
$ (((Int, b), (Int, Int)) -> ((Int, b), (Int, Int)) -> Ordering)
-> [((Int, b), (Int, Int))] -> [((Int, b), (Int, Int))]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy ((((Int, b), (Int, Int)) -> ((Int, b), (Int, Int)) -> Ordering)
-> ((Int, b), (Int, Int)) -> ((Int, b), (Int, Int)) -> Ordering
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Int, b), (Int, Int)) -> ((Int, b), (Int, Int)) -> Ordering
forall a b a a a a.
(Ord a, Ord b) =>
((a, a), (a, b)) -> ((a, a), (a, b)) -> Ordering
cmpNormal) [((Int, b), (Int, Int))]
newps'
        in (((Int, b), (Int, Int)) -> ((Int, b), (Int, Int)) -> Ordering)
-> [((Int, b), (Int, Int))]
-> [((Int, b), (Int, Int))]
-> [((Int, b), (Int, Int))]
forall a. (a -> a -> Ordering) -> [a] -> [a] -> [a]
mergeBy ((((Int, b), (Int, Int)) -> ((Int, b), (Int, Int)) -> Ordering)
-> ((Int, b), (Int, Int)) -> ((Int, b), (Int, Int)) -> Ordering
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Int, b), (Int, Int)) -> ((Int, b), (Int, Int)) -> Ordering
forall a b c a a.
(Ord a, Ord b, Ord c, Num a) =>
((a, b), (a, c)) -> ((a, b), (a, c)) -> Ordering
cmpSug) [((Int, b), (Int, Int))]
ps' [((Int, b), (Int, Int))]
newps''
        where
            discard1 :: [((a, a), (Int, b))]
-> [((a, a), (Int, b))] -> [((a, a), (Int, b))]
discard1 ls :: [((a, a), (Int, b))]
ls (r :: ((a, a), (Int, b))
r@((_sik :: a
_sik,tik :: a
tik),(i :: Int
i,_k :: b
_k)):rs :: [((a, a), (Int, b))]
rs) =
                if Vect b b -> b
forall b c. Vect b c -> c
lm ([Vect b b]
gs[Vect b b] -> Int -> Vect b b
forall a. [a] -> Int -> a
!Int
i) b -> b -> Bool
forall m. Monomial m => m -> m -> Bool
`mcoprime` Vect b b -> b
forall b c. Vect b c -> c
lm Vect b b
gk
                -- then discard [l | l@((_,tjk),_) <- ls, tjk /= tik] [r | r@((_,tjk),_) <- ls, tjk /= tik]
                then [((a, a), (Int, b))]
-> [((a, a), (Int, b))] -> [((a, a), (Int, b))]
discard1 ((((a, a), (Int, b)) -> Bool)
-> [((a, a), (Int, b))] -> [((a, a), (Int, b))]
forall a. (a -> Bool) -> [a] -> [a]
filter (\((_,tjk :: a
tjk),_) -> a
tjk a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
tik) [((a, a), (Int, b))]
ls) ((((a, a), (Int, b)) -> Bool)
-> [((a, a), (Int, b))] -> [((a, a), (Int, b))]
forall a. (a -> Bool) -> [a] -> [a]
filter (\((_,tjk :: a
tjk),_) -> a
tjk a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
tik) [((a, a), (Int, b))]
rs)
                else [((a, a), (Int, b))]
-> [((a, a), (Int, b))] -> [((a, a), (Int, b))]
discard1 (((a, a), (Int, b))
r((a, a), (Int, b)) -> [((a, a), (Int, b))] -> [((a, a), (Int, b))]
forall a. a -> [a] -> [a]
:[((a, a), (Int, b))]
ls) [((a, a), (Int, b))]
rs
            discard1 ls :: [((a, a), (Int, b))]
ls [] = [((a, a), (Int, b))]
ls
            discard2 :: [((a, m), (a, a))] -> [((a, m), (a, a))] -> [((a, m), (a, a))]
discard2 ls :: [((a, m), (a, a))]
ls (r :: ((a, m), (a, a))
r@((_sik :: a
_sik,tik :: m
tik),(i :: a
i,k :: a
k)):rs :: [((a, m), (a, a))]
rs) = [((a, m), (a, a))] -> [((a, m), (a, a))] -> [((a, m), (a, a))]
discard2 (((a, m), (a, a))
r((a, m), (a, a)) -> [((a, m), (a, a))] -> [((a, m), (a, a))]
forall a. a -> [a] -> [a]
:[((a, m), (a, a))]
ls) ([((a, m), (a, a))] -> [((a, m), (a, a))])
-> [((a, m), (a, a))] -> [((a, m), (a, a))]
forall a b. (a -> b) -> a -> b
$ (((a, m), (a, a)) -> Bool)
-> [((a, m), (a, a))] -> [((a, m), (a, a))]
forall a. (a -> Bool) -> [a] -> [a]
filter (\((_sjk :: a
_sjk,tjk :: m
tjk),(j :: a
j,k' :: a
k')) -> Bool -> Bool
not (a
k a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
k' Bool -> Bool -> Bool
&& m
tik m -> m -> Bool
forall m. Monomial m => m -> m -> Bool
`mdivides` m
tjk)) [((a, m), (a, a))]
rs
            discard2 ls :: [((a, m), (a, a))]
ls [] = [((a, m), (a, a))]
ls
-- The two calls to toMonic are designed to prevent coefficient explosion, but it is unproven that they are effective

-- sugar of sPoly f g, where h = lcm (lt f) (lt g)
-- this is an upper bound on deg (sPoly f g)
sugar :: Vect b m -> Vect b m -> m -> Int
sugar f :: Vect b m
f g :: Vect b m
g h :: m
h = m -> Int
forall m. Monomial m => m -> Int
mdeg m
h Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (Vect b m -> Int
forall m k. Monomial m => Vect k m -> Int
deg Vect b m
f Int -> Int -> Int
forall a. Num a => a -> a -> a
- m -> Int
forall m. Monomial m => m -> Int
mdeg (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
f)) (Vect b m -> Int
forall m k. Monomial m => Vect k m -> Int
deg Vect b m
g Int -> Int -> Int
forall a. Num a => a -> a -> a
- m -> Int
forall m. Monomial m => m -> Int
mdeg (Vect b m -> m
forall b c. Vect b c -> c
lm Vect b m
g))

cmpNormal :: ((a, a), (a, b)) -> ((a, a), (a, b)) -> Ordering
cmpNormal ((s1 :: a
s1,t1 :: a
t1),(i1 :: a
i1,j1 :: b
j1)) ((s2 :: a
s2,t2 :: a
t2),(i2 :: a
i2,j2 :: b
j2)) = (a, b) -> (a, b) -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (a
t1,b
j1) (a
t2,b
j2)

cmpSug :: ((a, b), (a, c)) -> ((a, b), (a, c)) -> Ordering
cmpSug ((s1 :: a
s1,t1 :: b
t1),(i1 :: a
i1,j1 :: c
j1)) ((s2 :: a
s2,t2 :: b
t2),(i2 :: a
i2,j2 :: c
j2)) = (a, b, c) -> (a, b, c) -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (-a
s1,b
t1,c
j1) (-a
s2,b
t2,c
j2)


{-
merge (t:ts) (u:us) =
    if t <= u
    then t : merge ts (u:us)
    else u : merge (t:ts) us
merge ts us = ts ++ us -- one of them is null
-}

-- OPERATIONS ON IDEALS

memberGB :: Vect k m -> [Vect k m] -> Bool
memberGB f :: Vect k m
f gs :: [Vect k m]
gs = Vect k m
f Vect k m -> [Vect k m] -> Vect k m
forall k m.
(Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
Vect k m -> [Vect k m] -> Vect k m
%% [Vect k m]
gs Vect k m -> Vect k m -> Bool
forall a. Eq a => a -> a -> Bool
== 0

-- |@memberI f gs@ returns whether f is in the ideal generated by gs
memberI :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
    Vect k m -> [Vect k m] -> Bool
memberI :: Vect k m -> [Vect k m] -> Bool
memberI f :: Vect k m
f gs :: [Vect k m]
gs = Vect k m -> [Vect k m] -> Bool
forall m k.
(Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
Vect k m -> [Vect k m] -> Bool
memberGB Vect k m
f ([Vect k m] -> [Vect k m]
forall k m.
(Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
[Vect k m] -> [Vect k m]
gb [Vect k m]
gs)

-- Cox et al, p181
-- |Given ideals I and J, their sum is defined as I+J = {f+g | f \<- I, g \<- J}.
--
-- If fs and gs are generators for I and J, then @sumI fs gs@ returns generators for I+J.
--
-- The geometric interpretation is that the variety of the sum is the intersection of the varieties,
-- ie V(I+J) = V(I) intersect V(J)
sumI :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
     [Vect k m] -> [Vect k m] -> [Vect k m]
sumI :: [Vect k m] -> [Vect k m] -> [Vect k m]
sumI fs :: [Vect k m]
fs gs :: [Vect k m]
gs = [Vect k m] -> [Vect k m]
forall k m.
(Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
[Vect k m] -> [Vect k m]
gb ([Vect k m]
fs [Vect k m] -> [Vect k m] -> [Vect k m]
forall a. [a] -> [a] -> [a]
++ [Vect k m]
gs)

-- Cox et al, p183
-- |Given ideals I and J, their product I.J is the ideal generated by all products {f.g | f \<- I, g \<- J}.
--
-- If fs and gs are generators for I and J, then @productI fs gs@ returns generators for I.J.
--
-- The geometric interpretation is that the variety of the product is the union of the varieties,
-- ie V(I.J) = V(I) union V(J)
productI :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
     [Vect k m] -> [Vect k m] -> [Vect k m]
productI :: [Vect k m] -> [Vect k m] -> [Vect k m]
productI fs :: [Vect k m]
fs gs :: [Vect k m]
gs = [Vect k m] -> [Vect k m]
forall k m.
(Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
[Vect k m] -> [Vect k m]
gb [Vect k m
f Vect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
* Vect k m
g | Vect k m
f <- [Vect k m]
fs, Vect k m
g <- [Vect k m]
gs]

-- Cox et al, p185-6
-- |The intersection of ideals I and J is the set of all polynomials which belong to both I and J.
--
-- If fs and gs are generators for I and J, then @intersectI fs gs@ returns generators for the intersection of I and J
--
-- The geometric interpretation is that the variety of the intersection is the union of the varieties,
-- ie V(I intersect J) = V(I) union V(J).
--
-- The reason for prefering the intersection over the product is that the intersection of radical ideals is radical,
-- whereas the product need not be.
intersectI :: (Fractional k, Ord k, Monomial m, Ord m) =>
     [Vect k m] -> [Vect k m] -> [Vect k m]
intersectI :: [Vect k m] -> [Vect k m] -> [Vect k m]
intersectI fs :: [Vect k m]
fs gs :: [Vect k m]
gs =
    let t :: f (Elim2 (Glex String) b)
t = f (Glex String) -> f (Elim2 (Glex String) b)
forall (f :: * -> *) b a.
(Functor f, Mon b) =>
f a -> f (Elim2 a b)
toElimFst (f (Glex String) -> f (Elim2 (Glex String) b))
-> f (Glex String) -> f (Elim2 (Glex String) b)
forall a b. (a -> b) -> a -> b
$ Glex String -> f (Glex String)
forall (m :: * -> *) a. Monad m => a -> m a
return (Glex String -> f (Glex String)) -> Glex String -> f (Glex String)
forall a b. (a -> b) -> a -> b
$ (String -> Glex String
forall (m :: * -> *) v. MonomialConstructor m => v -> m v
mvar "t" :: Glex String)
        hs :: [Vect k (Elim2 (Glex String) m)]
hs = (Vect k m -> Vect k (Elim2 (Glex String) m))
-> [Vect k m] -> [Vect k (Elim2 (Glex String) m)]
forall a b. (a -> b) -> [a] -> [b]
map ((Vect k (Elim2 (Glex String) m)
forall (f :: * -> *) b.
(Mon b, Monad f) =>
f (Elim2 (Glex String) b)
t Vect k (Elim2 (Glex String) m)
-> Vect k (Elim2 (Glex String) m) -> Vect k (Elim2 (Glex String) m)
forall a. Num a => a -> a -> a
*) (Vect k (Elim2 (Glex String) m) -> Vect k (Elim2 (Glex String) m))
-> (Vect k m -> Vect k (Elim2 (Glex String) m))
-> Vect k m
-> Vect k (Elim2 (Glex String) m)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vect k m -> Vect k (Elim2 (Glex String) m)
forall (f :: * -> *) a b.
(Functor f, Mon a) =>
f b -> f (Elim2 a b)
toElimSnd) [Vect k m]
fs [Vect k (Elim2 (Glex String) m)]
-> [Vect k (Elim2 (Glex String) m)]
-> [Vect k (Elim2 (Glex String) m)]
forall a. [a] -> [a] -> [a]
++ (Vect k m -> Vect k (Elim2 (Glex String) m))
-> [Vect k m] -> [Vect k (Elim2 (Glex String) m)]
forall a b. (a -> b) -> [a] -> [b]
map (((1Vect k (Elim2 (Glex String) m)
-> Vect k (Elim2 (Glex String) m) -> Vect k (Elim2 (Glex String) m)
forall a. Num a => a -> a -> a
-Vect k (Elim2 (Glex String) m)
forall (f :: * -> *) b.
(Mon b, Monad f) =>
f (Elim2 (Glex String) b)
t) Vect k (Elim2 (Glex String) m)
-> Vect k (Elim2 (Glex String) m) -> Vect k (Elim2 (Glex String) m)
forall a. Num a => a -> a -> a
*) (Vect k (Elim2 (Glex String) m) -> Vect k (Elim2 (Glex String) m))
-> (Vect k m -> Vect k (Elim2 (Glex String) m))
-> Vect k m
-> Vect k (Elim2 (Glex String) m)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vect k m -> Vect k (Elim2 (Glex String) m)
forall (f :: * -> *) a b.
(Functor f, Mon a) =>
f b -> f (Elim2 a b)
toElimSnd) [Vect k m]
gs
    in [Vect k (Elim2 (Glex String) m)] -> [Vect k m]
forall b a b.
(Fractional b, Monomial a, Monomial b, Ord b, Ord a, Ord b) =>
[Vect b (Elim2 a b)] -> [Vect b b]
eliminateFst [Vect k (Elim2 (Glex String) m)]
hs

toElimFst :: f a -> f (Elim2 a b)
toElimFst = (a -> Elim2 a b) -> f a -> f (Elim2 a b)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\m :: a
m -> a -> b -> Elim2 a b
forall a b. a -> b -> Elim2 a b
Elim2 a
m b
forall m. Mon m => m
munit)
toElimSnd :: f b -> f (Elim2 a b)
toElimSnd = (b -> Elim2 a b) -> f b -> f (Elim2 a b)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\m :: b
m -> a -> b -> Elim2 a b
forall a b. a -> b -> Elim2 a b
Elim2 a
forall m. Mon m => m
munit b
m)
isElimFst :: Vect b (Elim2 a b) -> Bool
isElimFst = (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
forall m. Mon m => m
munit) (a -> Bool)
-> (Vect b (Elim2 a b) -> a) -> Vect b (Elim2 a b) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (\(Elim2 m :: a
m _) -> a
m) (Elim2 a b -> a)
-> (Vect b (Elim2 a b) -> Elim2 a b) -> Vect b (Elim2 a b) -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vect b (Elim2 a b) -> Elim2 a b
forall b c. Vect b c -> c
lm
fromElimSnd :: f (Elim2 a b) -> f b
fromElimSnd = (Elim2 a b -> b) -> f (Elim2 a b) -> f b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\(Elim2 _ m :: b
m) -> b
m)
eliminateFst :: [Vect b (Elim2 a b)] -> [Vect b b]
eliminateFst = (Vect b (Elim2 a b) -> Vect b b)
-> [Vect b (Elim2 a b)] -> [Vect b b]
forall a b. (a -> b) -> [a] -> [b]
map Vect b (Elim2 a b) -> Vect b b
forall (f :: * -> *) a b. Functor f => f (Elim2 a b) -> f b
fromElimSnd ([Vect b (Elim2 a b)] -> [Vect b b])
-> ([Vect b (Elim2 a b)] -> [Vect b (Elim2 a b)])
-> [Vect b (Elim2 a b)]
-> [Vect b b]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vect b (Elim2 a b) -> Bool)
-> [Vect b (Elim2 a b)] -> [Vect b (Elim2 a b)]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile Vect b (Elim2 a b) -> Bool
forall a b b. (Eq a, Mon a) => Vect b (Elim2 a b) -> Bool
isElimFst ([Vect b (Elim2 a b)] -> [Vect b (Elim2 a b)])
-> ([Vect b (Elim2 a b)] -> [Vect b (Elim2 a b)])
-> [Vect b (Elim2 a b)]
-> [Vect b (Elim2 a b)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Vect b (Elim2 a b)] -> [Vect b (Elim2 a b)]
forall k m.
(Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
[Vect k m] -> [Vect k m]
gb


-- Cox et al, p193-4
-- |Given ideals I and J, their quotient is defined as I:J = {f | f \<- R, f.g is in I for all g in J}.
--
-- If fs and gs are generators for I and J, then @quotientI fs gs@ returns generators for I:J.
--
-- The ideal quotient is the algebraic analogue of the Zariski closure of a difference of varieties.
-- V(I:J) contains the Zariski closure of V(I)-V(J), with equality if k is algebraically closed and I is a radical ideal.
quotientI :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
    [Vect k m] -> [Vect k m] -> [Vect k m]
quotientI :: [Vect k m] -> [Vect k m] -> [Vect k m]
quotientI _ [] = [1]
quotientI fs :: [Vect k m]
fs gs :: [Vect k m]
gs = ([Vect k m] -> [Vect k m] -> [Vect k m])
-> [[Vect k m]] -> [Vect k m]
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 [Vect k m] -> [Vect k m] -> [Vect k m]
forall k m.
(Fractional k, Ord k, Monomial m, Ord m) =>
[Vect k m] -> [Vect k m] -> [Vect k m]
intersectI ([[Vect k m]] -> [Vect k m]) -> [[Vect k m]] -> [Vect k m]
forall a b. (a -> b) -> a -> b
$ (Vect k m -> [Vect k m]) -> [Vect k m] -> [[Vect k m]]
forall a b. (a -> b) -> [a] -> [b]
map ([Vect k m] -> Vect k m -> [Vect k m]
forall m k.
(Monomial m, Fractional k, Algebra k m, Ord m, Ord k) =>
[Vect k m] -> Vect k m -> [Vect k m]
quotientP [Vect k m]
fs) [Vect k m]
gs
-- quotientI fs gs = foldl intersectI [1] $ map (quotientP fs) gs

quotientP :: [Vect k m] -> Vect k m -> [Vect k m]
quotientP fs :: [Vect k m]
fs g :: Vect k m
g = (Vect k m -> Vect k m) -> [Vect k m] -> [Vect k m]
forall a b. (a -> b) -> [a] -> [b]
map ( Vect k m -> Vect k m -> Vect k m
forall m b.
(Monomial m, Fractional b, Ord m, Algebra b m, Eq b) =>
Vect b m -> Vect b m -> Vect b m
// Vect k m
g ) ([Vect k m] -> [Vect k m]) -> [Vect k m] -> [Vect k m]
forall a b. (a -> b) -> a -> b
$ [Vect k m] -> [Vect k m] -> [Vect k m]
forall k m.
(Fractional k, Ord k, Monomial m, Ord m) =>
[Vect k m] -> [Vect k m] -> [Vect k m]
intersectI [Vect k m]
fs [Vect k m
g]
    where h :: Vect b m
h // :: Vect b m -> Vect b m -> Vect b m
// g :: Vect b m
g = let ([u :: Vect b m
u],_) = Vect b m -> [Vect b m] -> ([Vect b m], Vect b m)
forall b m.
(Monomial m, Fractional b, Eq b, Ord m, Algebra b m) =>
Vect b m -> [Vect b m] -> ([Vect b m], Vect b m)
quotRemMP Vect b m
h [Vect b m
g] in Vect b m
u

-- |@eliminate vs gs@ returns the elimination ideal obtained from the ideal generated by gs by eliminating the variables vs.
eliminate :: (Eq k, Fractional k, Ord k, MonomialConstructor m, Monomial (m v), Ord (m v)) =>
    [Vect k (m v)] -> [Vect k (m v)] -> [Vect k (m v)]
eliminate :: [Vect k (m v)] -> [Vect k (m v)] -> [Vect k (m v)]
eliminate vs :: [Vect k (m v)]
vs gs :: [Vect k (m v)]
gs = let subs :: v -> Vect k (Elim2 (m v) (m v))
subs = [Vect k (m v)] -> v -> Vect k (Elim2 (m v) (m v))
forall k (m :: * -> *) v.
(Eq k, Num k, MonomialConstructor m, Eq (m v), Mon (m v)) =>
[Vect k (m v)] -> v -> Vect k (Elim2 (m v) (m v))
subFst [Vect k (m v)]
vs in [Vect k (Elim2 (m v) (m v))] -> [Vect k (m v)]
forall b a b.
(Fractional b, Monomial a, Monomial b, Ord b, Ord a, Ord b) =>
[Vect b (Elim2 a b)] -> [Vect b b]
eliminateFst [Vect k (m v)
g Vect k (m v)
-> (v -> Vect k (Elim2 (m v) (m v))) -> Vect k (Elim2 (m v) (m v))
forall k (m :: * -> *) a v.
(Eq k, Num k, MonomialConstructor m, Ord a, Show a, Algebra k a) =>
Vect k (m v) -> (v -> Vect k a) -> Vect k a
`bind` v -> Vect k (Elim2 (m v) (m v))
subs | Vect k (m v)
g <- [Vect k (m v)]
gs]
    where subFst :: (Eq k, Num k, MonomialConstructor m, Eq (m v), Mon (m v)) =>
              [Vect k (m v)] -> v -> Vect k (Elim2 (m v) (m v))
          subFst :: [Vect k (m v)] -> v -> Vect k (Elim2 (m v) (m v))
subFst vs :: [Vect k (m v)]
vs = (\v :: v
v -> let v' :: Vect k (m v)
v' = v -> Vect k (m v)
forall k (m :: * -> *) v.
(Num k, MonomialConstructor m) =>
v -> Vect k (m v)
var v
v in if Vect k (m v)
forall k (m :: * -> *).
(Num k, MonomialConstructor m) =>
Vect k (m v)
v' Vect k (m v) -> [Vect k (m v)] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Vect k (m v)]
vs then Vect k (m v) -> Vect k (Elim2 (m v) (m v))
forall (f :: * -> *) b a.
(Functor f, Mon b) =>
f a -> f (Elim2 a b)
toElimFst Vect k (m v)
forall k (m :: * -> *).
(Num k, MonomialConstructor m) =>
Vect k (m v)
v' else Vect k (m v) -> Vect k (Elim2 (m v) (m v))
forall (f :: * -> *) a b.
(Functor f, Mon a) =>
f b -> f (Elim2 a b)
toElimSnd Vect k (m v)
forall k (m :: * -> *).
(Num k, MonomialConstructor m) =>
Vect k (m v)
v')

{-
-- !! NOT WORKING
-- |@elimExcept vs gs@ returns the elimination ideal obtained from the ideal generated by gs by eliminating all variables except vs.
elimExcept :: (Fractional k, Ord k, MonomialConstructor m, Monomial (m v), Ord (m v)) =>
    [Vect k (m v)] -> [Vect k (m v)] -> [Vect k (m v)]
elimExcept vs gs = let subs = subSnd vs in eliminateFst [g `bind` subs | g <- gs]
    where subSnd :: (Num k, MonomialConstructor m, Eq (m v), Mon (m v)) =>
              [Vect k (m v)] -> v -> Vect k (Elim2 (m v) (m v))
          subSnd vs = (\v -> let v' = var v in if v' `elem` vs then toElimSnd v' else toElimFst v')
-}

-- MONOMIAL BASES FOR QUOTIENT ALGEBRAS

-- basis for the polynomial ring in variables vs
mbasis :: [a] -> [a]
mbasis vs :: [a]
vs = [a] -> [a]
mbasis' [1]
    where mbasis' :: [a] -> [a]
mbasis' ms :: [a]
ms = [a]
ms [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a] -> [a]
mbasis' ([a] -> [a]
forall a. Ord a => [a] -> [a]
toSet [a
va -> a -> a
forall a. Num a => a -> a -> a
*a
m | a
v <- [a]
vs, a
m <- [a]
ms])

-- |Given variables vs, and a Groebner basis gs, @mbasisQA vs gs@ returns a monomial basis for the quotient algebra k[vs]/\<gs\>.
-- For example, @mbasisQA [x,y] [x^2+y^2-1]@ returns a monomial basis for k[x,y]/\<x^2+y^2-1\>.
-- In general, the monomial basis is likely to be infinite.
mbasisQA :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
     [Vect k m] -> [Vect k m] -> [Vect k m]
mbasisQA :: [Vect k m] -> [Vect k m] -> [Vect k m]
mbasisQA vs :: [Vect k m]
vs gs :: [Vect k m]
gs = [Vect k m] -> [Vect k m]
mbasisQA' [1]
    where mbasisQA' :: [Vect k m] -> [Vect k m]
mbasisQA' [] = []  -- the quotient algebra is finite-dimensional
          mbasisQA' ms :: [Vect k m]
ms = [Vect k m]
ms [Vect k m] -> [Vect k m] -> [Vect k m]
forall a. [a] -> [a] -> [a]
++ [Vect k m] -> [Vect k m]
mbasisQA' ([Vect k m] -> [Vect k m]
forall a. Ord a => [a] -> [a]
toSet [Vect k m
f | Vect k m
v <- [Vect k m]
vs, Vect k m
m <- [Vect k m]
ms, let f :: Vect k m
f = Vect k m
vVect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
*Vect k m
m, Vect k m
f Vect k m -> [Vect k m] -> Vect k m
forall k m.
(Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
Vect k m -> [Vect k m] -> Vect k m
%% [Vect k m]
gs Vect k m -> Vect k m -> Bool
forall a. Eq a => a -> a -> Bool
== Vect k m
f])

-- |Given an ideal I, the leading term ideal lt(I) consists of the leading terms of all elements of I.
-- If I is generated by gs, then @ltIdeal gs@ returns generators for lt(I).
ltIdeal :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
    [Vect k m] -> [Vect k m]
ltIdeal :: [Vect k m] -> [Vect k m]
ltIdeal gs :: [Vect k m]
gs = (Vect k m -> Vect k m) -> [Vect k m] -> [Vect k m]
forall a b. (a -> b) -> [a] -> [b]
map (m -> Vect k m
forall (m :: * -> *) a. Monad m => a -> m a
return (m -> Vect k m) -> (Vect k m -> m) -> Vect k m -> Vect k m
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vect k m -> m
forall b c. Vect b c -> c
lm) ([Vect k m] -> [Vect k m]) -> [Vect k m] -> [Vect k m]
forall a b. (a -> b) -> a -> b
$ [Vect k m] -> [Vect k m]
forall k m.
(Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
[Vect k m] -> [Vect k m]
gb [Vect k m]
gs

-- number of monomials of degree i in n variables
numMonomials :: a -> a -> Integer
numMonomials n :: a
n i :: a
i = a -> Integer
forall a. Integral a => a -> Integer
toInteger (a
ia -> a -> a
forall a. Num a => a -> a -> a
+a
na -> a -> a
forall a. Num a => a -> a -> a
-1) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`choose` a -> Integer
forall a. Integral a => a -> Integer
toInteger (a
na -> a -> a
forall a. Num a => a -> a -> a
-1)

-- |Given variables vs, and a homogeneous ideal gs, @hilbertFunQA vs gs@ returns the Hilbert function for the quotient algebra k[vs]/\<gs\>.
-- Given an integer i, the Hilbert function returns the number of degree i monomials in a basis for k[vs]/\<gs\>.
-- For a homogeneous ideal, this number is independent of the monomial ordering used
-- (even though the elements of the monomial basis themselves are dependent on the ordering).
--
-- If the ideal I is not homogeneous, then R/I is not graded, and the Hilbert function is not well-defined.
-- Specifically, the number of degree i monomials in a basis is likely to depend on which monomial ordering you use.
hilbertFunQA :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
    [Vect k m] -> [Vect k m] -> Int -> Integer
hilbertFunQA :: [Vect k m] -> [Vect k m] -> Int -> Integer
hilbertFunQA vs :: [Vect k m]
vs gs :: [Vect k m]
gs i :: Int
i = [Vect k m] -> Int -> Integer
forall m k.
(Monomial m, Fractional k, Algebra k m, Ord m, Ord k) =>
[Vect k m] -> Int -> Integer
hilbertFunQA' ([Vect k m] -> [Vect k m]
forall k m.
(Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
[Vect k m] -> [Vect k m]
ltIdeal [Vect k m]
gs) Int
i
    where n :: Int
n = [Vect k m] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Vect k m]
vs
          hilbertFunQA' :: [Vect k m] -> Int -> Integer
hilbertFunQA' _ i :: Int
i | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< 0 = 0
          hilbertFunQA' (m :: Vect k m
m:ms :: [Vect k m]
ms) i :: Int
i = [Vect k m] -> Int -> Integer
hilbertFunQA' [Vect k m]
ms Int
i Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- [Vect k m] -> Int -> Integer
hilbertFunQA' ([Vect k m]
ms [Vect k m] -> Vect k m -> [Vect k m]
forall m k.
(Monomial m, Fractional k, Algebra k m, Ord m, Ord k) =>
[Vect k m] -> Vect k m -> [Vect k m]
`quotientP` Vect k m
m) (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Vect k m -> Int
forall m k. Monomial m => Vect k m -> Int
deg Vect k m
m)
          hilbertFunQA' [] i :: Int
i = Int -> Int -> Integer
forall a. Integral a => a -> a -> Integer
numMonomials Int
n Int
i
-- For example, consider k[x,y]/<x-y^2>
-- Under Lex ordering, the monomial basis is 1,y,y^2,y^3,...
-- Under Glex ordering, the monomial basis is 1,x,y,x^2,xy,x^3,x^2y,...
-- So the Hilbert function is not well-defined.
-- Note though that this function does still correctly return the number of degree i monomials for the given monomial ordering

-- naive implementation which simply counts monomials
hilbertSeriesQA1 :: [Vect k m] -> [Vect k m] -> [Int]
hilbertSeriesQA1 vs :: [Vect k m]
vs gs :: [Vect k m]
gs = [Vect k m] -> [Int]
hilbertSeriesQA1' [1]
    where hilbertSeriesQA1' :: [Vect k m] -> [Int]
hilbertSeriesQA1' [] = [] -- repeat 0
          hilbertSeriesQA1' ms :: [Vect k m]
ms = [Vect k m] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Vect k m]
ms Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Vect k m] -> [Int]
hilbertSeriesQA1' ([Vect k m] -> [Vect k m]
forall a. Ord a => [a] -> [a]
toSet [Vect k m
f | Vect k m
v <- [Vect k m]
vs, Vect k m
m <- [Vect k m]
ms, let f :: Vect k m
f = Vect k m
vVect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
*Vect k m
m, Vect k m
f Vect k m -> [Vect k m] -> Vect k m
forall k m.
(Eq k, Fractional k, Monomial m, Ord m, Algebra k m) =>
Vect k m -> [Vect k m] -> Vect k m
%% [Vect k m]
gs Vect k m -> Vect k m -> Bool
forall a. Eq a => a -> a -> Bool
== Vect k m
f])

-- Eisenbud p325, p357 / Schenck p56
-- This can be made more efficient by choosing which m to recurse on
-- |Given variables vs, and a homogeneous ideal gs, @hilbertSeriesQA vs gs@ returns the Hilbert series for the quotient algebra k[vs]/\<gs\>.
-- The Hilbert series should be interpreted as a formal power series where the coefficient of t^i is the Hilbert function evaluated at i.
-- That is, the i'th element in the series is the number of degree i monomials in a basis for k[vs]/\<gs\>.
hilbertSeriesQA :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
    [Vect k m] -> [Vect k m] -> [Integer]
hilbertSeriesQA :: [Vect k m] -> [Vect k m] -> [Integer]
hilbertSeriesQA vs :: [Vect k m]
vs gs :: [Vect k m]
gs = [Vect k m] -> [Integer]
forall m k.
(Monomial m, Fractional k, Ord k, Ord m, Algebra k m) =>
[Vect k m] -> [Integer]
hilbertSeriesQA' ([Vect k m] -> [Integer]) -> [Vect k m] -> [Integer]
forall a b. (a -> b) -> a -> b
$ [Vect k m] -> [Vect k m]
forall k m.
(Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
[Vect k m] -> [Vect k m]
ltIdeal [Vect k m]
gs
    where hilbertSeriesQA' :: [Vect k m] -> [Integer]
hilbertSeriesQA' (m :: Vect k m
m:ms :: [Vect k m]
ms) = [Vect k m] -> [Integer]
hilbertSeriesQA' [Vect k m]
ms [Integer] -> [Integer] -> [Integer]
forall b. Num b => [b] -> [b] -> [b]
<-> (Int -> Integer -> [Integer]
forall a. Int -> a -> [a]
replicate (Vect k m -> Int
forall m k. Monomial m => Vect k m -> Int
deg Vect k m
m) 0 [Integer] -> [Integer] -> [Integer]
forall a. [a] -> [a] -> [a]
++ [Vect k m] -> [Integer]
hilbertSeriesQA' ([Vect k m]
ms [Vect k m] -> [Vect k m] -> [Vect k m]
forall k m.
(Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
[Vect k m] -> [Vect k m] -> [Vect k m]
`quotientI` [Vect k m
m]))
          hilbertSeriesQA' [] = [Int -> Int -> Integer
forall a. Integral a => a -> a -> Integer
numMonomials Int
n Int
i | Int
i <- [0..] ]
          n :: Int
n = [Vect k m] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Vect k m]
vs
          (a :: b
a:as :: [b]
as) <-> :: [b] -> [b] -> [b]
<-> (b :: b
b:bs :: [b]
bs) = (b
ab -> b -> b
forall a. Num a => a -> a -> a
-b
b) b -> [b] -> [b]
forall a. a -> [a] -> [a]
: ([b]
as [b] -> [b] -> [b]
<-> [b]
bs)
          as :: [b]
as <-> [] = [b]
as
          [] <-> bs :: [b]
bs = (b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map b -> b
forall a. Num a => a -> a
negate [b]
bs

-- |In the case where every variable v occurs in some generator g of the homogeneous ideal (the usual case),
-- then the vs can be inferred from the gs.
-- @hilbertSeriesQA' gs@ returns the Hilbert series for the quotient algebra k[vs]/\<gs\>.
hilbertSeriesQA' :: (Fractional k, Ord k, MonomialConstructor m, Ord (m v), Monomial (m v), Algebra k (m v)) =>
    [Vect k (m v)] -> [Integer]
hilbertSeriesQA' :: [Vect k (m v)] -> [Integer]
hilbertSeriesQA' gs :: [Vect k (m v)]
gs = [Vect k (m v)] -> [Vect k (m v)] -> [Integer]
forall k m.
(Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
[Vect k m] -> [Vect k m] -> [Integer]
hilbertSeriesQA [Vect k (m v)]
vs [Vect k (m v)]
gs where vs :: [Vect k (m v)]
vs = [Vect k (m v)] -> [Vect k (m v)]
forall a. Ord a => [a] -> [a]
toSet ((Vect k (m v) -> [Vect k (m v)])
-> [Vect k (m v)] -> [Vect k (m v)]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap Vect k (m v) -> [Vect k (m v)]
forall k (m :: * -> *) v.
(Num k, Ord k, MonomialConstructor m, Ord (m v)) =>
Vect k (m v) -> [Vect k (m v)]
vars [Vect k (m v)]
gs)

-- |For i \>\> 0, the Hilbert function becomes a polynomial in i, called the Hilbert polynomial.
hilbertPolyQA :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
     [Vect k m] -> [Vect k m] -> GlexPoly Q String
hilbertPolyQA :: [Vect k m] -> [Vect k m] -> GlexPoly Q String
hilbertPolyQA vs :: [Vect k m]
vs gs :: [Vect k m]
gs = [Vect k m] -> GlexPoly Q String -> GlexPoly Q String
forall m a k.
(Monomial m, Fractional a, Fractional k, Algebra k m, Ord m,
 Ord k) =>
[Vect k m] -> a -> a
hilbertPolyQA' ([Vect k m] -> [Vect k m]
forall k m.
(Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
[Vect k m] -> [Vect k m]
ltIdeal [Vect k m]
gs) GlexPoly Q String
i
    where n :: Integer
n = Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int -> Integer) -> Int -> Integer
forall a b. (a -> b) -> a -> b
$ [Vect k m] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Vect k m]
vs
          i :: GlexPoly Q String
i = String -> GlexPoly Q String
forall v. v -> GlexPoly Q v
glexvar "i"
          hilbertPolyQA' :: [Vect k m] -> a -> a
hilbertPolyQA' [] x :: a
x = [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ a
x a -> a -> a
forall a. Num a => a -> a -> a
+ Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
j | Integer
j <- [1..Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-1] ] a -> a -> a
forall a. Fractional a => a -> a -> a
/ (Integer -> a
forall a. Num a => Integer -> a
fromInteger (Integer -> a) -> Integer -> a
forall a b. (a -> b) -> a -> b
$ [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [1..Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-1])
          hilbertPolyQA' (m :: Vect k m
m:ms :: [Vect k m]
ms) x :: a
x = [Vect k m] -> a -> a
hilbertPolyQA' [Vect k m]
ms a
x a -> a -> a
forall a. Num a => a -> a -> a
- [Vect k m] -> a -> a
hilbertPolyQA' ([Vect k m]
ms [Vect k m] -> Vect k m -> [Vect k m]
forall m k.
(Monomial m, Fractional k, Algebra k m, Ord m, Ord k) =>
[Vect k m] -> Vect k m -> [Vect k m]
`quotientP` Vect k m
m) (a
x a -> a -> a
forall a. Num a => a -> a -> a
- Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Vect k m -> Int
forall m k. Monomial m => Vect k m -> Int
deg Vect k m
m))

hilbertPolyQA' :: (Fractional k, Ord k, MonomialConstructor m, Ord (m v), Monomial (m v), Algebra k (m v)) =>
    [Vect k (m v)] -> GlexPoly Q String
hilbertPolyQA' :: [Vect k (m v)] -> GlexPoly Q String
hilbertPolyQA' gs :: [Vect k (m v)]
gs = [Vect k (m v)] -> [Vect k (m v)] -> GlexPoly Q String
forall k m.
(Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
[Vect k m] -> [Vect k m] -> GlexPoly Q String
hilbertPolyQA [Vect k (m v)]
vs [Vect k (m v)]
gs where vs :: [Vect k (m v)]
vs = [Vect k (m v)] -> [Vect k (m v)]
forall a. Ord a => [a] -> [a]
toSet ((Vect k (m v) -> [Vect k (m v)])
-> [Vect k (m v)] -> [Vect k (m v)]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap Vect k (m v) -> [Vect k (m v)]
forall k (m :: * -> *) v.
(Num k, Ord k, MonomialConstructor m, Ord (m v)) =>
Vect k (m v) -> [Vect k (m v)]
vars [Vect k (m v)]
gs)

-- The dimension of a variety
dim :: [Vect k m] -> [Vect k m] -> Int
dim vs :: [Vect k m]
vs gs :: [Vect k m]
gs = 1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ GlexPoly Q String -> Int
forall m k. Monomial m => Vect k m -> Int
deg ([Vect k m] -> [Vect k m] -> GlexPoly Q String
forall k m.
(Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>
[Vect k m] -> [Vect k m] -> GlexPoly Q String
hilbertPolyQA [Vect k m]
vs [Vect k m]
gs)

dim' :: [Vect k (m v)] -> Int
dim' gs :: [Vect k (m v)]
gs = 1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ GlexPoly Q String -> Int
forall m k. Monomial m => Vect k m -> Int
deg ([Vect k (m v)] -> GlexPoly Q String
forall k (m :: * -> *) v.
(Fractional k, Ord k, MonomialConstructor m, Ord (m v),
 Monomial (m v), Algebra k (m v)) =>
[Vect k (m v)] -> GlexPoly Q String
hilbertPolyQA' [Vect k (m v)]
gs)