module Numeric.SpecFunctions.Extra (
bd0
, chooseExact
, logChooseFast
) where
import Numeric.MathFunctions.Constants (m_NaN)
import Numeric.SpecFunctions.Internal (chooseExact,logChooseFast)
bd0 :: Double
-> Double
-> 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'
| otherwise -> loop (j+1) (ej*vv) s'