{-# LANGUAGE ForeignFunctionInterface #-}
module Data.Number.Erf(Erf(..), InvErf(..)) where
import Foreign.C
foreign import ccall "erf" c_erf :: CDouble -> CDouble
foreign import ccall "erfc" c_erfc :: CDouble -> CDouble
foreign import ccall "erff" c_erff :: CFloat -> CFloat
foreign import ccall "erfcf" c_erfcf :: CFloat -> CFloat
class (Floating a) => Erf a where
erf :: a -> a
erfc :: a -> a
erfcx :: a -> a
normcdf :: a -> a
erf x :: a
x = 1 a -> a -> a
forall a. Num a => a -> a -> a
- a -> a
forall a. Erf a => a -> a
erfc a
x
erfc x :: a
x = 2 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Erf a => a -> a
normcdf (-a
x a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
sqrt 2)
erfcx x :: a
x = a -> a
forall a. Floating a => a -> a
exp (a
xa -> a -> a
forall a. Num a => a -> a -> a
*a
x) a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Erf a => a -> a
erfc a
x
normcdf x :: a
x = a -> a
forall a. Erf a => a -> a
erfc(-a
x a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Floating a => a -> a
sqrt 2) a -> a -> a
forall a. Fractional a => a -> a -> a
/ 2
instance Erf Double where
erf :: Double -> Double
erf = CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (CDouble -> Double) -> (Double -> CDouble) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CDouble -> CDouble
c_erf (CDouble -> CDouble) -> (Double -> CDouble) -> Double -> CDouble
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac
erfc :: Double -> Double
erfc = CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (CDouble -> Double) -> (Double -> CDouble) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CDouble -> CDouble
c_erfc (CDouble -> CDouble) -> (Double -> CDouble) -> Double -> CDouble
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac
instance Erf Float where
erf :: Float -> Float
erf = CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac (CFloat -> Float) -> (Float -> CFloat) -> Float -> Float
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CFloat -> CFloat
c_erff (CFloat -> CFloat) -> (Float -> CFloat) -> Float -> CFloat
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Float -> CFloat
forall a b. (Real a, Fractional b) => a -> b
realToFrac
erfc :: Float -> Float
erfc = CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac (CFloat -> Float) -> (Float -> CFloat) -> Float -> Float
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CFloat -> CFloat
c_erfcf (CFloat -> CFloat) -> (Float -> CFloat) -> Float -> CFloat
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Float -> CFloat
forall a b. (Real a, Fractional b) => a -> b
realToFrac
class (Floating a) => InvErf a where
inverf :: a -> a
inverfc :: a -> a
invnormcdf :: a -> a
inverf p :: a
p = a -> a
forall a. InvErf a => a -> a
inverfc (1 a -> a -> a
forall a. Num a => a -> a -> a
- a
p)
inverfc p :: a
p = - a -> a
forall a. InvErf a => a -> a
invnormcdf (a
pa -> a -> a
forall a. Fractional a => a -> a -> a
/2) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Floating a => a -> a
sqrt 2
instance InvErf Double where
invnormcdf :: Double -> Double
invnormcdf 0 = -1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/0
invnormcdf 1 = 1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/0
invnormcdf p :: Double
p =
let x :: Double
x = Double -> Double
forall a. (Ord a, Floating a) => a -> a
inorm Double
p
e :: Double
e = 0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Erf a => a -> a
erfc (-Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt 2) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p
u :: Double
u = Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp (Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ 2)
in Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ 2)
instance InvErf Float where
invnormcdf :: Float -> Float
invnormcdf = Float -> Float
forall a. (Ord a, Floating a) => a -> a
inorm
inorm :: (Ord a, Floating a) => a -> a
inorm :: a -> a
inorm p :: a
p =
let a1 :: a
a1 = -3.969683028665376e+01
a2 :: a
a2 = 2.209460984245205e+02
a3 :: a
a3 = -2.759285104469687e+02
a4 :: a
a4 = 1.383577518672690e+02
a5 :: a
a5 = -3.066479806614716e+01
a6 :: a
a6 = 2.506628277459239e+00
b1 :: a
b1 = -5.447609879822406e+01
b2 :: a
b2 = 1.615858368580409e+02
b3 :: a
b3 = -1.556989798598866e+02
b4 :: a
b4 = 6.680131188771972e+01
b5 :: a
b5 = -1.328068155288572e+01
c1 :: a
c1 = -7.784894002430293e-03
c2 :: a
c2 = -3.223964580411365e-01
c3 :: a
c3 = -2.400758277161838e+00
c4 :: a
c4 = -2.549732539343734e+00
c5 :: a
c5 = 4.374664141464968e+00
c6 :: a
c6 = 2.938163982698783e+00
d1 :: a
d1 = 7.784695709041462e-03
d2 :: a
d2 = 3.224671290700398e-01
d3 :: a
d3 = 2.445134137142996e+00
d4 :: a
d4 = 3.754408661907416e+00
pLow :: a
pLow = 0.02425
nan :: a
nan = 0a -> a -> a
forall a. Fractional a => a -> a -> a
/0
in if a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< 0 then
a
nan
else if a
p a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== 0 then
-1a -> a -> a
forall a. Fractional a => a -> a -> a
/0
else if a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
pLow then
let q :: a
q = a -> a
forall a. Floating a => a -> a
sqrt(-2a -> a -> a
forall a. Num a => a -> a -> a
*a -> a
forall a. Floating a => a -> a
log(a
p))
in (((((a
c1a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
c2)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
c3)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
c4)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
c5)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
c6) a -> a -> a
forall a. Fractional a => a -> a -> a
/
((((a
d1a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
d2)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
d3)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
d4)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+1)
else if a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< 1 a -> a -> a
forall a. Num a => a -> a -> a
- a
pLow then
let q :: a
q = a
p a -> a -> a
forall a. Num a => a -> a -> a
- 0.5
r :: a
r = a
qa -> a -> a
forall a. Num a => a -> a -> a
*a
q
in (((((a
a1a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
a2)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
a3)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
a4)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
a5)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
a6)a -> a -> a
forall a. Num a => a -> a -> a
*a
q a -> a -> a
forall a. Fractional a => a -> a -> a
/
(((((a
b1a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
b2)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
b3)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
b4)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
b5)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+1)
else if a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= 1 then
- a -> a
forall a. (Ord a, Floating a) => a -> a
inorm (1 a -> a -> a
forall a. Num a => a -> a -> a
- a
p)
else
a
nan