number.bas

实现M5加密算法的源程序

Attribute VB_Name = "Calculate"
Option Explicit

'函数性质缺省函数
Function f(x As Double) As Double
f = x * x + 3 * x + 1
End Function
Function g(x As Double) As Double
g = 2 * x + 3
End Function
Function h(x As Double) As Double
h = 2
End Function

Function TwoEq(a As Double, b As Double, e As Double) As Double
Dim s0 As Long, s1 As Long, s As Long, y As Double
y = f(b)
If Abs(y) < e Then TwoEq = b: Exit Function
s1 = Sgn(y)
y = f(a)
If Abs(y) < e Then TwoEq = a: Exit Function
s0 = Sgn(y)
If s1 = s0 Then Exit Function
Do While True
    TwoEq = (a + b) / 2
    y = f(TwoEq)
    If Abs(y) < e Then Exit Function
    If Sgn(y) = s0 Then a = TwoEq Else b = TwoEq
    If Abs(a - b) < e Then TwoEq = (a + b) / 2: Exit Function
Loop
End Function

Function Newton(x As Double, e As Double) As Double
Dim q As Double
Dim k As Long
Do
    q = f(x) / g(x)
    x = x - q
    If Abs(q) < e Then Newton = x: Exit Function
    k = k + 1
Loop Until k > 1000
MsgBox "Newton法失效"
End Function

Function Stevensen(x As Double, e As Double) As Double
Dim q As Double, t As Double, d As Double
Dim k As Long
Do
    t = f(x)
    d = f(x + t) - t
    If Abs(d) < e Then Stevensen = x: Exit Function
    q = t * t / d
    x = x - q
    If Abs(q) < e Then Stevensen = x: Exit Function
    k = k + 1
Loop Until k > 1000
MsgBox "斯梯芬森方法失效"
End Function

Function GenQieBi(x As Double, e As Double) As Double
Dim q As Double, t As Double, d As Double
Dim k As Long
Do
    t = f(x): d = g(x)
    q = t / d: d = d * d * d: d = d + d: t = t * t
    q = q + h(x) * t / d
    x = x - q
    If Abs(q) < e Then GenQieBi = x: Exit Function
    k = k + 1
Loop Until k > 200
MsgBox "切比雪夫迭代法失效"
End Function

Function Alpha(x As Double, e As Double) As Double
Dim a As Double, b As Double, x1 As Double, x2 As Double
Dim f1 As Double, f2 As Double, al As Double
al = (Sqr(5) - 1) / 2
ConvertAB x, 0.01, a, b
x1 = a + (1 - al) * (b - a): x2 = a + b - x1: f1 = f(x1): f2 = f(x2)
Do While True
    If Abs(b - a) < e Then Alpha = (a + b) / 2: Exit Function
    If f1 < f2 Then
        b = x2: x2 = x1: f2 = f1
        x1 = a + (1 - al) * (b - a): f1 = f(x1)
    ElseIf f1 > f2 Then
        a = x1: x1 = x2: f1 = f2
        x2 = a + al * (b - a): f2 = f(x2)
    Else
        a = x1: b = x2: x1 = a + (1 - al) * (b - a): x2 = a + b - x1
        f1 = f(x1): f2 = f(x2)
    End If
Loop
End Function

Function FenShu(x As Double, e As Double) As Double
Dim p() As Double
Dim total As Double
Dim n As Integer, k As Integer
Dim a As Double, b As Double, x1 As Double, x2 As Double
Dim f1 As Double, f2 As Double, al As Double
ReDim p(100)
ConvertAB x, 0.01, a, b
n = 1: total = (b - a) / e: p(0) = 1: p(1) = 1
Do
    n = n + 1
    p(n) = p(n - 1) + p(n - 2)
Loop Until p(n) >= total
x1 = a + (b - a) * p(n - 2) / p(n): x2 = a + (b - a) * p(n - 1) / p(n)
f1 = f(x1): f2 = f(x2)
Do
    If f1 < f2 Then
        b = x2: x2 = x1: f2 = f1: x1 = a + (b - a) * p(n - k - 3) / p(n - k - 1): f1 = f(x1)
    Else
        a = x1: x1 = x2: f1 = f2: x2 = a + (b - a) * p(n - k - 2) / p(n - k - 1): f2 = f(x2)
    End If
    k = k + 1
Loop Until k = n - 3
If f1 < f2 Then
    b = x2: x2 = x1: f2 = f1
Else
    a = x1
