{-|
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 commutative ring. 
See README for examples.
-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}

module Math.Algebra.Hspray
  ( 
  -- * Types
    Powers (..)
  , Spray
  , QSpray
  , Monomial
  -- * Basic sprays
  , lone
  , unitSpray
  , zeroSpray
  , constantSpray
  -- * Operations on sprays
  , (*^)
  , (.^)
  , (^+^)
  , (^-^)
  , (^*^)
  , (^**^)
  -- * 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
  , substituteSpray
  , composeSpray
  -- * Differentiation of a spray
  , derivSpray
  -- * Permutation of the variables of a spray
  , permuteVariables
  , swapVariables
  -- * Division of a spray
  , sprayDivision
  , sprayDivisionRemainder
  -- * Gröbner basis
  , groebner
  , reduceGroebnerBasis
  -- * Symmetric polynomials
  , esPolynomial
  , psPolynomial
  , isSymmetricSpray
  -- * Resultant and subresultants
  , resultant
  , resultant'
  , resultant1
  , subresultants
  , subresultants1
  -- * Greatest common divisor
  , gcdSpray
  -- * Miscellaneous
  , fromList
  , toList
  , fromRationalSpray
  , leadingTerm
  , isPolynomialOf
  , bombieriSpray
  ) 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 Algebra.ZeroTestable          as AlgZT
import qualified Data.Foldable                 as DF
import           Data.Function                  ( on )
import           Data.HashMap.Strict            ( HashMap )
import qualified Data.HashMap.Strict           as HM
import           Data.Hashable                  ( Hashable(hashWithSalt) )
import qualified Data.IntMap.Strict            as IM
import           Data.List                      ( sortBy
                                                , maximumBy 
                                                , (\\)
                                                , findIndices
                                                , elemIndices
                                                , nub
                                                , foldl1'
                                                )
import           Data.Matrix                    ( Matrix 
                                                , fromLists
                                                , minorMatrix
                                                , nrows
                                                , submatrix
                                                )
import qualified Data.Matrix                   as DM
import           Data.Maybe                     ( isJust
                                                , fromJust, fromMaybe
                                                )
import           Data.Ord                       ( comparing )
import qualified Data.Sequence                 as S
import           Data.Sequence                  ( (><)
                                                , Seq 
                                                , dropWhileR
                                                , (|>)
                                                , index
                                                , adjust
                                                , fromFunction
                                                )
