{-# OPTIONS_HADDOCK show-extensions #-}
{-# LANGUAGE ScopedTypeVariables #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Data.Polynomial.GroebnerBasis
-- Copyright   :  (c) Masahiro Sakai 2012-2013
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- Gröbner basis
--
-- References:
--
-- * Monomial order <http://en.wikipedia.org/wiki/Monomial_order>
--
-- * Gröbner basis <http://en.wikipedia.org/wiki/Gr%C3%B6bner_basis>
--
-- * グレブナー基底 <http://d.hatena.ne.jp/keyword/%A5%B0%A5%EC%A5%D6%A5%CA%A1%BC%B4%F0%C4%EC>
--
-- * Gröbner Bases and Buchberger’s Algorithm <http://math.rice.edu/~cbruun/vigre/vigreHW6.pdf>
--
-- * Docon <http://www.haskell.org/docon/>
--
-----------------------------------------------------------------------------

module ToySolver.Data.Polynomial.GroebnerBasis
  (
  -- * Options
    Options (..)
  , Strategy (..)

  -- * Gröbner basis computation
  , basis
  , basis'
  , spolynomial
  , reduceGBasis
  ) where

import Data.Default.Class
import qualified Data.Set as Set
import qualified Data.Heap as H -- http://hackage.haskell.org/package/heaps
import ToySolver.Data.Polynomial.Base (Polynomial, Monomial, MonomialOrder)
import qualified ToySolver.Data.Polynomial.Base as P

-- | Options for Gröbner Basis computation.
--
-- The default option can be obtained by 'def'.
data Options
  = Options
  { Options -> Strategy
optStrategy :: Strategy
  }

instance Default Options where
  def :: Options
def =
    Options
    { optStrategy :: Strategy
optStrategy = Strategy
NormalStrategy
    }

data Strategy
  = NormalStrategy
  | SugarStrategy  -- ^ sugar strategy (not implemented yet)
  deriving (Strategy -> Strategy -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Strategy -> Strategy -> Bool
$c/= :: Strategy -> Strategy -> Bool
== :: Strategy -> Strategy -> Bool
$c== :: Strategy -> Strategy -> Bool
Eq, Eq Strategy
Strategy -> Strategy -> Bool
Strategy -> Strategy -> Ordering
Strategy -> Strategy -> Strategy
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: Strategy -> Strategy -> Strategy
$cmin :: Strategy -> Strategy -> Strategy
max :: Strategy -> Strategy -> Strategy
$cmax :: Strategy -> Strategy -> Strategy
>= :: Strategy -> Strategy -> Bool
$c>= :: Strategy -> Strategy -> Bool
> :: Strategy -> Strategy -> Bool
$c> :: Strategy -> Strategy -> Bool
<= :: Strategy -> Strategy -> Bool
$c<= :: Strategy -> Strategy -> Bool
< :: Strategy -> Strategy -> Bool
$c< :: Strategy -> Strategy -> Bool
compare :: Strategy -> Strategy -> Ordering
$ccompare :: Strategy -> Strategy -> Ordering
Ord, Int -> Strategy -> ShowS
[Strategy] -> ShowS
Strategy -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Strategy] -> ShowS
$cshowList :: [Strategy] -> ShowS
show :: Strategy -> String
$cshow :: Strategy -> String
showsPrec :: Int -> Strategy -> ShowS
$cshowsPrec :: Int -> Strategy -> ShowS
Show, ReadPrec [Strategy]
ReadPrec Strategy
Int -> ReadS Strategy
ReadS [Strategy]
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Strategy]
$creadListPrec :: ReadPrec [Strategy]
readPrec :: ReadPrec Strategy
$creadPrec :: ReadPrec Strategy
readList :: ReadS [Strategy]
$creadList :: ReadS [Strategy]
readsPrec :: Int -> ReadS Strategy
$creadsPrec :: Int -> ReadS Strategy
Read, Strategy
forall a. a -> a -> Bounded a
maxBound :: Strategy
$cmaxBound :: Strategy
minBound :: Strategy
$cminBound :: Strategy
Bounded, Int -> Strategy
Strategy -> Int
Strategy -> [Strategy]
Strategy -> Strategy
Strategy -> Strategy -> [Strategy]
Strategy -> Strategy -> Strategy -> [Strategy]
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
enumFromThenTo :: Strategy -> Strategy -> Strategy -> [Strategy]
$cenumFromThenTo :: Strategy -> Strategy -> Strategy -> [Strategy]
enumFromTo :: Strategy -> Strategy -> [Strategy]
$cenumFromTo :: Strategy -> Strategy -> [Strategy]
enumFromThen :: Strategy -> Strategy -> [Strategy]
$cenumFromThen :: Strategy -> Strategy -> [Strategy]
enumFrom :: Strategy -> [Strategy]
$cenumFrom :: Strategy -> [Strategy]
fromEnum :: Strategy -> Int
$cfromEnum :: Strategy -> Int
toEnum :: Int -> Strategy
$ctoEnum :: Int -> Strategy
pred :: Strategy -> Strategy
$cpred :: Strategy -> Strategy
succ :: Strategy -> Strategy
$csucc :: Strategy -> Strategy
Enum)

