d4r27.frm

来自「常用数值算法集哈」· FRM 代码 · 共 155 行

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   2820
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   5490
   LinkTopic       =   "Form1"
   ScaleHeight     =   2820
   ScaleWidth      =   5490
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3840
      TabIndex        =   0
      Top             =   2040
      Width           =   1215
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Private Sub Command1_ClicK()
    'PROGRAM D4R27
    'Driver for routine BESJAN
    Dim ACTUAL(5), F(5)
    ACTUAL(0) = 0.540974
    ACTUAL(1) = 0.0917017
    ACTUAL(2) = 0.00923641
    ACTUAL(3) = 0.000662379
    ACTUAL(4) = 0.0000368921
    A = 0.5
    IH = 1
    X = 0.5
    Print
    Print "   N     A     IH        X          Actual       BESJAN"
    For I = 0 To 4
        Call BESJAN(X, A, 4, IH, F())
        Print Tab(4); Format$(I, "#0");
        Print Tab(9); Format$(A, "0.#");
        Print Tab(16); Format$(IH, "#");
        Print Tab(22); Format$(X, "0.000000");
        Print Tab(34); Format$(ACTUAL(I), ".000000E+00");
        Print Tab(48); Format$(F(I), ".000000E+00");
    Next I
End Sub
Sub BESJAN(X, A, NM, IH, F())
    If IH = 1 Then
        Call JAP(X, A, NM, F())
        Exit Sub
    Else
        Call JAP(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 JAP(X, A, NMAX, F())
    MMAX = 10
    EPS = 0.00005
    E = 4
    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#)) / Sum
    D1 = 2.3026 * E + 1.3863
    If NMAX > 0 Then
        R = T(0.5 * D1 / NMAX) * NMAX
    Else
        R = 0#
    End If
    S = T(0.73576 * D1 / X) * 1.3591 * X
    If R <= S Then
        NU = 1 + Int(S)
    Else
        NU = 1 + Int(R)
    End If
1   M = 0
    AL = 1#
    LI = Int(NU / 2)
2   M = M + 1
    AL = AL * (M + A) / (M + 1)
    If M < LI Then GoTo 2
    N = 2 * M
    R = 0#
    S = 0#
3   R = 1# / (2# * (A + N) / X - R)
    If Int(N / 2) <> N / 2# Then
        AM = 0#
    Else
        AL = AL * (N + 2) / (N + 2# * A)
        AM = AL * (N + A)
    End If
    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
        BBB = 1
        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
   Function GAMMLN(XX)
      Dim COF(6)
      COF(1) = 76.18009173
      COF(2) = -86.50532033
      COF(3) = 24.01409822
      COF(4) = -1.231739516
      COF(5) = 0.00120858003
      COF(6) = -0.00000536382
      STP = 2.50662827465
      HALF = 0.5
      ONE = 1#
      FPF = 5.5
      X = XX - ONE
      TMP = X + FPF
      TMP = (X + HALF) * Log(TMP) - TMP
      SER = ONE
      For J = 1 To 6
          X = X + ONE
          SER = SER + COF(J) / X
    Next J
      GAMMLN = TMP + Log(STP * SER)
   End Function