probability-polynomial 1.0.0.0 → 1.0.0.1

raw patch · 5 files changed

+137/−64 lines, 5 filesdep ~QuickCheckdep ~containersnew-uploaderPVP: minor bump suggested

API additions: PVP suggests at least a minor version bump

Dependency ranges changed: QuickCheck, containers

API changes (from Hackage documentation)

Files

• CHANGELOG.md +6/−0
• probability-polynomial.cabal +3/−1
• src/Numeric/Polynomial/Simple.hs +113/−45
• src/Numeric/Probability/Moments.hs +2/−1
• test/Numeric/Polynomial/SimpleSpec.hs +13/−17

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

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

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

@@ -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] \).

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