hspray 0.2.4.0 → 0.2.5.0

raw patch · 5 files changed

+651/−124 lines, 5 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• CHANGELOG.md +18/−1
• README.md +58/−0
• hspray.cabal +3/−1
• src/Math/Algebra/Hspray.hs +493/−116
• tests/Main.hs +79/−6

@@ -88,4 +88,21 @@ 
 * The function `sprayDivision` has been renamed to `sprayDivisionRemainder`.
 
-* New function `sprayDivision`, returning the quotient and the remainder of the division of two sprays.+* New function `sprayDivision`, returning the quotient and the remainder of the division of two sprays.
+
+
+## 0.2.5.0 - 2024-04-02
+
+* New function `resultant'` which computes the resultant of two sprays with coefficients in a field. 
+Thus it is less general than the function `resultant` but it is more efficient. 
+
+* Fixed a small mistake in `isSymmetricSpray` and `isPolynomialOf`: these functions didn't deal with 
+the constant term of the spray.
+
+* New function `psPolynomial` which computes the power sum polynomials.
+
+* A particular type of sprays, namely `SymbolicSpray a`, has been introduced. The coefficients of these 
+sprays are ratios of univariate polynomials with `a` coefficients. There is a specialization 
+`SymbolicQSpray` for the case when `a` is a type of rational numbers. The necessary instances have 
+been defined and there is the function `prettySymbolic(Q)Spray` to display such sprays. There are also 
+some functions to perform evaluation of such sprays.

