hspray 0.1.3.0 → 0.2.0.0

raw patch · 6 files changed

+588/−21 lines, 6 filesdep ~containersdep ~hashabledep ~numeric-preludePVP ok

version bump matches the API change (PVP)

Dependency ranges changed: containers, hashable, numeric-prelude, text, unordered-containers

API changes (from Hackage documentation)

Files

• CHANGELOG.md +14/−0
• README.md +46/−0
• hspray.cabal +10/−7
• src/Math/Algebra/Hspray.hs +458/−12
• tests/Approx.hs +10/−0
• tests/Main.hs +50/−2

@@ -32,3 +32,17 @@ 
 * Completed the README to show how to deal with symbolic coefficients.
 
+
+## 0.2.0.0 - 2024-03-14
+
+* New functions `prettySpray'` and `prettySprayXYZ`.
+
+* New function `substituteSpray`.
+
+* New function `sprayDivision`, to perform the division of a spray by a list of sprays.
+
+* New function `groebner`, to compute a Groebner basis of a list of sprays.
+
+* New function `isSymmetricSpray`, to check whether a spray is a symmetric polynomial.
+
+* New function `isPolynomialOf`, to check whether a spray can be expressed as a polynomial of a given list of sprays.

@@ -60,6 +60,24 @@ -- 8.0
 ```
 
+Partial evaluation:
+
+```haskell
+import Math.Algebra.Hspray
+import Data.Ratio
+x1 = lone 1 :: Spray Rational
+x2 = lone 2 :: Spray Rational
+x3 = lone 3 :: Spray Rational
+poly = x1^**^2 ^+^ x2 ^+^ x3 ^-^ unitSpray
+prettySpray' poly
+-- "((-1) % 1) + (1 % 1) x3 + (1 % 1) x2 + (1 % 1) x1^2"
+--
+-- substitute x1 -> 2 and x3 -> 3
+poly' = substituteSpray [Just 2, Nothing, Just 3] p
+prettySpray' poly'
+-- "(6 % 1) + (1 % 1) x2"
+```
+
 Differentiation:
 
 ```haskell
