{-| 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 :: Double
a1 = -Double
3.969683028665376e+01;
a2 ::  Double
a2 :: Double
a2 =  Double
2.209460984245205e+02;
a3 ::  Double
a3 :: Double
a3 = -Double
2.759285104469687e+02;
a4 ::  Double
a4 :: Double
a4 =  Double
1.383577518672690e+02;
a5 ::  Double
a5 :: Double
a5 = -Double
3.066479806614716e+01;
a6 ::  Double
a6 :: Double
a6 =  Double
2.506628277459239e+00;

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

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

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

-- Limits of the approximation regions (break-points)
xlow ::  Double
xlow :: Double
xlow = Double
0.02425;
xhigh ::  Double
xhigh :: Double
xhigh = Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xlow;

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

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


inverseRecovery :: Double -> Double
inverseRecovery :: Double -> Double
inverseRecovery Double
x
        | Bool
isCloseToZero = Double
0.0
        | Bool
isCloseToOne  = Double
1.0
        | Bool
otherwise     = Double
0.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0.0 -- NaN effectively
        where   isCloseToZero :: Bool
isCloseToZero   = Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
ulp
                isCloseToOne :: Bool
isCloseToOne    = Double
diff Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
tolerance Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Num a => a -> a
abs Double
x Bool -> Bool -> Bool
|| Double
diff Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
tolerance
                diff :: Double
diff            = Double -> Double
forall a. Num a => a -> a
abs (Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1.0)
                tolerance :: Double
tolerance       = Double
42Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ulp
{-# INLINE inverseRecovery #-}

{-# ANN inverseInLowerRegion "NoHerbie" #-}
inverseInLowerRegion ::  Double -> Double
inverseInLowerRegion :: Double -> Double
inverseInLowerRegion Double
x = (((((Double
c1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
c2)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
c3)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
c4)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
c5)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
c6) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((((Double
d1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
d2)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
d3)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
d4)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1.0)
        where   z :: Double
z      = Double -> Double
forall a. Floating a => a -> a
sqrt (-Double
2.0Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double -> Double
forall a. Floating a => a -> a
log Double
x)
{-# INLINE inverseInLowerRegion #-}

{-# ANN inverseInCentralRegion "NoHerbie" #-}
inverseInCentralRegion ::  Double -> Double
inverseInCentralRegion :: Double -> Double
inverseInCentralRegion Double
x = (((((Double
a1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a2)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a3)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a4)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a5)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
a6)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (((((Double
b1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
b2)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
b3)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
b4)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
b5)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1.0)
        where   r :: Double
r       = Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
z
                z :: Double
z       = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5
{-# INLINE inverseInCentralRegion #-}

{-# ANN inverseInHigherRegion "NoHerbie" #-}
inverseInHigherRegion ::  Double -> Double
inverseInHigherRegion :: Double -> Double
inverseInHigherRegion Double
x =  -(((((Double
c1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
c2)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
c3)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
c4)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
c5)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
c6) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((((Double
d1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
d2)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
d3)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
d4)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1.0)
        where   z :: Double
z       = Double -> Double
forall a. Floating a => a -> a
sqrt (-Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log (Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x))
{-# INLINE inverseInHigherRegion #-}