hspray 0.5.1.0 → 0.5.2.0

raw patch · 5 files changed

+579/−28 lines, 5 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• CHANGELOG.md +19/−1
• README.md +2/−1
• hspray.cabal +2/−2
• src/Math/Algebra/Hspray.hs +405/−21
• tests/Main.hs +151/−3

@@ -279,4 +279,22 @@ ratio of sprays), and `dropVariables`, to drop a given number of leading 
 variables in a function-like object. The `dropVariables` functions is very 
 unsafe: if a variable is dropped while it is involved, the result can be an 
-invalid function-like object.+invalid function-like object.
+
+
+## 0.5.2.0 - 2024-05-05
+
+* Function `polynomialSubresultants`, to compute the polynomial subresultants 
+of two sprays (while the `subresultants` function computes the principal 
+subresultants).
+
+* Function `sturmHabichtSequence`, to compute the Sturm-Habicht sequence of a 
+spray.
+
+* Function `principalSturmHabichtSequence`, to compute the principal 
+Sturm-Habicht sequence of a spray.
+
+* Functions `numberOfRealRoots`, `numberOfRealRootsInOpenInterval` and 
+`numberOfRealRootsInClosedInterval`, to compute the total number of real 
+roots of a suitable univariate spray or its number of real roots in a given 
+interval. This can be very slow if the degree of the spray is not small.

@@ -552,5 +552,6 @@ ## Other features
 
 Other features offered by the package include: resultant and subresultants of 
-two polynomials, and greatest common divisor of two polynomials with 
+two polynomials, Sturm-Habicht sequence of a polynomial, number of real roots 
+of a univariate polynomial, and greatest common divisor of two polynomials with 
 coefficients in a field.

@@ -1,7 +1,7 @@ name:                hspray
-version:             0.5.1.0
+version:             0.5.2.0
 synopsis:            Multivariate polynomials and fractions of multivariate polynomials.
-description:         Manipulation of multivariate polynomials over a commutative ring and fractions of multivariate polynomials over a commutative field, Gröbner bases, resultant and subresultants, and greatest common divisor. It is possible to deal with multivariate polynomials whose coefficients are fractions of multivariate polynomials, and they can be interpreted as parametric polynomials with symbolic parameters.
+description:         Manipulation of multivariate polynomials over a commutative ring and fractions of multivariate polynomials over a commutative field, Gröbner bases, resultant, subresultants, Sturm-Habicht sequence, and greatest common divisor. It is possible to deal with multivariate polynomials whose coefficients are fractions of multivariate polynomials, and they can be interpreted as parametric polynomials with symbolic parameters.
 homepage:            https://github.com/stla/hspray#readme
 license:             GPL-3
 license-file:        LICENSE

@@ -5,7 +5,7 @@ License     : GPL-3
 Maintainer  : laurent_step@outlook.fr
 
-Deals with multivariate polynomials on a commutative ring. 
+Deals with multivariate polynomials over a commutative ring. 
 See README for examples.
 -}
 {-# LANGUAGE BangPatterns #-}
@@ -182,6 +182,7 @@   -- * Queries on a spray
   , getCoefficient
   , getConstantTerm
+  , isZeroSpray
   , isConstantSpray
   , isHomogeneousSpray
   , allExponents
@@ -191,9 +192,10 @@   , substituteSpray
   , composeSpray
   , evalSpraySpray
-  -- * Division of a spray
+  -- * Division of sprays
   , sprayDivision
   , sprayDivisionRemainder
+  , pseudoDivision
   -- * Gröbner basis
   , groebnerBasis
   , reduceGroebnerBasis
@@ -201,12 +203,23 @@   , esPolynomial
   , psPolynomial
   , isSymmetricSpray
-  -- * Resultant and subresultants
+  -- * Resultant, subresultants, and Sturm-Habicht sequence
   , resultant
   , resultant'
   , resultant1
   , subresultants
   , subresultants1
+  , polynomialSubresultants
+  , sturmHabichtSequence
+  , principalSturmHabichtSequence
+  -- * Number of real roots of a univariate spray. These functions can be very
+  -- slow if the degree of the spray is not small. 
+  , numberOfRealRoots
+  , numberOfRealRoots'
+  , numberOfRealRootsInOpenInterval
+  , numberOfRealRootsInOpenInterval'
+  , numberOfRealRootsInClosedInterval
+  , numberOfRealRootsInClosedInterval'
   -- * Greatest common divisor
   , gcdSpray
   -- * Matrices
@@ -214,17 +227,19 @@   , detLaplace'
   , characteristicPolynomial
   -- * Miscellaneous
+  , sumOfSprays
+  , productOfSprays
   , (.^)
   , (/>)
   , (/^)
   , fromList
   , toList
   , fromRationalSpray
-  , leadingTerm
   , isPolynomialOf
   , bombieriSpray
   , collinearSprays
   ) where