import           Data.Text                      ( Text
                                                , append
                                                , cons
                                                , intercalate
                                                , pack
                                                , 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


-- Univariate polynomials -----------------------------------------------------

newtype A a = A a 
  deriving
    (Eq, AlgAdd.C, AlgRing.C, AlgField.C)

type Rational' = NumberRatio.Rational
type Q = A Rational'

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

growSequence' :: Int -> Seq Int -> Seq Int
growSequence' n s = growSequence s (S.length s) n

-- | append trailing zeros to get the same length
harmonize :: (Powers, Powers) -> (Powers, Powers)
harmonize (pows1, pows2) = (Powers e1' n, Powers e2' n)
 where
  e1            = exponents pows1
  e2            = exponents pows2
  n1            = nvariables pows1
  n2            = nvariables pows2
  (e1', e2', n) = if n1 < n2
    then (growSequence e1 n1 n2, e2, n2)
    else (e1, growSequence e2 n2 n1, n1)

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

-- | 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
  zero :: Spray a
  zero   = HM.empty
  negate :: Spray a -> Spray a
  negate = negateSpray

instance (AlgRing.C a, Eq a) => AlgMod.C a (Spray a) where
  (*>) :: a -> Spray a -> Spray a
  lambda *> p = scaleSpray lambda p

instance (AlgRing.C a, Eq a) => AlgRing.C (Spray a) where
  (*) :: Spray a -> Spray a -> Spray a
  p * q = multSprays p q
  one :: Spray a
  one   = lone 0

{- instance (AlgRing.C a, Eq a) => Num (Spray a) where
  p + q = addSprays p q
  negate = negateSpray
  p * q = multSprays p q
  fromInteger n = fromInteger n .^ AlgRing.one
  abs _ = error "Prelude.Num.abs: inappropriate abstraction"
  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."

infixr 7 *^
-- | Scale a spray by a scalar
(*^) :: (AlgRing.C a, Eq a) => a -> Spray a -> Spray a
(*^) lambda pol = lambda AlgMod.*> pol

infixr 7 .^
-- | Scale a spray by an integer
--
-- prop> 3 .^ p == p ^+^ p ^+^ p
(.^) :: (AlgAdd.C a, Eq a) => Int -> Spray a -> Spray a
(.^) k pol = if k >= 0
  then AlgAdd.sum (replicate k pol)
  else AlgAdd.negate $ AlgAdd.sum (replicate (-k) pol)

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

-- | 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')
   where
    i'     = i - 1
    expts  = exponents pows
    expt_i = expts `index` i'
    expts' = adjust (subtract 1) i' expts
    coef'  = AlgAdd.sum (replicate expt_i coef)
    pows'  = Powers expts' (nvariables pows) 

-- | Derivative of a spray
derivSpray 
  :: (AlgRing.C a, Eq a) 
  => Int     -- ^ index of the variable of differentiation (starting at 1)
  -> Spray a -- ^ the spray to be derivated
  -> Spray a -- ^ the derivated spray
derivSpray i p = if i >= 1 
  then cleanSpray $ HM.fromListWith (AlgAdd.+) monomials
  else error "derivSpray: invalid index."
 where
  p'        = HM.toList p
  monomials = [ derivMonomial i mp | mp <- p' ]

-- | Spray corresponding to the basic monomial x_n
--
-- >>> x :: lone 1 :: Spray Int
-- >>> y :: lone 2 :: Spray Int
-- >>> p = 2*^x^**^2 ^-^ 3*^y
-- >>> putStrLn $ prettySpray' p
-- (2) x1^2 + (-3) x2
--
-- prop> lone 0 == unitSpray
lone :: AlgRing.C a => Int -> Spray a
lone n = if n >= 0 
  then HM.singleton pows AlgRing.one
  else error "lone: invalid index."
 where
  pows = if n == 0
    then Powers S.empty 0
    else Powers (S.replicate (n - 1) AlgAdd.zero |> AlgRing.one) n

-- | The unit spray
--
-- prop> p ^*^ unitSpray == p
unitSpray :: AlgRing.C a => Spray a
unitSpray = lone 0

-- | The null spray
--
-- prop> p ^+^ zeroSpray == p
zeroSpray :: (Eq a, AlgAdd.C a) => Spray a
zeroSpray = AlgAdd.zero

-- | Constant spray
--
-- prop> constantSpray 3 == 3 *^ unitSpray
constantSpray :: (AlgRing.C a, Eq a) => a -> Spray a
constantSpray c = c *^ lone 0

-- | Get coefficient of a term of a spray 
--
-- >>> 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
-- >>> getCoefficient [3, 2, 0] p
-- 4
-- >>> getCoefficient [0, 4] p
-- 0
getCoefficient :: AlgAdd.C a => [Int] -> Spray a -> a
getCoefficient expnts spray = fromMaybe AlgAdd.zero (HM.lookup powers spray)
  where
    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 =
  if null powers then 0 else maximum (map nvariables powers)
  where
    powers = HM.keys spray

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

-- | Evaluates a spray
--
-- >>> x :: lone 1 :: Spray Int
-- >>> y :: lone 2 :: Spray Int
-- >>> p = 2*^x^**^2 ^-^ 3*^y
-- >>> evalSpray p [2, 1]
-- 5
evalSpray :: AlgRing.C a => Spray a -> [a] -> a
evalSpray p xyz = if length xyz >= numberOfVariables p
  then AlgAdd.sum $ map (evalMonomial xyz) (HM.toList p)
  else error "evalSpray: not enough values provided."

-- | spray from monomial
fromMonomial :: Monomial a -> Spray a
fromMonomial (pows, coeff) = HM.singleton pows coeff

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

-- | Substitutes some variables in a spray by some values
--
-- >>> x1 :: lone 1 :: Spray Int
-- >>> x2 :: lone 2 :: Spray Int
-- >>> x3 :: lone 3 :: Spray Int
-- >>> p = x1^**^2 ^-^ x2 ^+^ x3 ^-^ unitSpray
-- >>> p' = substituteSpray [Just 2, Nothing, Just 3] p
-- >>> putStrLn $ prettySpray' p'
-- (-1) x2 + (6) 
substituteSpray :: (Eq a, AlgRing.C a) => [Maybe a] -> Spray a -> Spray a
substituteSpray subs spray = if length subs == n 
  then spray'
  else error "substituteSpray: incorrect length of the substitutions list."
  where
    n         = numberOfVariables spray
    monomials = HM.toList spray
    spray'    = 
      foldl1' (^+^) (map (fromMonomial . substituteMonomial subs) monomials)

-- | Converts a spray with rational coefficients to a spray with double 
-- coefficients (useful for evaluation)
fromRationalSpray :: Spray Rational -> Spray Double
fromRationalSpray = HM.map fromRational

-- | Sustitutes the variables of a spray with some sprays 
-- (e.g. change of variables)
--
-- >>> x :: lone 1 :: Spray Int
-- >>> y :: lone 2 :: Spray Int
-- >>> z :: lone 3 :: Spray Int
-- >>> p = x ^+^ y
-- >>> 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 p = evalSpray (identify p)
  where 
    identify :: Spray a -> Spray (Spray a)
    identify = HM.map constantSpray

-- | Creates a spray from a list of terms
fromList :: (AlgRing.C a, Eq a) => [([Int], a)] -> Spray a
fromList x = cleanSpray $ HM.fromList $ map
  (\(expts, coef) -> (Powers (S.fromList expts) (length expts), coef)) x

-- | Permutes the variables of a spray
--
-- >>> f :: Spray Rational -> Spray Rational -> Spray Rational -> Spray Rational
-- >>> f p1 p2 p3 = p1^**^4 ^+^ (2*^p2^**^3) ^+^ (3*^p3^**^2) ^-^ (4*^unitSpray)
-- >>> x1 = lone 1 :: Spray Rational
-- >>> x2 = lone 2 :: Spray Rational
-- >>> x3 = lone 3 :: Spray Rational
-- >>> p = f x1 x2 x3
--
-- prop> permuteVariables [3, 1, 2] p == f x3 x1 x2
permuteVariables :: [Int] -> Spray a -> Spray a
permuteVariables permutation spray = 
  if n' >= n && isPermutation permutation  
    then spray'
    else error "permuteVariables: invalid permutation."
  where
    n  = numberOfVariables spray
    n' = maximum permutation
    isPermutation pmtn = minimum pmtn == 1 && length (nub pmtn) == n'
    intmap         = IM.fromList (zip permutation [1 .. n'])
    invpermutation = [intmap IM.! i | i <- [1 .. n']]
    permuteSeq x   = 
      S.mapWithIndex (\i _ -> x `index` (invpermutation !! i - 1)) x 
    (powers, coeffs) = unzip (HM.toList spray)
    expnts  = map exponents powers
    expnts' = map (permuteSeq . growSequence' n') expnts
    powers' = map (\exps -> simplifyPowers (Powers exps n')) expnts'
    spray'  = HM.fromList (zip powers' coeffs)

-- | Swaps two variables of a spray
-- 
-- prop> swapVariables (1, 3) p == permuteVariables [3, 2, 1] p
swapVariables :: (Int, Int) -> Spray a -> Spray a
swapVariables (i, j) spray = 
  if i>=1 && j>=1  
    then spray'
    else error "swapVariables: invalid indices."
  where
    n = maximum [numberOfVariables spray, i, j]
    f k | k == i    = j
        | k == j    = i
        | otherwise = k
    transposition = map f [1 .. n]
    permuteSeq x  = 
      S.mapWithIndex (\ii _ -> x `index` (transposition !! ii - 1)) x 
    (powers, coeffs) = unzip (HM.toList spray)
    expnts  = map exponents powers
    expnts' = map (permuteSeq . growSequence' n) expnts
    powers' = map (\exps -> simplifyPowers (Powers exps n)) expnts'
    spray'  = HM.fromList (zip powers' coeffs)


-- 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
  x      = " " ++ var ++ "^"
  string = intercalate (pack ", ") (map (pack . show) pows)

-- | Pretty form of a spray
--
-- >>> x :: lone 1 :: Spray Int
-- >>> y :: lone 2 :: Spray Int
-- >>> z :: lone 3 :: Spray Int
-- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3
-- >>> putStrLn $ prettySpray show "x" p
-- (2) * x^(1) + (3) * x^(0, 2) + (-4) * x^(0, 0, 3)
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 (flip compare `on` fexpts) (HM.toList p))
  fexpts term     = exponents $ fst term
  stringTerm term = append
    (snoc (snoc (cons '(' $ snoc stringCoef ')') ' ') '*')
    (prettyPowers var pows)
   where
    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 shown as "x1x3^2"
--
-- >>> x :: lone 1 :: Spray Int
-- >>> y :: lone 2 :: Spray Int
-- >>> z :: lone 3 :: Spray Int
-- >>> 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' 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.length pows == 0
    stringCoef   = pack $ show (snd term)
    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 "prettyPowersXYZ: 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
--
-- >>> x :: lone 1 :: Spray Int
-- >>> y :: lone 2 :: Spray Int
-- >>> z :: lone 3 :: Spray Int
-- >>> p = 2*^x ^+^ 3*^y^**^2 ^-^ 4*^z^**^3
-- >>> putStrLn $ prettySprayXYZ p
-- (2) X + (3) Y^2 + (-4) Z^3
prettySprayXYZ :: (Show a) => Spray a -> String
prettySprayXYZ 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'' (prettyPowersXYZ pows)
   where
    pows         = exponents (fst term)
    constant     = S.null pows
    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

-- | Spray as a list
toList :: Spray a -> [([Int], a)]
toList p = HM.toList $ HM.mapKeys (DF.toList . exponents) p

-- | Bombieri spray (for internal usage in the \'scubature\' library)
bombieriSpray :: AlgAdd.C a => Spray a -> Spray a
bombieriSpray = HM.mapWithKey f
 where
  f pows          = times (pfactorial $ exponents pows)
  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 (uncurry (<=)) (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
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." 
    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

-- | Division of a spray by a spray
sprayDivision :: forall a. (Eq a, AlgField.C a) 
  => Spray a            -- ^ dividend 
  -> Spray a            -- ^ divisor
  -> (Spray a, Spray a) -- ^ (quotient, remainder)
sprayDivision sprayA sprayB =
  if sprayB == AlgAdd.zero 
    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 ltp !p !q r !i !divoccured
      | divoccured = (p, q, r)
      | i == 1     = (p ^-^ ltpspray, q, r ^+^ ltpspray)
      | otherwise  = go ltp newp newq r 1 newdivoccured
        where
          ltpspray      = fromMonomial ltp
          ltB           = leadingTerm sprayB
          newdivoccured = divides ltB ltp
          (newp, newq)  = if newdivoccured
            then (p ^-^ (qtnt ^*^ sprayB), q ^+^ qtnt)
            else (p, q)
            where
              qtnt = fromMonomial $ quotient ltp ltB
    ogo :: Spray a -> Spray a -> Spray a -> (Spray a, Spray a)
    ogo !p !q !r 
      | p == AlgAdd.zero = (q, r)
      | otherwise        = ogo p' q' r'
        where
          (p', q', r') = go (leadingTerm p) p q r 0 False


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

-- combinations of two among n
combn2 :: Int -> Int -> HashMap Int (Int, Int)
combn2 n s = HM.fromList (zip range0 (zip row1 row2)) 
  where
    range0 = [0 .. n-2]
    range1 = [1 .. n-1]
    row1   = drop s $ concatMap (\i -> [0 .. i-1]) range1 
    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 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   = sprayDivisionRemainder' 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 []

-- | Reduces 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 = sprayDivisionRemainder (ngbasis !! i) rest
      where
        rest = [ngbasis !! k | k <- [0 .. n-1] \\ [i]]

-- | Groebner basis (always minimal and possibly reduced)
--
-- prop> groebner sprays True == reduceGroebnerBasis (groebner sprays False)
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 -------------------------------------------

-- | combinations of k elements among a list
combinationsOf :: Int -> [a] -> [[a]]
combinationsOf _ []        = error "combinationsOf: should not happen."
combinationsOf 1 as        = map pure as
combinationsOf k as@(_:xs) = 
  run (l-1) (k-1) as $ combinationsOf (k-1) xs
  where
    l = length as
    run :: Int -> Int -> [a] -> [[a]] -> [[a]]
    run n i ys cs 
      | n == i    = map (ys ++) cs
      | otherwise = map (q:) cs ++ run (n-1) i qs (drop dc cs)
      where
        f :: [a] -> (a, [a])
        f []     = error "combinationsOf: should not happen."
        f (b:bs) = (b, bs)
        (q, qs)  = f (take (n-i+1) ys)
        dc       = product [(n-k+1) .. (n-1)] `div` product [1 .. i-1]

-- | generates all permutations of a binary sequence
permutationsBinarySequence :: Int -> Int -> [Seq Int]
permutationsBinarySequence nzeros nones = 
  let n = nzeros + nones in 
    map (binarySequence n) (combinationsOf nones [0 .. n-1])
  where
    binarySequence :: Int -> [Int] -> Seq Int
    binarySequence n combo = fromFunction n f 
      where
        f :: Int -> Int
        f i = fromEnum (i `elem` combo)

-- | Elementary symmetric polynomial
--
-- >>> putStrLn $ prettySpray' (esPolynomial 3 2)
-- (1) x1x2 + (1) x1x3 + (1) x2x3
esPolynomial 
  :: (AlgRing.C a, Eq a) 
  => Int -- ^ number of variables
  -> Int -- ^ index
  -> Spray a
esPolynomial n k
  | 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 
  where
    n = numberOfVariables spray
    indices = [1 .. n]
    gPolys = map (\i -> esPolynomial n i ^-^ lone (n + i)) indices
    gbasis  = groebner0 gPolys
    spray'  = spray ^-^ (constantSpray (getConstantTerm spray))
    g       = sprayDivisionRemainder 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); 
-- this polynomial is returned if this is true
--
-- >>> x = lone 1 :: Spray Rational
-- >>> y = lone 2 :: Spray Rational
-- >>> p1 = x ^+^ y
-- >>> p2 = x ^-^ y
-- >>> 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 spray sprays = result 
  where
    n = numberOfVariables spray
    n' = maximum $ map numberOfVariables sprays
    result
      | n > n'    = (False, Nothing)
      | 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
          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


-- resultant ------------------------------------------------------------------

-- | sylvester matrix
sylvesterMatrix :: AlgAdd.C a => [a] -> [a] -> Matrix a
sylvesterMatrix x y = fromLists (xrows ++ yrows) 
  where
    m = length x - 1
    n = length y - 1
    xrows = [replicate i AlgAdd.zero ++ x ++ replicate (n-i-1) AlgAdd.zero 
             | i <- [0 .. n-1]]
    yrows = [replicate i AlgAdd.zero ++ y ++ replicate (m-i-1) AlgAdd.zero 
             | i <- [0 .. m-1]]

-- | "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: 
                                 -- 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
    n = length y - 1
    s = m + n - 2*k
    xrows = [replicate i AlgAdd.zero ++ x ++ replicate (n-i-1) AlgAdd.zero 
             | i <- [0 .. n-1-k]]
    yrows = [replicate i AlgAdd.zero ++ y ++ replicate (m-i-1) AlgAdd.zero 
             | i <- [0 .. m-1-k]]

-- | 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]]
  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
sprayCoefficients :: (Eq a, AlgRing.C a) => Spray a -> [Spray a]
sprayCoefficients spray = 
  if n == 0 
    then [constantTerm]
    else reverse sprays
  where
    n = numberOfVariables spray 
    (powers, coeffs) = unzip (HM.toList spray)
    expnts           = map exponents powers
    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'
    expnts''           = map (S.deleteAt 0) expnts'
    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 constantTerm imap
    permutation = [2 .. n] ++ [1]
    sprays = [
        permuteVariables permutation (fromMaybe AlgAdd.zero (IM.lookup i imap')) 
        | i <- [0 .. maximum xpows]
      ]

-- | Resultant of two /univariate/ sprays
resultant1 :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> a
resultant1 p q = 
  if n <= 1 
    then detLaplace $ sylvesterMatrix pcoeffs qcoeffs
    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))
    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 
      then [p0]
      else [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) p) 
            | i <- [maxp, maxp-1 .. 1]] ++ [p0]
      where
        maxp = maximum pexpnts
    qcoeffs = if null qexpnts 
      then [q0]
      else [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) q) 
            | i <- [maxq, maxq-1 .. 1]] ++ [q0]
      where
        maxq = maximum qexpnts