spolynomial
  :: (Eq k, Fractional k, Ord v)
  => MonomialOrder v -> Polynomial k v -> Polynomial k v -> Polynomial k v
spolynomial :: forall k v.
(Eq k, Fractional k, Ord v) =>
MonomialOrder v
-> Polynomial k v -> Polynomial k v -> Polynomial k v
spolynomial MonomialOrder v
cmp Polynomial k v
f Polynomial k v
g =
      forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
P.fromTerm ((k
1,Monomial v
xs) forall k v.
(Fractional k, Ord v) =>
Term k v -> Term k v -> Term k v
`P.tdiv` Term k v
lt1) forall a. Num a => a -> a -> a
* Polynomial k v
f
    forall a. Num a => a -> a -> a
- forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
P.fromTerm ((k
1,Monomial v
xs) forall k v.
(Fractional k, Ord v) =>
Term k v -> Term k v -> Term k v
`P.tdiv` Term k v
lt2) forall a. Num a => a -> a -> a
* Polynomial k v
g
  where
    xs :: Monomial v
xs = forall v. Ord v => Monomial v -> Monomial v -> Monomial v
P.mlcm Monomial v
xs1 Monomial v
xs2
    lt1 :: Term k v
lt1@(k
c1, Monomial v
xs1) = forall k v. Num k => MonomialOrder v -> Polynomial k v -> Term k v
P.lt MonomialOrder v
cmp Polynomial k v
f
    lt2 :: Term k v
lt2@(k
c2, Monomial v
xs2) = forall k v. Num k => MonomialOrder v -> Polynomial k v -> Term k v
P.lt MonomialOrder v
cmp Polynomial k v
g

basis
  :: forall k v. (Eq k, Fractional k, Ord k, Ord v)
  => MonomialOrder v
  -> [Polynomial k v]
  -> [Polynomial k v]
basis :: forall k v.
(Eq k, Fractional k, Ord k, Ord v) =>
MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]
basis = forall k v.
(Eq k, Fractional k, Ord k, Ord v) =>
Options -> MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]
basis' forall a. Default a => a
def

basis'
  :: forall k v. (Eq k, Fractional k, Ord k, Ord v)
  => Options
  -> MonomialOrder v
  -> [Polynomial k v]
  -> [Polynomial k v]
basis' :: forall k v.
(Eq k, Fractional k, Ord k, Ord v) =>
Options -> MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]
basis' Options
opt MonomialOrder v
cmp [Polynomial k v]
fs =
  forall k v.
