module ZkFold.Base.Algebra.Polynomials.Multivariate.Groebner where

import           Data.Bool                                               (bool)
import           Data.List                                               (sortBy)
import           GHC.Natural                                             (Natural)
import           Prelude                                                 hiding (Num (..), drop, lcm, length, sum, take,
                                                                          (!!), (/))

import           ZkFold.Base.Algebra.Basic.Class
import           ZkFold.Base.Algebra.Polynomials.Multivariate.Monomial
import           ZkFold.Base.Algebra.Polynomials.Multivariate.Polynomial
import           ZkFold.Prelude                                          (length, (!!))

reducable :: Polynomial c i j  => Poly c i j -> Poly c i j -> Bool
reducable :: forall c i j. Polynomial c i j => Poly c i j -> Poly c i j -> Bool
reducable Poly c i j
l Poly c i j
r = Mono i j -> Mono i j -> Bool
forall i j. Monomial i j => Mono i j -> Mono i j -> Bool
dividable ((c, Mono i j) -> Mono i j
forall a b. (a, b) -> b
snd ((c, Mono i j) -> Mono i j) -> (c, Mono i j) -> Mono i j
forall a b. (a -> b) -> a -> b
$ Poly c i j -> (c, Mono i j)
forall c i j. Polynomial c i j => Poly c i j -> (c, Mono i j)
lt Poly c i j
l) ((c, Mono i j) -> Mono i j
forall a b. (a, b) -> b
snd ((c, Mono i j) -> Mono i j) -> (c, Mono i j) -> Mono i j
forall a b. (a -> b) -> a -> b
$ Poly c i j -> (c, Mono i j)
forall c i j. Polynomial c i j => Poly c i j -> (c, Mono i j)
lt Poly c i j
r)

-- TODO: refactor reduction methods so that they never applied to non-reducible polynomials

reduce ::
    forall c i j . (Ring j, Polynomial c i j)
    =>Poly c i j
    -> Poly c i j
    -> Poly c i j
reduce :: forall c i j.
(Ring j, Polynomial c i j) =>
Poly c i j -> Poly c i j -> Poly c i j
reduce Poly c i j
l Poly c i j
r =
    let (c
cl, Mono i j
ml) = Poly c i j -> (c, Mono i j)
forall c i j. Polynomial c i j => Poly c i j -> (c, Mono i j)
lt Poly c i j
l
        (c
cr, Mono i j
mr) = Poly c i j -> (c, Mono i j)
forall c i j. Polynomial c i j => Poly c i j -> (c, Mono i j)
lt Poly c i j
r
    in Poly c i j
l Poly c i j -> Poly c i j -> Poly c i j
forall a. AdditiveGroup a => a -> a -> a
-  (c, Mono i j) -> Poly c i j -> Poly c i j
forall c i j.
Polynomial c i j =>
(c, Mono i j) -> Poly c i j -> Poly c i j
scaleM (c
cl c -> c -> c
forall a. Field a => a -> a -> a
// c
cr, Mono i j
ml Mono i j -> Mono i j -> Mono i j
forall a. MultiplicativeGroup a => a -> a -> a
/ Mono i j
mr) Poly c i j
r

reduceMany ::
       forall c i j . (Ring j, Polynomial c i j)
    => Poly c i j
    -> [Poly c i j]
    -> Poly c i j
reduceMany :: forall c i j.
(Ring j, Polynomial c i j) =>
Poly c i j -> [Poly c i j] -> Poly c i j
reduceMany Poly c i j
h [Poly c i j]
fs = if Bool
reduced then Poly c i j -> [Poly c i j] -> Poly c i j
forall c i j.
(Ring j, Polynomial c i j) =>
Poly c i j -> [Poly c i j] -> Poly c i j
reduceMany Poly c i j
h' [Poly c i j]
fs else Poly c i j
h'
  where
    (Poly c i j
h', Bool
reduced) = Poly c i j -> [Poly c i j] -> Bool -> (Poly c i j, Bool)
reduceStep Poly c i j
h [Poly c i j]
fs Bool
False
    reduceStep :: Poly c i j -> [Poly c i j] -> Bool -> (Poly c i j, Bool)
reduceStep Poly c i j
p (Poly c i j
q:[Poly c i j]
qs) Bool
r
      | Poly c i j -> Bool
forall c i j. Poly c i j -> Bool
zeroP Poly c i j
p   = (Poly c i j
h, Bool
r)
      | Bool
otherwise =
        if Poly c i j -> Poly c i j -> Bool
