equation.bas

实现M5加密算法的源程序

        x(i) = s
    Next i
    k = Abs(x(1))
    For i = 2 To n
        t = Abs(x(i))
        If t > k Then k = t
    Next i
    MultipleMi = k
    For i = 1 To n
        y(i) = x(i) / k
    Next i
Next z
End Function

'反幂法由初始向量y()求出特征向量方向x(),及更精确的特征值TeZheng
'其中用到了带主元选取的Gauss消元法
Function MultipleFan(a() As Double, x() As Double, jama As Double, n As Integer, num As Integer) As Boolean
Dim i As Integer, z As Integer, j As Integer
Dim r As Double, t As Double, l As Double, y() As Double, Sig() As Integer
Dim k As Integer, m() As Integer, q() As Integer, IsSign As Boolean
ReDim y(n), m(n), q(n), Sig(n)
For i = 1 To n
    a(i, i) = a(i, i) - jama
Next i
For i = 1 To n - 1
    k = 0
    For j = 1 To n
        If q(j) = 0 Then
            If Abs(a(j, i)) > Abs(a(k, i)) Then k = j
        End If
    Next j
    If Abs(a(k, i)) < 1E-16 Then Exit Function
    q(k) = i: m(i) = k
    For j = 1 To n
        If q(j) = 0 Then
            l = a(j, i) / a(k, i)
            a(j, i) = l
            For z = i + 1 To n
                a(j, z) = a(j, z) - l * a(k, z)
            Next z
        End If
    Next j
Next i
For i = 1 To n
    If q(i) = 0 Then m(n) = i: q(i) = n: Exit For
Next i
For z = 1 To num
    If z = num Then
        For i = 1 To n
            Sig(i) = Sgn(x(i))
        Next i
    End If
    For i = 1 To n - 1
        For j = 1 To n
            k = m(i)
            If q(j) > i Then
                x(j) = x(j) - a(j, i) * x(k)
            End If
        Next j
    Next i
    k = m(n)
    y(n) = x(k) / a(k, n)
    For i = n - 1 To 1 Step -1
        k = m(i): t = 0
        For j = n To i + 1 Step -1
            t = t + y(j) * a(k, j)
        Next j
        y(i) = (x(k) - t) / a(k, i)
    Next i
    r = Abs(y(1))
    For i = 2 To n
        t = Abs(y(i))
        If t > r Then r = t
    Next i
    If z = num Then
        For i = 1 To n
            If Abs(Sgn(y(i)) - Sig(n)) = 2 Then IsSign = True: Exit For
        Next i
    End If
    For i = 1 To n
        x(i) = y(i) / r
    Next i
Next z
If Not IsSign Then r = -r
jama = jama - 1 / r
MultipleFan = True
End Function

'求特征值
'需对称的Jaccobi法
'有问题
Sub Jaccobi(a() As Double, n As Integer, t As Integer)
Dim z As Integer, i As Integer, j As Integer
Dim r As Double, s(1) As Double
For i = 1 To n - 1
    For j = i + 1 To n
        r = r + a(i, j) * a(i, j)
    Next j
Next i
r = Sqr(r + r)
For z = 1 To t
    For i = 1 To n - 1
        For j = i + 1 To n
            If Abs(a(i, j)) > r Then JaDo a(), i, j, n, False, s()
        Next j
    Next i
    r = r / (n + n)
Next z
End Sub
'返回特征向量s()
'有问题
Sub JaccobiEx(a() As Double, n As Integer, t As Integer, s() As Double)
Dim z As Integer, i As Integer, j As Integer
Dim r As Double
For i = 1 To n
    s(i, i) = 1
Next i
For i = 1 To n - 1
    For j = i + 1 To n
        r = r + a(i, j) * a(i, j)
    Next j
Next i
r = Sqr(r + r)
For z = 1 To t
    For i = 1 To n - 1
        For j = i + 1 To n
            If Abs(a(i, j)) > r Then JaDo a(), i, j, n, True, s()
        Next j
    Next i
    r = r / (n + n)
Next z
End Sub
Sub JaDo(a() As Double, p As Integer, q As Integer, n As Integer, IsVector As Boolean, s() As Double)
Dim t As Double, m As Double, w As Double, bp As Double, bq As Double
Dim k As Integer, i As Integer
t = -a(p, q): m = (a(p, p) - a(q, q)) / 2: k = Sgn(m)
If k = 0 Then
    t = Sqr(2) / 2: m = t