End If
x1 = (b - a) / 4 + a: f1 = f(x1)
If f1 < f2 Then
    FenShu = (a + x2) / 2
Else
    FenShu = (b + x2) / 2
End If
End Function

Sub ConvertAB(x As Double, e As Double, a As Double, b As Double)
Dim t As Double, q As Double
Dim fx As Double, ft As Double, fq As Double
t = x + e: fx = f(x): ft = f(t)
If ft <= fx Then
    Do While True
        e = e + e
        q = t + e: fq = f(q)
        If ft <= fq Then
            a = x: b = q: Exit Sub
        Else
            x = t: t = q: fx = ft: ft = fq
        End If
    Loop
Else
    Do While True
        e = e + e
        q = x - e: fq = f(q)
        If fq >= fx Then
            a = q: b = t: Exit Sub
        Else
            t = x: x = q: ft = fx: fx = fq
        End If
    Loop
End If
End Sub

'抛物线法
'要求:x1<x0<x2,f(x1)>f(x0)<f(x2)
'未调试
Function PaoWu(x0 As Double, x1 As Double, x2 As Double, e As Double) As Double
Dim f0 As Double, f1 As Double, f2 As Double, ft As Double
Dim xt As Double, t As Double
f0 = f(x1): f1 = f(x1): f2 = f(x2)
Do
    ft = f(xt)
    If Abs(x2 - x1) < e Then PaoWu = xt: Exit Function
    If f0 > ft Then
        If x0 > xt Then
            x2 = x0: f2 = f0: x0 = xt: f0 = ft
        ElseIf x0 < xt Then
            x1 = x0: f1 = f0: x0 = xt: f0 = ft
        End If
    ElseIf f0 < ft Then
        If x0 > xt Then
            x1 = xt: f1 = ft
        ElseIf x0 < xt Then
            x2 = xt: f2 = ft
        End If
    Else
        If x0 < xt Then
            x1 = x0: x0 = (x0 + xt) / 2: x2 = xt: f1 = f0: f2 = ft: f0 = f(x0)
        ElseIf x0 > xt Then
            x1 = xt: x2 = x0: x0 = (x0 + xt) / 2: f1 = ft: f2 = f0: f0 = f(x0)
        Else
            t = f((x1 + x0) / 2)
            If t < f0 Then
                x0 = (x1 + x0) / 2: f0 = t: x2 = x0: f2 = ft
            ElseIf t > f0 Then
                x1 = (x2 + x0) / 2: f1 = t
            End If
        End If
    End If
Loop While True
End Function

Function Three(a As Double, b As Double, e As Double) As Double
Dim u As Double, v As Double, q As Double, m As Double, n As Double, s As Double, z As Double, x As Double, w As Double
u = g(b): v = g(a): m = f(b): n = f(a)
Do
    s = 3 * (m - n) / (b - a): z = s - u - v
    w = Sqr(z * z - u * v)
    x = a + (b - a) * (1 - (u + w + z) / (u - v + w + w))
    q = g(x)
    If (Abs(a - b) < e Or q = 0) Then
        Three = x
        Exit Function
    End If
    If q > 0 Then
        b = x
    Else
        a = x
    End If
Loop While True
End Function

'数值微分
'Larrange公式(一阶)
Function LarrangeS(x As Double, h As Double, num As Long) As Double
Dim i As Long, j As Long
Dim tm As Double, tn As Double, m As Double, n As Double
Dim k As Double, c As Double
Dim bl As Boolean
c = num / (num + 1): j = num + 2: k = c
bl = True
For i = 1 To num
    tm = f(x + h * i) * c: tn = f(x - h * i) * c
    If bl Then
        m = m + tm: n = n + tn
    Else
        m = m + tn: n = n + tm
    End If
    bl = Not bl
    k = k * (num - i) / j: c = k / (j - num): j = j + 1
Next i
LarrangeS = (m - n) / h
End Function

'Larrange公式(二阶)
Function LarrangeD(x As Double, h As Double, num As Long) As Double
Dim i As Long, j As Long
Dim tm As Double, tn As Double, m As Double, n As Double
Dim k As Double, c As Double, s As Double, sm As Double, sn As Double
Dim bl As Boolean
c = num / (num + 1): c = c + c: j = num + 2: k = c
bl = True
For i = 1 To num
    tm = f(x + h * i) * c: tn = f(x - h * i) * c
    If bl Then
        m = m + tm + tn
        sm = sm + c
    Else
        n = n + tm + tn
        sn = sn + c
    End If
    bl = Not bl
    k = k * (num - i) / j: c = k / ((j - num) * (j - num)): j = j + 1
