Imports System.Math

Module Acklam

   Public Function SNormInv(ByVal p As Double) As Double

      Dim q As Double
      Dim r As Double

      'Coefficients in rational approximations.

      Const A1 As Double = -39.696830286653757
      Const A2 As Double = 220.9460984245205
      Const A3 As Double = -275.92851044696869
      Const A4 As Double = 138.357751867269
      Const A5 As Double = -30.66479806614716
      Const A6 As Double = 2.5066282774592392

      Const B1 As Double = -54.476098798224058
      Const B2 As Double = 161.58583685804089
      Const B3 As Double = -155.69897985988661
      Const B4 As Double = 66.80131188771972
      Const B5 As Double = -13.280681552885721

      Const C1 As Double = -0.0077848940024302926
      Const C2 As Double = -0.32239645804113648
      Const C3 As Double = -2.4007582771618381
      Const C4 As Double = -2.5497325393437338
      Const C5 As Double = 4.3746641414649678
      Const C6 As Double = 2.9381639826987831

      Const D1 As Double = 0.0077846957090414622
      Const D2 As Double = 0.32246712907003983
      Const D3 As Double = 2.445134137142996
      Const D4 As Double = 3.7544086619074162

      'Define break-points.

      Const P_LOW As Double = 0.02425
      Const P_HIGH As Double = 1 - P_LOW

      If p > 0 AndAlso p < P_LOW Then

         'Rational approximation for lower region.

         q = Sqrt(-2 * Log(p))

         Return (((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6) / _
            ((((D1 * q + D2) * q + D3) * q + D4) * q + 1)

      ElseIf p >= P_LOW AndAlso p <= P_HIGH Then

         'Rational approximation for central region.

         q = p - 0.5
         r = q * q

         Return (((((A1 * r + A2) * r + A3) * r + A4) * r + A5) * r + A6) * q / _
           (((((B1 * r + B2) * r + B3) * r + B4) * r + B5) * r + 1)

      ElseIf p > P_HIGH AndAlso p < 1 Then

         'Rational approximation for upper region.

         q = Sqrt(-2 * Log(1 - p))

         Return -(((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6) / _
                ((((D1 * q + D2) * q + D3) * q + D4) * q + 1)

      Else

         Throw New ArgumentOutOfRangeException()

      End If

   End Function

End Module