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;
xlow :: Double
xlow = 0.02425;
xhigh :: Double
xhigh = 1.0 - xlow;
ulp :: Double
ulp = 2.220446049250313E-16
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
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 #-}