math-functions 0.1.6.0 → 0.1.7.0

raw patch · 10 files changed

+246/−66 lines, 10 filesdep +primitivedep −ieee754dep ~QuickCheckdep ~erfPVP ok

version bump matches the API change (PVP)

Dependencies added: primitive

Dependencies removed: ieee754

Dependency ranges changed: QuickCheck, erf

API changes (from Hackage documentation)

Files

• Numeric/MathFunctions/Comparison.hs +114/−0
• Numeric/SpecFunctions/Internal.hs +7/−4
• Numeric/Sum.hs +3/−1
• changelog.md +17/−0
• math-functions.cabal +9/−3
• tests/Tests/Chebyshev.hs +28/−20
• tests/Tests/Comparison.hs +37/−0
• tests/Tests/Helpers.hs +3/−19
• tests/Tests/SpecFunctions.hs +24/−18
• tests/tests.hs +4/−1

@@ -0,0 +1,114 @@+-- |+-- Module    : Numeric.MathFunctions.Comparison+-- Copyright : (c) 2011 Bryan O'Sullivan+-- License   : BSD3+--+-- Maintainer  : bos@serpentine.com+-- Stability   : experimental+-- Portability : portable+--+-- Functions for approximate comparison of floating point numbers.+--+-- Approximate floating point comparison, based on Bruce Dawson's+-- \"Comparing floating point numbers\":+-- <http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm>+module Numeric.MathFunctions.Comparison+    ( -- * Relative erros+      relativeError+    , eqRelErr+      -- * Ulps-based comparison+    , addUlps+    , ulpDistance+    , within+    ) where++import Control.Monad.ST (runST)+import Data.Primitive.ByteArray (newByteArray, readByteArray, writeByteArray)+import Data.Word (Word64)+import Data.Int (Int64)++++----------------------------------------------------------------+-- Ulps-based comparison+----------------------------------------------------------------++-- |+-- Calculate relative error of two numbers:+--+-- > |a - b| / max |a| |b|+--+-- It lies in [0,1) interval for numbers with same sign and (1,2] for+-- numbers with different sign. If both arguments are zero or negative+-- zero function returns 0. If at least one argument is transfinite it+-- returns NaN+relativeError :: Double -> Double -> Double+relativeError a b+  | a == 0 && b == 0 = 0+  | otherwise        = abs (a - b) / max (abs a) (abs b)++-- | Check that relative error between two numbers @a@ and @b@. If+-- 'relativeError' returns NaN it returns @False@.+eqRelErr :: Double -- ^ @eps@ relative error should be in [0,1) range+         -> Double -- ^ @a@+         -> Double -- ^ @b@+         -> Bool+eqRelErr eps a b = relativeError a b < eps++++----------------------------------------------------------------+-- Ulps-based comparison+----------------------------------------------------------------++-- |+-- Add N ULPs (units of least precision) to @Double@ number.+addUlps :: Int -> Double -> Double+addUlps n a = runST $ do+  buf <- newByteArray 8+  ai0 <- writeByteArray buf 0 a >> readByteArray buf 0+  -- Convert to ulps number represented as Int64+  let big     = 0x8000000000000000+      order :: Word64 -> Int64+      order i | i < big   = fromIntegral i+              | otherwise = fromIntegral $ maxBound - (i - big)+      unorder :: Int64 -> Word64+      unorder i | i >= 0    = fromIntegral i+                | otherwise = big + (maxBound - (fromIntegral i))+  let ai0' = unorder $ order ai0 + fromIntegral n+  writeByteArray buf 0 ai0' >> readByteArray buf 0++-- |+-- Measure distance between two @Double@s in ULPs (units of least+-- precision).+ulpDistance :: Double+            -> Double+            -> Word64+ulpDistance a b = runST $ do+  buf <- newByteArray 8+  ai0 <- writeByteArray buf 0 a >> readByteArray buf 0+  bi0 <- writeByteArray buf 0 b >> readByteArray buf 0+  -- IEEE754 floats use most significant bit as sign bit (not+  -- 2-complement) and we need to rearrange representations of float+  -- number so that they could be compared lexicographically as+  -- Word64.+  let big     = 0x8000000000000000+      order i | i < big   = i + big+              | otherwise = maxBound - i+      ai = order ai0+      bi = order bi0+      d  | ai > bi   = ai - bi+         | otherwise = bi - ai+  return $! d++-- | Compare two 'Double' values for approximate equality, using+-- Dawson's method.+--+-- The required accuracy is specified in ULPs (units of least+-- precision).  If the two numbers differ by the given number of ULPs+-- or less, this function returns @True@.+within :: Int                   -- ^ Number of ULPs of accuracy desired.+       -> Double -> Double -> Bool+within ulps a b+  | ulps < 0  = False+  | otherwise = ulpDistance a b <= fromIntegral ulps