-- | Subresultants of two /univariate/ sprays
subresultants1 :: (Eq a, AlgRing.C a) => Spray a -> Spray a -> [a]
subresultants1 p q = if n <= 1 
  then map (detLaplace . sylvesterMatrix' pcoeffs qcoeffs) [0 .. min d e - 1]
  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))
    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 
      then [p0]
      else [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) p) 
            | i <- [maxp, maxp-1 .. 1]] ++ [p0]
      where
        maxp = maximum pexpnts
    qcoeffs = if null qexpnts 
      then [q0]
      else [fromMaybe AlgAdd.zero (HM.lookup (Powers (S.singleton i) 1) q) 
            | i <- [maxq, maxq-1 .. 1]] ++ [q0]
      where
        maxq = maximum qexpnts
    d = length pcoeffs
    e = length qcoeffs

-- | Resultant of two sprays
resultant :: (Eq a, AlgRing.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 p q = 
  if var >= 1 && var <= n 
    then permuteVariables permutation' det
    else error "resultant: invalid variable index."
  where
    n = max (numberOfVariables p) (numberOfVariables q)
    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 $ 
          sylvesterMatrix (sprayCoefficients p') (sprayCoefficients q')

-- | Subresultants of two sprays
subresultants :: (Eq a, AlgRing.C a) 
  => Int     -- ^ indicator of the variable with respect to which the subresultants are desired (e.g. 1 for x)
  -> Spray a 
  -> Spray a 
  -> [Spray a]
subresultants var p q 
  | var < 1 = error "subresultants: invalid variable index."
  | var > n = error "subresultants: too large variable index."
  | otherwise = map (permute' . detLaplace . sylvesterMatrix' pcoeffs qcoeffs) 
                    [0 .. min d e - 1]
  where
    pcoeffs = sprayCoefficients p'
    qcoeffs = sprayCoefficients q'
    d = length pcoeffs
    e = length qcoeffs
    n = max (numberOfVariables p) (numberOfVariables q)
    permutation = var : [1 .. var-1] ++ [var+1 .. n]
    permute     = permuteVariables permutation
    p' = permute p 
    q' = permute q 
    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
sprayCoefficients' :: (Eq a, AlgRing.C a) => Int -> Spray a -> [Spray a]
sprayCoefficients' n spray 
  | numberOfVariables spray /= n = [spray]
  | n == 0                       = [constantSpray constantTerm]
  | otherwise                    = sprays 
  where
    permutation = [2 .. n] ++ [1]
    spray'      = permuteVariables permutation spray
    (powers, coeffs) = unzip (HM.toList spray')
    expnts           = map exponents powers
    constantTerm = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) spray')
    (expnts', coeffs') = 
      unzip $ filter (\(s,_) -> (not . S.null) s) (zip expnts coeffs)
    xpows = map (`index` 0) expnts'
    expnts'' = map (S.deleteAt 0) expnts'
    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
    deg    = maximum xpows
    sprays = [
        fromMaybe AlgAdd.zero (IM.lookup i imap')
        | i <- [deg, deg-1 .. 0]
      ]

-- | 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 -- (should not happen)
        else 0
  | numberOfVariables spray /= n = 0
  | otherwise                    = maximum xpows
    where
      permutation = [2 .. n] ++ [1]
      spray'      = permuteVariables permutation spray
      expnts      = map exponents $ HM.keys spray'
      expnts'     = filter (not . S.null) expnts
      xpows       = map (`index` 0) expnts'

-- | 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 n spray 
  | n == 0                       = (
                                    if constantTerm == AlgAdd.zero 
                                      then minBound -- (should not happen)
                                      else 0, 
                                    constantSpray constantTerm
                                   )
  | numberOfVariables spray /= n = (0, spray)
  | otherwise                    = (deg, leadingCoeff)
  where
    permutation  = [2 .. n] ++ [1]
    spray'       = permuteVariables permutation spray
    (powers, coeffs) = unzip (HM.toList spray')
    expnts           = map exponents powers
    constantTerm = fromMaybe AlgAdd.zero (HM.lookup (Powers S.empty 0) spray')
    (expnts', coeffs') = 
      unzip $ filter (\(s,_) -> not $ S.null s) (zip expnts coeffs)
    xpows = map (`index` 0) expnts'
    deg   = maximum xpows
    is    = elemIndices deg xpows
    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'')

-- | pseudo-division of two sprays, assuming degA >= degB >= 0
pseudoDivision :: (Eq a, AlgRing.C a)
  => Int                           -- ^ number of variables
  -> Spray a                       -- ^ A
  -> Spray a                       -- ^ B
  -> (Spray a, (Spray a, Spray a)) -- ^ (c, (Q, R)) such that c^*^A = B^*^Q ^+^ R
pseudoDivision n sprayA sprayB 
  | degB == minBound = error "pseudoDivision: pseudo-division by 0."
  | degA < degB      = error "pseudoDivision: degree(A) < degree(B)."
  | otherwise        = (ellB ^**^ delta , go sprayA zeroSpray delta)
  where
    degA         = degree n sprayA
    (degB, ellB) = degreeAndLeadingCoefficient n sprayB
    delta        = degA - degB + 1
    go sprayR sprayQ e = 
      if degR < degB || sprayR == zeroSpray
        then (q ^*^ sprayQ, q ^*^ sprayR)
        else go (ellB ^*^ sprayR ^-^ sprayS ^*^ sprayB) 
                (ellB ^*^ sprayQ ^+^ sprayS) 
                (e - 1)
      where
        (degR, ellR) = degreeAndLeadingCoefficient n sprayR
        q            = ellB ^**^ e
        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
gcdKX1dotsXn n sprayA sprayB
  | n == 0              = constantSpray $ gcdKX0 sprayA sprayB
  | degB > degA         = gcdKX1dotsXn n sprayB sprayA 
  | sprayB == zeroSpray = sprayA
  | otherwise           = go sprayA' sprayB' unitSpray unitSpray
  where
    gcdKX0 :: Spray a -> Spray a -> a
    gcdKX0 = const $ const AlgRing.one 
    n' = max (numberOfVariables sprayA) (numberOfVariables sprayB)
    degA = degree n' sprayA
    degB = degree n' sprayB
    gcdKX1dotsXm = gcdKX1dotsXn (n-1)
    content :: Spray a -> Spray a
    content spray = foldl1' gcdKX1dotsXm (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
    reduceSpray :: Spray a -> Spray a
    reduceSpray spray = exactDivisionBy cntnt spray 
      where
        coeffs = sprayCoefficients' n' spray
        cntnt  = foldl1' gcdKX1dotsXm coeffs
    contA   = content sprayA
    contB   = content sprayB
    d       = gcdKX1dotsXm contA contB 
    sprayA' = exactDivisionBy contA sprayA 
    sprayB' = exactDivisionBy contB sprayB 
    go :: Spray a -> Spray a -> Spray a -> Spray a -> Spray a
    go sprayA'' sprayB'' g h 
      | sprayR == zeroSpray           = d ^*^ reduceSpray sprayB''
      | numberOfVariables sprayR == 0 = d
      | otherwise = go sprayB'' 
                       (exactDivisionBy (g ^*^ h^**^delta) sprayR)
                       ellA''
                       (exactDivisionBy (h^**^delta) (h ^*^ g^**^delta))
        where
          (_, (_, sprayR)) = pseudoDivision n' sprayA'' sprayB''
          (degA'', ellA'') = degreeAndLeadingCoefficient n' sprayA''
          degB''           = degree n' sprayB'' 
          delta            = degA'' - degB''

-- | Greatest common divisor of two sprays with coefficients in a field
gcdSpray :: forall a. (Eq a, AlgField.C a) => Spray a -> Spray a -> Spray a
gcdSpray sprayA sprayB = gcdKX1dotsXn n sprayA sprayB 
  where
    n = max (numberOfVariables sprayA) (numberOfVariables sprayB)