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

{-# OPTIONS_GHC -fglasgow-exts #-}

module Math.Algebra.Commutative.GBasis where

import Data.List
import qualified Data.Map as M

import Math.Algebra.Commutative.Monomial
import Math.Algebra.Commutative.MPoly

-- 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 f g = let h = lcmT (lt f) (lt g)
            in h `divT` lt f .* f - h `divT` lt g .* g
-- The point about the s-poly is that it cancels out the leading terms of the two polys, exposing their second terms

isGB fs = all (\h -> h %% fs == 0) (pairWith sPoly fs)


-- Cox p87
gb1 fs = gb' fs (pairWith sPoly fs) where
    gb' gs (h:hs) = let h' = h %% gs in
                        if h' == 0 then gb' gs hs else gb' (h':gs) (hs ++ map (sPoly h') gs)
    gb' gs [] = gs

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


-- Cox p89-90
reduce gs = reduce' [] gs where
    reduce' gs' (g:gs) | g' == 0   = reduce' gs' gs
                       | otherwise = reduce' (g':gs') gs
                       where g' = g %% (gs'++gs)
    reduce' gs' [] = reverse $ sort $ map toMonic gs'
-- the reverse means that 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 fs = reduce $ gb' fs (pairs [1..s]) s where
    s = length fs
    gb' gs ((i,j):ps) t =
        let fi = gs!i; fj = gs!j in
        if coprimeM (lm fi) (lm fj) || criterion fi fj
        -- if lcmM (lm fi) (lm fj) == lm fi * lm fj || criterion fi fj
        then gb' gs ps t
        else let h = sPoly fi fj %% gs in
             if h == 0 then gb' gs ps t else gb' (gs++[h]) (ps ++ [ (i,t+1) | i <- [1..t] ]) (t+1)
        where
            criterion fi fj = let l = lcmM (lm fi) (lm fj) in any (test l) [1..t]
            test l k = k `notElem` [i,j]
                    && ordpair i k `notElem` ps
                    && ordpair j k `notElem` ps
                    && lm (gs!k) `dividesM` l
    gb' gs [] _ = gs


pairs (x:xs) = map (\y->(x,y)) xs ++ pairs xs
pairs [] = []

xs ! i = xs !! (i-1) -- in other words, index the list from 1 not 0

ordpair x y | x < y     = (x,y)
            | otherwise = (y,x)

-- version of gb2 where we eliminate pairs as they're created, rather than as they're processed
gb2b fs = reduce $ gb1' [] fs [] 0 where
    gb1' gs (f:fs) ps t = gb1' (gs ++ [f]) fs ps' (t+1) where
        ps' = updatePairs gs ps f t
    gb1' ls [] ps t = gb2' ls ps t
    gb2' gs ((i,j):ps) t =
        let h = sPoly (gs!i) (gs!j) %% gs in
        if h == 0
        then gb2' gs ps t
        else let ps' = updatePairs gs ((i,j):ps) h t in gb2' (gs++[h]) ps' (t+1)
    gb2' gs [] _ = gs
    updatePairs gs ps f t =
        [p | p@(i,j) <- ps,
             not (lm f `dividesM` lcmM (lm (gs!i)) (lm (gs!j)))] 
     ++ [ (i,t+1) | (gi,i) <- zip gs [1..t],
                    not (coprimeM (lm gi) (lm f)),
                    not (criterion (lcmM (lm gi) (lm f)) i) ]
        where criterion l i = any (`dividesM` l) [lm gk | (gk,k) <- zip gs [1..t], k /= i, ordpair i k `notElem` 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")
gb3b fs =
    let fs' = sort $ filter (/=0) fs
    in reduce $ gb1' [] fs' [] 0 where
    gb1' gs (f:fs) ps t = gb1' (gs ++ [f]) fs ps' (t+1) where
        ps' = updatePairs gs ps f t
    gb1' ls [] ps t = gb2' ls ps t
    gb2' gs ((l,(i,j)):ps) t =
        let h = sPoly (gs!i) (gs!j) %% gs in
        if h == 0
        then gb2' gs ps t
        else let ps' = updatePairs gs ((l,(i,j)):ps) h t in gb2' (gs++[h]) ps' (t+1)
    gb2' gs [] _ = gs
    updatePairs :: (Ord (Monomial ord), Fractional r) => [MPoly ord r] -> [(Monomial ord, (Int,Int))] -> (MPoly ord r) -> Int -> [(Monomial ord, (Int,Int))]
    updatePairs gs ps f t =
        mergeBy cmpFst
        [p | p@(l,(i,j)) <- ps,
             not (lm f `dividesM` l)]
        $ sortBy cmpFst
        [ (l,(i,t+1)) | (gi,i) <- zip gs [1..t], l <- [lcmM (lm gi) (lm f)],
                    not (coprimeM (lm gi) (lm f)),
                    not (criterion l i) ]
        where criterion l i = any (`dividesM` l) [lm gk | (gk,k) <- zip gs [1..t], k /= i, ordpair i k `notElem` map snd ps]

cmpFst (a,_) (b,_) = compare a b

mergeBy cmp (t:ts) (u:us) =
    case cmp t u of
    LT -> t : mergeBy cmp ts (u:us)
    EQ -> t : mergeBy cmp ts (u:us)
    GT -> u : mergeBy cmp (t:ts) us
mergeBy _ ts us = ts ++ us -- one of them is null


-- naive implementation of "sugar strategy"
gb4b fs =
    let fs' = sort $ filter (/=0) fs
    in reduce $ gb1' [] fs' [] 0 where
    gb1' gs (f:fs) ps t = gb1' (gs ++ [f]) fs ps' (t+1) where
        ps' = updatePairs gs ps f t
    gb1' ls [] ps t = gb2' ls ps t
    gb2' gs ((sl,(i,j)):ps) t =
        let h = sPoly (gs!i) (gs!j) %% gs in
        if h == 0
        then gb2' gs ps t
        else let ps' = updatePairs gs ((sl,(i,j)):ps) h t in gb2' (gs++[h]) ps' (t+1)
    gb2' gs [] _ = gs
    updatePairs :: (Ord (Monomial ord), Fractional r) =>
                   [MPoly ord r] -> [((Int,Monomial ord), (Int,Int))] -> (MPoly ord r) -> Int -> [((Int,Monomial ord), (Int,Int))]
    updatePairs gs ps f t =
        mergeBy cmpFst
        [p | p@((s,l),(i,j)) <- ps,
             not (lm f `dividesM` l)]
        $ sortBy cmpFst
        [ ((s,l),(i,t+1)) | (gi,i) <- zip gs [1..t], l <- [lcmM (lm gi) (lm f)], s <- [sugar gi f l],
                    not (coprimeM (lm gi) (lm f)),
                    not (criterion l i) ]
        where criterion l i = any (`dividesM` l) [lm gk | (gk,k) <- zip gs [1..t], k /= i, ordpair i k `notElem` map snd ps]


-- 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

-- |Given a list of polynomials over a field, return a Groebner basis for the ideal generated by the polynomials
gb :: (Ord (Monomial ord), Fractional k, Ord k) =>
     [MPoly ord k] -> [MPoly ord k]
gb fs =
    -- let fs' = sort $ filter (/=0) fs
    let fs' = sort $ map toMonic $ filter (/=0) fs
    in reduce $ gb1' [] fs' [] 0 where
    gb1' gs (f:fs) ps t = gb1' (gs ++ [f]) fs ps' (t+1)
        where ps' = updatePairs gs ps f (t+1)
    gb1' ls [] ps t = gb2' ls ps t
    gb2' gs (p@(_,(i,j)):ps) t =
        if h == 0
        then gb2' gs ps t
        else gb2' (gs++[h]) ps' (t+1)
        where h = toMonic $ sPoly (gs!i) (gs!j) %% gs
              ps' = updatePairs gs (p:ps) h (t+1)
    gb2' gs [] _ = gs
    updatePairs :: (Ord (Monomial ord), Fractional r) =>
                   [MPoly ord r] -> [((Int,Monomial ord), (Int,Int))] -> (MPoly ord r) -> Int -> [((Int,Monomial ord), (Int,Int))]
    updatePairs gs ps gk k =
        let newps = [let l = lcmM (lm gi) (lm gk) in ((sugar gi gk l, l), (i,k)) | (gi,i) <- zip gs [1..k-1] ]
            ps' = [p | p@((sij,tij),(i,j)) <- ps, ((sik,tik),_) <- [newps ! i], ((sjk,tjk),_) <- [newps ! j],
                       not ( (tik `properlyDividesM` tij) && (tjk `properlyDividesM` tij) ) ] -- sloppy variant
            newps' = discard1 [] newps
            newps'' = sortBy cmpSug $ discard2 [] $ sortBy cmpNormal newps'
        in mergeBy cmpSug ps' newps''
        where
            discard1 ls (r@((_sik,tik),(i,_k)):rs) =
                if lm (gs!i) `coprimeM` lm gk
                -- then discard [l | l@((_,tjk),_) <- ls, tjk /= tik] [r | r@((_,tjk),_) <- ls, tjk /= tik]
                then discard1 (filter (\((_,tjk),_) -> tjk /= tik) ls) (filter (\((_,tjk),_) -> tjk /= tik) rs)
                else discard1 (r:ls) rs
            discard1 ls [] = ls
            discard2 ls (r@((_sik,tik),(i,k)):rs) = discard2 (r:ls) $ filter (\((_sjk,tjk),(j,k')) -> not (k == k' && tik `dividesM` tjk)) rs
            discard2 ls [] = ls
-- The type annotation on updatePairs appears to be required
-- 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 f g h = degM h + max (deg f - degM (lm f)) (deg g - degM (lm g))

cmpNormal ((s1,t1),(i1,j1)) ((s2,t2),(i2,j2)) = compare (t1,j1) (t2,j2)

cmpSug ((s1,t1),(i1,j1)) ((s2,t2),(i2,j2)) = compare (s1,t1,j1) (s2,t2,j2)


-- earlier version of gb3b
gb3 fs =
    let gs = sort $ filter (/=0) fs
        ps = sortBy cmpFst $ pairWith (\(i,fi) (j,fj) -> (lcmM (lm fi) (lm fj), (i,j)) ) $ zip [1..] gs
    in reduce $ gb' gs ps s
    where
        s = length fs
        gb' :: (Ord (Monomial ord), Fractional r) => [MPoly ord r] -> [(Monomial ord, (Int,Int))] -> Int -> [MPoly ord r]
        gb' gs ((l,(i,j)):ps) t =
            let fi = gs!i; fj = gs!j in
            if coprimeM (lm fi) (lm fj) || any (test l) [1..t]
            then gb' gs ps t
            else let h = sPoly fi fj %% gs in
                 if h == 0
                 then gb' gs ps t
                 else let ps' = mergeBy cmpFst ps $ sortBy cmpFst $ zip [lcmM (lm h) (lm fi) | fi <- gs] [(i,t+1) | i <- [1..t]]
                 -- else let ps' = mergeBy cmpFst ps $ zip [lcmM (lm h) (lm fi) | fi <- gs] [(i,t+1) | i <- [1..t]]
                      in gb' (gs++[h]) ps' (t+1)
            where
                test l k = k `notElem` [i,j]
                        && ordpair i k `notElem` map snd ps
                        && ordpair j k `notElem` map snd ps
                        && lm (gs!k) `dividesM` l
        gb' gs [] _ = gs
-- Note that the type annotation on gb' appears to be required. I think this is a bug in the type inference algorithm



-- earlier version of gb4b
gb4 fs =
    let gs = sort $ filter (/=0) fs
        ps = sortBy cmpFst $ pairWith (\(i,fi) (j,fj) -> let l = lcmM (lm fi) (lm fj) in ((sugar fi fj l, l), (i,j)) ) $ zip [1..] gs
    in reduce $ gb' gs ps s
    where
        s = length fs
        gb' :: (Ord (Monomial ord), Fractional r) => [MPoly ord r] -> [((Int,Monomial ord), (Int,Int))] -> Int -> [MPoly ord r]
        gb' gs (((s,l),(i,j)):ps) t =
            let fi = gs!i; fj = gs!j in
            if coprimeM (lm fi) (lm fj) || any (test l) [1..t]
            then gb' gs ps t
            else let h = sPoly fi fj %% gs in
                 if h == 0
                 then gb' gs ps t
                 else let ps' = mergeBy cmpFst ps $ sortBy cmpFst $ zip [let l = lcmM (lm fi) (lm h) in (sugar fi h l, l) | fi <- gs] [(i,t+1) | i <- [1..t]]
                      in gb' (gs++[h]) ps' (t+1)
            where
                test l k = k `notElem` [i,j]
                        && ordpair i k `notElem` map snd ps
                        && ordpair j k `notElem` map snd ps
                        && lm (gs!k) `dividesM` l
        gb' gs [] _ = gs
-- Note that the type annotation on gb' appears to be required. I think this is a bug in the type inference algorithm



{-
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

-- Cox et al, p181
-- Geometric interpretation: V(I+J) = V(I) `intersect` V(J)
sumI fs gs = gb $ fs ++ gs

-- Cox et al, p183
-- Geometric interpretation: V(I.J) = V(I) `union` V(J)
productI fs gs = gb [f * g | f <- fs, g <- gs]