forall c i j. Polynomial c i j => Poly c i j -> Poly c i j -> Bool
reducable Poly c i j
p Poly c i j
q
          then (Poly c i j -> Poly c i j -> Poly c i j
forall c i j.
(Ring j, Polynomial c i j) =>
Poly c i j -> Poly c i j -> Poly c i j
reduce Poly c i j
p Poly c i j
q, Bool
True)
          else Poly c i j -> [Poly c i j] -> Bool -> (Poly c i j, Bool)
reduceStep Poly c i j
p [Poly c i j]
qs Bool
r
    reduceStep Poly c i j
p [] Bool
r = (Poly c i j
p, Bool
r)

fullReduceMany ::
       forall c i j . (Ring j, Polynomial c i j)
    => Poly c i j
    -> [Poly c i j]
    -> Poly c i j
fullReduceMany :: forall c i j.
(Ring j, Polynomial c i j) =>
Poly c i j -> [Poly c i j] -> Poly c i j
fullReduceMany Poly c i j
h [Poly c i j]
fs =
    let h' :: Poly c i j
h' = Poly c i j -> [Poly c i j] -> Poly c i j
forall c i j.
(Ring j, Polynomial c i j) =>
Poly c i j -> [Poly c i j] -> Poly c i j
reduceMany Poly c i j
h [Poly c i j]
fs
    in case Poly c i j
h' of
        P []         -> Poly c i j
h'
        P ((c
c, Mono i j
m):[(c, Mono i j)]
_) -> [(c, Mono i j)] -> Poly c i j
forall c i j. [(c, Mono i j)] -> Poly c i j
P [(c
c, Mono i j
m)] Poly c i j -> Poly c i j -> Poly c i j
forall a. AdditiveSemigroup a => a -> a -> a
+ Poly c i j -> [Poly c i j] -> Poly c i j
forall c i j.
(Ring j, Polynomial c i j) =>
Poly c i j -> [Poly c i j] -> Poly c i j
fullReduceMany (Poly c i j
h' Poly c i j -> Poly c i j -> Poly c i j
forall a. AdditiveGroup a => a -> a -> a
- [(c, Mono i j)] -> Poly c i j
forall c i j. [(c, Mono i j)] -> Poly c i j
P [(c
c, Mono i j
m)]) [Poly c i j]
fs

systemReduce ::
       forall c i j . (Ring j, Polynomial c i j)
    => [Poly c i j]
    -> [Poly c i j]
systemReduce :: forall c i j.
(Ring j, Polynomial c i j) =>
[Poly c i j] -> [Poly c i j]
systemReduce = (Poly c i j -> [Poly c i j] -> [Poly c i j])
-> [Poly c i j] -> [Poly c i j] -> [Poly c i j]
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: Type -> Type) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Poly c i j -> [Poly c i j] -> [Poly c i j]
forall {c} {i} {j}.
(Field c, Ord i, Ord j, Eq c, Ring j) =>
Poly c i j -> [Poly c i j] -> [Poly c i j]
f []
    where
        f :: Poly c i j -> [Poly c i j] -> [Poly c i j]
f Poly c i j
p [Poly c i j]
ps =
            let p' :: Poly c i j
p' = Poly c i j -> [Poly c i j] -> Poly c i j
forall c i j.
(Ring j, Polynomial c i j) =>
Poly c i j -> [Poly c i j] -> Poly c i j
fullReduceMany Poly c i j
p [Poly c i j]
ps
            in [Poly c i j] -> [Poly c i j] -> Bool -> [Poly c i j]
