equation.bas

实现M5加密算法的源程序

Attribute VB_Name = "Equation"
Option Explicit

'尚需调试
Sub ReadMatrix(FileName As String, a() As Double, Cols As Integer, Rows As Integer)
Dim i As Integer, j As Integer
Open FileName For Input As #1
For i = 1 To Cols
    For j = 1 To Rows
        Input #1, a(i, j)
    Next j
Next i
End Sub

Function OutPutMatrix(a() As Double, n As Integer) As String
Dim i As Integer, j As Integer
Dim res As String
For i = 1 To n
    res = res & a(i, 1)
    For j = 2 To n
        res = res & "," & a(i, j)
    Next j
    res = res & vbCrLf
Next i
OutPutMatrix = res
End Function

Sub EquMulti(a() As Double, b() As Double, c() As Double, m As Integer, n As Integer, l As Integer)
Dim i As Integer, j As Integer, z As Integer
Dim s As Double
For i = 1 To m
    For j = 1 To l
        For z = 1 To n
            s = s + a(i, z) * b(z, j)
        Next z
        c(i, j) = s: s = 0
    Next j
Next i
End Sub

'Gauss消元法(不含主元选取)
Sub EquGauss(a() As Double, y() As Double, x() As Double, n As Integer)
Dim i As Integer, j As Integer, z As Integer
Dim l As Double
For i = 1 To n - 1
    For j = i + 1 To n
        l = a(j, i) / a(i, i)
        y(j) = y(j) - l * y(i)
        For z = i + 1 To n
            a(j, z) = a(j, z) - l * a(i, z)
        Next z
    Next j
Next i
x(n) = y(n) / a(n, n)
For i = n - 1 To 1 Step -1
    l = 0
    For j = n To i + 1 Step -1
        l = l + x(j) * a(i, j)
    Next j
    x(i) = (y(i) - l) / a(i, i)
Next i
End Sub

