{-| The inverse normal cumulative distribution is a non-linear function for which no closed-form solution exists. The function is continuous, monotonically increasing, infinitely differentiable, and maps the open interval (0,1) to the whole real line. By <a href="http://home.online.no/~pjacklam/notes/invnorm/"> An algorithm for computing the inverse normal cumulative distribution function</a>
-}

module QuantLib.Math.InverseNormal
        ( inverseNormal
        ) where


a1 ::  Double
a1 = -3.969683028665376e+01;
a2 ::  Double
a2 =  2.209460984245205e+02;
a3 ::  Double
a3 = -2.759285104469687e+02;
a4 ::  Double
a4 =  1.383577518672690e+02;
a5 ::  Double
a5 = -3.066479806614716e+01;
a6 ::  Double
a6 =  2.506628277459239e+00;

b1 ::  Double
b1 = -5.447609879822406e+01;
b2 ::  Double
b2 =  1.615858368580409e+02;
b3 ::  Double
b3 = -1.556989798598866e+02;
b4 ::  Double
b4 =  6.680131188771972e+01;
b5 ::  Double
b5 = -1.328068155288572e+01;

c1 ::  Double
c1 = -7.784894002430293e-03;
c2 ::  Double
c2 = -3.223964580411365e-01;
c3 ::  Double
c3 = -2.400758277161838e+00;
c4 ::  Double
c4 = -2.549732539343734e+00;
c5 ::  Double
c5 =  4.374664141464968e+00;
c6 ::  Double
c6 =  2.938163982698783e+00;

d1 ::  Double
d1 =  7.784695709041462e-03;
d2 ::  Double
d2 =  3.224671290700398e-01;
d3 ::  Double
d3 =  2.445134137142996e+00;
d4 ::  Double
d4 =  3.754408661907416e+00;

-- Limits of the approximation regions (break-points)
xlow ::  Double
xlow = 0.02425;
xhigh ::  Double
xhigh = 1.0 - xlow;

-- Precision of Double at 1.0 point
ulp ::  Double
ulp = 2.220446049250313E-16

-- | Computes the inverse cumulative standard normal distribution N(0, 1)
inverseNormal ::  Double -> Double
inverseNormal x
        | x < xlow      = inverseInLowerRegion z
        | x <= xhigh    = inverseInCentralRegion z
        | otherwise     = inverseInHigherRegion z
        where   z | x < 0.0 || x >1.0   = inverseRecovery x
                  | otherwise           = x


inverseRecovery :: Double -> Double
inverseRecovery x
        | isCloseToZero = 0.0
        | isCloseToOne  = 1.0
        | otherwise     = 0.0/0.0 -- NaN effectively
        where   isCloseToZero   = abs x < ulp
                isCloseToOne    = diff <= tolerance * abs x || diff <= tolerance
                diff            = abs (x-1.0)
                tolerance       = 42*ulp
{-# INLINE inverseRecovery #-}

{-# ANN inverseInLowerRegion "NoHerbie" #-}
inverseInLowerRegion ::  Double -> Double
inverseInLowerRegion x = (((((c1*z+c2)*z+c3)*z+c4)*z+c5)*z+c6) / ((((d1*z+d2)*z+d3)*z+d4)*z+1.0)
        where   z      = sqrt (-2.0*log x)
{-# INLINE inverseInLowerRegion #-}

{-# ANN inverseInCentralRegion "NoHerbie" #-}
inverseInCentralRegion ::  Double -> Double
inverseInCentralRegion x = (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*z / (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1.0)
        where   r       = z*z
                z       = x - 0.5
{-# INLINE inverseInCentralRegion #-}

{-# ANN inverseInHigherRegion "NoHerbie" #-}
inverseInHigherRegion ::  Double -> Double
inverseInHigherRegion x =  -(((((c1*z+c2)*z+c3)*z+c4)*z+c5)*z+c6) / ((((d1*z+d2)*z+d3)*z+d4)*z+1.0)
        where   z       = sqrt (-2.0 * log (1.0 - x))
{-# INLINE inverseInHigherRegion #-}