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