Else
    w = t / Sqr(t * t + m * m)
    w = w * k
    t = Sqr(1 - w * w) + 1
    t = w / Sqr(t + t)
    m = Sqr(1 - t * t)
    'w = m / Sqr(t * t + m * m)
    'm = Sqr((Abs(w) + 1) / 2)
    't = Sqr(1 - m * m)
    'If Sgn(t) <> Sgn(m) Then t = -t
End If
For i = q + 1 To n
    bp = a(p, i) * m - a(q, i) * t
    bq = a(q, i) * m + a(p, i) * t
    a(p, i) = bp: a(q, i) = bq
Next i
For i = p + 1 To q - 1
    a(p, i) = a(p, i) * m - a(i, q) * t
Next i
For i = 1 To p - 1
    bp = a(i, p) * m - a(i, q) * t
    bq = a(i, q) * m + a(i, p) * t
    a(i, p) = bp: a(i, q) = bq
Next i
For i = p + 1 To q - 1
    a(i, q) = a(i, q) * m + a(p, i) * t
Next i
bp = a(p, p) * m * m + a(q, q) * t * t - 2 * a(p, q) * t * m
bq = a(p, p) * t * t + a(q, q) * m * m + 2 * a(p, q) * t * m
a(p, p) = bp: a(q, q) = bq
a(p, q) = 0
If IsVector Then
    For i = 1 To n
        bp = s(i, p) * m - s(i, q) * t
        bq = s(i, p) * t + s(i, q) * m
        s(i, p) = bp: s(i, q) = bq
    Next i
End If
End Sub

'吉文斯-HouseHolder法
'三对角化
'y(n)返回所有特征值
'只用到了a()的下三角部分
Sub HouseHolder(a() As Double, y() As Double, n As Integer, e As Double)
Dim i As Integer, j As Integer, z As Integer
Dim t As Double, b() As Double, c() As Double, s() As Double, p() As Double
Dim af As Double, k As Double
ReDim b(n), c(n), p(n), s(n, n)
For i = 1 To n - 2
    t = 0
    For j = i + 1 To n
        t = t + a(j, i) * a(j, i)
    Next j
    af = t
    t = Sqr(t)
    k = a(i + 1, i)
    s(i, i + 1) = Sgn(k) * t + k
    If Sgn(k) = 0 Then s(i, i + 1) = t
    For j = i + 2 To n
        s(i, j) = a(j, i)
    Next j
    af = 1 / (af + Abs(k) * t)
    For j = i To n
        t = 0
        For z = i + 1 To n
            t = t + a(z, j) * s(i, z)
        Next z
        p(j) = af * t
    Next j
    t = 0
    For j = i + 1 To n
        t = t + s(i, j) * p(j)
    Next j
    t = t * af / 2
    For j = i + 1 To n
        p(j) = p(j) - t * s(i, j)
    Next j
    For j = i + 1 To n
        For z = j To n
            a(z, j) = a(z, j) - s(i, j) * p(z) - s(i, z) * p(j)
        Next z
    Next j
    t = s(i, i) * p(i): c(i) = a(i, i) - t - t
    b(i) = k - s(i, i) * p(i + 1) - s(i, i + 1) * p(i)
    s(n - 1, i) = af
Next i
b(n - 1) = a(n, n - 1)
c(n - 1) = a(n - 1, n - 1): c(n) = a(n, n)
GetJama b(), c(), y(), e, n
End Sub
Sub GetJama(b() As Double, c() As Double, y() As Double, e As Double, n As Integer)
Dim t As Double, min As Double, max As Double, k As Double
Dim u As Double, v As Double, l As Double
Dim i As Integer
t = Abs(b(1))
min = c(1) - t
max = c(1) + t
For i = 2 To n
    t = Abs(b(i)) + Abs(b(i - 1))
    k = c(i) - t: If k < min Then min = k
    k = c(i) + t: If k > max Then max = k
Next i
t = min: u = max
For i = 1 To n
    Do
        v = (l + u) / 2
        If GetJamaNumber(v, b(), c(), n) >= i Then
            l = v
        Else
            u = v
        End If
    Loop Until Abs(u - l) < e
    y(i) = (l + u) / 2
    l = min
