besian.txt

常用的数值算法的VB程序

Sub BESIAN(X, A, NM, IH, F())
    If IH = 1 Then
        Call IAP(X, A, NM, F())
        Exit Sub
    Else
        Call IAP(X, A, 1, F())
        F(1) = 2# * A * F(0) / X + F(1)
        For I = 1 To NM - 1
            F(I + 1) = 2# * (A - I) * F(I) / X + F(I - 1)
        Next I
        Exit Sub
    End If
End Sub
Sub IAP(X, A, NMAX, F())
    MMAX = 10
    E = 4
    EPS = 0.5 * 10 ^ (-E)
    Dim FA(10), RR(10)
    For N = 0 To NMAX
        FA(N) = 0#
    Next N
    Sum = Exp(GAMMLN(1# + A))
    Sum = Exp(A * Log(X / 2#)) * Exp(X) / Sum
    D1 = 2.3026 * E + 1.3863
    If NMAX > 0 Then
        R = T(0.5 * D1 / NMAX) * NMAX
    Else
        R = 0#
    End If
    If X < D1 Then
        S = T(0.73576 * (D1 - X) / X) * 1.3591 * X
    Else
        S = 1.3591 * X
    End If
    If R <= S Then
        NU = 1 + Int(S)
    Else
        NU = 1 + Int(R)
    End If
1   N = 0
    AL = 1#
2   N = N + 1
    AL = AL * (N + 2# * A) / (N + 1#)
    If N < NU Then GoTo 2
    R = 0#
    S = 0#
3   R = 1# / (2# * (A + N) / X + R)
    AL = AL * (N + 1#) / (N + 2# * A)
    AM = 2# * (N + A) * AL
    S = R * (AM + S)
    If N <= NMAX Then RR(N - 1) = R
    N = N - 1
    If N >= 1 Then GoTo 3
    F(0) = Sum / (1# + S)
    For N = 0 To NMAX - 1
        F(N + 1) = RR(N) * F(N)
    Next N
    For N = 0 To NMAX
        If Abs((F(N) - FA(N)) / F(N)) > EPS Then
            For M = 0 To NMAX
                FA(M) = F(M)
            Next M
            NU = NU + 5
            GoTo 1
        End If
    Next N
End Sub
Function T(Y)
    If Y <= 10# Then
        AAA = 0.000057941 * Y - 0.00176148 * Y + 0.0208645
        T = ((AAA * Y - 0.129013) * Y + 0.85777) * Y + 1.0125
    Else
        Z = Log(Y) - 0.775
        P = (0.775 - Log(Z)) / (1# + Z)
        T = Y / (P + 1#) / Z
    End If
End Function