laguer.txt

常用的数值算法的VB程序

Sub LAGUER(A(), M, X(), EPS, POLISH%)
    Dim ZERO(2), B(2), D(2), F(2), G(2), H(2)
    Dim G2(2), SQ(2), GP(2), GM(2), DX(2), X1(2)
    ZERO(1) = 0#
    ZERO(2) = 0#
    EPSS = 0.00000006
    MAXIT = 100
    DXOLD = CABS(X(1), X(2))
    For ITER = 1 To MAXIT
        B(1) = A(1, M + 1)
        B(2) = A(2, M + 1)
        ERQ = CABS(B(1), X(2))
        D(1) = ZERO(1)
        D(2) = ZERO(2)
        F(1) = ZERO(1)
        F(2) = ZERO(2)
        ABX = CABS(X(1), X(2))
        For J = M To 1 Step -1
            DUM = X(1) * F(1) - X(2) * F(2) + D(1)
            F(2) = X(2) * F(1) + X(1) * F(2) + D(2)
            F(1) = DUM
            DUM = X(1) * D(1) - X(2) * D(2) + B(1)
            D(2) = X(2) * D(1) + X(1) * D(2) + B(2)
            D(1) = DUM
            DUM = X(1) * B(1) - X(2) * B(2) + A(1, J)
            B(2) = X(2) * B(1) + X(1) * B(2) + A(2, J)
            B(1) = DUM
            ERQ = CABS(B(1), B(2)) + ABX * ERQ
        Next J
        ERQ = EPSS * ERQ
        If CABS(B(1), B(2)) <= ERQ Then
            Erase X1, DX, GM, GP, SQ, G2, H, F, D, B, ZERO
            Exit Sub
        Else
            G(1) = CDIV1(D(1), D(2), B(1), B(2))
            G(2) = CDIV2(D(1), D(2), B(1), B(2))
            G2(1) = G(1) * G(1) - G(2) * G(2)
            G2(2) = 2 * G(1) * G(2)
            H(1) = G2(1) - 2 * CDIV1(F(1), F(2), B(1), B(2))
            H(2) = G2(2) - 2 * CDIV2(F(1), F(2), B(1), B(2))
            DUM1 = (M - 1) * (M * H(1) - G2(1))
            DUM2 = (M - 1) * (M * H(2) - G2(2))
            SQ(1) = CSQR1(DUM1, DUM2)
            SQ(2) = CSQR2(DUM1, DUM2)
            GP(1) = G(1) + SQ(1)
            GP(2) = G(2) + SQ(2)
            GM(1) = G(1) - SQ(1)
            GM(2) = G(2) - SQ(2)
            If CABS(GP(1), GP(2)) < CABS(GM(1), GM(2)) Then
                GP(1) = GM(1)
                GP(2) = GM(2)
            End If
            DX(1) = CDIV1(M, 0, GP(1), GP(2))
            DX(2) = CDIV2(M, 0, GP(1), GP(2))
        End If
        X1(1) = X(1) - DX(1)
        X1(2) = X(2) - DX(2)
        If X(1) = X1(1) And X(2) = X1(2) Then
            Erase X1, DX, GM, GP, SQ, G2, H, F, D, B, ZERO
            Exit Sub
        End If
        X(1) = X1(1)
        X(2) = X1(2)
        CDX = CABS(DX(1), DX(2))
        DXOLD = CDX
        If Not POLISH% Then
            If CDX <= EPS * CABS(X(1), X(2)) Then
                Erase X1, DX, GM, GP, SQ, G2, H, F, D, B, ZERO
                Exit Sub
            End If
        End If
    Next ITER
    Print "too many iterations"
End Sub
Function CABS(A1, A2)
    X = Abs(A1)
    Y = Abs(A2)
    If X = 0 Then
        CABS = Y
    ElseIf Y = 0 Then
        CABS = X
    ElseIf X > Y Then
        CABS = X * Sqr(1 + Sqr(Y / X))
    Else
        CABS = Y * Sqr(1 + Sqr(X / Y))
    End If
End Function
Function CDIV1(A1, A2, B1, B2)
    If Abs(B1) >= Abs(B2) Then
        R = B2 / B1
        DEN = B1 + R * B2
        CDIV1 = (A1 + A2 * R) / DEN
    Else
        R = B1 / B2
        DEN = B2 + R * B1
        CDIV1 = (A1 * R + A2) / DEN
    End If
End Function
Function CDIV2(A1, A2, B1, B2)
    If Abs(B1) >= Abs(B2) Then
        R = B2 / B1
        DEN = B1 + R * B2
        CDIV2 = (A2 - A1 * R) / DEN
    Else
        R = B1 / B2
        DEN = B2 + R * B1
        CDIV2 = (A2 * R - A1) / DEN
    End If
End Function
Function CSQR1(X, Y)
    If X = 0 And Y = 0 Then
        U = 0
    Else
        If Abs(X) >= Abs(Y) Then
            W = Sqr(Abs(X)) * Sqr(0.5 * (1 + Sqr(1 + Sqr(Abs(Y / X)))))
        Else
            R = Abs(X / Y)
            W = Sqr(Abs(Y)) * Sqr(0.5 * (R + Sqr(1 + Sqr(R))))
        End If
        If X >= 0 Then
            U = W
            V = Y / (2 * U)
        Else
            If Y >= 0 Then
            V = W
        Else
            V = -W
        End If
        U = Y / (2 * V)
        End If
    End If
    CSQR1 = U
End Function
Function CSQR2(X, Y)
    If X = 0 And Y = 0 Then
        V = 0
    Else
        If Abs(X) >= Abs(Y) Then
            W = Sqr(Abs(X)) * Sqr(0.5 * (1 + Sqr(1 + Sqr(Abs(Y / X)))))
        Else
            R = Abs(X / Y)
            W = Sqr(Abs(Y)) * Sqr(0.5 * (R + Sqr(1 + Sqr(R))))
        End If
        If X >= 0 Then
            U = W
            V = Y / (2 * U)
        Else
            If Y >= 0 Then
            V = W
        Else
            V = -W
        End If
        U = Y / (2 * V)
        End If
    End If
    CSQR2 = V
End Function