Next i
End Sub
Function GetJamaNumber(x As Double, b() As Double, c() As Double, n As Integer) As Integer
Dim ps As Double, pe As Double, pn As Double
Dim sq As Integer, sn As Integer, m As Integer
Dim i As Integer
ps = 1: pe = c(1) - x: sq = Sgn(pe)
If sq <> -1 Then m = m + 1
If sq = 0 Then sq = 1
For i = 2 To n
    pn = (c(i) - x) * pe - b(i - 1) * b(i - 1) * ps
    sn = Sgn(pn)
    If sn = 0 Then
        m = m + 1
    ElseIf sn = sq Then
        m = m + 1
    Else
        sq = sn
    End If
    ps = pe: pe = pn
Next i
GetJamaNumber = m
End Function
'特殊的逆幂法
Function MultipleFanSpe(b() As Double, c() As Double, x() As Double, jama As Double, n As Integer, num As Integer) As Boolean
Dim i As Integer, z As Integer, j As Integer
Dim r As Double, t As Double, l As Double, y() As Double, Sig() As Integer, p() As Double, q() As Double
Dim IsSign As Boolean
ReDim y(n), Sig(n), p(n), q(n)
For i = 1 To n
    c(i) = c(i) - jama
Next i
p(1) = c(1)
For i = 2 To n
    l = b(i - 1) / p(i - 1)
    p(i) = c(i) - l * b(i - 1)
    q(i) = l
Next i
For z = 1 To num
    If z = num Then
        For i = 1 To n
            Sig(i) = Sgn(x(i))
        Next i
    End If
    For i = 2 To n
        x(i) = x(i) - q(i) * x(i - 1)
    Next i
    y(n) = x(n) / p(n)
    For i = n - 1 To 1 Step -1
        y(i) = (x(i) - b(i) * y(i + 1)) / p(i)
    Next i
    r = Abs(y(1))
    For i = 2 To n
        t = Abs(y(i))
        If t > r Then r = t
    Next i
    If z = num Then
        For i = 1 To n
            If Abs(Sgn(y(i)) - Sig(n)) = 2 Then IsSign = True: Exit For
        Next i
    End If
    For i = 1 To n
        x(i) = y(i) / r
    Next i
Next z
If Not IsSign Then r = -r
jama = jama - 1 / r
MultipleFanSpe = True
End Function
'转化特征向量(未调试)
Sub ConvertArr(z() As Double, b() As Double, c() As Double, s() As Double, n As Integer)
Dim i As Integer, j As Integer
Dim af As Double, t As Double, k As Double
For i = n - 2 To 1 Step -1
    af = s(n - 1, i): k = 0
    For j = i + 1 To n
        k = s(i, j) * z(j) + k
    Next j
    For j = i + 1 To n
        z(j) = z(j) - af * s(i, j) * k
    Next j
Next i
End Sub

'QR方法(可能还有问题)
'带原点位移('后为相应语句）
'未加入化上Hessenberg矩阵
Sub QR(a() As Double, n As Integer, total As Integer)
Dim i As Integer, j As Integer, num As Integer
Dim p() As Double, b() As Double, t As Double, s As Double, r As Double, d As Double
Dim m() As Double
Dim k As Integer, z As Integer
ReDim p(n), b(n * (n - 1) / 2), m(n)
For num = 1 To total
    d = a(n, n): a(n, n) = 0
    For i = 1 To n - 1
        a(i, i) = a(i, i) - d
    Next i
    'd = a(n, n): a(n, n) = 0
    'For i = 1 To n - 1
        'a(i, i) = a(i, i) - d
    'Next i
    For i = 1 To n - 1
        s = 0
        For j = 1 To n
            t = a(j, i)
            s = t * t + s
            p(j) = t
        Next j
        m(i) = s
        For j = i + 1 To n
            s = 0
            For z = 1 To n
                s = s + p(z) * a(z, j)
            Next z
            r = s / m(i)
            For z = 1 To n
                a(z, j) = a(z, j) - r * p(z)
            Next z
            k = k + 1
            b(k) = r
        Next j
    Next i
    s = 0
    For i = 1 To n
        s = s + a(i, n) * a(i, n)
    Next i
    m(n) = s
    For j = 1 To n
        k = 1
        For i = 1 To n
            s = 0
            For z = i + 1 To n
                s = s + b(k) * a(z, j)
                k = k + 1
            Next z
            a(i, j) = a(i, j) + s
        Next i
    Next j
    For i = 1 To n
        a(i, i) = a(i, i) + d
        'a(i, i) = a(i, i) + d
        For j = i + 1 To n
            t = Sqr(m(i) / m(j))
            a(i, j) = a(i, j) * t
            a(j, i) = a(j, i) / t
        Next j
    Next i
    k = 0
Next num
End Sub