(Eq k, Ord k, Fractional k, Ord v) =>
MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]
reduceGBasis MonomialOrder v
cmp forall a b. (a -> b) -> a -> b
$ [Polynomial k v] -> Heap (Item k v) -> [Polynomial k v]
go [Polynomial k v]
fs (forall a. Ord a => [a] -> Heap a
H.fromList [forall k v.
(Eq k, Num k, Ord v) =>
MonomialOrder v -> Polynomial k v -> Polynomial k v -> Item k v
item MonomialOrder v
cmp Polynomial k v
fi Polynomial k v
fj | (Polynomial k v
fi,Polynomial k v
fj) <- forall a. [a] -> [(a, a)]
pairs [Polynomial k v]
fs, forall {k} {k}.
(Num k, Num k) =>
Polynomial k v -> Polynomial k v -> Bool
checkGCD Polynomial k v
fi Polynomial k v
fj])
  where
    go :: [Polynomial k v] -> H.Heap (Item k v) -> [Polynomial k v]
    go :: [Polynomial k v] -> Heap (Item k v) -> [Polynomial k v]
go [Polynomial k v]
gs Heap (Item k v)
h | forall a. Heap a -> Bool
H.null Heap (Item k v)
h = [Polynomial k v]
gs
    go [Polynomial k v]
gs Heap (Item k v)
h
      | Polynomial k v
r forall a. Eq a => a -> a -> Bool
== Polynomial k v
0    = [Polynomial k v] -> Heap (Item k v) -> [Polynomial k v]
go [Polynomial k v]
gs Heap (Item k v)
h'
      | Bool
otherwise = [Polynomial k v] -> Heap (Item k v) -> [Polynomial k v]
go (Polynomial k v
rforall a. a -> [a] -> [a]
:[Polynomial k v]
gs) (forall a. Heap a -> Heap a -> Heap a
H.union Heap (Item k v)
h' (forall a. Ord a => [a] -> Heap a
H.fromList [forall k v.
(Eq k, Num k, Ord v) =>
MonomialOrder v -> Polynomial k v -> Polynomial k v -> Item k v
item MonomialOrder v
cmp Polynomial k v
r Polynomial k v
g | Polynomial k v
g <- [Polynomial k v]
gs, forall {k} {k}.
(Num k, Num k) =>
Polynomial k v -> Polynomial k v -> Bool
checkGCD  Polynomial k v
r Polynomial k v
g]))
      where
        Just (Item k v
i, Heap (Item k v)
h') = forall a. Heap a -> Maybe (a, Heap a)
H.viewMin Heap (Item k v)
h
        fi :: Polynomial k v
fi = forall k v. Item k v -> Polynomial k v
iFst Item k v
i
        fj :: Polynomial k v
fj = forall k v. Item k v -> Polynomial k v
iSnd Item k v
i
        spoly :: Polynomial k v
spoly = forall k v.
(Eq k, Fractional k, Ord v) =>
MonomialOrder v
-> Polynomial k v -> Polynomial k v -> Polynomial k v
spolynomial MonomialOrder v
cmp Polynomial k v
fi Polynomial k v
fj
        r :: Polynomial k v
r = forall k v.
(Eq k, Fractional k, Ord v) =>
MonomialOrder v
-> Polynomial k v -> [Polynomial k v] -> Polynomial k v
P.reduce MonomialOrder v
cmp Polynomial k v
spoly [Polynomial k v]
gs

    -- gcdが1となる組は選ばなくて良い
    checkGCD :: Polynomial k v -> Polynomial k v -> Bool
checkGCD Polynomial k v
fi Polynomial k v
fj = Bool -> Bool
not forall a b. (a -> b) -> a -> b
$ forall v. Ord v => Monomial v -> Monomial v -> Bool
P.mcoprime (forall k v.
Num k =>
MonomialOrder v -> Polynomial k v -> Monomial v
P.lm MonomialOrder v
cmp Polynomial k v
fi) (forall k v.
Num k =>
MonomialOrder v -> Polynomial k v -> Monomial v
P.lm MonomialOrder v
cmp Polynomial k v
fj)

reduceGBasis
  :: forall k v. (Eq k, Ord k, Fractional k, Ord v)
  => MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]
