diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,3 +6,9 @@ * Polynomials * Finite, signed measures on the number line * Probability measures + +## 1.0.0.1 - 2025-01-17 + +* Minor version update + * Improved inmplementation of findRoot to handle repeated roots + * Added test to check repeated roots handled correctly diff --git a/probability-polynomial.cabal b/probability-polynomial.cabal --- a/probability-polynomial.cabal +++ b/probability-polynomial.cabal @@ -5,7 +5,7 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change -version: 1.0.0.0 +version: 1.0.0.1 synopsis: Probability distributions via piecewise polynomials description: Package for manipulating finite probability distributions. @@ -30,6 +30,8 @@ tested-with: , GHC == 9.10.1 + , GHC == 9.6.6 + , GHC == 8.10.7 common warnings ghc-options: -Wall diff --git a/src/Numeric/Polynomial/Simple.hs b/src/Numeric/Polynomial/Simple.hs --- a/src/Numeric/Polynomial/Simple.hs +++ b/src/Numeric/Polynomial/Simple.hs @@ -1,5 +1,5 @@ {-# LANGUAGE DeriveGeneric #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeFamilies #-} @@ -39,6 +39,7 @@ , countRoots , isMonotonicallyIncreasingOn , root + , squareFreeFactorisation ) where import Control.DeepSeq @@ -188,8 +189,8 @@ > eval :: Poly a -> a -> a -} instance Num a => Fun.Function (Poly a) where - type instance Domain (Poly a) = a - type instance Codomain (Poly a) = a + type Domain (Poly a) = a + type Codomain (Poly a) = a eval = eval {-| @@ -227,13 +228,13 @@ -} display :: (Ord a, Eq a, Num a) => Poly a -> (a, a) -> a -> [(a, a)] display p (l, u) s - | s == 0 = map evalPoint [l, u] - | otherwise = map evalPoint (l : go (l + s)) + | s == 0 = map evalPoint [l, u] + | otherwise = map evalPoint (l : go (l + s)) where evalPoint x = (x, eval p x) go x - | x >= u = [u] -- always include the last point - | otherwise = x : go (x + s) + | x >= u = [u] -- always include the last point + | otherwise = x : go (x + s) {-| Linear polymonial connecting the points @(x1, y1)@ and @(x2, y2)@, assuming that @x1 ≠ x2@. @@ -536,13 +537,14 @@ remainder = snd $ euclidianDivision pIminusOne pI go _ = error "reversedSturmSequence: impossible" --- | Check whether a polynomial is monotonically increasing on --- a given interval. +{-| Check whether a polynomial is monotonically increasing on +a given interval. +-} isMonotonicallyIncreasingOn - :: (Fractional a, Eq a, Ord a) => Poly a -> (a,a) -> Bool -isMonotonicallyIncreasingOn p (x1,x2) = + :: (Fractional a, Eq a, Ord a) => Poly a -> (a, a) -> Bool +isMonotonicallyIncreasingOn p (x1, x2) = eval p x1 <= eval p x2 - && countRoots (x1, x2, differentiate p) == 0 + && countRoots (x1, x2, differentiate p) == 0 {-| Measure whether or not a polynomial is consistently above or below zero, @@ -566,49 +568,71 @@ upper = eval p u {-| -Find the root of a polynomial in a given interval, -assuming that there is exactly one root in the given interval. -This precondition has to be checked through other means, -e.g. 'countRoots'. +Find a root of a polynomial in a given interval. -We find the root by repeatedly halving the interval in which the root must lie +Return 'Nothing' if the polynomial does not have a root in the given interval. + +We find the root by first forming the square-free factorisation of the polynomial, +to eliminate repeated roots. One of the factors may have a root in the interval, +so we count roots for each factor until we find the one with a root in the interval. +Then we use the bisection method to find the root, +repeatedly halving the interval in which the root must lie until its width is less than the specified precision. Constant and linear polynomials, @degree p <= 1@, are treated as special cases. -} findRoot - :: (Fractional a, Eq a, Num a, Ord a) => a -> (a, a) -> Poly a -> Maybe a -findRoot precision (l, u) p - -- if the polynomial is zero, the whole interval is a root, so return the basepoint - | degp < 0 = Just l - -- if the poly is a non-zero constant, no root is present - | degp == 0 = Nothing - -- if the polynomial has degree 1, can calculate the root exactly - | degp == 1 = Just (-(head ps / last ps)) -- p0 + p1x = 0 => x = -p0/p1 - | precision <= 0 = error "Invalid precision value" - | otherwise = Just (halveInterval precision l u pl pu) + :: forall a. (Fractional a, Eq a, Num a, Ord a) => a -> (a, a) -> Poly a -> Maybe a +findRoot precision (lower, upper) poly = + if null rootFactors + then Nothing + else getRoot precision (lower, upper) (head rootFactors) where - Poly ps = p - degp = degree p - pu = eval p u - pl = eval p l - halveInterval eps x y px py - -- when the interval is small enough, stop: - -- the root is in this interval, so take the mid point - | width <= eps = mid - -- choose the lower half, - -- as the polynomial has different signs at the ends - | px * pmid < 0 = halveInterval eps x mid px pmid - -- choose the upper half - | otherwise = halveInterval eps mid y pmid py + rootFactors = + filter (\x -> countRoots (lower, upper, x) /= 0) (squareFreeFactorisation poly) + -- getRoot :: forall a. (Fractional a, Eq a, Num a, Ord a) => a -> (a, a) -> Poly a -> Maybe a + getRoot eps (l, u) p + -- if the polynomial is zero, the whole interval is a root, so return the basepoint + | degp < 0 = Just l + -- if the poly is a non-zero constant, no root is present + | degp == 0 = Nothing + -- if the polynomial has degree 1, can calculate the root exactly + | degp == 1 = Just (-(head ps / last ps)) -- p0 + p1x = 0 => x = -p0/p1 + | eps <= 0 = error "Invalid precision value" + | otherwise = bisect eps (l, u) (eval p l, eval p u) p where - width = y - x - mid = x + width / 2 - pmid = eval p mid + ps = toCoefficients p + degp = degree p + {- We bisect the interval exploiting the Intermediate Value Theorem: + if a polynomial has different signs at the ends of an interval, it must be zero somewhere in the interval. + If there is no change of sign, use countRoots to find which side of the interval the root is on. + -} + bisect + :: (Fractional a, Eq a, Num a, Ord a) => a -> (a, a) -> (a, a) -> Poly a -> Maybe a + bisect e (x, y) (px, py) p' + -- if we already have a root, choose it + | px == 0 = Just x + | py == 0 = Just y + | pmid == 0 = Just mid + -- when the interval is small enough, stop: + -- the root is in this interval, so take the mid point + | width <= e = Just mid + -- choose the lower half, if the polynomial has different signs at the ends + | signum px /= signum pmid = bisect e (x, mid) (px, pmid) p' + -- choose the upper half, if the polynomial has different signs at the ends + | signum py /= signum pmid = bisect e (mid, y) (pmid, py) p' + -- no sign change found, so we resort to counting roots + | countRoots (x, mid, p') > 0 = bisect e (x, mid) (px, pmid) p' + | countRoots (mid, y, p') > 0 = bisect e (mid, y) (pmid, py) p' + | otherwise = Nothing + where + width = y - x + mid = x + width / 2 + pmid = eval p' mid -{-| Otherwise we have a polynomial: +{-| We are seeking the point at which a polynomial has a specific value.: subtract the value we are looking for so that we seek a zero crossing -} -root +root -- TODO: this should probably be called something else such as 'findCrossing' :: (Ord a, Num a, Eq a, Fractional a) => a -> a @@ -616,3 +640,47 @@ -> Poly a -> Maybe a root e x (l, u) p = findRoot e (l, u) (p - constant x) + +-- | Greatest monic common divisor of two polynomials. +gcdPoly + :: forall a. (Fractional a, Eq a, Num a, Ord a) => Poly a -> Poly a -> Poly a +gcdPoly a b = if b == zero then a else makeMonic (gcdPoly b (polyRemainder a b)) + where + makeMonic :: Poly a -> Poly a + makeMonic (Poly as) = scale (1 / last as) (Poly as) + polyRemainder :: Poly a -> Poly a -> Poly a + polyRemainder x y = snd (euclidianDivision x y) + +{-| +We compute the square-free factorisation of a polynomial using Yun's algorithm. +Yun, David Y.Y. (1976). "On square-free decomposition algorithms". +SYMSAC '76 Proceedings of the third ACM Symposium on Symbolic and Algebraic Computation. +Association for Computing Machinery. pp. 26–35. doi:10.1145/800205.806320. ISBN 978-1-4503-7790-4. S2CID 12861227. +https://dl.acm.org/doi/10.1145/800205.806320 +G <- gcd (P, P') +C1 <- P / G +D1 <- P' / G - C1' +until Ci = 1 do + Pi <- gcd (Ci, Di) + Ci+1 <- Ci/Pi + Di+1 <- Di / Ai - Ci+1' +-} +squareFreeFactorisation + :: (Fractional a, Eq a, Num a, Ord a) => Poly a -> [Poly a] +squareFreeFactorisation p = + -- if p has degree <= 1 it can have no factors but itself + if degree p <= 1 then [p] else go c1 d1 + where + diffP = differentiate p + g0 = gcdPoly p diffP + c1 = p `divide` g0 + d1 = (diffP `divide` g0) - differentiate c1 + divide x y = fst (euclidianDivision x y) + go c d + | c == constant 1 = [] -- terminate the recursion + | a' == constant 1 = go c' d' -- skip over the constant polynomial + | otherwise = a' : go c' d' + where + a' = gcdPoly c d + c' = c `divide` a' + d' = (d `divide` a') - differentiate c' diff --git a/src/Numeric/Probability/Moments.hs b/src/Numeric/Probability/Moments.hs --- a/src/Numeric/Probability/Moments.hs +++ b/src/Numeric/Probability/Moments.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE DeriveFunctor #-} {-# LANGUAGE NamedFieldPuns #-} {-| @@ -35,7 +36,7 @@ -- -- The kurtosis is bounded below: \( \kappa \geq \gamma_1^2 + 1 \). } - deriving (Eq, Show) + deriving (Eq, Show, Functor) -- | Compute the 'Moments' of a probability distribution given -- the expectation values of the first four powers \( m_k = E[X^k] \). diff --git a/test/Numeric/Polynomial/SimpleSpec.hs b/test/Numeric/Polynomial/SimpleSpec.hs --- a/test/Numeric/Polynomial/SimpleSpec.hs +++ b/test/Numeric/Polynomial/SimpleSpec.hs @@ -45,22 +45,18 @@ ) import Test.Hspec ( Spec - , before_ , describe , it - , pendingWith ) import Test.QuickCheck ( Arbitrary , Gen , NonNegative (..) , Positive (..) - , Property , (===) , (==>) , (.&&.) , arbitrary - , counterexample , forAll , frequency , listOf @@ -74,9 +70,6 @@ {----------------------------------------------------------------------------- Tests ------------------------------------------------------------------------------} -xit' :: String -> String -> Property -> Spec -xit' reason label = before_ (pendingWith reason) . it label - spec :: Spec spec = do describe "constant" $ do @@ -217,7 +210,7 @@ in property $ abs (x2' - x2) <= epsilon - xit' "bug" "cubic polynomial, midpoint" $ property $ mapSize (`div` 5) $ + it "cubic polynomial, midpoint" $ property $ mapSize (`div` 5) $ \(x1 :: Rational) (Positive dx3) -> let xx = scaleX (constant 1) :: Poly Rational x2 = (x1 + x3) / 2 @@ -228,15 +221,18 @@ epsilon = (x3-x1)/(1000*1000*50) Just x2' = root epsilon 0 (l, u) p in - id - $ counterexample ("interval = " <> show (l,u)) - $ counterexample ("countRoots = " <> show (countRoots (l, u, p))) - $ counterexample ("expected root = " <> show x2) - $ counterexample ("eval polynomial at expected root = " <> show (eval p x2)) - $ counterexample ("epsilon = " <> show epsilon) - $ counterexample ("found root = " <> show x2') - $ counterexample ("root within range of other root " <> show (abs (x2' - x3) <= 20*epsilon)) - $ property $ abs (x2' - x2) <= epsilon + abs (x2' - x2) <= epsilon + + it "high-degree polynomial, repeated root" $ property $ mapSize (`div` 5) $ + \(Positive (x1 :: Rational)) (Positive (n :: Integer)) -> + let xx = scaleX (constant 1) :: Poly Rational + p = xx*(xx - constant x1)^n + l = x1 - 100 * epsilon + u = x1 + 100 * epsilon + epsilon = x1/(1000*1000*50) + Just x2' = root epsilon 0 (l, u) p + in + abs (x2' - x1) <= epsilon describe "isMonotonicallyIncreasingOn" $ it "quadratic polynomial" $ property $