{-# 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

-- |Error function related functions.
--
-- The derivative of 'erf' is @\ x -> 2 / sqrt pi * exp (x^2)@,
-- and this uniquely determines 'erf' by @erf 0 = 0@.
--
-- Minimal complete definition is 'erfc' or 'normcdf'.
class (Floating a) => Erf a where
    erf :: a -> a
    erfc :: a -> a       -- ^@erfc x = 1 - erf x@
    erfcx :: a -> a      -- ^@erfcx x = exp (x*x) * erfc x@
    normcdf :: a -> a    -- ^@normcdf x = erfc(-x / sqrt 2) / 2@

    -- All the functions are inter-related, here's some defaults.
    erf a
x = a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a -> a
forall a. Erf a => a -> a
erfc a
x
    erfc a
x = a
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 a
2)
    erfcx 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 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 a
2) a -> a -> a
forall a. Fractional a => 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

-- |Inverse error functions, e.g., @inverf . erf = id@ and @erf . inverf = id@ assuming
-- the appropriate codomain for 'inverf'.
-- Note that the accuracy may drop radically for extreme arguments.
class (Floating a) => InvErf a where
    inverf :: a -> a
    inverfc :: a -> a
--    inverfcx :: a -> a
    invnormcdf :: a -> a

    inverf a
p = a -> a
forall a. InvErf a => a -> a
inverfc (a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a
p)
    inverfc a
p = - a -> a
forall a. InvErf a => a -> a
invnormcdf (a
pa -> a -> a
forall a. Fractional a => a -> a -> a
/a
2) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Floating a => a -> a
sqrt a
2

instance InvErf Double where
    invnormcdf :: Double -> Double
invnormcdf Double
0 = -Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0
    invnormcdf Double
1 = Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0
    invnormcdf Double
p =
        -- Do one iteration with Halley's root finder to get a more accurate result.
        let x :: Double
x = Double -> Double
forall a. (Ord a, Floating a) => a -> a
inorm Double
p
            e :: Double
e = Double
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 Double
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 (Double
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
/ Double
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
/ (Double
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
/ Double
2)

instance InvErf Float where
    invnormcdf :: Float -> Float
invnormcdf = Float -> Float
forall a. (Ord a, Floating a) => a -> a
inorm

-- Taken from http://home.online.no/~pjacklam/notes/invnorm/
-- Accurate to about 1e-9.
inorm :: (Ord a, Floating a) => a -> a
inorm :: 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
0a -> a -> a
forall a. Fractional a => a -> a -> a
/a
0

    in  if a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 then
            a
nan
        else if a
p a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then
            -a
1a -> a -> a
forall a. Fractional a => 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(-a
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
+a
1)
        else if a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
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
- 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
+a
1)
        else if a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
1 then
            - a -> a
forall a. (Ord a, Floating a) => a -> a
inorm (a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a
p)
        else
            a
nan