@@ -72,6 +90,34 @@ prettySpray show "X" $ derivSpray 1 poly
 -- "(2.0) * X^(0, 1, 1) + (6.0) * X^(1)"
 ```
+
+## Grôbner bases
+
+As of version 2.0.0, it is possible to compute a Gröbner basis.
+
+```haskell
+import Math.Algebra.Hspray
+import Data.Ratio
+-- define the elementary monomials
+o = lone 0 :: Spray Rational
+x = lone 1 :: Spray Rational
+y = lone 2 :: Spray Rational
+z = lone 3 :: Spray Rational
+-- define three polynomials
+p1 = x^**^2 ^+^ y ^+^ z ^-^ o -- X² + Y + Z - 1
+p2 = x ^+^ y^**^2 ^+^ z ^-^ o -- X + Y² + Z - 1
+p3 = x ^+^ y ^+^ z^**^2 ^-^ o -- X + Y + Z² - 1
+-- compute the reduced Gröbner basis
+gbasis = groebner [p1, p2, p3] True
+-- show result
+prettyResult = map prettySprayXYZ gbasis
+mapM_ print prettyResult
+-- "((-1) % 1) + (1 % 1) Z^2 + (1 % 1) Y + (1 % 1) X"
+-- "(1 % 1) Z + ((-1) % 1) Z^2 + ((-1) % 1) Y + (1 % 1) Y^2"
+-- "((-1) % 2) Z^2 + (1 % 2) Z^4 + (1 % 1) YZ^2"
+-- "((-1) % 1) Z^2 + (4 % 1) Z^3 + ((-4) % 1) Z^4 + (1 % 1) Z^6"
+```
+
 
 ## Easier usage 
 

@@ -1,7 +1,7 @@ name:                hspray
-version:             0.1.3.0
+version:             0.2.0.0
 synopsis:            Multivariate polynomials.
-description:         Manipulation of multivariate polynomials on a ring.
+description:         Manipulation of multivariate polynomials on a ring and Gröbner bases.
 homepage:            https://github.com/stla/hspray#readme
 license:             GPL-3
 license-file:        LICENSE
@@ -18,15 +18,17 @@   hs-source-dirs:      src
   exposed-modules:     Math.Algebra.Hspray
   build-depends:       base >= 4.7 && < 5
-                     , containers >= 0.6.4.1
-                     , hashable >= 1.3.4.0
-                     , unordered-containers >= 0.2.17.0
-                     , numeric-prelude >= 0.4.4
-                     , text >= 1.2.5.0
+                     , containers >= 0.6.4.1 && < 0.8
+                     , hashable >= 1.3.4.0 && < 1.5
+                     , unordered-containers >= 0.2.17.0 && < 0.3
+                     , numeric-prelude >= 0.4.4 && < 0.5
+                     , text >= 1.2.5.0 && < 2.2
   other-extensions:    FlexibleInstances
                      , FlexibleContexts
                      , MultiParamTypeClasses
                      , InstanceSigs
+                     , ScopedTypeVariables
+                     , BangPatterns
   default-language:    Haskell2010
   ghc-options:         -Wall
                        -Wcompat
@@ -42,6 +44,7 @@   type:                 exitcode-stdio-1.0
   main-is:              Main.hs
   hs-source-dirs:       tests/
+  other-modules:        Approx
   Build-Depends:        base >= 4.7 && < 5
                       , tasty
                       , tasty-hunit

@@ -1,11 +1,23 @@+{-|
+Module      : Math.Algebra.Hspray
+Description : Multivariate polynomials on a ring.
+Copyright   : (c) Stéphane Laurent, 2023
+License     : GPL-3
+Maintainer  : laurent_step@outlook.fr
+
+Deals with multivariate polynomials on a ring. See README for examples.
+-}
+{-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE InstanceSigs #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 
 module Math.Algebra.Hspray
   ( Powers (..)
   , Spray
+  , Monomial
   , lone
   , unitSpray
   , constantSpray
@@ -19,12 +31,23 @@   , (^*^)
   , (^**^)
   , evalSpray
+  , substituteSpray
+  , fromRationalSpray
   , prettySpray
+  , prettySpray'
+  , prettySprayXYZ
   , composeSpray
   , bombieriSpray
   , derivSpray
+  , leadingTerm
+  , sprayDivision
+  , groebner
+  , esPolynomial
+  , isSymmetricSpray
+  , isPolynomialOf
   ) where
 import qualified Algebra.Additive              as AlgAdd
+import qualified Algebra.Field                 as AlgField
 import qualified Algebra.Module                as AlgMod
 import qualified Algebra.Ring                  as AlgRing
 import qualified Data.Foldable                 as DF
@@ -32,12 +55,21 @@ import           Data.HashMap.Strict            ( HashMap )
 import qualified Data.HashMap.Strict           as HM
 import           Data.Hashable                  ( Hashable(hashWithSalt) )
-import           Data.List                      ( sortBy )
+import           Data.List                      ( sortBy
+                                                , maximumBy 
+                                                , (\\)
+                                                , findIndices
+                                                )
+import           Data.Maybe                     ( isJust
+                                                , fromJust 
+                                                )
+import           Data.Ord                       ( comparing )
 import qualified Data.Sequence                 as S
 import           Data.Sequence                  ( (><)
-                                                , Seq
+                                                , Seq (Empty, (:<|))
                                                 , dropWhileR
                                                 , (|>)
+                                                , (<|)
                                                 , index
                                                 , adjust
                                                 )
@@ -66,9 +98,11 @@   }
   deriving Show
 
+-- | append trailing zeros
 growSequence :: Seq Int -> Int -> Int -> Seq Int
 growSequence s m n = s >< t where t = S.replicate (n - m) 0
 
+-- | append trailing zeros to get the same length
 harmonize :: (Powers, Powers) -> (Powers, Powers)
 harmonize (pows1, pows2) = (Powers e1' n, Powers e2' n)
  where
@@ -91,6 +125,8 @@ 
 type Spray a = HashMap Powers a
 
+type Monomial a = (Powers, a)
+
 instance (AlgAdd.C a, Eq a) => AlgAdd.C (Spray a) where
   p + q = addSprays p q
   zero   = HM.empty
@@ -129,27 +165,34 @@   then AlgAdd.sum (replicate k pol)
   else AlgAdd.negate $ AlgAdd.sum (replicate (-k) pol)
 
+-- | drop trailing zeros
 simplifyPowers :: Powers -> Powers
 simplifyPowers pows = Powers s (S.length s)
   where s = dropWhileR (== 0) (exponents pows)
 
+-- | drop trailing zeros in the powers of a spray
 simplifySpray :: Spray a -> Spray a
 simplifySpray = HM.mapKeys simplifyPowers
 
+-- | simplify powers and remove zero terms
 cleanSpray :: (AlgAdd.C a, Eq a) => Spray a -> Spray a
 cleanSpray p = HM.filter (/= AlgAdd.zero) (simplifySpray p)
 
+-- | addition of two sprays
 addSprays :: (AlgAdd.C a, Eq a) => Spray a -> Spray a -> Spray a
 addSprays p q = cleanSpray $ HM.foldlWithKey' f p q
   where f s powers coef = HM.insertWith (AlgAdd.+) powers coef s
 
+-- | opposite spray
 negateSpray :: AlgAdd.C a => Spray a -> Spray a
 negateSpray = HM.map AlgAdd.negate
 
+-- | scale a spray by a scalar
 scaleSpray :: (AlgRing.C a, Eq a) => a -> Spray a -> Spray a
 scaleSpray lambda p = cleanSpray $ HM.map (lambda AlgRing.*) p
 
-derivMonomial :: AlgRing.C a => Int -> (Powers, a) -> (Powers, a) 
+-- | derivative of a monomial
+derivMonomial :: AlgRing.C a => Int -> Monomial a -> Monomial a 
 derivMonomial i (pows, coef) = if i' >= S.length expts 
   then (Powers S.empty 0, AlgAdd.zero)
   else (pows', coef')
@@ -161,20 +204,26 @@     coef'  = AlgAdd.sum (replicate expt_i coef)
     pows'  = Powers expts' (nvariables pows) 
 
-derivSpray :: (AlgRing.C a, Eq a) => Int -> Spray a -> Spray a
+-- | Derivative of a spray
+derivSpray 
+  :: (AlgRing.C a, Eq a) 
+  => Int     -- ^ index of the variable of differentiation
+  -> Spray a -- ^ the spray
+  -> Spray a
 derivSpray i p = cleanSpray $ HM.fromListWith (AlgAdd.+) monomials
  where
   p'        = HM.toList p
   monomials = [ derivMonomial i mp | mp <- p' ]
 
-
-multMonomial :: AlgRing.C a => (Powers, a) -> (Powers, a) -> (Powers, a)
+-- | multiply two monomials
+multMonomial :: AlgRing.C a => Monomial a -> Monomial a -> Monomial a
 multMonomial (pows1, coef1) (pows2, coef2) = (pows, coef1 AlgRing.* coef2)
  where
   (pows1', pows2') = harmonize (pows1, pows2)
   expts            = S.zipWith (+) (exponents pows1') (exponents pows2')
   pows             = Powers expts (nvariables pows1')
 
+-- | multiply two sprays
 multSprays :: (AlgRing.C a, Eq a) => Spray a -> Spray a -> Spray a
 multSprays p q = cleanSpray $ HM.fromListWith (AlgAdd.+) prods
  where
@@ -182,7 +231,7 @@   q'    = HM.toList q
   prods = [ multMonomial mp mq | mp <- p', mq <- q' ]
 
--- | Spray corresponding to polynomial x_n
+-- | Spray corresponding to the basic monomial x_n
 lone :: AlgRing.C a => Int -> Spray a
 lone n = HM.singleton pows AlgRing.one
  where
@@ -198,15 +247,55 @@ constantSpray :: (AlgRing.C a, Eq a) => a -> Spray a
 constantSpray c = c *^ lone 0
 
-evalMonomial :: AlgRing.C a => [a] -> (Powers, a) -> a
-evalMonomial xyz (powers, coeff) = coeff
-  AlgRing.* AlgRing.product (zipWith (AlgRing.^) xyz pows)
+-- | evaluates a monomial
+evalMonomial :: AlgRing.C a => [a] -> Monomial a -> a
+evalMonomial xyz (powers, coeff) = 
+  coeff AlgRing.* AlgRing.product (zipWith (AlgRing.^) xyz pows)
   where pows = DF.toList (fromIntegral <$> exponents powers)
 
 -- | Evaluate a spray
 evalSpray :: AlgRing.C a => Spray a -> [a] -> a
 evalSpray p xyz = AlgAdd.sum $ map (evalMonomial xyz) (HM.toList p)
 
+-- | spray from monomial
+fromMonomial :: Monomial a -> Spray a
+fromMonomial (pows, coeff) = HM.singleton pows coeff
+
+-- | number of variables in a spray
+numberOfVariables :: Spray a -> Int
+numberOfVariables spray = maximum (map nvariables powers)
+  where
+    powers = HM.keys spray
+
+-- | substitute some variables in a monomial
+substituteMonomial :: AlgRing.C a => [Maybe a] -> Monomial a -> Monomial a
+substituteMonomial subs (powers, coeff) = (powers'', coeff')
+  where
+    pows = exponents powers
+    n = nvariables powers
+    indices = findIndices isJust (take n subs)
+    pows' = [fromIntegral (pows `index` i) | i <- indices]
+    xyz = [fromJust (subs !! i) | i <- indices]
+    coeff' = coeff AlgRing.* AlgRing.product (zipWith (AlgRing.^) xyz pows')
+    f i a = if i `elem` indices then 0 else a
+    pows'' = S.mapWithIndex f pows
+    powers'' = simplifyPowers $ Powers pows'' n
+
+-- | Substitute some variables in a spray
+substituteSpray :: (Eq a, AlgRing.C a) => [Maybe a] -> Spray a -> Spray a
+substituteSpray subs spray = if length subs == n 
+  then spray'
+  else error "incorrect length of the substitutions list"
+  where
+    n = numberOfVariables spray
+    monomials = HM.toList spray
+    spray' = foldl1 (^+^) (map (fromMonomial . substituteMonomial subs) monomials)
+
+-- | Convert a spray with rational coefficients to a spray with double coefficients
+fromRationalSpray :: Spray Rational -> Spray Double
+fromRationalSpray = HM.map fromRational
+
+-- | helper for `composeSpray`
 identify :: (AlgRing.C a, Eq a) => Spray a -> Spray (Spray a)
 identify = HM.map constantSpray
 
@@ -219,6 +308,10 @@ fromList x = cleanSpray $ HM.fromList $ map
   (\(expts, coef) -> (Powers (S.fromList expts) (length expts), coef)) x
 
+
+-- pretty stuff ---------------------------------------------------------------
+
+-- | prettyPowers "x" [0, 2, 1] = x^(0, 2, 1)
 prettyPowers :: String -> [Int] -> Text
 prettyPowers var pows = append (pack x) (cons '(' $ snoc string ')')
  where
@@ -226,7 +319,11 @@   string = intercalate (pack ", ") (map (pack . show) pows)
 
 -- | Pretty form of a spray
-prettySpray :: (a -> String) -> String -> Spray a -> String
+prettySpray 
+  :: (a -> String) -- ^ function mapping a coefficient to a string, typically 'show'
+  -> String        -- ^ a string denoting the variable, e.g. "x"
+  -> Spray a       -- ^ the spray
+  -> String
 prettySpray prettyCoef var p = unpack $ intercalate (pack " + ") stringTerms
  where
   stringTerms = map stringTerm (sortBy (compare `on` fexpts) (HM.toList p))
@@ -238,6 +335,65 @@     pows       = DF.toList $ exponents (fst term)
     stringCoef = pack $ prettyCoef (snd term)
 
+-- | prettyPowers' [0, 2, 1] = "x2^2x3"
+prettyPowers' :: Seq Int -> Text
+prettyPowers' pows = pack x1x2x3
+ where
+  n = S.length pows
+  f i p 
+    | p == 0 = ""
+    | p == 1 = "x" ++ show i
+    | otherwise = "x" ++ show i ++ "^" ++ show p
+  x1x2x3 = concatMap (\i -> f i (pows `index` (i-1))) [1 .. n]
+
+-- | Pretty form of a spray, with monomials showed as "x1x3^2"
+prettySpray' :: (Show a) => Spray a -> String
+prettySpray' spray = unpack $ intercalate (pack " + ") terms
+ where
+  terms = map stringTerm (sortBy (compare `on` fexpts) (HM.toList spray))
+  fexpts term = exponents $ fst term
+  stringTerm term = append stringCoef'' (prettyPowers' pows)
+   where
+    pows       = exponents (fst term)
+    constant   = S.length pows == 0
+    stringCoef = pack $ show (snd term)
+    stringCoef' = cons '(' $ snoc stringCoef ')'
+    stringCoef'' = if constant then stringCoef' else snoc stringCoef' ' '
+
+-- | prettyPowersXYZ [1, 2, 1] = XY^2Z
+prettyPowersXYZ :: Seq Int -> Text
+prettyPowersXYZ pows = if n <= 3 
+  then pack xyz
+  else error "there is more than three variables"
+ where
+  n = S.length pows
+  gpows = growSequence pows n 3
+  f letter p 
+    | p == 0 = ""
+    | p == 1 = letter
+    | otherwise = letter ++ "^" ++ show p
+  x = f "X" (gpows `index` 0)
+  y = f "Y" (gpows `index` 1)
+  z = f "Z" (gpows `index` 2)
+  xyz = x ++ y ++ z
+
+-- | Pretty form of a spray having at more three variables
+prettySprayXYZ :: (Show a) => Spray a -> String
+prettySprayXYZ spray = unpack $ intercalate (pack " + ") terms
+ where
+  terms = map stringTerm (sortBy (compare `on` fexpts) (HM.toList spray))
+  fexpts term = exponents $ fst term
+  stringTerm term = append stringCoef'' (prettyPowersXYZ pows)
+   where
+    pows         = exponents (fst term)
+    constant     = S.length pows == 0
+    stringCoef   = pack $ show (snd term)
+    stringCoef'  = cons '(' $ snoc stringCoef ')'
+    stringCoef'' = if constant then stringCoef' else snoc stringCoef' ' '
+
+
+-- misc -----------------------------------------------------------------------
+
 -- | Terms of a spray
 sprayTerms :: Spray a -> HashMap (Seq Int) a
 sprayTerms = HM.mapKeys exponents
@@ -246,7 +402,7 @@ toList :: Spray a -> [([Int], a)]
 toList p = HM.toList $ HM.mapKeys (DF.toList . exponents) p
 
--- | Bombieri spray
+-- | Bombieri spray (for internal usage in the 'scubature' library)
 bombieriSpray :: AlgAdd.C a => Spray a -> Spray a
 bombieriSpray = HM.mapWithKey f
  where
@@ -254,3 +410,293 @@   pfactorial pows = product $ DF.toList $ factorial <$> S.filter (/= 0) pows
   factorial n     = product [1 .. n]
   times k x       = AlgAdd.sum (replicate k x)
+
+
+-- division stuff -------------------------------------------------------------
+
+-- | index of the maximum of a list
+maxIndex :: Ord a => [a] -> Int
+maxIndex = fst . maximumBy (comparing snd) . zip [0..]
+
+-- | Leading term of a spray 
+leadingTerm :: Spray a -> Monomial a
+leadingTerm p = (biggest, p HM.! biggest) 
+  where
+    powers = HM.keys p
+    i = maxIndex $ map exponents powers
+    biggest = powers !! i
+
+-- | whether a monomial divides another monomial
+divides :: Monomial a -> Monomial a -> Bool
+divides (powsP, _) (powsQ, _) = S.length expntsP <= S.length expntsQ && lower
+  where
+    expntsP = exponents powsP
+    expntsQ = exponents powsQ
+    lower = DF.all (\(x, y) -> x <= y) (S.zip expntsP expntsQ)
+
+-- | quotient of monomial Q by monomial p, assuming P divides Q
+quotient :: AlgField.C a => Monomial a -> Monomial a -> Monomial a
+quotient (powsQ, coeffQ) (powsP, coeffP) = (pows, coeff)
+  where
+    (powsP', powsQ') = harmonize (powsP, powsQ)
+    expntsP = exponents powsP'
+    expntsQ = exponents powsQ'
+    expnts = S.zipWith (-) expntsQ expntsP
+    n = nvariables powsP'
+    pows = Powers expnts n
+    coeff = coeffQ AlgField./ coeffP
+
+-- | Remainder of the division of a spray by a list of divisors, 
+-- using the lexicographic ordering of the monomials
+sprayDivision :: forall a. (Eq a, AlgField.C a) => Spray a -> [Spray a] -> Spray a
+sprayDivision p qs = 
+  if n == 0 then error "the list of divisors is empty" else snd $ ogo p AlgAdd.zero
+  where
+    n = length qs
+    qsltqs = zip qs (map leadingTerm qs)
+    g :: Monomial a -> Spray a -> Spray a -> (Spray a, Spray a)
+    g lts s r = (s ^-^ ltsspray, r ^+^ ltsspray)
+      where
+        ltsspray = fromMonomial lts 
+    go :: Monomial a -> Spray a -> Spray a -> Int -> Bool -> (Spray a, Spray a)
+    go lts !s r !i !divoccured
+      | divoccured = (s, r)
+      | i == n = g lts s r 
+      | otherwise = go lts news r (i+1) newdivoccured
+        where
+          (q, ltq) = qsltqs !! i
+          newdivoccured = divides ltq lts
+          news = if newdivoccured
+            then s ^-^ (fromMonomial (quotient lts ltq) ^*^ q)
+            else s
+    ogo :: Spray a -> Spray a -> (Spray a, Spray a)
+    ogo !s !r 
+      | s == AlgAdd.zero = (s, r)
+      | otherwise = ogo s' r'
+        where
+          (s', r') = go (leadingTerm s) s r 0 False
+
+
+-- Groebner stuff -------------------------------------------------------------
+
+-- | slight modification of `sprayDivision` to speed up groebner00
+sprayDivision' :: forall a. (Eq a, AlgField.C a) => Spray a -> HashMap Int (Spray a, Monomial a) -> Spray a
+sprayDivision' p qsltqs = snd $ ogo p AlgAdd.zero
+  where
+    n = HM.size qsltqs
+    g :: Monomial a -> Spray a -> Spray a -> (Spray a, Spray a)
+    g lts s r = (s ^-^ ltsspray, r ^+^ ltsspray)
+      where
+        ltsspray = fromMonomial lts 
+    go :: Monomial a -> Spray a -> Spray a -> Int -> Bool -> (Spray a, Spray a)
+    go lts !s r !i !divoccured
+      | divoccured = (s, r)
+      | i == n = g lts s r 
+      | otherwise = go lts news r (i+1) newdivoccured
+        where
+          (q, ltq) = qsltqs HM.! i
+          newdivoccured = divides ltq lts
+          news = if newdivoccured
+            then s ^-^ (fromMonomial (quotient lts ltq) ^*^ q)
+            else s
+    ogo :: Spray a -> Spray a -> (Spray a, Spray a)
+    ogo !s !r 
+      | s == AlgAdd.zero = (s, r)
+      | otherwise = ogo s' r'
+        where
+          (s', r') = go (leadingTerm s) s r 0 False
+
+combn2 :: Int -> Int -> HashMap Int (Int, Int)
+combn2 n s = HM.fromList (zip [0 .. n-2] (zip row1 row2)) 
+  where
+    row1 = drop s $ concatMap (\i -> [0 .. (i-1)]) [1 .. n-1]
+    row2 = drop s $ concatMap (\i -> replicate i i) [1 .. n-1]
+
+sPolynomial :: (Eq a, AlgField.C a) => (Spray a, Monomial a) -> (Spray a, Monomial a) -> Spray a
+sPolynomial pltp qltq = wp ^*^ p ^-^ wq ^*^ q
+  where
+    p = fst pltp
+    q = fst qltq
+    (lpowsP, lcoefP) = snd pltp
+    (lpowsQ, lcoefQ) = snd qltq
+    (lpowsP', lpowsQ') = harmonize (lpowsP, lpowsQ)
+    lexpntsP = exponents lpowsP'
+    lexpntsQ = exponents lpowsQ'
+    gamma = S.zipWith max lexpntsP lexpntsQ
+    betaP = S.zipWith (-) gamma lexpntsP
+    betaQ = S.zipWith (-) gamma lexpntsQ
+    n = nvariables lpowsP'
+    wp = fromMonomial (Powers betaP n, AlgField.recip lcoefP)
+    wq = fromMonomial (Powers betaQ n, AlgField.recip lcoefQ)
+
+-- | groebner basis, not minimal and not reduced
+groebner00 :: forall a. (Eq a, AlgField.C a) => [Spray a] -> [Spray a]
+groebner00 sprays = go 0 j0 combins0 spraysMap
+  where
+    j0 = length sprays
+    combins0 = combn2 j0 0
+    ltsprays = map leadingTerm sprays
+    spraysltsprays = zip sprays ltsprays 
+    spraysMap = HM.fromList (zip [0 .. j0-1] spraysltsprays)
+    go :: Int -> Int -> HashMap Int (Int, Int) -> HashMap Int (Spray a, Monomial a) -> [Spray a]
+    go !i !j !combins !gpolysMap
+      | i == length combins = map fst (HM.elems gpolysMap)
+      | otherwise = go i' j' combins' gpolysMap'
+        where
+          (k, l) = combins HM.! i
+          sfg = sPolynomial (gpolysMap HM.! k) (gpolysMap HM.! l)
+          sbarfg = sprayDivision' sfg gpolysMap
+          ltsbarfg = leadingTerm sbarfg
+          (i', j', gpolysMap', combins') = if sbarfg == AlgAdd.zero
+            then
+              (i+1, j, gpolysMap, combins)
+            else
+              ( 0
+              , j+1
+              , HM.insert j (sbarfg, ltsbarfg) gpolysMap
+              , combn2 (j+1) (i+1)
+              )
+
+-- | groebner basis, minimal but not reduced
+groebner0 :: forall a. (Eq a, AlgField.C a) => [Spray a] -> [Spray a]
+groebner0 sprays = 
+  if n <= 1 then sprays else [basis00 !! k | k <- [0 .. n-1] \\ discard]
+  where
+    n = length basis00
+    basis00 = groebner00 sprays
+    go :: Int -> [Int] -> [Int]
+    go !i toRemove
+      | i == n = toRemove
+      | otherwise = go (i+1) toRemove'
+        where
+          ltf = leadingTerm (basis00 !! i)
+          toDrop = toRemove ++ [i]
+          igo :: Int -> Bool
+          igo !j 
+            | j == n = False
+            | j `elem` toDrop = igo (j+1)
+            | otherwise = ok || igo (j+1)
+              where 
+                ok = divides (leadingTerm (basis00 !! j)) ltf
+          toRemove' = if igo 0 then toDrop else toRemove
+    discard = go 0 []
+
+-- | reduce a Groebner basis
+reduceGroebnerBasis :: forall a. (Eq a, AlgField.C a) => [Spray a] -> [Spray a]
+reduceGroebnerBasis gbasis = 
+  if length gbasis >= 2 then map reduction [0 .. n-1] else ngbasis
+  where
+    normalize :: Spray a -> Spray a
+    normalize spray = AlgField.recip coef *^ spray
+      where
+        (_, coef) = leadingTerm spray
+    ngbasis = map normalize gbasis
+    n = length ngbasis
+    reduction :: Int -> Spray a
+    reduction i = sprayDivision (ngbasis !! i) rest
+      where
+        rest = [ngbasis !! k | k <- [0 .. n-1] \\ [i]]
+
+-- | Groebner basis (always minimal and possibly reduced)
+groebner 
+  :: forall a. (Eq a, AlgField.C a) 
+  => [Spray a] -- ^ list of sprays 
+  -> Bool      -- ^ whether to return the reduced basis
+  -> [Spray a]
+groebner sprays reduced = 
+  if reduced then reduceGroebnerBasis gbasis0 else map normalize gbasis0
+  where
+    gbasis0 = groebner0 sprays
+    normalize :: Spray a -> Spray a
+    normalize spray = AlgField.recip coef *^ spray
+      where
+        (_, coef) = leadingTerm spray
+
+
+-- elementary symmetric polynomials -------------------------------------------
+
+-- | generate all permutations of a binary sequence
+permutationsBinarySequence :: Int -> Int -> [Seq Int] 
+permutationsBinarySequence nzeros nones = unfold1 next z 
+  where
+    z = (><) (S.replicate nzeros False) (S.replicate nones True)
+    unfold1 :: (Seq Bool -> Maybe (Seq Bool)) -> Seq Bool -> [Seq Int]
+    unfold1 f x = case f x of 
+      Nothing -> [x'] 
+      Just y  -> x' : unfold1 f y 
+      where
+        x' = fmap fromEnum x
+    next :: Seq Bool -> Maybe (Seq Bool)
+    next xs = case findj (S.reverse xs, S.empty) of 
+      Nothing -> Nothing
+      Just ( l:<|ls , rs ) -> Just $ inc l ls (S.reverse rs, S.empty) 
+      Just ( Empty , _ ) -> error "permutationsBinarySequence: should not happen"
+    findj :: (Seq Bool, Seq Bool) -> Maybe (Seq Bool, Seq Bool)   
+    findj ( xxs@(x:<|xs), yys@(_:<|_) ) = if x 
+      then findj ( xs, True <| yys )
+      else Just ( xxs, yys )
+    findj ( x:<|xs, Empty ) = findj ( xs, S.singleton x )
+    findj ( Empty , _ ) = Nothing
+    inc :: Bool -> Seq Bool -> (Seq Bool, Seq Bool) -> Seq Bool
+    inc !u us ( x:<|xs , yys ) = if u
+      then inc True us ( xs , x <| yys ) 
+      else (><) (S.reverse (True <| us)) ((><) (S.reverse (u <| yys)) xs)
+    inc _ _ ( Empty , _ ) = error "permutationsBinarySequence: should not happen"
+
+-- | Elementary symmetric polynomial
+esPolynomial 
+  :: (AlgRing.C a, Eq a) 
+  => Int -- ^ number of variables
+  -> Int -- ^ index
+  -> Spray a
+esPolynomial n k
+  | k <= 0 || n <= 0 = error "both arguments must be positive integers"
+  | k > n = AlgAdd.zero
+  | otherwise = simplifySpray spray
+  where
+    perms = permutationsBinarySequence (n-k) k
+    spray = HM.fromList $ map (\expts -> (Powers expts n, AlgRing.one)) perms
+
+-- | Whether a spray is a symmetric polynomial
+isSymmetricSpray :: forall a. (AlgField.C a, Eq a) => Spray a -> Bool
+isSymmetricSpray spray = check1 && check2 
+  where
+    n = numberOfVariables spray
+    indices = [1 .. n]
+    esPolys = map (\i -> esPolynomial n i :: Spray a) indices
+    yPolys = map (\i -> lone (n + i) :: Spray a) indices
+    gPolys = zipWith (^-^) esPolys yPolys
+    gbasis = groebner0 gPolys
+    g = sprayDivision spray gbasis
+    gpowers = HM.keys g
+    check1 = minimum (map nvariables gpowers) > n
+    expnts = map exponents gpowers
+    check2 = DF.all (DF.all (0 ==)) (map (S.take n) expnts) 
+
+-- | Whether a spray can be written as a polynomial of a given list of sprays;
+-- the sprays in the list must belong to the same polynomial ring as the spray
+isPolynomialOf :: forall a. (AlgField.C a, Eq a) => Spray a -> [Spray a] -> (Bool, Maybe (Spray a))
+isPolynomialOf spray sprays = result 
+  where
+    n = numberOfVariables spray
+    n' = maximum $ map numberOfVariables sprays
+    result
+      | n > n' = (False, Nothing)
+      | n < n' = error "not enough variables in the spray" 
+      | otherwise = (checks, poly)
+        where
+          m = length sprays
+          yPolys = map (\i -> lone (n + i) :: Spray a) [1 .. m]
+          gPolys = zipWith (^-^) sprays yPolys
+          gbasis0 = groebner0 gPolys
+          g = sprayDivision spray gbasis0
+          gpowers = HM.keys g
+          check1 = minimum (map nvariables gpowers) > n
+          expnts = map exponents gpowers
+          check2 = DF.all (DF.all (0 ==)) (map (S.take n) expnts)
+          checks = check1 && check2
+          poly = if checks
+            then Just $ dropXis g
+            else Nothing
+          dropXis = HM.mapKeys f
+          f (Powers expnnts _) = Powers (S.drop n expnnts) n

@@ -0,0 +1,10 @@+module Approx (assertApproxEqual) where
+import           Test.Tasty.HUnit ( Assertion, assertEqual )
+
+-- round x to n digits
+approx :: Int -> Double -> Double
+approx n x = fromInteger (round $ x * (10^n)) / (10.0^^n)
+
+assertApproxEqual :: String -> Int -> Double -> Double -> Assertion
+assertApproxEqual prefix n x1 x2 = 
+  assertEqual prefix (approx n x1) (approx n x2)

@@ -1,22 +1,32 @@ module Main where
+import           Approx                         ( assertApproxEqual )
 import           Data.Ratio                     ( (%) )
 import           Math.Algebra.Hspray            ( Spray,
                                                   (^+^),
+                                                  (^-^),
                                                   (^*^),
                                                   (^**^),
                                                   (*^),
                                                   lone,
                                                   unitSpray,
                                                   evalSpray,
+                                                  substituteSpray,
                                                   composeSpray,
                                                   fromList,
                                                   toList,
                                                   bombieriSpray,
-                                                  derivSpray )
+                                                  derivSpray,
+                                                  groebner,
+                                                  fromRationalSpray,
+                                                  esPolynomial,
+                                                  isSymmetricSpray, 
+                                                  isPolynomialOf 
+                                                )
 import           Test.Tasty                     ( defaultMain
                                                 , testGroup
                                                 )
 import           Test.Tasty.HUnit               ( assertEqual
+                                                , assertBool
                                                 , testCase
                                                 )
 
@@ -68,6 +78,44 @@         p1' = derivSpray 1 p1
         p2' = derivSpray 1 p2
         q'  = derivSpray 1 q
-      assertEqual "" q' ((p1' ^*^ p2) ^+^ (p1 ^*^ p2'))
+      assertEqual "" q' ((p1' ^*^ p2) ^+^ (p1 ^*^ p2')),
 
+    testCase "groebner" $ do
+      let
+        x = lone 1 :: Spray Rational
+        y = lone 2 :: Spray Rational
+        z = lone 3 :: Spray Rational
+        p1 = x^**^2 ^+^ y ^+^ z ^-^ unitSpray
+        p2 = x ^+^ y^**^2 ^+^ z ^-^ unitSpray
+        p3 = x ^+^ y ^+^ z^**^2 ^-^ unitSpray
+        g = groebner [p1, p2, p3] True
+        xyz = [sqrt 2 - 1, sqrt 2 - 1, sqrt 2 - 1]
+        gxyz = map ((`evalSpray` xyz) . fromRationalSpray) g
+        sumAbsValues = sum $ map abs gxyz
+      assertApproxEqual "" 8 sumAbsValues 0,
+
+    testCase "symmetric polynomials" $ do
+      let
+        e2 = esPolynomial 4 2 :: Spray Rational
+        e3 = esPolynomial 4 3 :: Spray Rational
+        p = e2^**^2 ^+^ (2*^ e3)
+      assertBool "" (isSymmetricSpray p),
+
+    testCase "isPolynomialOf" $ do
+      let
+        x = lone 1 :: Spray Rational
+        y = lone 2 :: Spray Rational
+        p1 = x ^+^ y
+        p2 = x ^-^ y
+        p = p1 ^*^ p2
+      assertEqual "" (isPolynomialOf p [p1, p2]) (True, Just $ x ^*^ y),
+
+    testCase "substituteSpray" $ do
+      let
+        x1 = lone 1 :: Spray Rational
+        x2 = lone 2 :: Spray Rational
+        x3 = lone 3 :: Spray Rational
+        p = x1^**^2 ^+^ x2 ^+^ x3 ^-^ unitSpray
+        p' = substituteSpray [Just 2, Nothing, Just 3] p
+      assertEqual "" p' (x2 ^+^ (6*^ unitSpray))
   ]