packages feed

math-functions 0.1.1.1 → 0.1.1.2

raw patch · 4 files changed

+112/−41 lines, 4 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

Numeric/SpecFunctions.hs view
@@ -1,7 +1,7 @@-{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE BangPatterns, ScopedTypeVariables #-} -- | -- Module    : Numeric.SpecFunctions--- Copyright : (c) 2009, 2011 Bryan O'Sullivan+-- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan -- License   : BSD3 -- -- Maintainer  : bos@serpentine.com@@ -71,7 +71,7 @@                   ((((x + r2_8) * x + r2_7) * x + r2_6) * x + r2_5)     | x < 12    = ((((r3_4 * x + r3_3) * x + r3_2) * x + r3_1) * x + r3_0) /                   ((((x + r3_8) * x + r3_7) * x + r3_6) * x + r3_5)-    | x > 5.1e5 = k+    | x > 3e6   = k     | otherwise = k + x1 *                   ((r4_2 * x2 + r4_1) * x2 + r4_0) /                   ((x2 + r4_4) * x2 + r4_3)@@ -168,17 +168,18 @@                 -> Double       -- ^ /x/                 -> Double incompleteGamma p x-    | x < 0 || p <= 0 = m_pos_inf-    | x == 0          = 0-    | p >= 1000       = norm (3 * sqrt p * ((x/p) ** (1/3) + 1/(9*p) - 1))-    | x >= 1e8        = 1-    | x <= 1 || x < p = let a = p * log x - x - logGamma (p + 1)-                            g = a + log (pearson p 1 1)-                        in if g > limit then exp g else 0-    | otherwise       = let g = p * log x - x - logGamma p + log cf-                        in if g > limit then 1 - exp g else 1+    | isNaN p || isNaN x = m_NaN+    | x < 0 || p <= 0    = m_pos_inf+    | x == 0             = 0+    | p >= 1000          = norm (3 * sqrt p * ((x/p) ** (1/3) + 1/(9*p) - 1))+    | x >= 1e8           = 1+    | x <= 1 || x < p    = let a = p * log x - x - logGamma (p + 1)+                               g = a + log (pearson p 1 1)+                           in if g > limit then exp g else 0+    | otherwise          = let g = p * log x - x - logGamma p + log cf+                           in if g > limit then 1 - exp g else 1   where-    norm a = erfc (- a / m_sqrt_2)+    norm a = 0.5 * erfc (- a / m_sqrt_2)     pearson !a !c !g         | c' <= tolerance = g'         | otherwise       = pearson a' c' g'@@ -231,24 +232,25 @@     -- Solve equation γ(a,x) = p using Halley method     loop :: Int -> Double -> Double     loop i x-      | i >= 12   = x-      | otherwise =-         let -           -- Value of γ(a,x) - p-           f    = incompleteGamma a x - p-           -- dγ(a,x)/dx-           f'   | a > 1     = afac * exp( -(x - a1) + a1 * (log x - lna1))-                | otherwise = exp( -x + a1 * log x - gln)-           u    = f / f'-           -- Halley correction to Newton-Rapson step-           corr = u * (a1 / x - 1)-           dx   = u / (1 - 0.5 * min 1.0 corr)-           -- New approximation to x-           x'   | x < dx    = 0.5 * x -- Do not go below 0-                | otherwise = x - dx-         in if abs dx < eps * x'-            then x'-            else loop (i+1) x'+      | i >= 12           = x'+      -- For small s derivative becomes approximately 1/x*exp(-x) and+      -- skyrockets for small x. If it happens correct answer is 0.+      | isInfinite f'     = 0+      | abs dx < eps * x' = x'+      | otherwise         = loop (i + 1) x'+      where+        -- Value of γ(a,x) - p+        f    = incompleteGamma a x - p+        -- dγ(a,x)/dx+        f'   | a > 1     = afac * exp( -(x - a1) + a1 * (log x - lna1))+             | otherwise = exp( -x + a1 * log x - gln)+        u    = f / f'+        -- Halley correction to Newton-Rapson step+        corr = u * (a1 / x - 1)+        dx   = u / (1 - 0.5 * min 1.0 corr)+        -- New approximation to x+        x'   | x < dx    = 0.5 * x -- Do not go below 0+             | otherwise = x - dx     -- Calculate inital guess for root     guess       -- @@ -297,7 +299,7 @@       c   = logGammaCorrection q - logGammaCorrection pq  -- | Regularized incomplete beta function. Uses algorithm AS63 by---   Majumder abd Bhattachrjee.+-- Majumder and Bhattachrjee. incompleteBeta :: Double -- ^ /p/ > 0                -> Double -- ^ /q/ > 0                -> Double -- ^ /x/, must lie in [0,1] range@@ -312,22 +314,22 @@                 -> Double -- ^ /x/, must lie in [0,1] range                 -> Double incompleteBeta_ beta p q x-  | p <= 0 || q <= 0 = error "p <= 0 || q <= 0"-  | x <  0 || x >  1 = error "x <  0 || x >  1"-  | x == 0 || x == 1 = x+  | p <= 0 || q <= 0            = error "p <= 0 || q <= 0"+  | x <  0 || x >  1 || isNaN x = error "x out of [0,1] range"+  | x == 0 || x == 1            = x   | p >= (p+q) * x   = incompleteBetaWorker beta p q x   | otherwise        = 1 - incompleteBetaWorker beta q p (1 - x)  -- Worker for incomplete beta function. It is separate function to -- avoid confusion with parameter during parameter swapping incompleteBetaWorker :: Double -> Double -> Double -> Double -> Double-incompleteBetaWorker beta p q x = loop (p+q) (truncate $ q + cx * (p+q) :: Int) 1 1 1+incompleteBetaWorker beta p q x = loop (p+q) (truncate $ q + cx * (p+q)) 1 1 1   where     -- Constants     eps = 1e-15     cx  = 1 - x     -- Loop-    loop psq ns ai term betain+    loop !psq (ns :: Int) ai term betain       | done      = betain' * exp( p * log x + (q - 1) * log cx - beta) / p       | otherwise = loop psq' (ns - 1) (ai + 1) term' betain'       where
math-functions.cabal view
@@ -1,5 +1,5 @@ name:           math-functions-version:        0.1.1.1+version:        0.1.1.2 cabal-version:  >= 1.8 license:        BSD3 license-file:   LICENSE@@ -31,10 +31,12 @@     Numeric.MathFunctions.Constants  test-suite tests+  buildable:      False   type:           exitcode-stdio-1.0   hs-source-dirs: tests   main-is:        tests.hs   other-modules:+    Tests.Helpers     Tests.Chebyshev     Tests.SpecFunctions     Tests.SpecFunctions.Tables
tests/Tests/SpecFunctions.hs view
@@ -21,9 +21,11 @@   [ testProperty "Γ(x+1) = x·Γ(x) logGamma"  $ gammaReccurence logGamma  3e-8   , testProperty "Γ(x+1) = x·Γ(x) logGammaL" $ gammaReccurence logGammaL 2e-13   , testProperty "γ(1,x) = 1 - exp(-x)"      $ incompleteGammaAt1Check+  , testProperty "0 <= γ <= 1"               $ incompleteGammaInRange   , testProperty "γ - increases"             $       \s x y -> s > 0 && x > 0 && y > 0 ==> monotonicallyIncreases (incompleteGamma s) x y   , testProperty "invIncompleteGamma = γ^-1" $ invIGammaIsInverse+  , testProperty "0 <= I[B] <= 1"            $ incompleteBetaInRange   , testProperty "invIncompleteBeta  = B^-1" $ invIBetaIsInverse     -- Unit tests   , testAssertion "Factorial is expected to be precise at 1e-15 level"@@ -57,6 +59,9 @@   , 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..300], k <- [0..n]]+    ----------------------------------------------------------------+    -- Self tests+  , testProperty "Self-test: 0 <= range01 <= 1" $ \x -> let f = range01 x in f <= 1 && f >= 0   ]  ----------------------------------------------------------------@@ -71,19 +76,25 @@       g1 = logG x       g2 = logG (x+1) +-- γ(s,x) is in [0,1] range+incompleteGammaInRange :: Double -> Double -> Property+incompleteGammaInRange (abs -> s) (abs -> x) =+  x >= 0 && s > 0  ==> let i = incompleteGamma s x in i >= 0 && i <= 1  -- γ(1,x) = 1 - exp(-x) -- Since Γ(1) = 1 normalization doesn't make any difference incompleteGammaAt1Check :: Double -> Property-incompleteGammaAt1Check x =+incompleteGammaAt1Check (abs -> x) =   x > 0 ==> (incompleteGamma 1 x + exp(-x)) ≈ 1   where     (≈) = eq 1e-13  -- invIncompleteGamma is inverse of incompleteGamma invIGammaIsInverse :: Double -> Double -> Property-invIGammaIsInverse (abs -> a) (abs . snd . properFraction -> p) =-  a > 0 && p > 0 && p < 1  ==> ( printTestCase ("x  = " ++ show x )+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'))                                $ abs (p - p') <= 1e-12@@ -92,6 +103,11 @@     x  = invIncompleteGamma a p     p' = incompleteGamma    a x +-- B(s,x) is in [0,1] range+incompleteBetaInRange :: Double -> Double -> Double -> Property+incompleteBetaInRange (abs -> p) (abs -> q) (range01 -> x) =+  p > 0 && q > 0  ==> let i = incompleteBeta p q x in i >= 0 && i <= 1+ -- invIncompleteBeta is inverse of incompleteBeta invIBetaIsInverse :: Double -> Double -> Double -> Property invIBetaIsInverse (abs -> p) (abs -> q) (abs . snd . properFraction -> x) =@@ -125,3 +141,7 @@ -- Exact albeit slow implementation of choose choose' :: Integer -> Integer -> Integer choose' n k = factorial' n `div` (factorial' k * factorial' (n-k))++-- Truncate double to [0,1]+range01 :: Double -> Double+range01 = abs . snd . properFraction
+ tests/Tests/SpecFunctions/Tables.hs view
@@ -0,0 +1,47 @@+module Tests.SpecFunctions.Tables where++tableLogGamma :: [(Double,Double)]+tableLogGamma =+  [(0.000001250000000, 13.592366285131769033)+  , (0.000068200000000, 9.5930266308318756785)+  , (0.000246000000000, 8.3100370767447966358)+  , (0.000880000000000, 7.03508133735248542)+  , (0.003120000000000, 5.768129358365567505)+  , (0.026700000000000, 3.6082588918892977148)+  , (0.077700000000000, 2.5148371858768232556)+  , (0.234000000000000, 1.3579557559432759994)+  , (0.860000000000000, 0.098146578027685615897)+  , (1.340000000000000, -0.11404757557207759189)+  , (1.890000000000000, -0.0425116422978701336)+  , (2.450000000000000, 0.25014296569217625565)+  , (3.650000000000000, 1.3701041997380685178)+  , (4.560000000000000, 2.5375143317949580002)+  , (6.660000000000000, 5.9515377269550207018)+  , (8.250000000000000, 9.0331869196051233217)+  , (11.300000000000001, 15.814180681373947834)+  , (25.600000000000001, 56.711261598328121636)+  , (50.399999999999999, 146.12815158702164808)+  , (123.299999999999997, 468.85500075897556371)+  , (487.399999999999977, 2526.9846647543727158)+  , (853.399999999999977, 4903.9359135978220365)+  , (2923.300000000000182, 20402.93198938705973)+  , (8764.299999999999272, 70798.268343590112636)+  , (12630.000000000000000, 106641.77264982508495)+  , (34500.000000000000000, 325976.34838781820145)+  , (82340.000000000000000, 849629.79603036714252)+  , (234800.000000000000000, 2668846.4390507959761)+  , (834300.000000000000000, 10540830.912557534873)+  , (1230000.000000000000000, 16017699.322315014899)+  ]+tableIncompleteBeta :: [(Double,Double,Double,Double)]+tableIncompleteBeta =+  [(2.000000000000000, 3.000000000000000, 0.030000000000000, 0.0051864299999999996862)+  , (2.000000000000000, 3.000000000000000, 0.230000000000000, 0.22845923000000001313)+  , (2.000000000000000, 3.000000000000000, 0.760000000000000, 0.95465728000000005249)+  , (4.000000000000000, 2.300000000000000, 0.890000000000000, 0.93829812158347802864)+  , (1.000000000000000, 1.000000000000000, 0.550000000000000, 0.55000000000000004441)+  , (0.300000000000000, 12.199999999999999, 0.110000000000000, 0.95063000053947077639)+  , (13.100000000000000, 9.800000000000001, 0.120000000000000, 1.3483109941962659385e-07)+  , (13.100000000000000, 9.800000000000001, 0.420000000000000, 0.071321857831804780226)+  , (13.100000000000000, 9.800000000000001, 0.920000000000000, 0.99999578339197081611)+  ]