bnldev.txt

常用的数值算法的VB程序

Function BNLDEV(PP, N, IDUM&)
    PI = 3.141592654
    If PP <= 0.5 Then
        P = PP
    Else
        P = 1# - PP
    End If
    AM = N * P
    If N < 25 Then
        BNL = 0#
        For J = 1 To N
            If RAN1(IDUM&) < P Then BNL = BNL + 1#
        Next J
    ElseIf AM < 1# Then
        G = Exp(-AM)
        T = 1#
        For J = 0 To N
            T = T * RAN1(IDUM&)
            If T < G Then Exit For
        Next J
        If T >= G Then J = N
        BNL = J
    Else
        If N <> NOLD Then
            EN = N
            OLDG = GAMMLN(EN + 1#)
            NOLD = N
        End If
        If P <> POLD Then
            PC = 1# - P
            PLOG = Log(P)
            PCLOG = Log(PC)
            POLD = P
        End If
        SQ = Sqr(2# * AM * PC)
        Do
          Do
            Y = Tan(PI * RAN1(IDUM&))
            EM = SQ * Y + AM
          Loop While (EM < 0) Or (EM >= EN + 1#)
          EM = Int(EM)
          T = EN - EM
          T = Exp(OLDG - GAMMLN(EM + 1#) - GAMMLN(T + 1#) + EM * PLOG + T * PCLOG)
          T = 1.2 * SQ * (1# + Y ^ 2) * T
        Loop While RAN1(IDUM&) > T
        BNL = EM
    End If
    If P <> PP Then BNL = N - BNL
    BNLDEV = BNL
End Function