Next i
s = sm - sn
s = s + s
LarrangeD = (m - n - f(x) * s) / (h * h)
End Function

'外推法
Sub WeiOut(x As Double, h As Double)
Dim a(16, 16) As Double
Dim i As Integer, m As Long, n As Integer, j As Integer
n = 16
a(1, 1) = cha(x, h)
For i = 2 To n
    h = h / 2: a(1, i) = cha(x, h)
    'If Len(CStr(a(1, i))) < 10 Then Exit Sub
    m = 1
    For j = 2 To i
        m = m * 4
        a(j, i) = (m * a(j - 1, i) - a(j - 1, i - 1)) / (m - 1)
    Next j
    Maths.tval = Maths.tval & a(i, i) & vbCrLf
Next i
End Sub
Function cha(x As Double, h As Double) As Double
cha = (f(x + h / 2) - f(x - h / 2)) / h
End Function

'Simpson公式(g为一阶导数)
'需给出f0-fn,y()范围1至N-1
Sub Simpson(f() As Double, y() As Double, g0 As Double, gn As Double, h As Double, n As Long)
Dim i As Long, l As Double, z() As Double, w() As Double
ReDim z(n), w(n)
z(1) = 4: w(1) = 3 * (f(2) - f(0)) / h - g0
For i = 2 To n - 1
    l = 1 / z(i - 1)
    z(i) = 4 - l
    w(i) = 3 * (f(i + 1) - f(i - 1)) / h - l * w(i - 1)
Next i
y(n - 1) = (w(n - 1) - gn) / z(n - 1)
For i = n - 2 To 1 Step -1
    y(i) = (w(i) - y(i + 1)) / z(i)
Next i
End Sub

Function Cotes(a As Double, b As Double, n As Long) As Double
Dim s As Double, t As Double, h As Double, h1 As Double
Dim x As Double, x1 As Double
Dim i As Long, j As Long, k(6) As Double
h = (b - a) / n: h1 = h / 6: x = a: x1 = x
k(0) = 41: k(1) = 216: k(2) = 27: k(3) = 272: k(4) = 27: k(5) = 216: k(6) = 41
For i = 0 To 6
    s = 0
    For j = 1 To n
        s = s + fun(x1)
        x1 = x1 + h
    Next j
    t = t + s * k(i)
    x = x + h1: x1 = x
Next i
Cotes = h * t / 840
End Function

Function QieBiJi(a As Double, b As Double, n As Long) As Double
Dim h As Double, s As Double, x As Double, w(4) As Double
Dim i As Long, j As Long
h = (b - a) / n: x = a + h / 2
Open DataPath1 & ".dat" For Binary As #1
For i = 1 To 4
    Get #1, , w(i)
    w(i) = w(i) * h / 2
Next i
Close #1
For i = 1 To n
    For j = 1 To 4
        s = s + fun(x + w(j)) + fun(x - w(j))
    Next j
    s = s + fun(x)
    x = x + h
Next i
QieBiJi = s * h / 9
End Function

Function Romberg(a As Double, b As Double, n As Integer) As Double
Dim k As Long, num As Long
Dim q(16, 16) As Double
Dim s As Double, h As Double, t As Double, v As Double, x As Double
Dim o As Long, total As Long, i As Long
h = (b - a) / 2: t = h * (fun(a) + fun(b)): k = 2
If n = 1 Then v = (gun(a) - gun(b)) * h * h / 3 Else v = 0
q(1, 1) = t + v
If n = 0 Then num = 17 Else num = 16
total = 1
Do
    'For i = a + h To b - h Step h + h
        's = s + fun(i)
    'Next i
    x = a + h
    For i = 1 To total
        s = s + fun(x)
        x = x + h + h
    Next i
    t = t / 2 + h * s: s = 0: h = h / 2: v = v / 4: q(1, k) = t + v
    If n = 0 Then o = 1 Else o = 4
    For i = 2 To k
        o = o + o: o = o + o
        q(i, k) = (q(i - 1, k) * o - q(i - 1, k - 1)) / (o - 1)
    Next i
    Maths.tval = Maths.tval & q(k, k) & vbCrLf
    k = k + 1: total = total + total
Loop Until k = num
End Function

