diff --git a/Numeric/MathFunctions/Constants.hs b/Numeric/MathFunctions/Constants.hs new file mode 100644 --- /dev/null +++ b/Numeric/MathFunctions/Constants.hs @@ -0,0 +1,85 @@ +-- | +-- Module : Numeric.MathFunctions.Constants +-- Copyright : (c) 2009, 2011 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Constant values common to much numeric code. + +module Numeric.MathFunctions.Constants + ( + m_epsilon + , m_huge + , m_tiny + , m_1_sqrt_2 + , m_2_sqrt_pi + , m_ln_sqrt_2_pi + , m_max_exp + , m_sqrt_2 + , m_sqrt_2_pi + , m_pos_inf + , m_neg_inf + , m_NaN + ) where + +-- | A very large number. +m_huge :: Double +m_huge = 1.7976931348623157e308 +{-# INLINE m_huge #-} + +m_tiny :: Double +m_tiny = 2.2250738585072014e-308 +{-# INLINE m_tiny #-} + +-- | The largest 'Int' /x/ such that 2**(/x/-1) is approximately +-- representable as a 'Double'. +m_max_exp :: Int +m_max_exp = 1024 + +-- | @sqrt 2@ +m_sqrt_2 :: Double +m_sqrt_2 = 1.4142135623730950488016887242096980785696718753769480731766 +{-# INLINE m_sqrt_2 #-} + +-- | @sqrt (2 * pi)@ +m_sqrt_2_pi :: Double +m_sqrt_2_pi = 2.5066282746310005024157652848110452530069867406099383166299 +{-# INLINE m_sqrt_2_pi #-} + +-- | @2 / sqrt pi@ +m_2_sqrt_pi :: Double +m_2_sqrt_pi = 1.1283791670955125738961589031215451716881012586579977136881 +{-# INLINE m_2_sqrt_pi #-} + +-- | @1 / sqrt 2@ +m_1_sqrt_2 :: Double +m_1_sqrt_2 = 0.7071067811865475244008443621048490392848359376884740365883 +{-# INLINE m_1_sqrt_2 #-} + +-- | The smallest 'Double' ε such that 1 + ε ≠ 1. +m_epsilon :: Double +m_epsilon = encodeFloat (signif+1) expo - 1.0 + where (signif,expo) = decodeFloat (1.0::Double) + +-- | @log(sqrt((2*pi))@ +m_ln_sqrt_2_pi :: Double +m_ln_sqrt_2_pi = 0.9189385332046727417803297364056176398613974736377834128171 +{-# INLINE m_ln_sqrt_2_pi #-} + +-- | Positive infinity. +m_pos_inf :: Double +m_pos_inf = 1/0 +{-# INLINE m_pos_inf #-} + +-- | Negative infinity. +m_neg_inf :: Double +m_neg_inf = -1/0 +{-# INLINE m_neg_inf #-} + +-- | Not a number. +m_NaN :: Double +m_NaN = 0/0 +{-# INLINE m_NaN #-} diff --git a/Numeric/Polynomial/Chebyshev.hs b/Numeric/Polynomial/Chebyshev.hs --- a/Numeric/Polynomial/Chebyshev.hs +++ b/Numeric/Polynomial/Chebyshev.hs @@ -1,4 +1,14 @@ {-# LANGUAGE FlexibleContexts #-} +-- | +-- Module : Numeric.Polynomial.Chebyshev +-- Copyright : (c) 2009, 2011 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Chebyshev polynomials. module Numeric.Polynomial.Chebyshev ( -- * Chebyshev polinomials -- $chebyshev diff --git a/Numeric/SpecFunctions.hs b/Numeric/SpecFunctions.hs --- a/Numeric/SpecFunctions.hs +++ b/Numeric/SpecFunctions.hs @@ -1,4 +1,14 @@ {-# LANGUAGE BangPatterns #-} +-- | +-- Module : Numeric.SpecFunctions +-- Copyright : (c) 2009, 2011 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Special functions and factorials. module Numeric.SpecFunctions ( -- * Gamma function logGamma @@ -16,6 +26,7 @@ -- * Factorial , factorial , logFactorial + , stirlingError -- * Combinatorics , choose -- * References @@ -28,7 +39,9 @@ import Data.Number.Erf (erfc) import qualified Data.Vector.Unboxed as U -import Numeric.Polynomial.Chebyshev (chebyshevBroucke) +import Numeric.Polynomial.Chebyshev (chebyshevBroucke) +import Numeric.MathFunctions.Constants (m_epsilon, m_sqrt_2_pi, m_ln_sqrt_2_pi, + m_NaN, m_neg_inf, m_pos_inf, m_sqrt_2) @@ -470,6 +483,43 @@ z = ((-(5.95238095238e-4 * y) + 7.936500793651e-4) * y - 2.7777777777778e-3) * y + 8.3333333333333e-2 +-- | Calculate the error term of the Stirling approximation. This is +-- only defined for non-negative values. +-- +-- > stirlingError @n@ = @log(n!) - log(sqrt(2*pi*n)*(n/e)^n) +stirlingError :: Double -> Double +stirlingError n + | n <= 15.0 = case properFraction (n+n) of + (i,0) -> sfe `U.unsafeIndex` i + _ -> logGamma (n+1.0) - (n+0.5) * log n + n - + m_ln_sqrt_2_pi + | n > 500 = (s0-s1/nn)/n + | n > 80 = (s0-(s1-s2/nn)/nn)/n + | n > 35 = (s0-(s1-(s2-s3/nn)/nn)/nn)/n + | otherwise = (s0-(s1-(s2-(s3-s4/nn)/nn)/nn)/nn)/n + where + nn = n*n + s0 = 0.083333333333333333333 -- 1/12 + s1 = 0.00277777777777777777778 -- 1/360 + s2 = 0.00079365079365079365079365 -- 1/1260 + s3 = 0.000595238095238095238095238 -- 1/1680 + s4 = 0.0008417508417508417508417508 -- 1/1188 + sfe = U.fromList [ 0.0, + 0.1534264097200273452913848, 0.0810614667953272582196702, + 0.0548141210519176538961390, 0.0413406959554092940938221, + 0.03316287351993628748511048, 0.02767792568499833914878929, + 0.02374616365629749597132920, 0.02079067210376509311152277, + 0.01848845053267318523077934, 0.01664469118982119216319487, + 0.01513497322191737887351255, 0.01387612882307074799874573, + 0.01281046524292022692424986, 0.01189670994589177009505572, + 0.01110455975820691732662991, 0.010411265261972096497478567, + 0.009799416126158803298389475, 0.009255462182712732917728637, + 0.008768700134139385462952823, 0.008330563433362871256469318, + 0.007934114564314020547248100, 0.007573675487951840794972024, + 0.007244554301320383179543912, 0.006942840107209529865664152, + 0.006665247032707682442354394, 0.006408994188004207068439631, + 0.006171712263039457647532867, 0.005951370112758847735624416, + 0.005746216513010115682023589, 0.005554733551962801371038690 ] ---------------------------------------------------------------- @@ -507,46 +557,6 @@ max64 = fromIntegral (maxBound :: Int64) round64 x = round x :: Int64 - ----------------------------------------------------------------- --- Constants ----------------------------------------------------------------- - --- | @sqrt 2@ -m_sqrt_2 :: Double -m_sqrt_2 = 1.4142135623730950488016887242096980785696718753769480731766 -{-# INLINE m_sqrt_2 #-} - --- | @sqrt (2 * pi)@ -m_sqrt_2_pi :: Double -m_sqrt_2_pi = 2.5066282746310005024157652848110452530069867406099383166299 -{-# INLINE m_sqrt_2_pi #-} - - --- | The smallest 'Double' ε such that 1 + ε ≠ 1. -m_epsilon :: Double -m_epsilon = encodeFloat (signif+1) expo - 1.0 - where (signif,expo) = decodeFloat (1.0::Double) - --- | @log(sqrt((2*pi))@ -m_ln_sqrt_2_pi :: Double -m_ln_sqrt_2_pi = 0.9189385332046727417803297364056176398613974736377834128171 -{-# INLINE m_ln_sqrt_2_pi #-} - --- | Positive infinity. -m_pos_inf :: Double -m_pos_inf = 1/0 -{-# INLINE m_pos_inf #-} - --- | Negative infinity. -m_neg_inf :: Double -m_neg_inf = -1/0 -{-# INLINE m_neg_inf #-} - --- | Not a number. -m_NaN :: Double -m_NaN = 0/0 -{-# INLINE m_NaN #-} diff --git a/Numeric/SpecFunctions/Extra.hs b/Numeric/SpecFunctions/Extra.hs new file mode 100644 --- /dev/null +++ b/Numeric/SpecFunctions/Extra.hs @@ -0,0 +1,33 @@ +-- | +-- Module : Numeric.SpecFunctions.Extra +-- Copyright : (c) 2009, 2011 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Less common mathematical functions. +module Numeric.SpecFunctions.Extra ( + bd0 + ) where + +import Numeric.MathFunctions.Constants (m_NaN) + +-- | Evaluate the deviance term @x log(x/np) + np - x@. +bd0 :: Double -- ^ @x@ + -> Double -- ^ @np@ + -> Double +bd0 x np + | isInfinite x || isInfinite np || np == 0 = m_NaN + | abs x_np >= 0.1*(x+np) = x * log (x/np) - x_np + | otherwise = loop 1 (ej0*vv) s0 + where + x_np = x - np + v = x_np / (x+np) + s0 = x_np * v + ej0 = 2*x*v + vv = v*v + loop j ej s = case s + ej/(2*j+1) of + s' | s' == s -> s' -- FIXME: Comparing Doubles for equality! + | otherwise -> loop (j+1) (ej*vv) s' diff --git a/math-functions.cabal b/math-functions.cabal --- a/math-functions.cabal +++ b/math-functions.cabal @@ -1,5 +1,5 @@ name: math-functions -version: 0.1.0.0 +version: 0.1.1.0 cabal-version: >= 1.8 license: BSD3 license-file: LICENSE @@ -24,7 +24,9 @@ erf >= 2 exposed-modules: Numeric.SpecFunctions + Numeric.SpecFunctions.Extra Numeric.Polynomial.Chebyshev + Numeric.MathFunctions.Constants test-suite tests type: exitcode-stdio-1.0