@@ -639,11 +639,14 @@ logFactorial n     | n <  0    = error "Numeric.SpecFunctions.logFactorial: negative input"     | n <= 14   = log $ factorial $ fromIntegral n-    | otherwise = (x - 0.5) * log x - x + 9.1893853320467e-1 + z / x+    -- N.B. Γ(n+1) = n!+    --+    -- We use here asymptotic series for gamma function. See+    -- http://mathworld.wolfram.com/StirlingsSeries.html+    | otherwise = (x - 0.5) * log x - x+                + m_ln_sqrt_2_pi+                + evaluateOddPolynomialL (1/x) [1/12, -1/360, 1/1260, -1/1680]     where x = fromIntegral n + 1-          y = 1 / (x * x)-          z = ((-(5.95238095238e-4 * y) + 7.936500793651e-4) * y --               2.7777777777778e-3) * y + 8.3333333333333e-2 {-# SPECIALIZE logFactorial :: Int -> Double #-}  -- | Calculate the error term of the Stirling approximation.  This is

@@ -54,8 +54,10 @@ import Data.Bits (shiftR) import Data.Data (Typeable, Data) import Data.Vector.Generic (Vector(..), foldl')-import Data.Vector.Generic.Mutable (MVector(..)) import Data.Vector.Unboxed.Deriving (derivingUnbox)+-- Needed for GHC 7.2 & 7.4 to derive Unbox instances+import Data.Vector.Generic.Mutable (MVector(..))+ import qualified Data.Foldable as F import qualified Data.Vector as V import qualified Data.Vector.Generic as G

@@ -1,3 +1,16 @@+Changes in 0.1.7.0++  * Module `statistics: Statistics.Function.Comparison` moved to+    `Numeric.MathFunctions.Comparison`. Old implementation if `within` compared+    negative numbers incorrectly.++  * `addUlps` and `ulpDistance` added to `Numeric.MathFunctions.Comparison`.++  * `relativeError` and `eqRelErr` added to `Numeric.MathFunctions.Comparison`.++  * Precision of `logFactorial` is slightly improved.++ Changes in 0.1.6.0    * `logChoose` added for calculation of logarithm of binomial coefficient@@ -6,20 +19,24 @@    * `sinc` added + Changes in 0.1.5.3    * Fix for test suite on 32bit platform + Changes in 0.1.5    * Numeric.Sum: new module adds accurate floating point summation. + Changes in 0.1.4    * logFactorial type is genberalized. It accepts any `Integral` type    * Evaluation of polynomials using Horner's method where coefficients     are store in lists added+  Changes in 0.1.3 

@@ -1,5 +1,5 @@ name:           math-functions-version:        0.1.6.0+version:        0.1.7.0 cabal-version:  >= 1.8 license:        BSD3 license-file:   LICENSE@@ -31,9 +31,11 @@                         deepseq,                         erf >= 2,                         vector >= 0.7,+                        primitive,                         vector-th-unbox   exposed-modules:     Numeric.MathFunctions.Constants+    Numeric.MathFunctions.Comparison     Numeric.Polynomial     Numeric.Polynomial.Chebyshev     Numeric.SpecFunctions@@ -53,16 +55,20 @@   other-modules:     Tests.Helpers     Tests.Chebyshev+    Tests.Comparison     Tests.SpecFunctions     Tests.SpecFunctions.Tables     Tests.Sum   build-depends:     math-functions,     base >=3 && <5,+    deepseq,+    primitive,     vector >= 0.7,-    ieee754 >= 0.7.3,+    vector-th-unbox,+    erf,     HUnit      >= 1.2,-    QuickCheck >= 2.4,+    QuickCheck >= 2.7,     test-framework,     test-framework-hunit,     test-framework-quickcheck2

@@ -7,45 +7,53 @@ import Data.Vector.Unboxed                  (fromList) import Test.Framework import Test.Framework.Providers.QuickCheck2-import Test.QuickCheck                      (Arbitrary(..),printTestCase,Property)+import Test.QuickCheck                      (Arbitrary(..),counterexample,Property)  import Tests.Helpers import Numeric.Polynomial.Chebyshev+import Numeric.MathFunctions.Comparison   tests :: Test tests = testGroup "Chebyshev polynomials"   [ testProperty "Chebyshev 0" $ \a0 (Ch x) ->       testCheb [a0] x-  -- XXX FIXME DISABLED due to failure-  -- , testProperty "Chebyshev 1" $ \a0 a1 (Ch x) ->-  --   testCheb [a0,a1] x-  -- , testProperty "Chebyshev 2" $ \a0 a1 a2 (Ch x) ->-  --   testCheb [a0,a1,a2] x-  -- , testProperty "Chebyshev 3" $ \a0 a1 a2 a3 (Ch x) ->-  --   testCheb [a0,a1,a2,a3] x-  -- , testProperty "Chebyshev 4" $ \a0 a1 a2 a3 a4 (Ch x) ->-  --   testCheb [a0,a1,a2,a3,a4] x-  -- , testProperty "Broucke" $ testBroucke+  , testProperty "Chebyshev 1" $ \a0 a1 (Ch x) ->+      testCheb [a0,a1] x+  , testProperty "Chebyshev 2" $ \a0 a1 a2 (Ch x) ->+      testCheb [a0,a1,a2] x+  , testProperty "Chebyshev 3" $ \a0 a1 a2 a3 (Ch x) ->+      testCheb [a0,a1,a2,a3] x+  , testProperty "Chebyshev 4" $ \a0 a1 a2 a3 a4 (Ch x) ->+      testCheb [a0,a1,a2,a3,a4] x+  , testProperty "Broucke" $ testBroucke   ]   where -testBroucke :: Ch -> [Double] -> Bool-testBroucke _      []     = True-testBroucke (Ch x) (c:cs) = let c1 = chebyshev        x (fromList $ c : cs)-                                cb = chebyshevBroucke x (fromList $ c*2 : cs)-                            in eq 1e-15 c1 cb+testBroucke :: Ch -> [Double] -> Property+testBroucke _      []     = counterexample "" True+testBroucke (Ch x) (c:cs)+  = counterexample (">>> Chebyshev  = " ++ show c1)+  $ counterexample (">>> Brouke     = " ++ show cb)+  $ counterexample (">>> rel.err.   = " ++ show (relativeError c1 cb))+  $ counterexample (">>> diff. ulps = " ++ show (ulpDistance   c1 cb))+  $ within 64 c1 cb+  where+    c1 = chebyshev        x (fromList $ c : cs)+    cb = chebyshevBroucke x (fromList $ c*2 : cs) + testCheb :: [Double] -> Double -> Property testCheb as x-  = printTestCase (">>> Exact   = " ++ show exact)-  $ printTestCase (">>> Numeric = " ++ show num  )-  $ printTestCase (">>> rel.err.= " ++ show err  )+  = counterexample (">>> Exact      = " ++ show exact)+  $ counterexample (">>> Numeric    = " ++ show num  )+  $ counterexample (">>> rel.err.   = " ++ show err  )+  $ counterexample (">>> diff. ulps = " ++ show (ulpDistance num exact))   $ eq 1e-12 num exact   where     exact = evalCheb as x     num   = chebyshev x (fromList as)-    err   = abs (num - exact) / abs exact+    err   = relativeError num exact  evalCheb :: [Double] -> Double -> Double evalCheb as x

@@ -0,0 +1,37 @@+-- |+-- Tests for approximate comparison+module Tests.Comparison (tests) where++import Test.Framework+import Test.Framework.Providers.QuickCheck2++import Tests.Helpers++import Numeric.MathFunctions.Comparison+import Numeric.MathFunctions.Constants (m_epsilon)++tests :: Test+tests = testGroup "Comparison"+  [ testProperty "addUlps 0"       $ \x   -> x == addUlps 0 x+  , testProperty "addUlps sym"     $ \i x -> x == (addUlps (-i) . addUlps i) x+  , testProperty "ulpDistance==0"  $ \x   -> ulpDistance x x == 0+  , testProperty "ulpDistance sym" $ \x y -> ulpDistance x y == ulpDistance y x+  , testProperty "ulpDistance/addUlps" $ \x i -> ulpDistance x (addUlps i x) == fromIntegral (abs i)+    -- Test that code is correct for m_epsilon+  , testAssertion "eps distance" $ ulpDistance 1 (1+m_epsilon) == 1+  , testAssertion "eps add"      $ addUlps 1 1 == 1 + m_epsilon+    --+  , testProperty  "relativeError sym"     $ \x y -> relativeError x y == relativeError y x+  , testAssertion "relativeError inf   1" $ isNaN $ relativeError inf 1+  , testAssertion "relativeError 1   inf" $ isNaN $ relativeError 1 inf+  , testAssertion "relativeError -inf  1" $ isNaN $ relativeError (-inf) 1+  , testAssertion "relativeError 1  -inf" $ isNaN $ relativeError 1 (-inf)+  , testAssertion "relativeError inf inf" $ isNaN $ relativeError inf inf+  , testAssertion "relativeError inf-inf" $ isNaN $ relativeError inf (-inf)+  , testAssertion "relativeError   1 Nan" $ isNaN $ relativeError 1 nan+  , testAssertion "relativeError NaN   1" $ isNaN $ relativeError nan 1+  , testAssertion "relativeError NaN Nan" $ isNaN $ relativeError nan nan+  ]+  where+    inf = 1/0+    nan = 0/0

@@ -7,7 +7,6 @@   , eqC     -- * Generic QC tests   , monotonicallyIncreases-  , monotonicallyIncreasesIEEE     -- * HUnit helpers   , testAssertion   , testEquality@@ -16,12 +15,11 @@ import Data.Complex import Data.Typeable -import qualified Numeric.IEEE    as IEEE- import qualified Test.HUnit      as HU import Test.Framework import Test.Framework.Providers.HUnit +import Numeric.MathFunctions.Comparison   @@ -43,9 +41,7 @@ --   which are almost zero. eq :: Double                    -- ^ Relative error    -> Double -> Double -> Bool-eq eps a b -  | a == 0 && b == 0 = True-  | otherwise        = abs (a - b) <= eps * max (abs a) (abs b)+eq = eqRelErr  -- | Approximate equality for 'Complex Double' eqC :: Double                   -- ^ Relative error@@ -69,19 +65,7 @@ monotonicallyIncreases :: (Ord a, Ord b) => (a -> b) -> a -> a -> Bool monotonicallyIncreases f x1 x2 = f (min x1 x2) <= f (max x1 x2) --- Check that function is nondecreasing taking rounding errors into--- account.------ In fact funstion is allowed to decrease less than one ulp in order--- to guard againist problems with excess precision. On x86 FPU works--- with 80-bit numbers but doubles are 64-bit so rounding happens--- whenever values are moved from registers to memory-monotonicallyIncreasesIEEE :: (Ord a, IEEE.IEEE b)  => (a -> b) -> a -> a -> Bool-monotonicallyIncreasesIEEE f x1 x2 =-  y1 <= y2 || (y1 - y2) < y2 * IEEE.epsilon-  where-    y1 = f (min x1 x2)-    y2 = f (max x1 x2)+  ---------------------------------------------------------------- -- HUnit helpers

@@ -7,13 +7,14 @@ import qualified Data.Vector as V import           Data.Vector   ((!)) -import Test.QuickCheck  hiding (choose)+import Test.QuickCheck  hiding (choose,within) import Test.Framework import Test.Framework.Providers.QuickCheck2  import Tests.Helpers import Tests.SpecFunctions.Tables import Numeric.SpecFunctions+import Numeric.MathFunctions.Comparison (within)   tests :: Test@@ -79,6 +80,11 @@   , testAssertion "choose is expected to precise at 1e-12 level"       $ and [ eq 1e-12 (choose (fromIntegral n) (fromIntegral k)) (fromIntegral $ choose' n k)             | n <- [0..1000], k <- [0..n]]+  , testAssertion "logChoose == log . choose"+      $ and [ let n' = fromIntegral n+                  k' = fromIntegral k+              in within 2 (logChoose n' k') (log $ choose n' k')+            | n <- [0..1000], k <- [0..n]]     ----------------------------------------------------------------     -- Self tests   , testProperty "Self-test: 0 <= range01 <= 1" $ \x -> let f = range01 x in f <= 1 && f >= 0@@ -112,11 +118,11 @@ -- invIncompleteGamma is inverse of incompleteGamma invIGammaIsInverse :: Double -> Double -> Property invIGammaIsInverse (abs -> a) (range01 -> p) =-  a > 0 && p > 0 && p < 1  ==> ( printTestCase ("a  = " ++ show a )-                               $ printTestCase ("p  = " ++ show p )-                               $ printTestCase ("x  = " ++ show x )-                               $ printTestCase ("p' = " ++ show p')-                               $ printTestCase ("Δp = " ++ show (p - p'))+  a > 0 && p > 0 && p < 1  ==> ( counterexample ("a  = " ++ show a )+                               $ counterexample ("p  = " ++ show p )+                               $ counterexample ("x  = " ++ show x )+                               $ counterexample ("p' = " ++ show p')+                               $ counterexample ("Δp = " ++ show (p - p'))                                $ abs (p - p') <= 1e-12                                )   where@@ -126,9 +132,9 @@ -- invErfc is inverse of erfc invErfcIsInverse :: Double -> Property invErfcIsInverse ((*2) . range01 -> p)-  = printTestCase ("p  = " ++ show p )-  $ printTestCase ("x  = " ++ show x )-  $ printTestCase ("p' = " ++ show p')+  = counterexample ("p  = " ++ show p )+  $ counterexample ("x  = " ++ show x )+  $ counterexample ("p' = " ++ show p')   $ abs (p - p') <= 1e-14   where     x  = invErfc p@@ -137,9 +143,9 @@ -- invErf is inverse of erf invErfIsInverse :: Double -> Property invErfIsInverse a-  = printTestCase ("p  = " ++ show p )-  $ printTestCase ("x  = " ++ show x )-  $ printTestCase ("p' = " ++ show p')+  = counterexample ("p  = " ++ show p )+  $ counterexample ("x  = " ++ show x )+  $ counterexample ("p' = " ++ show p')   $ abs (p - p') <= 1e-14   where     x  = invErf p@@ -155,12 +161,12 @@ -- invIncompleteBeta is inverse of incompleteBeta invIBetaIsInverse :: Double -> Double -> Double -> Property invIBetaIsInverse (abs -> p) (abs -> q) (range01 -> x) =-  p > 0 && q > 0  ==> ( printTestCase ("p   = " ++ show p )-                      $ printTestCase ("q   = " ++ show q )-                      $ printTestCase ("x   = " ++ show x )-                      $ printTestCase ("x'  = " ++ show x')-                      $ printTestCase ("a   = " ++ show a)-                      $ printTestCase ("err = " ++ (show $ abs $ (x - x') / x))+  p > 0 && q > 0  ==> ( counterexample ("p   = " ++ show p )+                      $ counterexample ("q   = " ++ show q )+                      $ counterexample ("x   = " ++ show x )+                      $ counterexample ("x'  = " ++ show x')+                      $ counterexample ("a   = " ++ show a)+                      $ counterexample ("err = " ++ (show $ abs $ (x - x') / x))                       $ abs (x - x') <= 1e-12                       )   where

@@ -2,9 +2,12 @@ import qualified Tests.SpecFunctions import qualified Tests.Chebyshev import qualified Tests.Sum+import qualified Tests.Comparison  main :: IO () main = defaultMain [ Tests.SpecFunctions.tests-                   , Tests.Chebyshev.tests+                   -- FIXME: tests for chebyshev polynomials fail intermittently+                   -- , Tests.Chebyshev.tests                    , Tests.Sum.tests+                   , Tests.Comparison.tests                    ]