Excel函数NormSDist和NormSInv的VB实现

最近在写一个数值计算程序，为了省事，祭起了久违的VB6.0。算法已经在Excel中实现，但VB中缺了两个数学函数，标准正态累积分布NormSDist及其反函数NormSInv，一番搜索，折腾出近似算法如下，精度高于10E-7。

Function NormSDist(ByVal a As Double) As Double
    '网上找的多有bug，此算法修改诸多才通，原作者已不可考
    Const pi As Double = 3.14159265358979
    Dim p, b1, b2, b3, b4, b5 As Double
    
    p = 0.2316419
    b1 = 0.31938153
    b2 = -0.356563782
    b3 = 1.781477937
    b4 = -1.821255978
    b5 = 1.330274429

    x = Abs(a)
    t = 1 / (1 + p * x)
    
    NormSDist = 1 - (1 / (Sqr(2 * pi)) * Exp(-a ^ 2 / 2)) * (b1 * t + b2 * t ^ 2 + b3 * t ^ 3 + b4 * t ^ 4 + b5 * t ^ 5)
    
    If a < 0 Then
        NormSDist = 1 - NormSDist
    End If
    
End Function

Function NormSInv(ByVal p As Double) As Double
    ' It uses the algorithm of Peter J. Acklam to compute the inverse normal cumulative
    ' distribution. Refer to http://home.online.no/~pjacklam/notes/invnorm/index.html for
    ' a description of the algorithm.
    ' Adapted to VB by Christian d'Heureuse, http://www.source-code.biz

    Const a1 = -39.6968302866538, a2 = 220.946098424521, a3 = -275.928510446969
    Const a4 = 138.357751867269, a5 = -30.6647980661472, a6 = 2.50662827745924
    Const b1 = -54.4760987982241, b2 = 161.585836858041, b3 = -155.698979859887
    Const b4 = 66.8013118877197, b5 = -13.2806815528857, c1 = -7.78489400243029E-03
    Const c2 = -0.322396458041136, c3 = -2.40075827716184, c4 = -2.54973253934373
    Const c5 = 4.37466414146497, c6 = 2.93816398269878, d1 = 7.78469570904146E-03
    Const d2 = 0.32246712907004, d3 = 2.445134137143, d4 = 3.75440866190742
    Const p_low = 0.02425, p_high = 1 - p_low
   
    Dim q As Double, r As Double
   
    If p < 0 Or p > 1 Then
        Err.Raise vbObjectError, , "NormSInv: Argument out of range."
    ElseIf p < p_low Then
        q = Sqr(-2 * Log(p))
        NormSInv = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1)
    ElseIf p <= p_high Then
        q = p - 0.5: r = q * q
        NormSInv = (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q / (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1)
    Else
        q = Sqr(-2 * Log(1 - p))
        NormSInv = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1)
    End If

End Function

以上。
隔几年就要说不弹此调久矣的老狼