插值.bas

实现M5加密算法的源程序

Attribute VB_Name = "SimpleBase"
Option Explicit
'拟合数据数应小于32767

Function Omi2(x As Double) As Double
Dim p As Double
Dim a As Double, b As Double
p = Abs(x)
If p < 1.5 Then
    a = 1.5 - p
    a = a * a
End If
If p < 0.5 Then
    b = 0.5 - p
    b = b * b
    b = b + b + b
End If
Omi2 = (a - b) / 2
End Function

Function Omi3(x As Double) As Double
Dim p As Double, a As Double, b As Double
p = Abs(x)
If p < 2 Then
    a = 2 - p
    a = a * a * a
End If
If p < 1 Then
    b = 1 - p
    b = b * b * b
    b = b + b
    b = b + b
End If
Omi3 = (a - b) / 6
End Function

Function Omi3d(x As Double) As Double
Dim p As Double, a As Double, b As Double
p = Abs(x)
If p < 2 Then
    a = 2 - p
    a = a * a
End If
If p < 1 Then
    b = 1 - p
    b = b * b
    b = b + b
    b = b + b
End If
Omi3d = (a - b) / 2
If x > 0 Then Omi3d = -Omi3d
End Function

Function Omi3t(x As Double) As Double
Dim p As Double, a As Double, b As Double
p = Abs(x)
If p < 2 Then a = 2 - p
If p < 1 Then
    b = 1 - p
    b = b + b
    b = b + b
End If
Omi3t = a - b
End Function

'l正时为微分次数,负时为积分次数
Function Omi(k As Integer, l As Integer, ByVal x As Double) As Double
Dim i As Integer, j As Integer
Dim r As Double, p As Double
Dim c As Boolean
If l >= 0 Then
    If Abs(x) >= (k + 1) / 2 Then Exit Function
    If (l Mod 2) = 1 And x > 0 Then c = True
    x = -Abs(x)
End If
r = x + (k + 1) / 2
For i = 0 To k + 1
    If r > 0 Then
        p = 1
        For j = 1 To i
            p = p * (k + 2 - j) / j
        Next j
        If i Mod 2 Then p = -p
        For j = 1 To k - l
            p = p * r
        Next j
        For j = k - l To -1
            p = p / r
        Next j
        Omi = Omi + p
    Else
        Exit For
    End If
    r = r - 1
Next i
For i = k - l To 2 Step -1
    Omi = Omi / i
Next i
If c Then Omi = -Omi
End Function

'多节点基样条
Function OmiEx3(x As Double) As Double
If Abs(x) >= 3 Then Exit Function
Dim s As Double
s = Omi3(x + 0.5) + Omi3(x - 0.5)
s = s + s
s = s + s
s = s + s
s = 20 * Omi3(x) - s + Omi3(x - 1) + Omi3(x + 1)
OmiEx3 = s / 6
End Function

'保凸拟合(等距)
Function BaoTu(x() As Double, y() As Double, n As Long, xi As Double) As Double
Dim j As Long
Dim r As Double, w0 As Double, wn As Double, r0 As Double, rn As Double, h As Double
Dim a As Double, b As Double, q As Double, q1 As Double, q2 As Double, q3 As Double, xn As Double, s As Double
a = x(0): b = x(n)
h = (b - a) / n
q = (xi - a) / (b - a): q1 = k2(xi, a, b, q): q2 = k2(xi, a + h, b, q): xn = a + h + h
For j = 1 To n - 1
    q3 = k2(xi, xn, b, q)
    r = (y(j + 1) + y(j - 1) - y(j) - y(j)) / (h * h)
    If j = 1 Then r0 = r
    If j = n - 1 Then rn = r
    s = s + r * (q1 + q3 - q2 - q2)
    xn = xn + h: q1 = q2: q2 = q3
Next j
w0 = (k2(xi, a + h, b, q) - k2(xi, a, b, q)) / h - k1(xi, a, b, q)
wn = (k2(xi, b - h, b, q) - k2(xi, b, b, q)) / h
BaoTu = (y(0) * (b - xi) + y(n) * (xi - a)) / (b - a)
BaoTu = s / h + r0 * w0 + rn * wn + BaoTu
End Function
'保凸拟合(不等距)
Function BaoTuEx(x() As Double, y() As Double, n As Long, xi As Double) As Double
Dim j As Long
Dim r As Double, w0 As Double, wn As Double, r0 As Double, rn As Double, h As Double
Dim a As Double, b As Double, q As Double, q1 As Double, q2 As Double, q3 As Double, xn As Double, s As Double
a = x(0): b = x(n)
q = (xi - a) / (b - a): q1 = k2(xi, a, b, q): q2 = k2(xi, x(1), b, q)
For j = 1 To n - 1
    q3 = k2(xi, x(j + 1), b, q)
    r = ((y(j + 1) - y(j)) / (x(j + 1) - x(j)) - (y(j) - y(j - 1)) / (x(j) - x(j - 1))) / (x(j + 1) - x(j - 1))
    If j = 1 Then r0 = r
    If j = n - 1 Then rn = r
    r = r * ((q3 - q2) / (x(j + 1) - x(j)) - (q2 - q1) / (x(j) - x(j - 1))): s = s + r + r
    q1 = q2: q2 = q3
Next j
w0 = (k2(xi, x(1), b, q) - k2(xi, a, b, q)) / (x(1) - a) - k1(xi, a, b, q)
wn = (k2(xi, b - h, b, q) - k2(xi, b, b, q)) / (b - x(n - 1))
BaoTuEx = (y(0) * (b - xi) + y(n) * (xi - a)) / (b - a)
BaoTuEx = s + r0 * w0 + rn * wn + BaoTuEx
End Function
Function k1(xi As Double, t As Double, b As Double, q As Double) As Double
Dim a As Double
k1 = t - b: k1 = k1 * k1 * q
If xi > t Then
    a = xi - t
    a = a * a
    k1 = k1 - a
End If
k1 = k1 / 2
End Function
Function k2(xi As Double, t As Double, b As Double, q As Double) As Double
Dim a As Double
k2 = t - b
k2 = k2 * k2 * k2 * q
If xi > t Then
    a = xi - t
    a = a * a * a
    k2 = k2 + a
End If
k2 = k2 / 6
End Function

'max为函数最大值 , min为函数最小值
Sub FunDraw(x() As Double, y() As Double, frm As Form, max As Double, min As Double, n As Long)
Dim j As Double, s As Double
Dim i As Long
frm.Scale (x(0), max)-(x(n), min)
For j = 0 To n
    frm.PSet (x(j), BaoTu(x(), y(), n, x(j)))
Next j
'frm.DrawWidth = 2
'For i = 0 To UBound(x)
    'frm.PSet (x(i), y(i)), vbRed
'Next i
'frm.DrawWidth = 1
End Sub