Function Gauss(a As Double, b As Double, n As Integer) As Double
Dim i As Integer
Dim h As Double, d As Double, s As Double, w() As Double, x() As Double
ReDim x((n + 1) / 2), w((n + 1) / 2)
h = b - a
d = h / 2
Open DataPath & n & ".dat" For Binary As #1
For i = 1 To (n + 1) \ 2
    Get #1, , x(i)
    Get #1, , w(i)
Next i
Close #1
For i = 1 To (n + 1) \ 2
    s = s + w(i) * (fun(d * x(i) + d + a) + fun(-d * x(i) + d + a))
Next i
s = s * d
Gauss = s
End Function
Function fun(x As Double) As Double
fun = 1 / Log(x)
'fun = x * x
End Function
Function gun(x As Double) As Double
gun = -1 / (Log(x) * Log(x) * x)
'gun = x * Log(x) * Log(x)
End Function
'需要进一步调试
'舍入误差过大
Sub GetGen(l As Long)
Dim a() As Double, b() As Double, x() As Double
Dim m As Double, n As Double, k As Double
Dim i As Long
ReDim a(l), b(l), x(l + 1)
x(1) = 1 / 3: x(2) = 1
For i = 2 To l
    GetGenPre x(), i
Next i
a(l) = 1: m = l + l: n = l + l + l + l - 1
For i = 1 To l
    a(l - i) = a(l - i + 1) * m * (m - 1) / n
    a(l - i) = a(l - i) / (i + i)
    m = m - 2: n = n - 2
Next i
FToQ a(), b(), l
k = 1
For i = l + l + 1 To l + l + l + l - 1 Step 2
    k = k * i
Next i
For i = 4 To l + l Step 2
    k = k / i
Next i
For i = 1 To l
    a(i) = k * FSpe(b(), x(i), l)
    a(i) = a(i) * a(i)
    a(i) = a(i) * x(i) * (1 - x(i))
    a(i) = 2 / a(i)
    x(i) = Sqr(x(i))
    Maths.tval = Maths.tval & x(i) & " " & a(i) & Chr(13) & Chr(10)
Next i
If MsgBox("need save", vbOKCancel) = vbOK Then
Open "E:\范翔\素材\data" & l + l & ".dat" For Binary As #1
For i = 1 To l
    Put #1, , x(i)
    Put #1, , a(i)
Next i
Close #1
End If
End Sub
Sub GetGenPre(x() As Double, l As Long)
Dim a() As Double
Dim m As Double, n As Double, c As Double, d As Double
Dim i As Long
ReDim a(l)
a(l) = 1: m = l + l: n = m + m - 1
For i = 1 To l
    a(l - i) = a(l - i + 1) * m * (m - 1) / n
    a(l - i) = a(l - i) / (i + i)
    m = m - 2: n = n - 2
Next i
c = 0
For i = 1 To l
    d = TwoEqSpe(c, x(i), 2E-16, a(), l)
    c = x(i)
    x(i) = d
Next i
x(l + 1) = 1
End Sub
Function TwoEqSpe(ByVal a As Double, ByVal b As Double, e As Double, c() As Double, l As Long) As Double
Dim s0, s1, s As Double
s0 = Sgn(FSpe(c(), a, l))
s1 = Sgn(FSpe(c(), b, l))
Do While True
    TwoEqSpe = (a + b) / 2
    s = Sgn(FSpe(c(), TwoEqSpe, l))
    If s = 0 Then Exit Function
    If s = s0 Then a = TwoEqSpe Else b = TwoEqSpe
    If Abs(a - b) < e Then TwoEqSpe = (a + b) / 2: Exit Function
Loop
End Function
Function FSpe(a() As Double, x As Double, l As Long) As Double
Dim i As Long, bl As Boolean
Dim n As Double, m As Double, c As Double
bl = n Mod 2: c = 1
If bl Then n = a(0) Else m = a(0)
For i = 1 To l
    c = c * x
    bl = Not bl
    If bl Then n = n + a(i) * c Else m = m + a(i) * c
Next i
FSpe = m - n
End Function
Sub FToQ(a() As Double, b() As Double, l As Long)
Dim i As Long
For i = 1 To l
    b(i - 1) = a(i) * i
Next i
End Sub

'高维网格积分
'W2(n,h)(h为是整数)
Function GetW2(h() As Double, n As Long, s As Integer) As Double
Dim i As Long, j As Integer, d As Integer
Dim t As Double, k As Double, l As Double
For i = 1 To (n - 1) \ 2
    k = 1
    For j = 1 To s
        t = i * h(j) / n