forall a. a -> a -> Bool -> a
bool [Poly c i j]
ps (Poly c i j
p' Poly c i j -> [Poly c i j] -> [Poly c i j]
forall a. a -> [a] -> [a]
: [Poly c i j]
ps) (Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ Poly c i j -> Bool
forall c i j. Poly c i j -> Bool
zeroP Poly c i j
p')

data GroebnerParams c i j = GroebnerParams
    { forall c i j. GroebnerParams c i j -> Natural
groebnerMaxSteps      :: Natural
    , forall c i j. GroebnerParams c i j -> Poly c i j -> Bool
groebnerSPolySelector :: Poly c i j -> Bool
    }

makeSPoly ::
       (Ring j, Polynomial c i j)
    => Poly c i j
    -> Poly c i j
    -> Poly c i j
makeSPoly :: forall j c i.
(Ring j, Polynomial c i j) =>
Poly c i j -> Poly c i j -> Poly c i j
makeSPoly Poly c i j
l Poly c i j
r =
    let (c
cl, Mono i j
ml) = Poly c i j -> (c, Mono i j)
forall c i j. Polynomial c i j => Poly c i j -> (c, Mono i j)
lt Poly c i j
l
        (c
cr, Mono i j
mr) = Poly c i j -> (c, Mono i j)
forall c i j. Polynomial c i j => Poly c i j -> (c, Mono i j)
lt Poly c i j
r

        M Map i j
as = Mono i j -> Mono i j -> Mono i j
forall i j. Monomial i j => Mono i j -> Mono i j -> Mono i j
gcdM Mono i j
ml Mono i j
mr
        lcm :: Mono i j
lcm = Mono i j -> Mono i j -> Mono i j
forall i j. Monomial i j => Mono i j -> Mono i j -> Mono i j
lcmM Mono i j
ml Mono i j
mr

        ra :: Mono i j
ra  = Mono i j
lcm Mono i j -> Mono i j -> Mono i j
forall a. MultiplicativeGroup a => a -> a -> a
/ Mono i j
ml
        la :: Mono i j
la  = Mono i j
lcm Mono i j -> Mono i j -> Mono i j
forall a. MultiplicativeGroup a => a -> a -> a
/ Mono i j
mr

        l' :: Poly c i j
l'  = (c
cr, Mono i j
ra) (c, Mono i j) -> Poly c i j -> Poly c i j
forall c i j.
Polynomial c i j =>
(c, Mono i j) -> Poly c i j -> Poly c i j
`scaleM` Poly c i j
l
        r' :: Poly c i j
r'  = (c
cl, Mono i j
la) (c, Mono i j) -> Poly c i j -> Poly c i j
forall c i j.
Polynomial c i j =>
(c, Mono i j) -> Poly c i j -> Poly c i j
`scaleM` Poly c i j
r
    in if Map i j -> Bool
forall a. Map i a -> Bool
forall (t :: Type -> Type) a. Foldable t => t a -> Bool
null Map i j
as
        then Poly c i j
forall a. AdditiveMonoid a => a
zero
        else Poly c i j
r' Poly c i j -> Poly c i j -> Poly c i j
forall a. AdditiveGroup a => a -> a -> a
- Poly c i j
l'

groebnerStep ::
        (Ring j, Polynomial c i j)
    => GroebnerParams c i j
    -> [Poly c i j]
    -> [Poly c i j]
groebnerStep :: forall j c i.
(Ring j, Polynomial c i j) =>
GroebnerParams c i j -> [Poly c i j] -> [Poly c i j]
groebnerStep GroebnerParams{Natural
Poly c i j -> Bool
groebnerMaxSteps :: forall c i j. GroebnerParams c i j -> Natural
groebnerSPolySelector :: forall c i j. GroebnerParams c i j -> Poly c i j -> Bool
groebnerMaxSteps :: Natural
groebnerSPolySelector :: Poly c i j -> Bool
..} [Poly c i j]
ps =
    let n :: Natural
n = [Poly c i j] -> Natural
forall (t :: Type -> Type) a. Foldable t => t a -> Natural
length [Poly c i j]
ps
        inds :: [(Natural, Natural)]
inds = [(Natural
i, Natural
j) | Natural
i <- [Natural
0 .. Natural
nNatural -> Natural -> Natural
-!Natural
1], Natural
j <- [Natural
0 .. Natural
nNatural -> Natural -> Natural
-!Natural
1], Natural
i Natural -> Natural -> Bool
forall a. Ord a => a -> a -> Bool
< Natural
j]
        ss :: [Poly c i j]
ss   = ((Natural, Natural) -> Poly c i j)
-> [(Natural, Natural)] -> [Poly c i j]
forall a b. (a -> b) -> [a] -> [b]
map (\(Natural
i, Natural
j) -> Poly c i j -> Poly c i j -> Poly c i j
forall j c i.
(Ring j, Polynomial c i j) =>
Poly c i j -> Poly c i j -> Poly c i j
makeSPoly ([Poly c i j]
ps [Poly c i j] -> Natural -> Poly c i j
forall a. [a] -> Natural -> a
!! Natural
i) ([Poly c i j]
ps [Poly c i j] -> Natural -> Poly c i j
forall a. [a] -> Natural -> a
!! Natural
j) Poly c i j -> [Poly c i j] -> Poly c i j
forall c i j.
(Ring j, Polynomial c i j) =>
Poly c i j -> [Poly c i j] -> Poly c i j
`reduceMany` [Poly c i j]
ps) [(Natural, Natural)]
inds
        ss' :: [Poly c i j]
