{-# 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 a
x = a
1 forall a. Num a => a -> a -> a
- forall a. Erf a => a -> a
erfc a
x
erfc a
x = a
2 forall a. Num a => a -> a -> a
* forall a. Erf a => a -> a
normcdf (-a
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt a
2)
erfcx a
x = forall a. Floating a => a -> a
exp (a
xforall a. Num a => a -> a -> a
*a
x) forall a. Num a => a -> a -> a
* forall a. Erf a => a -> a
erfc a
x
normcdf a
x = forall a. Erf a => a -> a
erfc(-a
x forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt a
2) forall a. Fractional a => a -> a -> a
/ a
2
instance Erf Double where
erf :: Double -> Double
erf = forall a b. (Real a, Fractional b) => a -> b
realToFrac forall b c a. (b -> c) -> (a -> b) -> a -> c
. CDouble -> CDouble
c_erf forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (Real a, Fractional b) => a -> b
realToFrac
erfc :: Double -> Double
erfc = forall a b. (Real a, Fractional b) => a -> b
realToFrac forall b c a. (b -> c) -> (a -> b) -> a -> c
. CDouble -> CDouble
c_erfc forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (Real a, Fractional b) => a -> b
realToFrac
instance Erf Float where
erf :: Float -> Float
erf = forall a b. (Real a, Fractional b) => a -> b
realToFrac forall b c a. (b -> c) -> (a -> b) -> a -> c
. CFloat -> CFloat
c_erff forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (Real a, Fractional b) => a -> b
realToFrac
erfc :: Float -> Float
erfc = forall a b. (Real a, Fractional b) => a -> b
realToFrac forall b c a. (b -> c) -> (a -> b) -> a -> c
. CFloat -> CFloat
c_erfcf forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 a
p = forall a. InvErf a => a -> a
inverfc (a
1 forall a. Num a => a -> a -> a
- a
p)
inverfc a
p = - forall a. InvErf a => a -> a
invnormcdf (a
pforall a. Fractional a => a -> a -> a
/a
2) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt a
2
instance InvErf Double where
invnormcdf :: Double -> Double
invnormcdf Double
0 = -Double
1forall a. Fractional a => a -> a -> a
/Double
0
invnormcdf Double
1 = Double
1forall a. Fractional a => a -> a -> a
/Double
0
invnormcdf Double
p =
let x :: Double
x = forall a. (Ord a, Floating a) => a -> a
inorm Double
p
e :: Double
e = Double
0.5 forall a. Num a => a -> a -> a
* forall a. Erf a => a -> a
erfc (-Double
x forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt Double
2) forall a. Num a => a -> a -> a
- Double
p
u :: Double
u = Double
e forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt (Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
pi) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp (Double
xforall a. Num a => a -> a -> a
*Double
x forall a. Fractional a => a -> a -> a
/ Double
2)
in Double
x forall a. Num a => a -> a -> a
- Double
u forall a. Fractional a => a -> a -> a
/ (Double
1 forall a. Num a => a -> a -> a
+ Double
x forall a. Num a => a -> a -> a
* Double
u forall a. Fractional a => a -> a -> a
/ Double
2)
instance InvErf Float where
invnormcdf :: Float -> Float
invnormcdf = forall a. (Ord a, Floating a) => a -> a
inorm
inorm :: (Ord a, Floating a) => a -> a
inorm :: forall a. (Ord a, Floating a) => a -> a
inorm a
p =
let a1 :: a
a1 = -a
3.969683028665376e+01
a2 :: a
a2 = a
2.209460984245205e+02
a3 :: a
a3 = -a
2.759285104469687e+02
a4 :: a
a4 = a
1.383577518672690e+02
a5 :: a
a5 = -a
3.066479806614716e+01
a6 :: a
a6 = a
2.506628277459239e+00
b1 :: a
b1 = -a
5.447609879822406e+01
b2 :: a
b2 = a
1.615858368580409e+02
b3 :: a
b3 = -a
1.556989798598866e+02
b4 :: a
b4 = a
6.680131188771972e+01
b5 :: a
b5 = -a
1.328068155288572e+01
c1 :: a
c1 = -a
7.784894002430293e-03
c2 :: a
c2 = -a
3.223964580411365e-01
c3 :: a
c3 = -a
2.400758277161838e+00
c4 :: a
c4 = -a
2.549732539343734e+00
c5 :: a
c5 = a
4.374664141464968e+00
c6 :: a
c6 = a
2.938163982698783e+00
d1 :: a
d1 = a
7.784695709041462e-03
d2 :: a
d2 = a
3.224671290700398e-01
d3 :: a
d3 = a
2.445134137142996e+00
d4 :: a
d4 = a
3.754408661907416e+00
pLow :: a
pLow = a
0.02425
nan :: a
nan = a
0forall a. Fractional a => a -> a -> a
/a
0
in if a
p forall a. Ord a => a -> a -> Bool
< a
0 then
a
nan
else if a
p forall a. Eq a => a -> a -> Bool
== a
0 then
-a
1forall a. Fractional a => a -> a -> a
/a
0
else if a
p forall a. Ord a => a -> a -> Bool
< a
pLow then
let q :: a
q = forall a. Floating a => a -> a
sqrt(-a
2forall a. Num a => a -> a -> a
*forall a. Floating a => a -> a
log(a
p))
in (((((a
c1forall a. Num a => a -> a -> a
*a
qforall a. Num a => a -> a -> a
+a
c2)forall a. Num a => a -> a -> a
*a
qforall a. Num a => a -> a -> a
+a
c3)forall a. Num a => a -> a -> a
*a
qforall a. Num a => a -> a -> a
+a
c4)forall a. Num a => a -> a -> a
*a
qforall a. Num a => a -> a -> a
+a
c5)forall a. Num a => a -> a -> a
*a
qforall a. Num a => a -> a -> a
+a
c6) forall a. Fractional a => a -> a -> a
/
((((a
d1forall a. Num a => a -> a -> a
*a
qforall a. Num a => a -> a -> a
+a
d2)forall a. Num a => a -> a -> a
*a
qforall a. Num a => a -> a -> a
+a
d3)forall a. Num a => a -> a -> a
*a
qforall a. Num a => a -> a -> a
+a
d4)forall a. Num a => a -> a -> a
*a
qforall a. Num a => a -> a -> a
+a
1)
else if a
p forall a. Ord a => a -> a -> Bool
< a
1 forall a. Num a => a -> a -> a
- a
pLow then
let q :: a
q = a
p forall a. Num a => a -> a -> a
- a
0.5
r :: a
r = a
qforall a. Num a => a -> a -> a
*a
q
in (((((a
a1forall a. Num a => a -> a -> a
*a
rforall a. Num a => a -> a -> a
+a
a2)forall a. Num a => a -> a -> a
*a
rforall a. Num a => a -> a -> a
+a
a3)forall a. Num a => a -> a -> a
*a
rforall a. Num a => a -> a -> a
+a
a4)forall a. Num a => a -> a -> a
*a
rforall a. Num a => a -> a -> a
+a
a5)forall a. Num a => a -> a -> a
*a
rforall a. Num a => a -> a -> a
+a
a6)forall a. Num a => a -> a -> a
*a
q forall a. Fractional a => a -> a -> a
/
(((((a
b1forall a. Num a => a -> a -> a
*a
rforall a. Num a => a -> a -> a
+a
b2)forall a. Num a => a -> a -> a
*a
rforall a. Num a => a -> a -> a
+a
b3)forall a. Num a => a -> a -> a
*a
rforall a. Num a => a -> a -> a
+a
b4)forall a. Num a => a -> a -> a
*a
rforall a. Num a => a -> a -> a
+a
b5)forall a. Num a => a -> a -> a
*a
rforall a. Num a => a -> a -> a
+a
1)
else if a
p forall a. Ord a => a -> a -> Bool
<= a
1 then
- forall a. (Ord a, Floating a) => a -> a
inorm (a
1 forall a. Num a => a -> a -> a
- a
p)
else
a
nan