+import qualified Algebra.Absolute              as AlgAbs
 import qualified Algebra.Additive              as AlgAdd
 import qualified Algebra.Differential          as AlgDiff
 import qualified Algebra.Field                 as AlgField
@@ -247,7 +262,9 @@                                                 , foldl1'
                                                 , uncons
                                                 )
-import           Data.List.Extra                ( allSame )
+import           Data.List.Extra                ( allSame
+                                                , unsnoc 
+                                                )
 import           Data.Matrix                    ( Matrix 
                                                 , fromLists
                                                 , minorMatrix
@@ -696,7 +713,7 @@ -- | identify a `Polynomial a` to a `Spray a`, in order to apply the show spray 
 -- functions to the univariate polynomials
 polynomialToSpray :: forall a. (Eq a, AlgRing.C a) => Polynomial a -> Spray a
-polynomialToSpray pol = AlgAdd.sum terms
+polynomialToSpray pol = sumOfSprays terms
   where
     coeffs  = MathPol.coeffs pol
     indices = findIndices (/= A AlgAdd.zero) coeffs
@@ -705,7 +722,7 @@     terms = [get (coeffs !! i) *^ lone' 1 i | i <- indices]
 
 qPolynomialToQSpray :: QPolynomial -> QSpray
-qPolynomialToQSpray pol = AlgAdd.sum terms
+qPolynomialToQSpray pol = sumOfSprays terms
   where
     coeffs  = MathPol.coeffs pol
     indices = findIndices (/= A 0) coeffs
@@ -1348,7 +1365,7 @@ zeroSpray :: (Eq a, AlgAdd.C a) => Spray a
 zeroSpray = AlgAdd.zero
 
--- | whether the spray is zero
+-- | Whether a spray is the zero spray
 isZeroSpray :: Spray a -> Bool
 isZeroSpray = HM.null 
 
@@ -1868,6 +1885,14 @@ removeConstantTerm'' :: SafeSpray a -> SafeSpray a
 removeConstantTerm'' = HM.delete S.empty
 
+-- | Sum of sprays
+sumOfSprays :: (Eq a, AlgAdd.C a) => [Spray a] -> Spray a
+sumOfSprays = AlgAdd.sum
+
+-- | Product of sprays
+productOfSprays :: (Eq a, AlgRing.C a) => [Spray a] -> Spray a
+productOfSprays = AlgRing.product
+
 -- | ordered terms of a spray
 orderedTerms :: Spray a -> [Term a]
 orderedTerms spray = 
@@ -1919,7 +1944,7 @@ -- maxWithIndex :: Ord a => [a] -> (Int, a)
 -- maxWithIndex = maximumBy (comparing snd) . zip [0 .. ]
 
--- | Leading term of a spray 
+-- | leading term of a spray 
 leadingTerm :: Spray a -> Term a
 leadingTerm p = (biggest, p HM.! biggest) 
   where
@@ -1938,7 +1963,7 @@     expntsQ = exponents powsQ
     lower   = DF.all (uncurry (<=)) (S.zip expntsP expntsQ)
 
--- | quotient of term Q by term p, assuming P divides Q
+-- | quotient of term Q by term P, assuming P divides Q
 quotient :: AlgField.C a => Term a -> Term a -> Term a
 quotient (powsQ, coeffQ) (powsP, coeffP) = (pows, coeff)
   where
@@ -2363,7 +2388,7 @@     permutation = [2 .. n] ++ [1]
     deg = maximum xpows
     sprays = [
-      permuteVariables permutation (fromMaybe AlgAdd.zero (IM.lookup i imap'))
+      permuteVariables permutation (fromMaybe zeroSpray (IM.lookup i imap'))
       | i <- [deg, deg-1 .. 0]
       ]
 
@@ -2396,7 +2421,8 @@       where
         maxq = maximum qexpnts
 
--- | Subresultants of two /univariate/ sprays
+-- | Subresultants of two /univariate/ sprays. This function makes several calls
+-- to `detLaplace` and then it can be slow.
 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]
@@ -2445,7 +2471,9 @@     det = detLaplace $ 
           sylvesterMatrix (sprayCoefficients p') (sprayCoefficients q')
 
--- | Subresultants of two sprays
+-- | Subresultants of two sprays (the /principal/ subresultants, while the 
+-- `polynomialSubresultants` function returns the polynomial subresultants). 
+-- This function makes several calls to `detLaplace` and then it can be slow.
 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 
@@ -2523,7 +2551,7 @@           s' = if odd degp && odd degq 
             then AlgAdd.negate s 
             else s
-          (_, (_, r)) = pseudoDivision n p q
+          (_, (_, r)) = pseudoDivision' n p q
           p'             = q
           q'             = exactDivisionBy (g ^*^ h^**^delta) r
           (degp', ellp') = degreeAndLeadingCoefficient n p'
@@ -2532,11 +2560,357 @@           h'  = exactDivisionBy (h^**^delta) (h ^*^ g'^**^delta)
           h'' = exactDivisionBy (h'^**^degp') (h' ^*^ ellq'^**^degp')
 
+-- | Polynomial subresultants of two sprays (the `subresultants` function 
+-- computes the /principal/ subresultants).
+-- This function makes several calls to `detLaplace` and then it can be slow.
+polynomialSubresultants :: 
+    (Eq a, AlgRing.C a) 
+  => Int     -- ^ index of the variable with respect to which the subresultants will be computed (e.g. 2 for @y@)
+  -> Spray a 
+  -> Spray a 
+  -> [Spray a]
+polynomialSubresultants var p q 
+  | var < 1 || var > d               
+    = error "polynomialSubresultants: invalid variable index." 
+  | otherwise = [sres i | i <- [0 .. min n m]]
+  where
+    d = max (numberOfVariables p) (numberOfVariables q)
+    permutation  = [d-var+1 .. d] ++ [1 .. d-var]
+    permutation' = [var+1 .. d] ++ [1 .. var] 
+    p' = permuteVariables permutation p
+    q' = permuteVariables permutation q
+    pcoeffs = reverse $ sprayCoefficients' d p'
+    qcoeffs = reverse $ sprayCoefficients' d q'
+    pcoeff k = if k < 0 then zeroSpray else pcoeffs !! k
+    qcoeff k = if k < 0 then zeroSpray else qcoeffs !! k
+    n = length pcoeffs - 1
+    m = length qcoeffs - 1
+    x = loneTerm' d
+    prow i k = replicate k zeroSpray ++ 
+                 [pcoeff j | j <- [n, n-1 .. 2*i - m + k + 2]] ++ 
+                  [multSprayByTerm p' (x (m-i-1-k))]
+    prows i = [prow i k | k <- [0 .. m - i - 1]]
+    qrow i k = replicate k zeroSpray ++ 
+                 [qcoeff j | j <- [m, m-1 .. 2*i - n + k + 2]] ++
+                  [multSprayByTerm q' (x (n-i-1-k))]
+    qrows i = [qrow i k | k <- [0 .. n - i - 1]]
+    matrix i = fromLists (prows i ++ qrows i)
+    sres i = if n == m && i == n 
+      then q 
+      else permuteVariables permutation' $ detLaplace (matrix i)
 
+-- | Sturm-Habicht sequence of a spray. This function calls `polynomialSubresultants`
+-- and then it can be slow. 
+sturmHabichtSequence :: 
+    (Eq a, AlgRing.C a) 
+  => Int     -- ^ index of the variable with respect to which the Sturm-Habicht sequence will be computed (e.g. 2 for @y@)
+  -> Spray a 
+  -> [Spray a]
+sturmHabichtSequence var spray 
+  | var < 1 || var > d               
+    = error "sturmHabichtSequence: invalid variable index." 
+  | otherwise = signedSubresultants ++ [spray', spray]
+  where
+    d = numberOfVariables spray
+    spray' = derivative var spray
+    spraySubresultants = polynomialSubresultants var spray spray'
+    n = length spraySubresultants 
+    delta = [odd $ ((n - k - 1)*(n - k)) `div` 2 | k <- [0 .. n - 2]]
+    negateIf test = if test then AlgAdd.negate else id
+    signedSubresultants = zipWith negateIf delta spraySubresultants
+
+-- | Principal Sturm-Habicht sequence of a spray. This function calls `sturmHabicht` 
+-- sequence and then it can be slow.
+principalSturmHabichtSequence :: 
+    (Eq a, AlgRing.C a) 
+  => Int     -- ^ index of the variable with respect to which the Sturm-Habicht sequence will be computed (e.g. 2 for @y@)
+  -> Spray a 
+  -> [Spray a]
+principalSturmHabichtSequence var spray 
+  | var < 1 || var > d               
+    = error "principalSturmHabichtSequence: invalid variable index." 
+  | otherwise = [jcoeff j | j <- [0 .. length sHS - 1]]
+  where
+    d = numberOfVariables spray
+    sHS = sturmHabichtSequence var spray
+    permutation  = [d - var + 1 .. d] ++ [1 .. d - var]
+    permutation' = [var + 1 .. d] ++ [1 .. var] 
+    jcoeff j = if isZeroSpray sHSj || j >= length coeffs
+      then zeroSpray
+      else permuteVariables permutation' (coeffs !! j)
+      where
+        sHSj = sHS !! j
+        p = permuteVariables permutation sHSj
+        coeffs = reverse $ sprayCoefficients' d p -- what if numberOfVariables p /= d ?
+
+
+-- Number of real roots in interval -------------------------------------------
+
+_signVariations :: (a -> Char) -> [a] -> Int
+_signVariations signFunc as = v1 + v2 + 2*v3
+  where
+    count x xs = sum (map (fromEnum . (== x)) xs)
+    l = length as
+    signs = map signFunc as
+    chunks2 = [(signs !! i, signs !! (i+1)) | i <- [0 .. l-2]]
+    v1 = count ('+', '-') chunks2
+           + count ('-', '+') chunks2
+    chunks3 = [(signs !! i, signs !! (i+1), signs !! (i+2)) | i <- [0 .. l-3]]
+    chunks4 = [(signs !! i, signs !! (i+1), signs !! (i+2), signs !! (i+3)) 
+                | i <- [0 .. l-4]]
+    v2 = count ('-', '0', '+') chunks3 +
+           count ('+', '0', '-') chunks3 +
+             count ('+', '0', '0', '-') chunks4 +
+               count ('-', '0', '0', '+') chunks4
+    v3 = count ('+', '0', '0', '+') chunks4 + 
+           count ('-', '0', '0', '-') chunks4
+
+signVariations :: (Eq a, Num a) => [a] -> Int
+signVariations = _signVariations signFunc
+  where
+    signFunc a
+      | signum a == 0 = '0'
+      | signum a == 1 = '+'
+      | otherwise     = '-' 
+
+signVariations' :: (Eq a, AlgAbs.C a) => [a] -> Int
+signVariations' = _signVariations signFunc
+  where
+    signFunc a
+      | AlgAbs.signum a == AlgAdd.zero = '0'
+      | AlgAbs.signum a == AlgRing.one = '+'
+      | otherwise                      = '-' 
+
+_numberOfRealRootsInOpenInterval :: 
+  (AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> Maybe a -> Maybe a -> Int
+_numberOfRealRootsInOpenInterval signVariationsFunc spray alpha beta 
+  | isConstantSpray spray =
+    if isZeroSpray spray 
+      then error "numberOfRealRootsInOpenInterval: the spray is null."
+      else 0
+  | alphaIsJust && betaIsJust = 
+    if alpha' == beta'
+      then 0
+      else
+        if alpha > beta
+          then 
+            error "numberOfRealRootsInOpenInterval: the bounds are not ordered."
+          else 
+            if isZeroAtBeta then svDiff - 1 else svDiff
+  | alphaIsJust = svDiff
+  | betaIsJust  = if isZeroAtBeta then svDiff - 1 else svDiff
+  | otherwise   = svDiff 
+  where
+    alphaIsJust = isJust alpha
+    betaIsJust  = isJust beta
+    alpha' = if alphaIsJust then fromJust alpha else undefined
+    beta'  = if betaIsJust then fromJust beta else undefined
+    evaluateAtMinusInfinity p = if isConstantSpray p 
+      then getConstantTerm p
+      else let 
+        (deg, lCoeff) = first ((`S.index` 0) . exponents) (leadingTerm p) in
+        if even deg then lCoeff else AlgAdd.negate lCoeff
+    evaluateAtPlusInfinity p = if isConstantSpray p 
+      then getConstantTerm p
+      else (snd . leadingTerm) p
+    (ginit, glast) = 
+      fromJust $ unsnoc $ filter (not . isZeroSpray) (sturmHabichtSequence 1 spray)
+    eval1 = 
+      if alphaIsJust then evaluateAt [alpha'] else evaluateAtMinusInfinity
+    eval2 = 
+      if betaIsJust then evaluateAt [beta'] else evaluateAtPlusInfinity
+    sprayAtAlpha = eval1 glast
+    sprayAtBeta  = eval2 glast
+    isZeroAtBeta = sprayAtBeta == AlgAdd.zero
+    galpha = map eval1 ginit ++ [sprayAtAlpha]
+    gbeta  = map eval2 ginit ++ [sprayAtBeta]
+    svAtAlpha = signVariationsFunc galpha
+    svAtBeta  = signVariationsFunc gbeta
+    svDiff    = svAtAlpha - svAtBeta
+
+_numberOfRealRootsInClosedInterval :: 
+   (AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> Maybe a -> Maybe a -> Int
+_numberOfRealRootsInClosedInterval signVariationsFunc spray alpha beta = 
+  _numberOfRealRootsInOpenInterval signVariationsFunc spray alpha beta + toAdd
+  where
+    alphaIsJust = isJust alpha
+    betaIsJust  = isJust beta
+    alpha' = if alphaIsJust then fromJust alpha else undefined
+    beta'  = if betaIsJust then fromJust beta else undefined
+    isZeroAtAlpha = alphaIsJust && evaluateAt [alpha'] spray == AlgAdd.zero
+    isZeroAtBeta  = betaIsJust && evaluateAt [beta'] spray == AlgAdd.zero
+    toAdd = 
+      if alpha == beta
+        then fromEnum isZeroAtAlpha
+        else fromEnum isZeroAtAlpha + fromEnum isZeroAtBeta
+
+-- | Number of real roots of a spray in an open interval (that makes sense 
+-- only for a spray on a ring embeddable in the real numbers).
+numberOfRealRootsInOpenInterval :: 
+  (Num a, AlgRing.C a, Ord a) 
+  => Spray a -- ^ a spray
+  -> Maybe a -- ^ lower bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity
+  -> Maybe a -- ^ upper bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity
+  -> Int
+numberOfRealRootsInOpenInterval spray = 
+  if isUnivariate spray 
+    then _numberOfRealRootsInOpenInterval signVariations spray
+    else error "numberOfRealRootsInOpenInterval: the spray is not univariate."
+
+-- | Number of real roots of a spray in a closed interval (that makes sense 
+-- only for a spray on a ring embeddable in the real numbers). The roots are 
+-- not counted with their multiplicity.
+numberOfRealRootsInClosedInterval :: 
+  (Num a, AlgRing.C a, Ord a) 
+  => Spray a -- ^ a spray
+  -> Maybe a -- ^ lower bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity
+  -> Maybe a -- ^ upper bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity
+  -> Int
+numberOfRealRootsInClosedInterval spray = 
+  if isUnivariate spray 
+    then _numberOfRealRootsInClosedInterval signVariations spray
+    else error "numberOfRealRootsInClosedInterval: the spray is not univariate."
+ 
+-- | Number of real roots of a spray in an open interval (that makes sense 
+-- only for a spray on a ring embeddable in the real numbers).
+numberOfRealRootsInOpenInterval' :: 
+  (AlgAbs.C a, Ord a) 
+  => Spray a -- ^ a spray
+  -> Maybe a -- ^ lower bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity
+  -> Maybe a -- ^ upper bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity
+  -> Int
+numberOfRealRootsInOpenInterval' spray =
+  if isUnivariate spray 
+    then _numberOfRealRootsInOpenInterval signVariations' spray
+    else error "numberOfRealRootsInOpenInterval': the spray is not univariate."
+
+-- | Number of real roots of a spray in a closed interval (that makes sense 
+-- only for a spray on a ring embeddable in the real numbers). The roots are 
+-- not counted with their multiplicity.
+numberOfRealRootsInClosedInterval' :: 
+  (AlgAbs.C a, Ord a) 
+  => Spray a -- ^ a spray
+  -> Maybe a -- ^ lower bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity
+  -> Maybe a -- ^ upper bound of the interval; use @Just@ for a finite bound, and @Nothing@ for minus infinity
+  -> Int
+numberOfRealRootsInClosedInterval' spray =
+  if isUnivariate spray 
+    then _numberOfRealRootsInClosedInterval signVariations' spray
+    else error "numberOfRealRootsInClosedInterval': the spray is not univariate."
+
+
+-- Number of real roots -------------------------------------------------------
+
+runLengthEncoding :: Eq a => [a] -> [(a,Int)]
+runLengthEncoding = foldr code []
+  where 
+    code c []         = [(c,1)]
+    code c ((x,n):ts) 
+      | c == x        = (x,n+1):ts
+      | otherwise     = (c,1):(x,n):ts
+
+_signPermanencesAndVariations :: (a -> Char) -> [a] -> (Int, Int)
+_signPermanencesAndVariations signFunc as = (permanences, variations)
+  where
+    signs = map signFunc as
+    rle = runLengthEncoding signs
+    l = length rle
+    lengths = map snd rle
+    permanences = sum lengths - l
+    variations = l - 1
+
+signPermanencesAndVariations :: (Eq a, Num a) => [a] -> (Int, Int)
+signPermanencesAndVariations = _signPermanencesAndVariations signFunc
+  where
+    signFunc a = if signum a == 1 then '+' else '-'
+
+signPermanencesAndVariations' :: (Eq a, AlgAbs.C a) => [a] -> (Int, Int)
+signPermanencesAndVariations' = _signPermanencesAndVariations signFunc
+  where
+    signFunc a = if AlgAbs.signum a == AlgRing.one then '+' else '-'
+
+distributionsOfZeros :: (Eq a, AlgAdd.C a) => [a] -> ([Int], [Int], [Int])
+distributionsOfZeros as = (i_, k_, ik_)
+  where
+    symbol a = if a == AlgAdd.zero then '0' else '*'
+    symbols = map symbol as ++ ['0']
+    lengths = map snd (runLengthEncoding symbols)
+    l = length lengths
+    cumulativeLengths = scanl (+) 0 lengths
+    range = [0 .. l-1]
+    ik_ = [cumulativeLengths !! i | i <- range, even i]
+    i_ = [cumulativeLengths !! i - 1 | i <- range, odd i]
+    k_ = [lengths !! i | i <- [0 .. l-2], odd i] ++ [lengths !! (l-1) - 1]
+
+_blocksAndEpsilons :: (Eq a, AlgAdd.C a) => (a -> Int) -> [a] -> ([[a]], [Int])
+_blocksAndEpsilons signFunc as = (blocks, epsilons)
+  where
+    (i_, k_, ik_) = distributionsOfZeros as
+    t = length i_
+    blocks = [[as !! m | m <- [ik_ !! n .. i_ !! n]] | n <- [0 .. t-1]]
+    epsilon s = let ks = k_ !! s in if odd ks 
+      then
+        0
+      else
+        (if even (ks `div` 2) then 1 else -1) * 
+          signFunc (as !! (ik_ !! (s+1))) * signFunc (as !! (i_ !! s))
+    epsilons = [epsilon n | n <- [0 .. t-2]]
+
+blocksAndEpsilons :: (Eq a, AlgAdd.C a, Num a) => [a] -> ([[a]], [Int])
+blocksAndEpsilons = _blocksAndEpsilons signFunc 
+  where
+    signFunc a = if signum a == 1 then 1 else -1
+
+blocksAndEpsilons' :: (Eq a, AlgAbs.C a) => [a] -> ([[a]], [Int])
+blocksAndEpsilons' = _blocksAndEpsilons signFunc 
+  where
+    signFunc a = if AlgAbs.signum a == AlgRing.one then 1 else -1
+
+_C :: ([a] -> (Int, Int)) -> ([[a]], [Int]) -> Int
+_C signPermanencesAndVariationsFunc (blocks, epsilons) = sum pvs - sum epsilons
+  where
+    permanencesMinusVariations as = 
+      let (p, v) = signPermanencesAndVariationsFunc as in p - v 
+    pvs = map permanencesMinusVariations blocks
+
+-- | Number of real roots of a spray (that makes sense only for a spray on a 
+-- ring embeddable in the real numbers). The roots are not counted with their 
+-- multiplicity.
+numberOfRealRoots :: (Eq a, AlgRing.C a, Num a) => Spray a -> Int
+numberOfRealRoots spray = 
+  if isUnivariate spray 
+    then 
+      if isConstantSpray spray 
+        then if isZeroSpray spray
+          then error "numberOfRealRoots: the spray is null."
+          else 0
+        else _C signPermanencesAndVariations (blocksAndEpsilons $ reverse as)
+    else 
+      error "numberOfRealRoots: the spray is not univariate."
+  where
+    as = map getConstantTerm (principalSturmHabichtSequence 1 spray)
+
+-- | Number of real roots of a spray (that makes sense only for a spray on a 
+-- ring embeddable in the real numbers). The roots are not counted with their 
+-- multiplicity.
+numberOfRealRoots' :: (Eq a, AlgAbs.C a) => Spray a -> Int
+numberOfRealRoots' spray = 
+  if isUnivariate spray 
+    then 
+      if isConstantSpray spray 
+        then if isZeroSpray spray
+          then error "numberOfRealRoots: the spray is null."
+          else 0
+        else _C signPermanencesAndVariations' (blocksAndEpsilons' $ reverse as)
+    else 
+      error "numberOfRealRoots': the spray is not univariate."
+  where
+    as = map getConstantTerm (principalSturmHabichtSequence 1 spray)
+
+
 -- GCD stuff ------------------------------------------------------------------
 
--- | the coefficients of a spray as a univariate spray in x_n with 
--- spray coefficients
+-- | the spray coefficients of a spray as a univariate spray in x_n 
 sprayCoefficients' :: (Eq a, AlgRing.C a) => Int -> Spray a -> [Spray a]
 sprayCoefficients' n spray 
   | numberOfVariables spray /= n = [spray]
@@ -2555,7 +2929,7 @@     imap'  = IM.insertWith (^+^) 0 (constantSpray constantTerm) imap
     deg    = maximum xpows
     sprays = [
-        fromMaybe AlgAdd.zero (IM.lookup i imap')
+        fromMaybe zeroSpray (IM.lookup i imap')
         | i <- [deg, deg-1 .. 0]
       ]
 
@@ -2601,13 +2975,23 @@     coeffs'' = [coeffs' !! i | i <- is]
     leadingCoeff = sumTerms (zip powers'' coeffs'')
 
--- | Pseudo-division of two sprays, assuming degA >= degB >= 0
+-- | Pseudo-division of two sprays @A@ and @B@ such that @deg(A) >= deg(B)@
+-- where @deg@ is the degree with respect to the outermost variable.
 pseudoDivision :: (Eq a, AlgRing.C a)
+  => Spray a                       -- ^ A
+  -> Spray a                       -- ^ B
+  -> (Spray a, (Spray a, Spray a)) -- ^ (C, (Q, R)) such that C^*^A = B^*^Q ^+^ R
+pseudoDivision sprayA sprayB = 
+  pseudoDivision' (max (numberOfVariables sprayA) (numberOfVariables sprayB))
+                    sprayA sprayB
+
+-- | 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 
+  -> (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)
@@ -2679,7 +3063,7 @@                        ellB'' ---- ellA''
                        (exactDivisionBy (h^**^delta) (h ^*^ ellB''^**^delta)) --- g^**^delta
         where
-          (_, (_, sprayR)) = pseudoDivision n' sprayA'' sprayB''
+          (_, (_, sprayR)) = pseudoDivision' n' sprayA'' sprayB''
           -- (degA'', ellA'') = degreeAndLeadingCoefficient n' sprayA''
           -- degB''           = degree n' sprayB''
           degA''           = degree n' sprayA''

@@ -1,3 +1,4 @@+{-# LANGUAGE TupleSections #-}
 module Main (main) where
 import qualified Algebra.Additive               as AlgAdd
 import qualified Algebra.Module                 as AlgMod
@@ -43,6 +44,12 @@                                                   subresultants,
                                                   resultant1,
                                                   subresultants1,
+                                                  polynomialSubresultants,
+                                                  sturmHabichtSequence,
+                                                  principalSturmHabichtSequence,
+                                                  numberOfRealRootsInOpenInterval,
+                                                  numberOfRealRootsInClosedInterval,
+                                                  numberOfRealRoots, 
                                                   sprayDivision,
                                                   gcdSpray,
                                                   QSpray',
@@ -96,7 +103,8 @@                                                   qmonomial,
                                                   isHomogeneousSpray,
                                                   psPolynomial,
-                                                  prettyQSprayXYZ
+                                                  prettyQSprayXYZ,
+                                                  productOfSprays
                                                 )
 import           MathObj.Matrix                 ( fromRows )
 import qualified MathObj.Matrix                 as MathMatrix
@@ -804,9 +812,149 @@         test1 = rx == rx'
         test2 = ry == ry'
         test3 = rz == rz'
-      assertBool "" (test1 && test2 && test3),
+      assertBool "" (test1 && test2 && test3)
 
-    testCase "gcdSpray - univariate example" $ do
+    , testCase "polynomialSubresultants" $ do
+      let
+        x = qlone 1
+        y = qlone 2
+        p = x^**^3^*^y^**^3 ^+^ x^**^2 ^+^ 3*^x ^+^ unitSpray
+        q = x^**^4 ^+^ 2*^x^**^3 ^-^ x^*^y^**^4
+        srs1 = map prettyQSpray (polynomialSubresultants 1 p q)
+        srs2 = map prettyQSpray (polynomialSubresultants 2 p q)
+      assertEqual "" 
+        (srs1, srs2)
+        (
+          [
+            "y^21 + 15*y^14 - 11*y^11 + y^8 + 37*y^7 + 4*y^4 - 8*y^3 - 1"
+          , "x.y^17 + 14*x.y^10 - 11*x.y^7 + x.y^4 + 39*x.y^3 + 5*x + 2*y^10 - y^7 + 14*y^3 + 2"
+          , "-5*x^2.y^3 + x^2 - x.y^10 - 7*x.y^3 + 3*x - 2*y^3 + 1"
+          , "x^3.y^3 + x^2 + 3*x + 1"
+          ],
+          [
+            "x^24 + 6*x^23 + 12*x^22 + 8*x^21 - x^11 - 12*x^10 - 58*x^9 - 144*x^8 - 195*x^7 - 144*x^6 - 58*x^5 - 12*x^4 - x^3"
+          , "x^13 + 5*x^12 + 7*x^11 + 2*x^10 + x^9.y + 6*x^8.y + 11*x^7.y + 6*x^6.y + x^5.y"
+          , "x^10 + 2*x^9 + x^6.y + 3*x^5.y + x^4.y"
+          , "x^3.y^3 + x^2 + 3*x + 1"
+          ]
+        )
+
+    , testCase "principalSturmHabichtSequence" $ do
+      let
+        x = qlone 1
+        y = qlone 2
+        p = x^**^4 ^*^ y ^+^ y^**^2
+        psh1 = principalSturmHabichtSequence 1 p
+        psh2 = principalSturmHabichtSequence 2 p
+      assertEqual "" 
+        (psh1, psh2)
+        (
+          [
+              256 *^ y^**^10
+            , zeroSpray
+            , zeroSpray
+            , 4 *^ y
+            , y
+          ],
+          [
+              x^**^8
+            , constantSpray 2
+            , constantSpray 1
+          ]
+        )
+
+    , testCase "sturmHabichtSequence & principalSturmHabichtSequence" $ do
+      let
+        x = qlone 1
+        y = qlone 2
+        p = x^**^3^*^y^**^4 ^+^ x^**^2^*^y^**^2 ^+^ 3*^y ^-^ constantSpray 6
+        sh1 = map prettyQSpray (sturmHabichtSequence 1 p)
+        sh2 = map prettyQSpray (sturmHabichtSequence 2 p)
+        psh1 = map prettyQSpray (principalSturmHabichtSequence 1 p)
+        psh2 = map prettyQSpray (principalSturmHabichtSequence 2 p)
+      assertEqual "" 
+        ((sh1, psh1), (sh2, psh2))
+        (
+          (
+            [
+              "-243*y^14 + 972*y^13 - 972*y^12 - 12*y^11 + 24*y^10"
+            , "2*x.y^8 - 27*y^9 + 54*y^8"
+            , "3*x^2.y^4 + 2*x.y^2"
+            , "x^3.y^4 + x^2.y^2 + 3*y - 6"
+            ],
+            [
+              "-243*y^14 + 972*y^13 - 972*y^12 - 12*y^11 + 24*y^10"
+            , "2*y^8"
+            , "3*y^4"
+            , "y^4"
+            ]
+          ),
+          (
+            [
+              "-96*x^14 - 4608*x^13 - 55332*x^12 - 7776*x^11 - 2187*x^9"
+            , "-8*x^12.y - 192*x^11.y - 12*x^10 - 324*x^9.y + 864*x^9"
+            , "-8*x^8.y^2 - 36*x^6.y + 96*x^6"
+            , "4*x^3.y^3 + 2*x^2.y + 3"
+            , "x^3.y^4 + x^2.y^2 + 3*y - 6"
+            ],
+            [
+              "-96*x^14 - 4608*x^13 - 55332*x^12 - 7776*x^11 - 2187*x^9"
+            , "-8*x^12 - 192*x^11 - 324*x^9"
+            , "-8*x^8"
+            , "4*x^3"
+            , "x^3"
+            ]
+          )
+        )
+
+    , testCase "number of real roots in interval" $ do
+      let
+        x = qlone 1
+        factors = [x ^-^ constantSpray (toRational i) | i <- [1::Int .. 5]]
+        spray = productOfSprays factors
+        intervals = [
+            (Just 0, Just 9)
+          , (Just 1, Just 6)
+          , (Just 2, Just 3)
+          , (Just 0, Just 4)
+          , (Just $ 2 + (1%4), Just $ 3 - (1%4))
+          ]
+        nroots = 
+          map (uncurry (numberOfRealRootsInClosedInterval spray)) intervals
+        nroots' = 
+          map (uncurry (numberOfRealRootsInOpenInterval spray)) intervals
+        rightUnboundedIntervals = map ((, Nothing) . Just) [0, 1, 2, 2 + (1%4), 5]
+        nroots'' = map (uncurry (numberOfRealRootsInClosedInterval spray)) 
+                        rightUnboundedIntervals
+        nroots''' = map (uncurry (numberOfRealRootsInOpenInterval spray)) 
+                          rightUnboundedIntervals
+        leftUnboundedIntervals = map ((Nothing, ) . Just) [0, 1, 2, 2 + (1%4), 5]
+        nroots'''' = map (uncurry (numberOfRealRootsInClosedInterval spray)) 
+                        leftUnboundedIntervals
+        nroots''''' = map (uncurry (numberOfRealRootsInOpenInterval spray)) 
+                          leftUnboundedIntervals
+        ntotalroots = numberOfRealRootsInClosedInterval spray Nothing Nothing
+      assertEqual "" 
+        ( nroots, nroots'
+        , nroots'', nroots'''
+        , nroots'''', nroots'''''
+        , ntotalroots
+        ) 
+        ( [5, 5, 2, 4, 0], [5, 4, 0, 3, 0]
+        , [5, 5, 4, 3, 1], [5, 4, 3, 3, 0]
+        , [0, 1, 2, 2, 5], [0, 0, 1, 2, 4]
+        , 5
+        )
+
+    , testCase "number of real roots" $ do
+      let
+        x = qlone 1
+        factors = [x ^-^ constantSpray (toRational i) | i <- [1::Int .. 5]]
+        spray = productOfSprays factors
+        nroots = map numberOfRealRoots [spray, x^**^3, x^**^2 <+ 1]
+      assertEqual "" nroots [5, 1, 0]
+
+    , testCase "gcdSpray - univariate example" $ do
       let
         x = lone 1 :: Spray Rational
         sprayD = x^**^2 ^+^ unitSpray