ss'  = (Poly c i j -> Bool) -> [Poly c i j] -> [Poly c i j]
forall a. (a -> Bool) -> [a] -> [a]
filter (Bool -> Bool
not (Bool -> Bool) -> (Poly c i j -> Bool) -> Poly c i j -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Poly c i j -> Bool
forall c i j. Poly c i j -> Bool
zeroP) [Poly c i j]
ss
        ss'' :: [Poly c i j]
ss'' = (Poly c i j -> Bool) -> [Poly c i j] -> [Poly c i j]
forall a. (a -> Bool) -> [a] -> [a]
filter Poly c i j -> Bool
groebnerSPolySelector [Poly c i j]
ss'
        ps' :: [Poly c i j]
ps'  = (Poly c i j -> Poly c i j -> Ordering)
-> [Poly c i j] -> [Poly c i j]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy ((Poly c i j -> Poly c i j -> Ordering)
-> Poly c i j -> Poly c i j -> Ordering
forall a b c. (a -> b -> c) -> b -> a -> c
flip Poly c i j -> Poly c i j -> Ordering
forall a. Ord a => a -> a -> Ordering
compare) ([Poly c i j]
ss'' [Poly c i j] -> [Poly c i j] -> [Poly c i j]
forall a. [a] -> [a] -> [a]
++ [Poly c i j]
ps)
    in [Poly c i j] -> [Poly c i j]
forall c i j.
(Ring j, Polynomial c i j) =>
[Poly c i j] -> [Poly c i j]
systemReduce [Poly c i j]
ps'

groebner ::
       forall c i j . (Ring j, Polynomial c i j)
    => GroebnerParams c i j
    -> [Poly c i j]
    -> [Poly c i j]
groebner :: forall c i j.
(Ring j, Polynomial c i j) =>
GroebnerParams c i j -> [Poly c i j] -> [Poly c i j]
groebner GroebnerParams{Natural
Poly c i j -> Bool
groebnerMaxSteps :: forall c i j. GroebnerParams c i j -> Natural
groebnerSPolySelector :: forall c i j. GroebnerParams c i j -> Poly c i j -> Bool
groebnerMaxSteps :: Natural
groebnerSPolySelector :: Poly c i j -> Bool
..} [Poly c i j]
ps =
    let step :: [Poly c i j] -> [Poly c i j]
step = GroebnerParams c i j -> [Poly c i j] -> [Poly c i j]
forall j c i.
(Ring j, Polynomial c i j) =>
GroebnerParams c i j -> [Poly c i j] -> [Poly c i j]
groebnerStep GroebnerParams{Natural
Poly c i j -> Bool
groebnerMaxSteps :: Natural
groebnerSPolySelector :: Poly c i j -> Bool
groebnerMaxSteps :: Natural
groebnerSPolySelector :: Poly c i j -> Bool
..}
        go :: Natural -> [Poly c i j] -> [Poly c i j]
go Natural
0 [Poly c i j]
ps' = [Poly c i j]
ps'
        go Natural
n [Poly c i j]
ps' = Natural -> [Poly c i j] -> [Poly c i j]
go (Natural
n Natural -> Natural -> Natural
-! Natural
1) ([Poly c i j] -> [Poly c i j]) -> [Poly c i j] -> [Poly c i j]
forall a b. (a -> b) -> a -> b
$ [Poly c i j] -> [Poly c i j]
step [Poly c i j]
ps'
    in Natural -> [Poly c i j] -> [Poly c i j]
go Natural
groebnerMaxSteps [Poly c i j]
ps

verifyGroebner ::
       forall c i j . (Ring j, Polynomial c i j)
    => GroebnerParams c i j
    -> [Poly c i j]
    -> Poly c i j
    -> Bool
verifyGroebner :: forall c i j.
(Ring j, Polynomial c i j) =>
GroebnerParams c i j -> [Poly c i j] -> Poly c i j -> Bool
verifyGroebner GroebnerParams c i j
params [Poly c i j]
ps = Poly c i j -> Bool
forall c i j. Poly c i j -> Bool
zeroP (Poly c i j -> Bool)
-> (Poly c i j -> Poly c i j) -> Poly c i j -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Poly c i j -> [Poly c i j] -> Poly c i j
forall c i j.
(Ring j, Polynomial c i j) =>
Poly c i j -> [Poly c i j] -> Poly c i j
`fullReduceMany` GroebnerParams c i j -> [Poly c i j] -> [Poly c i j]
forall c i j.
(Ring j, Polynomial c i j) =>
GroebnerParams c i j -> [Poly c i j] -> [Poly c i j]
groebner GroebnerParams c i j
params [Poly c i j]
ps)