reduceGBasis :: forall k v.
(Eq k, Ord k, Fractional k, Ord v) =>
MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]
reduceGBasis MonomialOrder v
cmp [Polynomial k v]
ps = forall a. Set a -> [a]
Set.toList forall a b. (a -> b) -> a -> b
$ forall a. Ord a => [a] -> Set a
Set.fromList forall a b. (a -> b) -> a -> b
$ forall {k}.
(Fractional k, Eq k) =>
[Polynomial k v] -> [Polynomial k v] -> [Polynomial k v]
go [Polynomial k v]
ps []
  where
    go :: [Polynomial k v] -> [Polynomial k v] -> [Polynomial k v]
go [] [Polynomial k v]
qs = [Polynomial k v]
qs
    go (Polynomial k v
p:[Polynomial k v]
ps) [Polynomial k v]
qs
      | Polynomial k v
q forall a. Eq a => a -> a -> Bool
== Polynomial k v
0    = [Polynomial k v] -> [Polynomial k v] -> [Polynomial k v]
go [Polynomial k v]
ps [Polynomial k v]
qs
      | Bool
otherwise = [Polynomial k v] -> [Polynomial k v] -> [Polynomial k v]
go [Polynomial k v]
ps (forall r v.
(Eq r, Fractional r) =>
MonomialOrder v -> Polynomial r v -> Polynomial r v
P.toMonic MonomialOrder v
cmp Polynomial k v
q forall a. a -> [a] -> [a]
: [Polynomial k v]
qs)
      where
        q :: Polynomial k v
q = forall k v.
(Eq k, Fractional k, Ord v) =>
MonomialOrder v
-> Polynomial k v -> [Polynomial k v] -> Polynomial k v
P.reduce MonomialOrder v
cmp Polynomial k v
p ([Polynomial k v]
psforall a. [a] -> [a] -> [a]
++[Polynomial k v]
qs)

{--------------------------------------------------------------------
  Item
--------------------------------------------------------------------}

data Item k v
  = Item
  { forall k v. Item k v -> Polynomial k v
iFst :: Polynomial k v
  , forall k v. Item k v -> Polynomial k v
iSnd :: Polynomial k v
  , forall k v. Item k v -> MonomialOrder v
iCmp :: MonomialOrder v
  , forall k v. Item k v -> Monomial v
iLCM :: Monomial v
  }

item :: (Eq k, Num k, Ord v) => MonomialOrder v -> Polynomial k v -> Polynomial k v -> Item k v
item :: forall k v.
(Eq k, Num k, Ord v) =>
MonomialOrder v -> Polynomial k v -> Polynomial k v -> Item k v
item MonomialOrder v
cmp Polynomial k v
f Polynomial k v
g = forall k v.
Polynomial k v
-> Polynomial k v -> MonomialOrder v -> Monomial v -> Item k v
Item Polynomial k v
f Polynomial k v
g MonomialOrder v
cmp (forall v. Ord v => Monomial v -> Monomial v -> Monomial v
P.mlcm (forall k v.
Num k =>
MonomialOrder v -> Polynomial k v -> Monomial v
P.lm MonomialOrder v
cmp Polynomial k v
f) (forall k v.
Num k =>
MonomialOrder v -> Polynomial k v -> Monomial v
P.lm MonomialOrder v
cmp Polynomial k v
g))

instance Ord v => Ord (Item k v) where
  Item k v
a compare :: Item k v -> Item k v -> Ordering
`compare` Item k v
b = forall k v. Item k v -> MonomialOrder v
iCmp Item k v
a (forall k v. Item k v -> Monomial v
iLCM Item k v
a) (forall k v. Item k v -> Monomial v
iLCM Item k v
b)

instance Ord v => Eq (Item k v) where
  Item k v
a == :: Item k v -> Item k v -> Bool
== Item k v
b = forall a. Ord a => a -> a -> Ordering
compare Item k v
a Item k v
b forall a. Eq a => a -> a -> Bool
== Ordering
EQ

{--------------------------------------------------------------------
  Utilities
--------------------------------------------------------------------}

pairs :: [a] -> [(a,a)]
pairs :: forall a. [a] -> [(a, a)]
pairs [] = []
pairs (a
x:[a]
xs) = [(a
x,a
y) | a
y <- [a]
xs] forall a. [a] -> [a] -> [a]
++ forall a. [a] -> [(a, a)]
pairs [a]
xs