@@ -175,6 +175,64 @@ -- [[([0,1],2 % 3)],[([1],1 % 1)],[([1],1 % 1)]]
 ```
 
+## The `SymbolicSpray` type
+
+If you have only one symbolic coefficient, it is easier to deal with the sprays of type 
+`SymbolicSpray`. These are sprays whose coefficients are ratios of univariate polynomials, 
+so this allows more possibilities than a `Spray (Spray a)`. 
+Assume you want to deal with the polynomial `4/5 * a/(a² + a + 1) * (x² + y²) + 2a/3 * yz`. 
+Then you define it as follows:
+
+```haskell
+import           Prelude hiding ((*), (+), (-), (/), (^), (*>))
+import qualified Prelude as P
+import           Algebra.Additive              
+import           Algebra.Module            
+import           Algebra.Ring
+import           Algebra.Field                
+import           Math.Algebra.Hspray
+import           Number.Ratio       ( (%), T ( (:%) ) )
+x = lone 1 :: SymbolicQSpray 
+y = lone 2 :: SymbolicQSpray 
+z = lone 3 :: SymbolicQSpray 
+a = outerQVariable  
+sSpray 
+  = ((4%5) *. (a :% (a^2 + a + one))) *> (x^2 + y^2)  +  (constQPoly (2%3) * a) *> (y * z)
+putStrLn $ prettySymbolicQSpray "a" sSpray
+-- ([(4/5)a] / [(1) + a + a^2])*x1^2 + ([(4/5)a] / [(1) + a + a^2])*x2^2 + ((2/3)a)*x2x3
+```
+
+This pretty form of the symbolic qspray will be improved in a future version.
+
+There are three possible evaluations of this symbolic spray:
+
+```haskell
+-- substitute a value for 'a':
+putStrLn $ prettySpray' $ evalSymbolicSpray sSpray (6%5)
+-- (24 % 91) x1^2 + (24 % 91) x2^2 + (4 % 5) x2x3
+
+-- substitue a value for 'a' and some values for 'x1', 'x2', 'x3':
+evalSymbolicSpray' sSpray (6%5) [2, 3, 4%7]
+-- 24 % 5
+
+-- substitue some values for 'x1', 'x2', 'x3':
+putStrLn $ 
+  prettyRatioOfQPolynomials "a" $ evalSymbolicSpray'' sSpray [2, 3, 4%7]
+-- [(404/35)a + (8/7)a^2 + (8/7)a^3] / [(1) + a + a^2]
+```
+
+The nice point regarding these ratios of univariate polynomials is that they are automatically 
+"simplified". For example:
+
+```haskell
+polyFrac = (a^8 - one) :% (a - one)
+putStrLn $ prettyRatioOfQPolynomials "a" polyFrac
+-- (1) + a + a^2 + a^3 + a^4 + a^5 + a^6 + a^7
+```
+
+Maybe you prefer the fractional form, but it is nice to see that this ratio of 
+polynomials actually is a polynomial.
+
 
 ## Other features
 

@@ -1,5 +1,5 @@ name:                hspray
-version:             0.2.4.0
+version:             0.2.5.0
 synopsis:            Multivariate polynomials.
 description:         Manipulation of multivariate polynomials on a commutative ring, Gröbner basis, resultant and subresultants, and greatest common divisor.
 homepage:            https://github.com/stla/hspray#readme
@@ -29,6 +29,7 @@                      , MultiParamTypeClasses
                      , InstanceSigs
                      , ScopedTypeVariables
+                     , GeneralizedNewtypeDeriving
                      , BangPatterns
   default-language:    Haskell2010
   ghc-options:         -Wall
@@ -49,6 +50,7 @@   Build-Depends:        base >= 4.7 && < 5
                       , tasty >= 1.4 && < 1.6
                       , tasty-hunit >= 0.10 && < 0.11
+                      , numeric-prelude >= 0.4.4 && < 0.5
                       , hspray
   Default-Language:     Haskell2010
 

@@ -5,7 +5,8 @@ License     : GPL-3
 Maintainer  : laurent_step@outlook.fr
 
-Deals with multivariate polynomials on a commutative ring. See README for examples.
+Deals with multivariate polynomials on a commutative ring. 
+See README for examples.
 -}
 {-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE FlexibleInstances #-}
@@ -13,12 +14,14 @@ {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE InstanceSigs #-}
 {-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
 
 module Math.Algebra.Hspray
   ( 
   -- * Types
     Powers (..)
   , Spray
+  , QSpray
   , Monomial
   -- * Basic sprays
   , lone
@@ -35,9 +38,39 @@   -- * Showing a spray
   , prettySpray
   , prettySpray'
+  , prettySpray''
   , prettySprayXYZ
+  -- * Univariate polynomials 
+  , A (..)
+  , Rational'
+  , Q
+  , scalarQ
+  , Polynomial 
+  , RatioOfPolynomials
+  , QPolynomial 
+  , RatioOfQPolynomials
+  , prettyRatioOfPolynomials
+  , prettyRatioOfQPolynomials
+  , (*.)
+  , constPoly
+  , polyFromCoeffs
+  , outerVariable
+  , constQPoly
+  , qpolyFromCoeffs
+  , outerQVariable
+  , evalRatioOfPolynomials
+  -- * Symbolic sprays (with univariate polynomials coefficients)
+  , SymbolicSpray
+  , SymbolicQSpray
+  , prettySymbolicSpray
+  , prettySymbolicQSpray
+  , simplifySymbolicSpray
+  , evalSymbolicSpray
+  , evalSymbolicSpray'
+  , evalSymbolicSpray''
   -- * Queries on a spray
   , getCoefficient
+  , getConstantTerm
   , sprayTerms
   -- * Evaluation of a spray
   , evalSpray
@@ -56,9 +89,11 @@   , reduceGroebnerBasis
   -- * Symmetric polynomials
   , esPolynomial
+  , psPolynomial
   , isSymmetricSpray
   -- * Resultant and subresultants
   , resultant
+  , resultant'
   , resultant1
   , subresultants
   , subresultants1
@@ -76,6 +111,7 @@ import qualified Algebra.Field                 as AlgField
 import qualified Algebra.Module                as AlgMod
 import qualified Algebra.Ring                  as AlgRing
+import qualified Algebra.ZeroTestable          as AlgZT
 import qualified Data.Foldable                 as DF
 import           Data.Function                  ( on )
 import           Data.HashMap.Strict            ( HashMap )
@@ -118,23 +154,231 @@                                                 , snoc
                                                 , unpack
                                                 )
+import qualified MathObj.Polynomial            as MathPol
+import qualified Number.Ratio                  as NumberRatio
+-- import qualified Algebra.PrincipalIdealDomain  as AlgPID
+-- import qualified Algebra.Units  as AlgUnits
+-- import qualified Algebra.IntegralDomain  as AlgID
 
 
-infixr 7 *^, .^
+-- Univariate polynomials -----------------------------------------------------
 
-infixl 6 ^+^, ^-^
+newtype A a = A a 
+  deriving
+    (Eq, AlgAdd.C, AlgRing.C, AlgField.C)
 
-infixl 7 ^*^
+type Rational' = NumberRatio.Rational
+type Q = A Rational'
 
-infixr 8 ^**^
+-- | Identify a rational to a @A Rational'@ element
+scalarQ :: Rational' -> Q
+scalarQ = A 
 
+type Polynomial a         = MathPol.T (A a)
+type RatioOfPolynomials a = NumberRatio.T (Polynomial a)
+type QPolynomial          = Polynomial Rational'
+type RatioOfQPolynomials  = RatioOfPolynomials Rational'
 
+instance (Eq a, AlgField.C a) => AlgZT.C (A a) where
+  isZero :: A a -> Bool
+  isZero (A r) = r == AlgAdd.zero
+
+instance (Eq a, AlgField.C a) => AlgMod.C (A a) (RatioOfPolynomials a) where
+  (*>) :: A a -> RatioOfPolynomials a -> RatioOfPolynomials a
+  r *> rop = NumberRatio.scale (MathPol.const r) rop 
+
+instance (Eq a, AlgField.C a) => AlgMod.C (Polynomial a) (RatioOfPolynomials a) where
+  (*>) :: Polynomial a -> RatioOfPolynomials a -> RatioOfPolynomials a
+  p *> r = NumberRatio.scale p r 
+
+instance (Eq a, AlgField.C a) => AlgMod.C (Polynomial a) (SymbolicSpray a) where
+  (*>) :: Polynomial a -> SymbolicSpray a -> SymbolicSpray a
+  p *> r = constantSpray (p NumberRatio.:% AlgRing.one) ^*^ r
+
+infixr 7 *.
+-- | Scale a ratio of univariate polynomials by a scalar
+(*.) :: (Eq a, AlgField.C a) => a -> RatioOfPolynomials a -> RatioOfPolynomials a
+(*.) scalar rop = A scalar AlgMod.*> rop
+
+-- | Constant univariate polynomial
+constPoly :: a -> Polynomial a
+constPoly x = MathPol.const (A x)
+
+-- | Univariate polynomial from its coefficients (ordered by increasing degrees)
+polyFromCoeffs :: [a] -> Polynomial a
+polyFromCoeffs as = MathPol.fromCoeffs (map A as)
+
+-- | The variable of a univariate polynomial; it is called \"outer\" because this is the variable 
+-- occuring in the polynomial coefficients of a `SymbolicSpray` 
+outerVariable :: AlgRing.C a => Polynomial a
+outerVariable = polyFromCoeffs [AlgAdd.zero, AlgRing.one] 
+
+-- | Constant rational univariate polynomial
+-- 
+-- >>> import Number.Ratio ( (%) )
+-- >>> constQPoly (2 % 3)
+constQPoly :: Rational' -> QPolynomial
+constQPoly = constPoly
+
+-- | Rational univariate polynomial from coefficients
+-- 
+-- >>> import Number.Ratio ( (%) )
+-- >>> qpolyFromCoeffs [2 % 3, 5, 7 % 4]
+qpolyFromCoeffs :: [Rational'] -> QPolynomial
+qpolyFromCoeffs = polyFromCoeffs
+
+-- | The variable of a univariate qpolynomial; it is called \"outer\" because this is the variable 
+-- occuring in the polynomial coefficients of a `SymbolicQSpray` 
+--
+-- prop> outerQVariable == qpolyFromCoeffs [0, 1] 
+outerQVariable :: QPolynomial
+outerQVariable = qpolyFromCoeffs [0, 1] 
+
+-- show a ratio
+showQ :: (Eq a, Num a, Show a) => NumberRatio.T a -> String
+showQ q = if d == 1 
+  then show n 
+  else show n ++ "/" ++ show d
+  where
+    n = NumberRatio.numerator q
+    d = NumberRatio.denominator q 
+
+-- | helper function for prettyRatioOfPolynomials (and prettySymbolicSpray)
+showQpol :: forall a. (Eq a, AlgField.C a) 
+         => Polynomial a -> String -> (a -> String) -> Bool -> String
+showQpol pol variable showCoeff brackets = if brackets 
+  then '[' : polyString ++ "]"
+  else polyString
+  where
+    showCoeff' :: Int -> A a -> String
+    showCoeff' i (A coeff) = case i of 
+      0 -> '(' : showCoeff coeff ++ ")"
+      _ -> if coeff == AlgRing.one 
+        then "" 
+        else '(' : showCoeff coeff ++ ")"
+    coeffs   = MathPol.coeffs pol
+    nonzeros = findIndices (/= A AlgAdd.zero) coeffs
+    terms    = map (pack . showTerm) nonzeros
+      where
+        showTerm i = case i of 
+          0 -> showCoeff' 0 (coeffs !! 0)
+          1 -> showCoeff' 1 (coeffs !! 1) ++ variable
+          _ -> showCoeff' i (coeffs !! i) ++ variable ++ "^" ++ show i
+    polyString = unpack (intercalate (pack " + ") terms)
+
+-- | helper function for prettyRatioOfPolynomials (and prettySymbolicSpray)
+showQpolysRatio :: forall a. (Eq a, AlgField.C a) 
+                   => String -> (a -> String) -> RatioOfPolynomials a -> String
+showQpolysRatio var showCoeff polysRatio = numeratorString ++ denominatorString
+  where
+    denominator       = NumberRatio.denominator polysRatio
+    brackets          = denominator /= MathPol.const (A AlgRing.one)
+    numeratorString   = showQpol (NumberRatio.numerator polysRatio) var showCoeff brackets
+    denominatorString = if not brackets
+      then ""
+      else " / " ++ showQpol denominator var showCoeff True
+
+-- | Pretty form of a ratio of univariate polynomials
+prettyRatioOfPolynomials :: (Eq a, AlgField.C a, Show a) 
+  => String               -- ^ a string to denote the variable, e.g. @"a"@
+  -> RatioOfPolynomials a 
+  -> String 
+prettyRatioOfPolynomials var = showQpolysRatio var show 
+
+-- | Pretty form of a ratio of univariate qpolynomials
+prettyRatioOfQPolynomials 
+  :: String               -- ^ a string to denote the variable, e.g. @"a"@ 
+  -> RatioOfQPolynomials 
+  -> String 
+prettyRatioOfQPolynomials var = showQpolysRatio var showQ
+
+-- | Evaluates a ratio of univariate polynomials
+evalRatioOfPolynomials :: AlgField.C a 
+  => a                    -- ^ the value at which the evaluation is desired
+  -> RatioOfPolynomials a 
+  -> a
+evalRatioOfPolynomials value polysRatio = 
+  resultNumerator AlgField./ resultDenominator
+  where
+    A resultNumerator   = MathPol.evaluate (NumberRatio.numerator polysRatio) (A value)
+    A resultDenominator = MathPol.evaluate (NumberRatio.denominator polysRatio) (A value)
+
+
+-- Symbolic sprays ------------------------------------------------------------
+
+type SymbolicSpray a = Spray (RatioOfPolynomials a)
+type SymbolicQSpray  = SymbolicSpray Rational'
+
+-- | Simplifies the coefficients (the ratio of univariate polynomials) of a 
+-- symbolic spray
+simplifySymbolicSpray :: 
+  (Eq a, AlgField.C a) => SymbolicSpray a -> SymbolicSpray a
+simplifySymbolicSpray = HM.map (AlgAdd.+ AlgAdd.zero)
+
+-- | Pretty form of a symbolic spray
+prettySymbolicSpray 
+  :: (Eq a, Show a, AlgField.C a) 
+  => String          -- ^ a string to denote the outer variable of the spray, e.g. @"a"@
+  -> SymbolicSpray a -- ^ a symbolic spray; note that this function does not simplify it
+  -> String 
+prettySymbolicSpray var = prettySpray'' (showQpolysRatio var show)
+
+-- | Pretty form of a symbolic qspray
+prettySymbolicQSpray 
+  :: String          -- ^ a string to denote the outer variable of the spray, e.g. @"a"@
+  -> SymbolicQSpray  -- ^ a symbolic qspray; note that this function does not simplify it
+  -> String 
+prettySymbolicQSpray var = prettySpray'' (showQpolysRatio var showQ)
+
+-- | Substitutes a value to the outer variable of a symbolic spray
+evalSymbolicSpray :: AlgField.C a => SymbolicSpray a -> a -> Spray a
+evalSymbolicSpray spray x = HM.map (evalRatioOfPolynomials x) spray 
+
+-- | Substitutes a value to the outer variable of a symbolic spray as well 
+-- as some values to the inner variables of this spray
+evalSymbolicSpray' :: AlgField.C a 
+  => SymbolicSpray a -- ^ symbolic spray to be evaluated
+  -> a               -- ^ a value for the outer variable
+  -> [a]             -- ^ some values for the inner variables 
+  -> a
+evalSymbolicSpray' spray x xs = if length xs >= numberOfVariables spray 
+  then evalSpray (evalSymbolicSpray spray x) xs
+  else error "evalSymbolicSpray': not enough values provided."
+
+-- helper function for evalSymbolicSpray''
+evalSymbolicMonomial :: (Eq a, AlgField.C a) 
+  => [a] -> Monomial (RatioOfPolynomials a) -> RatioOfPolynomials a
+evalSymbolicMonomial xs (powers, coeff) = 
+  AlgRing.product (zipWith (AlgRing.^) xs pows) *. coeff
+  where 
+    pows = DF.toList (fromIntegral <$> exponents powers)
+
+-- | Substitutes some values to the inner variables of a symbolic spray
+evalSymbolicSpray'' 
+  :: (Eq a, AlgField.C a) => SymbolicSpray a -> [a] -> RatioOfPolynomials a
+evalSymbolicSpray'' spray xs = if length xs >= numberOfVariables spray
+  then AlgAdd.sum $ map (evalSymbolicMonomial xs) (HM.toList spray)
+  else error "evalSymbolicSpray'': not enough values provided."
+
+
+-- Sprays ---------------------------------------------------------------------
+
 data Powers = Powers
   { exponents  :: Seq Int
   , nvariables :: Int
   }
   deriving Show
 
+instance Eq Powers where
+  (==) :: Powers -> Powers -> Bool
+  pows1 == pows2 = exponents pows1' == exponents pows2'
+    where 
+      (pows1', pows2') = harmonize (pows1, pows2)
+
+instance Hashable Powers where
+  hashWithSalt :: Int -> Powers -> Int
+  hashWithSalt k pows = hashWithSalt k (exponents pows, nvariables pows)
+
 -- | append trailing zeros
 growSequence :: Seq Int -> Int -> Int -> Seq Int
 growSequence s m n = s >< t where t = S.replicate (n - m) 0
@@ -154,20 +398,46 @@     then (growSequence e1 n1 n2, e2, n2)
     else (e1, growSequence e2 n2 n1, n1)
 
-instance Eq Powers where
-  (==) :: Powers -> Powers -> Bool
-  pows1 == pows2 = exponents pows1' == exponents pows2'
-    where 
-      (pows1', pows2') = harmonize (pows1, pows2)
-
-instance Hashable Powers where
-  hashWithSalt :: Int -> Powers -> Int
-  hashWithSalt k pows = hashWithSalt k (exponents pows, nvariables pows)
+-- | drop trailing zeros
+simplifyPowers :: Powers -> Powers
+simplifyPowers pows = Powers s (S.length s)
+  where 
+    s = dropWhileR (== 0) (exponents pows)
 
+type Monomial a = (Powers, a)
 type Spray a = HashMap Powers a
+type QSpray = Spray Rational'
 
-type Monomial a = (Powers, a)
+-- | 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
+
+-- | 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
+  p'    = HM.toList p
+  q'    = HM.toList q
+  prods = [ multMonomial mp mq | mp <- p', mq <- q' ]
+
 instance (AlgAdd.C a, Eq a) => AlgAdd.C (Spray a) where
   (+) :: Spray a -> Spray a -> Spray a
   p + q  = addSprays p q
@@ -195,29 +465,35 @@   signum _ = error "Prelude.Num.signum: inappropriate abstraction"
  -} 
 
+infixl 6 ^+^
 -- | Addition of two sprays
 (^+^) :: (AlgAdd.C a, Eq a) => Spray a -> Spray a -> Spray a
 (^+^) p q = p AlgAdd.+ q
 
+infixl 6 ^-^
 -- | Substraction of two sprays
 (^-^) :: (AlgAdd.C a, Eq a) => Spray a -> Spray a -> Spray a
 (^-^) p q = p AlgAdd.- q
 
+infixl 7 ^*^
 -- | Multiply two sprays
 (^*^) :: (AlgRing.C a, Eq a) => Spray a -> Spray a -> Spray a
 (^*^) p q = p AlgRing.* q
 
+infixr 8 ^**^
 -- | Power of a spray
 (^**^) :: (AlgRing.C a, Eq a) => Spray a -> Int -> Spray a
 (^**^) p n = if n >= 0 
   then AlgRing.product (replicate n p)
   else error "(^**^): negative power of a spray is not allowed."
 
--- | Scale spray by a scalar
+infixr 7 *^
+-- | Scale a spray by a scalar
 (*^) :: (AlgRing.C a, Eq a) => a -> Spray a -> Spray a
 (*^) lambda pol = lambda AlgMod.*> pol
 
--- | Scale spray by an integer
+infixr 7 .^
+-- | Scale a spray by an integer
 --
 -- prop> 3 .^ p == p ^+^ p ^+^ p
 (.^) :: (AlgAdd.C a, Eq a) => Int -> Spray a -> Spray a
@@ -225,12 +501,6 @@   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
@@ -239,20 +509,6 @@ 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
-
 -- | derivative of a monomial
 derivMonomial :: AlgRing.C a => Int -> Monomial a -> Monomial a 
 derivMonomial i (pows, coef) = if i' >= S.length expts 
@@ -279,22 +535,6 @@   p'        = HM.toList p
   monomials = [ derivMonomial i mp | mp <- p' ]
 
--- | 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
-  p'    = HM.toList p
-  q'    = HM.toList q
-  prods = [ multMonomial mp mq | mp <- p', mq <- q' ]
-
 -- | Spray corresponding to the basic monomial x_n
 --
 -- >>> x :: lone 1 :: Spray Int
@@ -347,6 +587,14 @@     expnts' = S.dropWhileR (== 0) (S.fromList expnts)
     powers  = Powers expnts' (S.length expnts')
 
+-- | Get the constant term of a spray
+--
+-- prop> getConstantTerm p == getCoefficient [] p 
+getConstantTerm :: AlgAdd.C a => Spray a -> a
+getConstantTerm spray = fromMaybe AlgAdd.zero (HM.lookup powers spray)
+  where
+    powers  = Powers S.empty 0
+
 -- | number of variables in a spray
 numberOfVariables :: Spray a -> Int
 numberOfVariables spray =
@@ -391,7 +639,7 @@     pows''   = S.mapWithIndex f pows
     powers'' = simplifyPowers $ Powers pows'' n
 
--- | Substitutes some variables in a spray
+-- | Substitutes some variables in a spray by some values
 --
 -- >>> x1 :: lone 1 :: Spray Int
 -- >>> x2 :: lone 2 :: Spray Int
@@ -410,12 +658,13 @@     spray'    = 
       foldl1' (^+^) (map (fromMonomial . substituteMonomial subs) monomials)
 
--- | Converts a spray with rational coefficients to a spray with double coefficients
--- (useful for evaluation)
+-- | Converts a spray with rational coefficients to a spray with double 
+-- coefficients (useful for evaluation)
 fromRationalSpray :: Spray Rational -> Spray Double
 fromRationalSpray = HM.map fromRational
 
--- | Composes a spray with a change of variables
+-- | Sustitutes the variables of a spray with some sprays 
+-- (e.g. change of variables)
 --
 -- >>> x :: lone 1 :: Spray Int
 -- >>> y :: lone 2 :: Spray Int
@@ -424,7 +673,8 @@ -- >>> q = composeSpray p [z, x ^+^ y ^+^ z]
 -- >>> putStrLn $ prettySprayXYZ q
 -- (1) X + (1) Y + (2) Z
-composeSpray :: forall a. (AlgRing.C a, Eq a) => Spray a -> [Spray a] -> Spray a
+composeSpray :: forall a. (AlgRing.C a, Eq a) 
+                => Spray a -> [Spray a] -> Spray a
 composeSpray p = evalSpray (identify p)
   where 
     identify :: Spray a -> Spray (Spray a)
@@ -532,7 +782,7 @@     | 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"
+-- | Pretty form of a spray, with monomials shown as "x1x3^2"
 --
 -- >>> x :: lone 1 :: Spray Int
 -- >>> y :: lone 2 :: Spray Int
@@ -540,7 +790,7 @@ -- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3
 -- >>> putStrLn $ prettySpray' p
 -- (2) x1 + (3) x2^2 + (-4) x3^3 
-prettySpray' :: (Show a) => Spray a -> String
+prettySpray' :: Show a => Spray a -> String
 prettySpray' spray = unpack $ intercalate (pack " + ") terms
  where
   terms           = map stringTerm 
@@ -554,11 +804,29 @@     stringCoef'  = cons '(' $ snoc stringCoef ')'
     stringCoef'' = if constant then stringCoef' else snoc stringCoef' ' '
 
+-- | Pretty form of a spray, with monomials shown as "x1x3^2", and with 
+-- a user-defined showing function for the coefficients
+--
+-- prop> prettySpray' p == prettySpray'' show p
+prettySpray'' :: (a -> String) -> Spray a -> String
+prettySpray'' showCoeff spray = unpack $ intercalate (pack " + ") terms
+ where
+  terms           = map stringTerm 
+                        (sortBy (flip compare `on` fexpts) (HM.toList spray))
+  fexpts term     = exponents $ fst term
+  stringTerm term = append stringCoef'' (prettyPowers' pows)
+   where
+    pows         = exponents (fst term)
+    constant     = S.null pows
+    stringCoef   = pack $ showCoeff (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"
+  else error "prettyPowersXYZ: there is more than three variables"
  where
   n     = S.length pows
   gpows = growSequence pows n 3
@@ -587,7 +855,7 @@   stringTerm term = append stringCoef'' (prettyPowersXYZ pows)
    where
     pows         = exponents (fst term)
-    constant     = S.length pows == 0
+    constant     = S.null pows
     stringCoef   = pack $ show (snd term)
     stringCoef'  = cons '(' $ snoc stringCoef ')'
     stringCoef'' = if constant then stringCoef' else snoc stringCoef' ' '
@@ -649,7 +917,8 @@ 
 -- | Remainder of the division of a spray by a list of divisors, 
 -- using the lexicographic ordering of the monomials
-sprayDivisionRemainder :: forall a. (Eq a, AlgField.C a) => Spray a -> [Spray a] -> Spray a
+sprayDivisionRemainder :: forall a. (Eq a, AlgField.C a) 
+                          => Spray a -> [Spray a] -> Spray a
 sprayDivisionRemainder p qs = 
   if n == 0 
     then error "sprayDivisionRemainder: the list of divisors is empty." 
@@ -689,7 +958,8 @@     then error "sprayDivision: division by zero."
     else ogo sprayA AlgAdd.zero AlgAdd.zero
   where
-    go :: Monomial a -> Spray a -> Spray a -> Spray a -> Int -> Bool -> (Spray a, Spray a, Spray a)
+    go :: Monomial a -> Spray a -> Spray a -> Spray a -> Int -> Bool 
+          -> (Spray a, Spray a, Spray a)
     go ltp !p !q r !i !divoccured
       | divoccured = (p, q, r)
       | i == 1     = (p ^-^ ltpspray, q, r ^+^ ltpspray)
@@ -714,8 +984,9 @@ -- Groebner stuff -------------------------------------------------------------
 
 -- | slight modification of `sprayDivisionRemainder` to speed up groebner00
-sprayDivisionRemainder' :: forall a. (Eq a, AlgField.C a) 
-                           => Spray a -> HashMap Int (Spray a, Monomial a) -> Spray a
+sprayDivisionRemainder' 
+  :: forall a. (Eq a, AlgField.C a) 
+  => Spray a -> HashMap Int (Spray a, Monomial a) -> Spray a
 sprayDivisionRemainder' p qsltqs = snd $ ogo p AlgAdd.zero
   where
     n = HM.size qsltqs
@@ -751,7 +1022,8 @@     row2   = drop s $ concatMap (\i -> replicate i i) range1
 
 -- the "S polynomial"
-sPolynomial :: (Eq a, AlgField.C a) => (Spray a, Monomial a) -> (Spray a, Monomial a) -> Spray a
+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
@@ -777,7 +1049,8 @@     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 :: 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'
@@ -840,7 +1113,7 @@ 
 -- | Groebner basis (always minimal and possibly reduced)
 --
--- prop> groebner ps True == reduceGroebnerBasis (groebner ps False)
+-- prop> groebner sprays True == reduceGroebnerBasis (groebner sprays False)
 groebner 
   :: forall a. (Eq a, AlgField.C a) 
   => [Spray a] -- ^ list of sprays 
@@ -899,13 +1172,34 @@   -> Int -- ^ index
   -> Spray a
 esPolynomial n k
-  | k <= 0 || n <= 0 = error "esPolynomial: both arguments must be positive integers."
-  | k > n            = AlgAdd.zero
-  | otherwise        = simplifySpray spray
+  | k < 0 || n < 0 
+    = error "esPolynomial: both arguments must be positive integers."
+  | k > n     = AlgAdd.zero
+  | k == 0    = unitSpray
+  | otherwise = simplifySpray spray
   where
     perms = permutationsBinarySequence (n-k) k
     spray = HM.fromList $ map (\expts -> (Powers expts n, AlgRing.one)) perms
 
+-- | Power sum polynomial
+psPolynomial 
+  :: forall a. (AlgRing.C a, Eq a) 
+  => Int -- ^ number of variables
+  -> Int -- ^ power
+  -> Spray a
+psPolynomial n k
+  | k < 0 || n < 0 
+    = error "psPolynomial: both arguments must be positive integers."
+  | k > n     = AlgAdd.zero
+  | k == 0    = n .^ unitSpray
+  | otherwise = spray
+  where
+    spray = HM.fromList $ map f [1 .. n]
+    f :: Int -> (Powers, a)
+    f j = (Powers expts j, AlgRing.one)
+      where
+        expts = S.replicate (j-1) 0 |> k
+
 -- | Whether a spray is a symmetric polynomial
 isSymmetricSpray :: forall a. (AlgField.C a, Eq a) => Spray a -> Bool
 isSymmetricSpray spray = check1 && check2 
@@ -914,7 +1208,8 @@     indices = [1 .. n]
     gPolys = map (\i -> esPolynomial n i ^-^ lone (n + i)) indices
     gbasis  = groebner0 gPolys
-    g       = sprayDivisionRemainder spray gbasis
+    spray'  = spray ^-^ (constantSpray (getConstantTerm spray))
+    g       = sprayDivisionRemainder spray' gbasis
     gpowers = HM.keys g
     check1  = minimum (map nvariables gpowers) > n
     expnts  = map exponents gpowers
@@ -931,28 +1226,31 @@ -- >>> p = p1 ^*^ p2
 -- 
 -- prop> isPolynomialOf p [p1, p2] == (True, Just $ x ^*^ y)
-isPolynomialOf :: forall a. (AlgField.C a, Eq a) => Spray a -> [Spray a] -> (Bool, Maybe (Spray a))
+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" 
+      | n < n'    = error "isPolynomialOf: 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       = sprayDivisionRemainder 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
+          m            = length sprays
+          yPolys       = map (\i -> lone (n + i) :: Spray a) [1 .. m]
+          gPolys       = zipWith (^-^) sprays yPolys
+          gbasis0      = groebner0 gPolys
+          constantTerm = constantSpray (getConstantTerm spray)
+          spray'       = spray ^-^ constantTerm
+          g            = sprayDivisionRemainder 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 ^+^ constantTerm
             else Nothing
           dropXis = HM.mapKeys f
           f (Powers expnnts _) = Powers (S.drop n expnnts) n
@@ -960,7 +1258,7 @@ 
 -- resultant ------------------------------------------------------------------
 
--- sylvester matrix
+-- | sylvester matrix
 sylvesterMatrix :: AlgAdd.C a => [a] -> [a] -> Matrix a
 sylvesterMatrix x y = fromLists (xrows ++ yrows) 
   where
@@ -971,11 +1269,13 @@     yrows = [replicate i AlgAdd.zero ++ y ++ replicate (m-i-1) AlgAdd.zero 
              | i <- [0 .. m-1]]
 
--- "truncated" Sylvester matrix
+-- | "truncated" Sylvester matrix
 sylvesterMatrix' :: AlgRing.C a => [a] -> [a] -> Int -> Matrix a
 sylvesterMatrix' x y k = if s == 0 
-  then fromLists [[AlgRing.one]] -- plays the role of the empty matrix: determinant=1 
-                                 -- (because the empty matrix is not allowed in the matrix package)
+  then fromLists [[AlgRing.one]] -- plays the role of the empty matrix: 
+                                 -- the point to get is determinant=1 
+                                 -- (because the empty matrix is not allowed
+                                 -- in the matrix package)
   else submatrix 1 s 1 s $ fromLists (xrows ++ yrows) 
   where
     m = length x - 1
@@ -986,29 +1286,33 @@     yrows = [replicate i AlgAdd.zero ++ y ++ replicate (m-i-1) AlgAdd.zero 
              | i <- [0 .. m-1-k]]
 
--- determinant
+-- | determinant of a matrix
 detLaplace :: forall a. (Eq a, AlgRing.C a) => Matrix a -> a
 detLaplace m = if nrows m == 1 
-  then m DM.! (1,1)
-  else suml1 [negateIf i (times (m DM.! (i,1)) (detLaplace (minorMatrix i 1 m))) 
-              | i <- [1 .. nrows m]]
+  then 
+    m DM.! (1,1)
+  else 
+    suml1 [negateIf i (times (m DM.! (i,1)) (detLaplace (minorMatrix i 1 m))) 
+           | i <- [1 .. nrows m]]
   where 
     suml1      = foldl1' (AlgAdd.+)
     negateIf i = if even i then AlgAdd.negate else id
     times :: a -> a -> a
     times x y = if x == AlgAdd.zero then AlgAdd.zero else x AlgRing.* y
 
--- the coefficients of a spray as a univariate spray in x with spray coefficients
+-- | the coefficients of a spray as a univariate spray in x with 
+-- spray coefficients
 sprayCoefficients :: (Eq a, AlgRing.C a) => Spray a -> [Spray a]
 sprayCoefficients spray = 
   if n == 0 
-    then [constantSpray constantTerm]
+    then [constantTerm]
     else reverse sprays
   where
     n = numberOfVariables spray 
     (powers, coeffs) = unzip (HM.toList spray)
     expnts           = map exponents powers
-    constantTerm = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) spray)
+    constantTerm = 
+      constantSpray $ fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) spray)
     (expnts', coeffs') = 
       unzip $ filter (\(s,_) -> S.length s > 0) (zip expnts coeffs)
     xpows              = map (`index` 0) expnts'
@@ -1016,7 +1320,7 @@     powers''           = map (\s -> Powers s (S.length s)) expnts''
     sprays''           = zipWith (curry fromMonomial) powers'' coeffs'
     imap               = IM.fromListWith (^+^) (zip xpows sprays'')
-    imap'              = IM.insertWith (^+^) 0 (constantSpray constantTerm) imap
+    imap'              = IM.insertWith (^+^) 0 constantTerm imap
     permutation = [2 .. n] ++ [1]
     sprays = [
         permuteVariables permutation (fromMaybe AlgAdd.zero (IM.lookup i imap')) 
@@ -1031,8 +1335,10 @@     else error "resultant1: the two sprays must be univariate."
   where
     n = max (numberOfVariables p) (numberOfVariables q)
-    pexpnts = map (`index` 0) $ filter (not . S.null) (map exponents (HM.keys p))
-    qexpnts = map (`index` 0) $ filter (not . S.null) (map exponents (HM.keys q))
+    pexpnts = 
+      map (`index` 0) $ filter (not . S.null) (map exponents (HM.keys p))
+    qexpnts = 
+      map (`index` 0) $ filter (not . S.null) (map exponents (HM.keys q))
     p0 = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) p)
     q0 = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) q)
     pcoeffs = if null pexpnts 
@@ -1055,8 +1361,10 @@   else error "subresultants1: the two sprays must be univariate."
   where
     n = max (numberOfVariables p) (numberOfVariables q)
-    pexpnts = map (`index` 0) $ filter (not . S.null) (map exponents (HM.keys p))
-    qexpnts = map (`index` 0) $ filter (not . S.null) (map exponents (HM.keys q))
+    pexpnts = 
+      map (`index` 0) $ filter (not . S.null) (map exponents (HM.keys p))
+    qexpnts = 
+      map (`index` 0) $ filter (not . S.null) (map exponents (HM.keys q))
     p0 = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) p)
     q0 = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) q)
     pcoeffs = if null pexpnts 
@@ -1086,8 +1394,8 @@     else error "resultant: invalid variable index."
   where
     n = max (numberOfVariables p) (numberOfVariables q)
-    permutation  = var : [1 .. var-1] ++ [var+1 .. n]
-    permutation' = [2 .. var] ++ (1 : [var+1 .. n])
+    permutation  = [n-var+2 .. n] ++ [1 .. n-var+1]
+    permutation' = [var .. n] ++ [1 .. var-1]
     p' = permuteVariables permutation p
     q' = permuteVariables permutation q
     det = detLaplace $ 
@@ -1095,7 +1403,7 @@ 
 -- | Subresultants of two sprays
 subresultants :: (Eq a, AlgRing.C a) 
-  => Int     -- ^ indicator of the variable with respect to which the resultant is desired (e.g. 1 for x)
+  => Int     -- ^ indicator of the variable with respect to which the subresultants are desired (e.g. 1 for x)
   -> Spray a 
   -> Spray a 
   -> [Spray a]
@@ -1117,10 +1425,74 @@     permutation' = [2 .. var] ++ (1 : [var+1 .. n])
     permute'     = permuteVariables permutation'
 
+-- | Resultant of two sprays with coefficients in a field; this function is more 
+-- efficient than the function `resultant`
+resultant' :: forall a. (Eq a, AlgField.C a) 
+  => Int     -- ^ indicator of the variable with respect to which the resultant is desired (e.g. 1 for x)
+  -> Spray a 
+  -> Spray a 
+  -> Spray a
+resultant' var sprayA sprayB 
+  | var < 1 || var > n                         
+    = error "resultant': invalid variable index." 
+  | sprayA == zeroSpray || sprayB == zeroSpray 
+    = zeroSpray
+  | otherwise 
+    = permuteVariables permutation' $ go unitSpray unitSpray s0 p0 q0
+  where
+    n = max (numberOfVariables sprayA) (numberOfVariables sprayB)
+    permutation  = [n-var+1 .. n] ++ [1 .. n-var]
+    permutation' = [var+1 .. n] ++ [1 .. var] 
+    sprayA' = permuteVariables permutation sprayA
+    sprayB' = permuteVariables permutation sprayB
+    degA = degree n sprayA'
+    degB = degree n sprayB'
+    content :: Spray a -> Spray a
+    content spray = foldl1' gcdSpray (sprayCoefficients' n spray)
+    exactDivisionBy :: Spray a -> Spray a -> Spray a
+    exactDivisionBy b a = 
+      if snd division == zeroSpray 
+        then fst division 
+        else error "exactDivisionBy: should not happen."
+      where
+        division = sprayDivision a b
+    contA = content sprayA'
+    contB = content sprayB'
+    sprayA'' = exactDivisionBy contA sprayA'
+    sprayB'' = exactDivisionBy contB sprayB'
+    t = contA^**^degB ^*^ contB^**^degA
+    s0 = if degA < degB && odd degA && odd degB 
+      then AlgAdd.negate unitSpray :: Spray a
+      else unitSpray
+    (p0, q0) = if degA >= degB
+      then (sprayA'', sprayB'')
+      else (sprayB'', sprayA'')
+    go :: Spray a -> Spray a -> Spray a -> Spray a -> Spray a -> Spray a
+    go g h s p q = 
+      if degq' == 0
+        then s' ^*^ t ^*^ h''
+        else go g' h' s' p' q'
+        where
+          degp           = degree n p
+          degq           = degree n q
+          delta          = degp - degq
+          s' = if odd degp && odd degq 
+            then AlgAdd.negate s 
+            else s
+          (_, (_, r)) = pseudoDivision n p q
+          p'             = q
+          q'             = exactDivisionBy (g ^*^ h^**^delta) r
+          (degp', ellp') = degreeAndLeadingCoefficient n p'
+          (degq', ellq') = degreeAndLeadingCoefficient n q'
+          g'  = ellp'
+          h'  = exactDivisionBy (h^**^delta) (h ^*^ g'^**^delta)
+          h'' = exactDivisionBy (h'^**^degp') (h' ^*^ ellq'^**^degp')
 
+
 -- GCD stuff ------------------------------------------------------------------
 
--- the coefficients of a spray as a univariate spray in x_n with spray coefficients
+-- | the coefficients of a spray as a univariate spray in x_n with 
+-- spray coefficients
 sprayCoefficients' :: (Eq a, AlgRing.C a) => Int -> Spray a -> [Spray a]
 sprayCoefficients' n spray 
   | numberOfVariables spray /= n = [spray]
@@ -1146,11 +1518,13 @@         | i <- [deg, deg-1 .. 0]
       ]
 
--- the degree of a spray as a univariate spray in x_n with spray coefficients
+-- | the degree of a spray as a univariate spray in x_n with spray coefficients
 degree :: (Eq a, AlgAdd.C a) => Int -> Spray a -> Int
 degree n spray 
   | numberOfVariables spray == 0 = 
-      if spray == zeroSpray then minBound else 0
+      if spray == zeroSpray 
+        then minBound -- (should not happen)
+        else 0
   | numberOfVariables spray /= n = 0
   | otherwise                    = maximum xpows
     where
@@ -1160,13 +1534,14 @@       expnts'     = filter (not . S.null) expnts
       xpows       = map (`index` 0) expnts'
 
--- the degree and the leading coefficient of a spray as a univariate spray 
+-- | the degree and the leading coefficient of a spray as a univariate spray 
 -- in x_n with spray coefficients
-degreeAndLeadingCoefficient :: (Eq a, AlgRing.C a) => Int -> Spray a -> (Int, Spray a)
+degreeAndLeadingCoefficient :: (Eq a, AlgRing.C a) 
+                                => Int -> Spray a -> (Int, Spray a)
 degreeAndLeadingCoefficient n spray 
   | n == 0                       = (
                                     if constantTerm == AlgAdd.zero 
-                                      then minBound 
+                                      then minBound -- (should not happen)
                                       else 0, 
                                     constantSpray constantTerm
                                    )
@@ -1186,9 +1561,10 @@     expnts'' = [S.deleteAt 0 (expnts' !! i) | i <- is]
     powers'' = map (\s -> Powers s (S.length s)) expnts''
     coeffs'' = [coeffs' !! i | i <- is]
-    leadingCoeff = foldl1' (^+^) (zipWith (curry fromMonomial) powers'' coeffs'')
+    leadingCoeff = 
+      foldl1' (^+^) (zipWith (curry fromMonomial) powers'' coeffs'')
 
--- pseudo-division of two sprays, assuming degA >= degB >= 0
+-- | pseudo-division of two sprays, assuming degA >= degB >= 0
 pseudoDivision :: (Eq a, AlgRing.C a)
   => Int                           -- ^ number of variables
   -> Spray a                       -- ^ A
@@ -1214,8 +1590,9 @@         sprayXn      = lone n 
         sprayS       = ellR ^*^ sprayXn ^**^ (degR - degB)
 
--- recursive GCD function
-gcdKX1dotsXn :: forall a. (Eq a, AlgField.C a) => Int -> Spray a -> Spray a -> Spray a
+-- | recursive GCD function
+gcdKX1dotsXn :: forall a. (Eq a, AlgField.C a) 
+                => Int -> Spray a -> Spray a -> Spray a
 gcdKX1dotsXn n sprayA sprayB
   | n == 0              = constantSpray $ gcdKX0 sprayA sprayB
   | degB > degA         = gcdKX1dotsXn n sprayB sprayA 

@@ -1,5 +1,8 @@ module Main where
+import qualified Algebra.Additive               as AlgAdd             
+import qualified Algebra.Ring                   as AlgRing      
 import           Approx                         ( assertApproxEqual )
+import           Data.Maybe                     ( fromJust )
 import           Data.Ratio                     ( (%) )
 import           Math.Algebra.Hspray            ( Spray,
                                                   (^+^),
@@ -12,6 +15,7 @@                                                   zeroSpray,
                                                   constantSpray,
                                                   getCoefficient,
+                                                  getConstantTerm,
                                                   evalSpray,
                                                   substituteSpray,
                                                   composeSpray,
@@ -24,15 +28,24 @@                                                   groebner,
                                                   fromRationalSpray,
                                                   esPolynomial,
-                                                  isSymmetricSpray, 
+                                                  psPolynomial,
+                                                  isSymmetricSpray,
                                                   isPolynomialOf,
                                                   resultant,
+                                                  resultant',
                                                   subresultants,
                                                   resultant1,
                                                   subresultants1,
                                                   sprayDivision,
-                                                  gcdSpray
+                                                  gcdSpray,
+                                                  QSpray,
+                                                  SymbolicQSpray,
+                                                  evalRatioOfPolynomials,
+                                                  evalSymbolicSpray',
+                                                  qpolyFromCoeffs,
+                                                  constQPoly, evalSymbolicSpray''
                                                 )
+import           Number.Ratio                   ( T ( (:%) ) )
 import           Test.Tasty                     ( defaultMain
                                                 , testGroup
                                                 )
@@ -78,6 +91,14 @@         p = 2 *^ (2 *^ (x^**^3 ^*^ y^**^2)) ^+^ 4 *^ z ^+^ 5 *^ unitSpray
       assertEqual "" (getCoefficient [3, 2, 0] p, getCoefficient [0, 4] p) (4, 0),
 
+    testCase "getConstantTerm" $ do
+      let
+        x = lone 1 :: Spray Int
+        y = lone 2 :: Spray Int
+        z = lone 3 :: Spray Int
+        p = 2 *^ (2 *^ (x^**^3 ^*^ y^**^2)) ^+^ 4 *^ z ^+^ 5 *^ unitSpray
+      assertEqual "" (getConstantTerm p) 5,
+
     testCase "fromList . toList = identity" $ do
       let
         x = lone 1 :: Spray Int
@@ -113,7 +134,7 @@         sumAbsValues = sum $ map abs gxyz
       assertApproxEqual "" 8 sumAbsValues 0,
 
-    testCase "symmetric polynomials" $ do
+    testCase "symmetric polynomial" $ do
       let
         e2 = esPolynomial 4 2 :: Spray Rational
         e3 = esPolynomial 4 3 :: Spray Rational
@@ -134,9 +155,20 @@         y = lone 2 :: Spray Rational
         p1 = x ^+^ y
         p2 = x ^-^ y
-        p = p1 ^*^ p2
-      assertEqual "" (isPolynomialOf p [p1, p2]) (True, Just $ x ^*^ y),
+        p = p1 ^*^ p2 ^+^ unitSpray
+      assertEqual "" (isPolynomialOf p [p1, p2]) (True, Just $ x ^*^ y ^+^ unitSpray),
 
+    testCase "power sum polynomials" $ do
+      let
+        x = lone 1 :: Spray Rational
+        y = lone 2 :: Spray Rational
+        symSpray = x^**^2 ^+^ y^**^2 ^+^ x ^+^ y
+        p1 = psPolynomial 2 1 :: Spray Rational
+        p2 = psPolynomial 2 2 :: Spray Rational
+        p = fromJust $ snd $ isPolynomialOf symSpray [p1, p2]
+        symSpray' = composeSpray p [p1, p2]
+      assertEqual "" symSpray symSpray',
+
     testCase "substituteSpray" $ do
       let
         x1 = lone 1 :: Spray Rational
@@ -233,6 +265,24 @@         r1 = resultant1 f g
       assertEqual "" r1 (getCoefficient [] r),
 
+    testCase "resultant and resultant' are in agreement" $ do
+      let
+        x = lone 1 :: Spray Rational
+        y = lone 2 :: Spray Rational
+        z = lone 3 :: Spray Rational
+        p = x^**^4 ^-^ x^**^3 ^+^ x^**^2 ^-^ 2*^ (x ^*^ y^**^2) ^+^ z^**^4 
+        q = x ^-^ (2*^ y^**^2) ^*^ z^**^2 ^*^ unitSpray
+        rx  = resultant 1 p q
+        rx' = resultant' 1 p q
+        ry  = resultant 2 p q
+        ry' = resultant' 2 p q
+        rz  = resultant 3 p q
+        rz' = resultant' 3 p q
+        test1 = rx == rx'
+        test2 = ry == ry'
+        test3 = rz == rz'
+      assertBool "" (test1 && test2 && test3),
+
     testCase "gcdSpray - univariate example" $ do
       let
         x = lone 1 :: Spray Rational
@@ -280,6 +330,29 @@         sprayA = sprayD^**^1 ^*^ (x^**^4  ^-^  x  ^+^  y   ^+^  x ^*^ y ^*^ z^**^2)
         sprayB = sprayD^**^1 ^*^ y ^*^ (2*^x  ^+^  unitSpray) ^*^ z
         g = gcdSpray sprayA sprayB
-      assertEqual "" g sprayD
+      assertEqual "" g sprayD,
+
+    testCase "evaluation of symbolic spray" $ do
+      let 
+        a    = qpolyFromCoeffs [0, 1]  
+        p    = a AlgRing.^ 2 AlgAdd.- constQPoly 4 
+        q1   = a AlgAdd.- constQPoly 3
+        q2   = a AlgAdd.- constQPoly 2
+        rop1 = p :% q1 
+        rop2 = p :% q2
+        f :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> Spray a -> (Spray a, Spray a)
+        f x y z = (x^**^2 ^+^ y^**^2, z)
+        g :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> Spray a -> (a, a, a) -> (a, a)
+        g px py pz (x, y, z) = (evalSpray f1 [x, y, AlgAdd.zero], evalSpray f2 [AlgAdd.zero, AlgAdd.zero, z])
+          where (f1, f2) = f px py pz
+        (r1, r2) = g (lone 1 :: QSpray) (lone 2) (lone 3) (2, 3, 4) 
+        r = evalRatioOfPolynomials 5 rop1 AlgRing.* r1  AlgAdd.+  evalRatioOfPolynomials 5 rop2 AlgRing.* r2
+        (f1', f2')  = f (lone 1 :: SymbolicQSpray) (lone 2) (lone 3)
+        symSpray  = rop1 *^ f1'  ^+^  rop2 *^ f2' 
+        r' = evalSymbolicSpray' symSpray 5 [2, 3, 4]
+        rop1' = evalSymbolicSpray'' f1' [2, 3]
+        rop2' = evalSymbolicSpray'' f2' [0, 0, 4]
+        r'' = evalRatioOfPolynomials 5 (rop1 AlgRing.* rop1' AlgAdd.+ rop2 AlgRing.* rop2')
+      assertEqual "" (r, r') (r', r'') 
 
   ]