number.bas

实现M5加密算法的源程序

        t = t - Int(t)
        t = t * t - t + 1 / 6
        t = 2 * PI * PI * t + 1
        k = k * t
    Next j
    l = l + k
Next i
l = (1 + PI * PI / 3) ^ s + l + l
If n Mod 2 = 0 Then
    d = 0
    For i = 1 To s
        If h(i) Mod 2 = 1 Then d = d + 1
    Next i
    l = l + (1 - PI * PI / 6) ^ d * (1 + PI * PI / 3) ^ (s - d)
End If
l = l / n - 1
GetW2 = l
End Function

'W2(n,h)(h为是实数)
Function GetW2Esp(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
    k = 1
    For j = 1 To s
        t = i * h(j)
        t = t - Int(t)
        t = t * t - t + 1 / 6
        t = 2 * PI * PI * t + 1
        k = k * t
    Next j
    l = l + (n - i) * k / n
Next i
l = (1 + PI * PI / 3) ^ s + l + l
l = l / n - 1
GetW2Esp = l
End Function

'W4(n,h)(h为是整数)
Function GetW4(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, temp As Double
temp = PI * PI: temp = temp * temp
For i = 1 To (n - 1) \ 2
    k = 1
    For j = 1 To s
        t = i * h(j) / n
        t = t - Int(t)
        t = BFour(t)
        t = 1 - (temp + temp) * t / 3
        k = k * t
    Next j
    l = l + k
Next i
l = (1 + temp / 45) ^ s + l + l
If n Mod 2 = 0 Then
    d = 0
    For i = 1 To s
        If h(i) Mod 2 = 1 Then d = d + 1
    Next i
    l = l + (1 - temp * 7 / 360) ^ d * (1 + temp / 45) ^ (s - d)
End If
l = l / n - 1
GetW4 = l
End Function
Function BFour(x As Double) As Double
Dim s As Double
s = x * x
BFour = s * s - s * x * 2 + s - 1 / 30
End Function

Sub SetW(DataName As String, ws As Integer)
Dim db As Database
Dim k As Recordset
Dim tbl As TableDef
Dim fld As Field
Dim h(10) As Double, s As Integer, n As Long, i As Integer
Set db = OpenDatabase(DataName)
h(1) = 1
Do
    s = InputBox("")
    If s = 0 Then db.Close: Exit Sub
    Set tbl = db.TableDefs(CStr(s))
    Set k = tbl.OpenRecordset
    Do Until k.EOF
        Set fld = k.Fields("n")
        n = fld.Value
        For i = 2 To s
            Set fld = k.Fields("h" & i)
            h(i) = fld.Value
        Next i
        k.Edit
        If ws = 2 Then
            Set fld = k.Fields("W2")
            fld.Value = GetW2(h(), n, s)
        Else
            Set fld = k.Fields("W4")
            fld.Value = GetW4(h(), n, s)
        End If
        k.Update
        k.MoveNext
    Loop
Loop While True
db.Close
End Sub

'带权积分
'带权w(x)=exp(x)的二次样条积分
'g0为a处一阶导数值,g1为b处一阶导数值
Function ExpD(a As Double, b As Double, g0 As Double, g1 As Double, n As Long) As Double
ExpD = ExpSimpleD(b, g1, n) - ExpSimpleD(a, g0, n)
End Function
Function ExpSimpleD(b As Double, g0 As Double, n As Long) As Double
Dim h As Double, c As Double, t As Double, x As Double
Dim i As Long
h = b / n
c = f(h) - f(0) - g0 * h: t = c
For i = 1 To n - 1
    x = x + h
    c = f(x + h) - f(x) * 2 + f(x - h) - c
    t = t + c * Exp(-x)
Next i
t = t / (h * h)
t = (1 - Exp(-h)) * t: t = t + t
t = t - Exp(-b) * ((1 + h) * f(b) - f(x) + c) / h + f(0) + g0
ExpSimpleD = t
End Function

'获取傅利叶级数系数a(k)
'即振荡积分
'尚未调试
Function FuShiA(m As Integer, n As Long, g0 As Double, gn As Double) As Double
Dim i As Long
Dim h As Double, x As Double, xe As Double, t As Double, f() As Double, c() As Double
ReDim f(n), c(-1 To n + 1)
h = (PI + PI) / n
For i = 0 To n
    f(i) = fuli(x)
    x = x + h
Next i
SimpleT1 f(), c(), g0, gn, h, n
xe = h * m: x = xe
For i = 1 To n
    t = t + Cos(x) * (f(i) - c(i))
    x = x + xe
Next i
x = Sin(xe / 2)
xe = m * m
x = x * x * 12 / (xe * xe * h * h * h)
t = x * t: t = t + t
t = (f(n) - f(0) + c(0) - c(n)) * x - t
t = (gn - g0) / xe + t
FuShiA = t
End Function

Function FuShiB(m As Integer, n As Long, g0 As Double, gn As Double) As Double
Dim i As Long
Dim h As Double, x As Double, xe As Double, t As Double, f() As Double, c() As Double
ReDim f(n), c(-1 To n + 1)
h = (PI + PI) / n
For i = 0 To n
    f(i) = fuli(x)
    x = x + h
Next i
SimpleT1 f(), c(), g0, gn, h, n
xe = h * m: x = xe
For i = 1 To n
    t = t + Sin(x) * (f(i) - c(i))
    x = x + xe
Next i
x = Sin(xe / 2)
x = x * x * 24 / (xe * xe * xe * m)
t = x * t
t = 6 * (f(n) - f(0) + c(0) - c(n)) * (1 - Sin(xe) / xe) / (xe * xe * m) - t
t = t - (f(n) - f(0)) / m
FuShiB = t
End Function
Function fuli(x As Double) As Double
fuli = x * Cos(x)
End Function

'常微分方程初值问题
'Euler折线法
Function Euler(a As Double, b As Double, y As Double, n As Long) As Double
Dim i As Long
Dim h As Double, x As Double
x = a: h = (b - a) / n
For i = 1 To n
    y = y + h * yu(x, y)
    x = x + h
Next i
Euler = y
End Function

'改进的Euler折线法
Function EulerEx(a As Double, b As Double, y As Double, n As Long, total As Long) As Double
Dim i As Long, j As Long
Dim h As Double, x As Double, t As Double, d As Double, s As Double
x = a: h = (b - a) / n
For i = 1 To n
    t = y: s = yu(x, y)
    y = y + h * s
    x = x + h
    For j = 1 To total
        d = yu(x, y)
        y = t + h * (s + d) / 2
    Next j
Next i
EulerEx = y
End Function

Function Runge(a As Double, b As Double, y As Double, n As Long) As Double
Dim k1 As Double, k2 As Double, k3 As Double, k4 As Double
Dim h As Double, x As Double
Dim i As Long
x = a: h = (b - a) / n
For i = 1 To n
    k1 = h * yu(x, y)
    k2 = h * yu(x + h / 2, y + k1 / 2)
    k3 = h * yu(x + h / 2, y + k2 / 2)
    k4 = h * yu(x + h, y + k3)
    y = y + (k1 + k2 + k2 + k3 + k3 + k4) / 6
    x = x + h
Next i
Runge = y
End Function

'未调试
Function Adams(a As Double, b As Double, y0 As Double, y1 As Double, y2 As Double, y3 As Double, n As Long, total As Long) As Double
Dim i As Long, j As Long
Dim at1 As Double, at2 As Double, at3 As Double, at4 As Double
Dim h As Double, x As Double, t As Double, d As Double, s As Double, y As Double
h = (b - a) / n
at1 = yu(a, y0): at2 = yu(a + h, y1): at3 = yu(a + h + h, y2): at4 = yu(a + h * 3, y3)
x = a + h * 4: y = y3
For i = 1 To n - 3
    t = y
    y = y + h * (55 * at4 - 59 * at3 + 37 * at2 - 9 * at1) / 24
    d = 19 * at4 - 5 * at3 + at2
    For j = 1 To total
        s = yu(x, y)
        y = t + h * (s * 9 + d) / 24
    Next j
    at1 = at2: at2 = at3: at3 = at4: at4 = yu(x, y)
    x = x + h
Next i
Adams = y
End Function

'未调试
Function Milne(a As Double, b As Double, y0 As Double, y1 As Double, y2 As Double, y3 As Double, n As Long, total As Long) As Double
Dim i As Long, j As Long
Dim at1 As Double, at2 As Double, at3 As Double
Dim h As Double, x As Double, t As Double, d As Double, s As Double, y As Double
h = (b - a) / n
at1 = yu(a + h, y1): at2 = yu(a + h + h, y2): at3 = yu(a + h * 3, y3)
x = a + h * 4
For i = 1 To n - 3
    y = y0 + h * (at1 + at1 - at2 + at3 + at3) * 4 / 3
    d = at2 + 4 * at3: t = y2
    For j = 1 To total
        s = yu(x, y)
        y = t + h * (s + d) / 3
    Next j
    at1 = at2: at2 = at3: at3 = yu(x, y)
    y0 = y1: y1 = y2: y2 = y3: y3 = y
    x = x + h
Next i
Milne = y
End Function

'未调试
Function HaMing(a As Double, b As Double, y0 As Double, y1 As Double, y2 As Double, y3 As Double, plus As Double, n As Long) As Double
Dim i As Long
Dim at1 As Double, at2 As Double, at3 As Double
Dim h As Double, x As Double, p As Double, c As Double, d As Double, m As Double, y As Double
h = (b - a) / n
at1 = yu(a + h, y1): at2 = yu(a + h + h, y2): at3 = yu(a + h * 3, y3)
x = a + h * 4: d = plus
For i = 1 To n - 3
    p = y1 + h * (at1 + at1 - at2 + at3 + at3) * 4 / 3
    m = p - d * 112 / 121
    c = (9 * y4 - y2 + 8 * h * (yu(x, m) + at3 + at3 - at2)) / 8
    d = p - c
    y = c + d * 9 / 121
    at1 = at2: at2 = at3: at3 = yu(x, y)
    y1 = y2: y2 = y3: y3 = y4: y4 = y
    x = x + h
Next i
HaMing = y
End Function

'未完成
Function Obrechkoff(a As Double, b As Double, y0 As Double, y1 As Double, y2 As Double, y3 As Double, plus As Double, n As Long) As Double
Dim i As Long
Dim at1 As Double, at2 As Double, at3 As Double
Dim h As Double, x As Double, p As Double, c As Double, d As Double, m As Double, y As Double
h = (b - a) / n
Obrechkoff = y
End Function

Function yu(x As Double, y As Double) As Double
yu = y - 2 * x / y
End Function

'快速傅立叶变换(FFT)
'n是2的方幂
'可能有问题
'C语言改写而成
Sub FFT(ar() As Double, ai() As Double, n As Long)
Dim wr() As Double, wi() As Double, xr() As Double, xi() As Double
ReDim wr(n / 2), wi(n / 2), xr(n), xi(n)
Dim tr As Double, ti As Double, tr1 As Double, ti1 As Double, tr2 As Double, ti2 As Double, mr1 As Double, mi1 As Double, mr2 As Double, mi2 As Double
Dim u As Long, m As Long, i As Long, j As Long, k As Long, l As Long, s As Long
Dim m1 As Long, m2 As Long, m3 As Long, d As Long
u = n - 1
Do While u > 0
    m = m + 1
    u = u \ 2
Loop
For i = 0 To n - 1
    xr(i) = ar(i)
    xi(i) = ai(i)
Next i
For i = 0 To n - 1
    l = 0: k = i
    For j = 1 To m
        l = l + l + k Mod 2
        k = k \ 2
    Next j
    If l > i Then
        tr = xr(i): ti = xi(i)
        xr(i) = xr(l): xi(i) = xi(l)
        xr(l) = tr: xi(l) = ti
    End If
Next i
For i = 0 To n \ 2 - 1
    wr(i) = Cos(PI2 / n * i)
    wi(i) = Sin(PI2 / n * i)
Next i
s = 1: d = n
For i = 1 To m
    d = d / 2
    For j = 0 To n - 1 Step s + s
        m3 = 0
        For k = 0 To s - 1
            m1 = j + k: m2 = m1 + s
            mr1 = xr(m2): mi1 = xi(m2)
            mr2 = wr(m3): mi2 = wi(m3)
            tr = mr1 * mr2 - mi1 * mi2
            ti = mr1 * mi2 + mr2 * mi1
            tr1 = xr(m1) + tr: ti1 = xi(m1) + ti
            tr2 = xr(m1) - tr: ti2 = xi(m1) - ti
            xr(m1) = tr1: xi(m1) = ti1
            xr(m2) = tr2: xi(m2) = ti2
            m3 = m3 + d
        Next k
    Next j
    s = s + s
Next i
For i = 0 To n - 1
    ar(i) = xr(i)
    ai(i) = xi(i)
Next i
Erase wr, wi, xi, xr
End Sub

'我编的
Sub FFT1(ar() As Double, ai() As Double, n As Long)
Dim wr() As Double, wi() As Double, xr() As Double, xi() As Double
ReDim wr(n / 2), wi(n / 2), xr(n), xi(n)
Dim tr As Double, ti As Double, tr1 As Double, ti1 As Double, tr2 As Double, ti2 As Double, mr1 As Double, mi1 As Double, mr2 As Double, mi2 As Double
Dim u As Long, m As Long, i As Long, j As Long, k As Long, s As Long, t As Long, z As Long, w As Long
u = n - 1
Do While u > 0
    m = m + 1
    u = u \ 2
Loop
For i = 0 To n - 1
    xr(i) = ar(i)
    xi(i) = ai(i)
Next i
For i = 0 To n \ 2 - 1
    wr(i) = Cos(PI2 / n * i)
    wi(i) = Sin(PI2 / n * i)
Next i
k = n
For i = 1 To m
    s = k: k = k / 2
    For j = 0 To n - 1 Step s
        w = j
        For z = 0 To k - 1
            t = con1(w, m, i)
            mr1 = xr(w + k): mi1 = xi(w + k)
            mr2 = wr(t): mi2 = wi(t)
            tr = mr1 * mr2 - mi1 * mi2
            ti = mr1 * mi2 + mr2 * mi1
            tr1 = xr(w) + tr: ti1 = xi(w) + ti
            tr2 = xr(w) - tr: ti2 = xi(w) - ti
            xr(w) = tr1: xi(w) = ti1
            xr(w + k) = tr2: xi(w + k) = ti2
            w = w + 1
        Next z
    Next j
Next i
For i = 0 To n - 1
    w = con2(i, m)
    ar(w) = xr(i)
    ai(w) = xi(i)
Next i
Erase wr, wi, xi, xr
End Sub
Function con1(ByVal x As Long, n As Long, l As Long) As Long
Dim b As Long, i As Long
For i = 1 To n - l + 1
    x = x \ 2
Next i
For i = 1 To l - 1
    b = b + b + x Mod 2
    x = x \ 2
Next i
For i = 1 To n - l
    b = b + b
Next i
con1 = b
End Function
Function con2(ByVal x As Long, m As Long) As Long
Dim b As Long, i As Long
For i = 1 To m
    b = b + b + x Mod 2
    x = x \ 2
Next i
con2 = b
End Function

'FFT的逆变换(只需将wi()数组变号即可,并将a()数组除以n)
Sub FFTNi(ar() As Double, ai() As Double, n As Long)
Dim wr() As Double, wi() As Double, xr() As Double, xi() As Double
ReDim wr(n / 2), wi(n / 2), xr(n), xi(n)
Dim tr As Double, ti As Double, tr1 As Double, ti1 As Double, tr2 As Double, ti2 As Double, mr1 As Double, mi1 As Double, mr2 As Double, mi2 As Double
Dim u As Long, m As Long, i As Long, j As Long, k As Long, l As Long, s As Long
Dim m1 As Long, m2 As Long, m3 As Long, d As Long
u = n - 1
Do While u > 0
    m = m + 1
    u = u \ 2
Loop
For i = 0 To n - 1
    xr(i) = ar(i)
    xi(i) = ai(i)
Next i
For i = 0 To n - 1
    l = 0: k = i
    For j = 1 To m
        l = l + l + k Mod 2
        k = k \ 2
    Next j
    If l > i Then
        tr = xr(i): ti = xi(i)
        xr(i) = xr(l): xi(i) = xi(l)
        xr(l) = tr: xi(l) = ti
    End If
Next i
For i = 0 To n \ 2 - 1
    wr(i) = Cos(PI2 / n * i)
    wi(i) = -Sin(PI2 / n * i)
Next i
s = 1: d = n
For i = 1 To m
    d = d / 2
    For j = 0 To n - 1 Step s + s
        m3 = 0
        For k = 0 To s - 1
            m1 = j + k: m2 = m1 + s
            mr1 = xr(m2): mi1 = xi(m2)
            mr2 = wr(m3 \ (s + s)): mi2 = wi(m3 \ (s + s))
            tr = mr1 * mr2 - mi1 * mi2
            ti = mr1 * mi2 + mr2 * mi1
            tr1 = xr(m1) + tr: ti1 = xi(m1) + ti
            tr2 = xr(m1) - tr: ti2 = xi(m1) - ti
            xr(m1) = tr1: xi(m1) = ti1
            xr(m2) = tr2: xi(m2) = ti2
            m3 = m3 + d
        Next k
    Next j
    s = s + s
Next i
For i = 0 To n - 1
    ar(i) = xr(i) / n
    ai(i) = xi(i) / n
Next i
Erase wr, wi, xi, xr
End Sub