{-| 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 An algorithm for computing the inverse normal cumulative distribution function -} 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 #-}