numberOfRealRootsInInterval for semi-infinite intervals

stla · May 4, 2024 · d553695 · d553695
1 parent 4174371
commit d553695
Show file tree

Hide file tree

Showing 3 changed files with 107 additions and 75 deletions.
diff --git a/scripts/signAnalysis2.hs b/scripts/signAnalysis2.hs
@@ -2,7 +2,7 @@ import qualified Algebra.Absolute as AlgAbs
 import qualified Algebra.Additive as AlgAdd
 import qualified Algebra.Ring     as AlgRing
 import Data.List.Extra (unsnoc)
-import Data.Maybe (fromJust)
+import Data.Maybe (fromJust, isJust, isNothing)
 import Math.Algebra.Hspray
 
 runLengthEncoding :: Eq a => [a] -> [(a,Int)]
@@ -33,7 +33,7 @@ signPermanencesAndVariations' = _signPermanencesAndVariations signFunc
   where
     signFunc a = if AlgAbs.signum a == AlgRing.one then '+' else '-'
 
-{- _signVariations :: Eq a => (a -> Char) -> [a] -> Int
+_signVariations :: Eq a => (a -> Char) -> [a] -> Int
 _signVariations signFunc as = v1 + v2 + 2*v3
   where
     count x xs = sum (map (fromEnum . (== x)) xs)
@@ -69,58 +69,74 @@ signVariations' = _signVariations signFunc
       | otherwise                      = '-' 
 
 _numberOfRealRootsInOpenInterval :: 
-  (Eq a, AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> (a, a) -> Int
-_numberOfRealRootsInOpenInterval signVariationsFunc spray (alpha, beta) 
-  | alpha == beta = 0
-  | isConstantSpray spray = if isZeroSpray spray 
-    then error "numberOfRealRoots: the spray is null."
-    else 0
-  | otherwise = if sprayAtBeta == AlgAdd.zero then svDiff - 1 else svDiff
+  (Eq a, AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> a -> Maybe a -> Int
+_numberOfRealRootsInOpenInterval signVariationsFunc spray alpha beta =
+  if isConstantSpray spray 
+    then 
+      if isZeroSpray spray 
+        then error "numberOfRealRootsInInterval: the spray is null."
+        else 0
+    else 
+      if betaIsJust
+        then 
+          if alpha == beta'
+            then 0
+            else
+              if alpha > beta' 
+                then 
+                  error "numberOfRealRootsInInterval: the bounds are not ordered."
+                else 
+                  if isZeroAtBeta then svDiff - 1 else svDiff
+        else 
+          svAtAlpha
   where
-    (alpha', beta') = if alpha < beta then (alpha, beta) else (beta, alpha)
+    betaIsJust = isJust beta
+    beta' = if betaIsJust then fromJust beta else undefined
     (ginit, glast) = 
       fromJust $ unsnoc $ filter (not . isZeroSpray) (sturmHabichtSequence 1 spray)
     sprayAtAlpha = evaluateAt [alpha] glast
-    sprayAtBeta  = evaluateAt [beta] glast
+    sprayAtBeta  = evaluateAt [beta'] glast
+    isZeroAtBeta = sprayAtBeta == AlgAdd.zero
     galpha = map (evaluateAt [alpha]) ginit ++ [sprayAtAlpha]
-    gbeta  = map (evaluateAt [beta]) ginit ++ [sprayAtBeta]
-    svalpha = signVariationsFunc galpha
-    svbeta  = signVariationsFunc gbeta
-    svDiff  = svalpha - svbeta
+    gbeta  = map (evaluateAt [beta']) ginit ++ [sprayAtBeta]
+    svAtAlpha = signVariationsFunc galpha
+    svAtBeta  = signVariationsFunc gbeta
+    svDiff    = svAtAlpha - svAtBeta
 
 _numberOfRealRootsInClosedInterval :: 
-  (Eq a, AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> (a, a) -> Int
-_numberOfRealRootsInClosedInterval signVariationsFunc spray (alpha, beta) = 
-  if alpha == beta 
-    then 
-      fromEnum (sprayAtAlpha == AlgAdd.zero)
-    else 
-      _numberOfRealRootsInOpenInterval signVariationsFunc spray (alpha, beta) +
-        fromEnum (sprayAtAlpha == AlgAdd.zero) + 
-          fromEnum (sprayAtBeta == AlgAdd.zero)
+  (Eq a, AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> a -> Maybe a -> Int
+_numberOfRealRootsInClosedInterval signVariationsFunc spray alpha beta = 
+  _numberOfRealRootsInOpenInterval signVariationsFunc spray alpha beta + toAdd
   where
-    sprayAtAlpha = evaluateAt [alpha] spray
-    sprayAtBeta  = evaluateAt [beta] spray
- 
+    betaIsJust = isJust beta
+    beta' = if betaIsJust then fromJust beta else undefined
+    isZeroAtAlpha = evaluateAt [alpha] spray == AlgAdd.zero
+    isZeroAtBeta  = evaluateAt [beta'] spray == AlgAdd.zero
+    toAdd = if betaIsJust
+      then if alpha == beta' 
+        then fromEnum isZeroAtAlpha
+        else fromEnum isZeroAtAlpha + fromEnum isZeroAtBeta
+      else fromEnum isZeroAtAlpha
+
 -- | 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 :: 
-  (Eq a, Num a, AlgRing.C a, Ord a) => Spray a -> (a, a) -> Int
-numberOfRealRootsInOpenInterval spray = 
+{- xnumberOfRealRootsInOpenInterval :: 
+  (Eq a, Num a, AlgRing.C a, Ord a) => Spray a -> (a, Maybe a) -> Int
+xnumberOfRealRootsInOpenInterval spray (alpha, beta) = 
   if isUnivariate spray 
-    then _numberOfRealRootsInOpenInterval signVariations spray
+    then _numberOfRealRootsInOpenInterval signVariations spray (alpha, beta)
     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).
-numberOfRealRootsInClosedInterval :: 
+{- xnumberOfRealRootsInClosedInterval :: 
   (Eq a, Num a, AlgRing.C a, Ord a) => Spray a -> (a, a) -> 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 
+ -} 
+{- -- | 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' :: 
   (Eq a, AlgAbs.C a, Ord a) => Spray a -> (a, a) -> Int
@@ -137,7 +153,6 @@ numberOfRealRootsInClosedInterval' spray =
   if isUnivariate spray 
     then _numberOfRealRootsInClosedInterval signVariations' spray
     else error "numberOfRealRootsInClosedInterval': the spray is not univariate."
-
  -}
 x = qlone 1
 factors = [x ^-^ constantSpray (toRational i) | i <- [1::Int .. 3]]
@@ -146,7 +161,7 @@ spray = AlgRing.product factors
 test = map (numberOfRealRootsInClosedInterval spray) [(0, 10), (0, 5), (2, 3)]
 
 
-distributionsOfZeros :: (Eq a, AlgAdd.C a) => [a] -> ([Int], [Int], [Int])
+{- 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 '*'
@@ -221,3 +236,4 @@ numberOfRealRoots' spray =
     as = map getConstantTerm (principalSturmHabichtSequence 1 spray)
 
 spray' = x^**^2 ^+^ unitSpray
+ -}
diff --git a/src/Math/Algebra/Hspray.hs b/src/Math/Algebra/Hspray.hs
@@ -2682,45 +2682,55 @@ signVariations' = _signVariations signFunc
       | otherwise                      = '-' 
 
 _numberOfRealRootsInOpenInterval :: 
-  (AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> (a, a) -> Int
-_numberOfRealRootsInOpenInterval signVariationsFunc spray (alpha, beta) 
-  | alpha == beta         = 0
-  | isConstantSpray spray = if isZeroSpray spray 
-    then error "numberOfRealRoots: the spray is null (hence the number of real roots is infinite)."
-    else 0
-  | otherwise             = if sprayAtBeta == AlgAdd.zero 
-      then svDiff - 1 
-      else svDiff
-  where
-    (alpha', beta') = if alpha < beta then (alpha, beta) else (beta, alpha)
-    (ginit, glast) = fromJust $ 
-      unsnoc $ filter (not . isZeroSpray) (sturmHabichtSequence 1 spray)
-    sprayAtAlpha = evaluateAt [alpha'] glast
+  (AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> a -> Maybe a -> Int
+_numberOfRealRootsInOpenInterval signVariationsFunc spray alpha beta 
+  | isConstantSpray spray =
+    if isZeroSpray spray 
+      then error "numberOfRealRootsInInterval: the spray is null."
+      else 0
+  | betaIsJust = 
+    if alpha == beta'
+      then 0
+      else
+        if alpha > beta' 
+          then 
+            error "numberOfRealRootsInInterval: the bounds are not ordered."
+          else 
+            if isZeroAtBeta then svDiff - 1 else svDiff
+  | otherwise = svAtAlpha
+  where
+    betaIsJust = isJust beta
+    beta' = if betaIsJust then fromJust beta else undefined
+    (ginit, glast) = 
+      fromJust $ unsnoc $ filter (not . isZeroSpray) (sturmHabichtSequence 1 spray)
+    sprayAtAlpha = evaluateAt [alpha] glast
     sprayAtBeta  = evaluateAt [beta'] glast
-    galpha = map (evaluateAt [alpha']) ginit ++ [sprayAtAlpha]
+    isZeroAtBeta = sprayAtBeta == AlgAdd.zero
+    galpha = map (evaluateAt [alpha]) ginit ++ [sprayAtAlpha]
     gbeta  = map (evaluateAt [beta']) ginit ++ [sprayAtBeta]
-    svalpha = signVariationsFunc galpha
-    svbeta  = signVariationsFunc gbeta
-    svDiff  = svalpha - svbeta
+    svAtAlpha = signVariationsFunc galpha
+    svAtBeta  = signVariationsFunc gbeta
+    svDiff    = svAtAlpha - svAtBeta
 
 _numberOfRealRootsInClosedInterval :: 
-  (AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> (a, a) -> Int
-_numberOfRealRootsInClosedInterval signVariationsFunc spray (alpha, beta) = 
-  if alpha == beta 
-    then 
-      fromEnum (sprayAtAlpha == AlgAdd.zero)
-    else 
-      _numberOfRealRootsInOpenInterval signVariationsFunc spray (alpha, beta) +
-        fromEnum (sprayAtAlpha == AlgAdd.zero) + 
-          fromEnum (sprayAtBeta == AlgAdd.zero)
-  where
-    sprayAtAlpha = evaluateAt [alpha] spray
-    sprayAtBeta  = evaluateAt [beta] spray
+  (AlgRing.C a, Ord a) => ([a] -> Int) -> Spray a -> a -> Maybe a -> Int
+_numberOfRealRootsInClosedInterval signVariationsFunc spray alpha beta = 
+  _numberOfRealRootsInOpenInterval signVariationsFunc spray alpha beta + toAdd
+  where
+    betaIsJust = isJust beta
+    beta' = if betaIsJust then fromJust beta else undefined
+    isZeroAtAlpha = evaluateAt [alpha] spray == AlgAdd.zero
+    isZeroAtBeta  = evaluateAt [beta'] spray == AlgAdd.zero
+    toAdd = if betaIsJust
+      then if alpha == beta' 
+        then fromEnum isZeroAtAlpha
+        else fromEnum isZeroAtAlpha + fromEnum isZeroAtBeta
+      else fromEnum isZeroAtAlpha
 
 -- | 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, a) -> Int
+  (Num a, AlgRing.C a, Ord a) => Spray a -> a -> Maybe a -> Int
 numberOfRealRootsInOpenInterval spray = 
   if isUnivariate spray 
     then _numberOfRealRootsInOpenInterval signVariations spray
@@ -2730,7 +2740,7 @@ numberOfRealRootsInOpenInterval spray =
 -- 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, a) -> Int
+  (Num a, AlgRing.C a, Ord a) => Spray a -> a -> Maybe a -> Int
 numberOfRealRootsInClosedInterval spray = 
   if isUnivariate spray 
     then _numberOfRealRootsInClosedInterval signVariations spray
@@ -2739,7 +2749,7 @@ numberOfRealRootsInClosedInterval spray =
 -- | 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, a) -> Int
+  (AlgAbs.C a, Ord a) => Spray a -> a -> Maybe a -> Int
 numberOfRealRootsInOpenInterval' spray =
   if isUnivariate spray 
     then _numberOfRealRootsInOpenInterval signVariations' spray
@@ -2749,7 +2759,7 @@ numberOfRealRootsInOpenInterval' spray =
 -- 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, a) -> Int
+  (AlgAbs.C a, Ord a) => Spray a -> a -> Maybe a -> Int
 numberOfRealRootsInClosedInterval' spray =
   if isUnivariate spray 
     then _numberOfRealRootsInClosedInterval signVariations' spray

diff --git a/tests/Main.hs b/tests/Main.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE TupleSections #-}
 module Main (main) where
 import qualified Algebra.Additive               as AlgAdd
 import qualified Algebra.Module                 as AlgMod
@@ -887,10 +888,15 @@ main = defaultMain $ testGroup
         x = qlone 1
         factors = [x ^-^ constantSpray (toRational i) | i <- [1::Int .. 5]]
         spray = productOfSprays factors
-        intervals = [(0, 9), (1, 6), (2, 3), (0, 4), (2 + (1%4), 3 - (1%4))]
-        nroots = map (numberOfRealRootsInClosedInterval spray) intervals
-        nroots' = map (numberOfRealRootsInOpenInterval spray) intervals
-      assertEqual "" (nroots, nroots') ([5, 5, 2, 4, 0], [5, 4, 0, 3, 0])
+        intervals = [(0, Just 9), (1, Just 6), (2, Just 3), (0, Just 4), (2 + (1%4), Just $ 3 - (1%4))]
+        nroots = map (uncurry (numberOfRealRootsInClosedInterval spray)) intervals
+        nroots' = map (uncurry (numberOfRealRootsInOpenInterval spray)) intervals
+        infiniteIntervals = map (, Nothing) [0, 1, 2, 2 + (1%4), 5]
+        nroots'' = map (uncurry (numberOfRealRootsInClosedInterval spray)) infiniteIntervals
+        nroots''' = map (uncurry (numberOfRealRootsInOpenInterval spray)) infiniteIntervals
+      assertEqual "" 
+        (nroots, nroots', nroots'', nroots''') 
+        ([5, 5, 2, 4, 0], [5, 4, 0, 3, 0], [5, 5, 4, 3, 1], [5, 4, 3, 3, 0])
 
     , testCase "number of real roots" $ do
       let