'Gauss消元法(含主元选取)
Function EquGaussEx(a() As Double, y() As Double, x() As Double, n As Integer) As Boolean
Dim i As Integer, j As Integer, z As Integer, k As Integer, t() As Integer, m() As Integer
ReDim t(n), m(n)
Dim w As Double, l As Double
For i = 1 To n - 1
    k = 0
    For j = 1 To n
        If t(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
    t(k) = i: m(i) = k
    For j = 1 To n
        If t(j) = 0 Then
            l = a(j, i) / a(k, i)
            y(j) = y(j) - l * y(k)
            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 t(i) = 0 Then k = i: Exit For
Next i
m(n) = k
x(n) = y(k) / a(k, n)
For i = n - 1 To 1 Step -1
    k = m(i): w = 0
    For j = n To i + 1 Step -1
        w = w + x(j) * a(k, j)
    Next j
    x(i) = (y(k) - w) / a(k, i)
Next i
EquGaussEx = True
End Function

'对称正定矩阵的Cholesky法
Sub Cholesky(a() As Double, y() As Double, x() As Double, n As Integer)
Dim m As Double, t As Double
Dim i As Integer, j As Integer, z As Integer
m = a(1, 1)
For i = 2 To n
    a(i, 1) = a(i, 1) / m
Next i
For i = 2 To n
    For j = i To n
        t = 0
        For z = 1 To i - 1
            t = t + a(i, z) * a(z, j)
        Next z
        a(i, j) = a(i, j) - t
    Next j
    m = a(i, i)
    t = 0
    For z = 1 To i - 1
        t = t + y(z) * a(i, z)
    Next z
    y(i) = y(i) - t
    For j = i + 1 To n
        a(j, i) = a(i, j) / m
    Next j
Next i
x(n) = y(n) / a(n, n)
For i = n - 1 To 1 Step -1
    t = 0
    For j = n To i + 1 Step -1
        t = t + x(j) * a(i, j)
    Next j
    x(i) = (y(i) - t) / a(i, i)
Next i
End Sub

'求矩阵逆(不含主元选取)
Sub GetEquNi(a() As Double, n As Integer)
Dim i As Integer, j As Integer, z As Integer
Dim d As Double
For i = 1 To n
    d = a(i, i)
    For j = 1 To n
        If j <> i Then a(j, i) = -a(j, i) / d
    Next j
    For j = 1 To n
        For z = 1 To n
            If j <> i And z <> i Then a(j, z) = a(j, z) + a(i, z) * a(j, i)
        Next z
    Next j
    For j = 1 To n
        If j <> i Then a(i, j) = a(i, j) / d
    Next j
Next i
End Sub

'求矩阵逆(含主元选取)
Function GetEquNiEx(a() As Double, b() As Double, n As Integer) As Boolean
Dim i As Integer, j As Integer, z As Integer, k As Integer, t() As Integer, m() As Integer
ReDim t(n), m(n)
Dim d As Double
For i = 1 To n
    k = 0
    For j = 1 To n
        If t(j) = 0 Then
            If Abs(a(j, i)) > Abs(a(k, i)) Then k = j
        End If
    Next j
    t(k) = i: m(i) = k: d = a(k, i)
    If Abs(d) < 1E-16 Then Exit Function
    a(k, i) = 1 / d
    For j = 1 To n
        If j <> k Then a(j, i) = -a(j, i) / d
    Next j
    For j = 1 To n
        For z = 1 To n
            If j <> k And z <> i Then a(j, z) = a(j, z) + a(k, z) * a(j, i)
        Next z
    Next j
    For j = 1 To n
        If j <> i Then a(k, j) = a(k, j) / d
    Next j
Next i
For i = 1 To n
    For j = 1 To n
       b(i, j) = a(m(i), t(j))
    Next j
Next i
GetEquNiEx = True
End Function

Function GetEquNiEx2(a() As Double, c() As Double, n As Integer) As Boolean
Dim i As Integer, j As Integer, z As Integer, k As Integer, t() As Integer, m() As Integer, b() As Double, x() As Double
ReDim t(n), m(n), b(n), x(n)
Dim w As Double, l As Double
For i = 1 To n - 1
    k = 0
    For j = 1 To n
        If t(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
    t(k) = i: m(i) = k
    For j = 1 To n
        If t(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 t(i) = 0 Then m(n) = i: t(i) = n: Exit For
Next i
For i = 1 To n
    b(i) = 1
    For j = t(i) To n - 1
        For z = 1 To n
            k = m(j)
            If t(z) > j Then
                b(z) = b(z) - a(z, j) * b(k)
            End If
        Next z
    Next j
    k = m(n)
    x(n) = b(k) / a(k, n)
    For j = n - 1 To 1 Step -1
        k = m(j): w = 0
        For z = n To j + 1 Step -1
            w = w + x(z) * a(k, z)
        Next z
        x(j) = (b(k) - w) / a(k, j)
    Next j
    For j = 1 To n
        c(j, i) = x(j): b(j) = 0
    Next j
Next i
GetEquNiEx2 = True
End Function

'解三对角方程
'解存在x()数组中
Sub EquThree(a() As Double, b() As Double, c() As Double, y() As Double, x() As Double, n As Integer)
Dim z() As Double, w() As Double
Dim i As Integer
Dim l As Double
ReDim z(n), w(n)
z(1) = b(1): w(1) = y(1)
For i = 2 To n
    l = a(i) / z(i - 1)
    z(i) = b(i) - l * c(i - 1)
    w(i) = y(i) - l * w(i - 1)
Next i
x(n) = w(n) / z(n)
For i = n - 1 To 1 Step -1
    x(i) = (w(i) - c(i) * x(i + 1)) / z(i)
Next i
End Sub

'Seidel迭代法
Sub Seidel(a() As Double, y() As Double, x() As Double, n As Integer, m As Integer)
Dim z As Integer, i As Integer, j As Integer
Dim s As Double
For z = 1 To m
    For i = 1 To n
        s = 0
        For j = 1 To n
            If j <> i Then s = s + a(i, j) * x(j)
        Next j
        x(i) = (y(i) - s) / a(i, i)
    Next i
Next z
End Sub
'Seidel法的收敛速度的下界,为-1时则未必收敛
Function SeidelSpeed(a() As Double, n As Integer) As Double
Dim i As Integer, j As Integer
Dim t As Double, r As Double, max As Double
For i = 1 To n
    t = 0
    For j = 1 To i - 1
        t = t + Abs(a(i, j))
    Next j
    r = t: t = 0
    For j = i + 1 To n
        t = t + Abs(a(i, j))
    Next j
    If t + r >= Abs(a(i, i)) Then SeidelSpeed = -1: Exit Function
    t = t / (Abs(a(i, i)) - r)
    If t > max Then max = t
Next i
SeidelSpeed = max
End Function

'正定对称的共轭斜量法
Sub GongE(a() As Double, y() As Double, x() As Double, n As Integer, t As Integer)
Dim i As Integer, j As Integer, z As Integer
Dim s As Double, k As Double, d As Double, p() As Double, r() As Double, m() As Double
ReDim p(n), r(n), m(n)
For i = 1 To n
    s = 0
    For j = 1 To n
        s = s + a(i, j) * x(j)
    Next j
    r(i) = y(i) - s: p(i) = r(i)
Next i
For z = 1 To t
k = Inner(r(), r(), n)
For i = 1 To n
    s = 0
    For j = 1 To n
        s = s + a(i, j) * p(j)
    Next j
    m(i) = s
Next i
d = k / Inner(m(), p(), n)
For i = 1 To n
    x(i) = x(i) + d * p(i)
    r(i) = r(i) - d * m(i)
Next i
d = Inner(r(), r(), n) / k
For i = 1 To n
    p(i) = d * p(i) + r(i)
Next i
Next z
End Sub
Function Inner(a() As Double, b() As Double, n As Integer) As Double
Dim i As Integer
Dim t As Double
For i = 1 To n
    t = a(i) * b(i) + t
Next i
Inner = t
End Function
'一般情况下的共轭斜量法
Sub GongEx(a() As Double, y() As Double, x() As Double, n As Integer, num As Integer)
Dim i As Integer, j As Integer, z As Integer
Dim s As Double, k As Double, t As Double, p() As Double, r() As Double, m() As Double
ReDim p(n), r(n), m(n)
For i = 1 To n
    s = 0
    For j = 1 To n
        s = s + a(i, j) * x(j)
    Next j
    r(i) = y(i) - s
Next i
For i = 1 To n
    s = 0
    For j = 1 To n
        s = s + a(j, i) * r(j)
    Next j
    p(i) = s: t = t + s * s
Next i
For z = 1 To num
    k = 0
    For i = 1 To n
        s = 0
        For j = 1 To n
            s = s + a(i, j) * p(j)
        Next j
        m(i) = s: k = k + s * s
    Next i
    k = t / k
    For i = 1 To n
        x(i) = x(i) + k * p(i)
        r(i) = r(i) - k * m(i)
    Next i
    k = 0
    For i = 1 To n
        s = 0
        For j = 1 To n
            s = s + a(j, i) * r(j)
        Next j
        m(i) = s: k = k + s * s
    Next i
    t = k / t
    For i = 1 To n
        p(i) = t * p(i) + m(i)
    Next i
    t = k
Next z
End Sub

'乘幂法
Function MultipleMi(a() As Double, x() As Double, n As Double, m As Double)
Dim i As Integer, j As Integer, z As Integer
Dim s As Double, t As Double, k As Double, y() As Double
ReDim y(n)
k = Abs(x(i))
For i = 2 To n
    t = Abs(x(i))
    If t > k Then k = t
Next i
For i = 1 To n
    x(i) = x(i) / k
    y(i) = x(i)
Next i
For z = 1 To m
    For i = 1 To n
        s = 0
        For j = 1 To n
            s = s + a(i, j) * y(j)
        Next j