回复 344楼 ysr2857
Private Sub Command1_Click()
Dim m, n
m = Trim(Text1): n = Trim(Text2)
ts = Timer
c = MbC4(Trim(m), Trim(n))
Text3 = c & "用时" & Timer - ts & "秒，有" & Len(c) & "位"
End Sub
Private Sub Command2_Click()
Text1 = ""
Text2 = ""
Text3 = ""
End Sub
Public Function MbC4(D1 As String, D2 As String) As String '快速乘法
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Currency, w1 As Currency, wlr As Currency, wl1 As Currency, tr As Currency, t1 As Currency
Dim pi As Currency, t As Currency, tr1 As Currency
Dim xr() As Currency, a As String
a = Trim(D1)
B = Trim(D2)
X = Len(a) \ 6: Y = Len(B) \ 6
a = String(Val(X * 6 + 6 - Len(a)), "0") & a
B = String(Val(Y * 6 + 6 - Len(B)), "0") & B
X = X + 1: Y = Y + 1
sb1 = X + Y
sb2 = Log(sb1) / Log(2)
If InStr(sb2, ".") = 0 Then
sb2 = sb2
Else
sb2 = Int(sb2) + 1
End If
sb = 2 ^ sb2
a = String(Val(sb) * 6 - Len(a), "0") & a
B = String(Val(sb) * 6 - Len(B), "0") & B
ReDim x_(1 To sb): ReDim y_(1 To sb)
For i1 = 1 To sb
x_(i1) = Mid(a, (sb - i1 + 1) * 6 - 5, 6): y_(i1) = Mid(B, (sb - i1 + 1) * 6 - 5, 6)
If Len(x_(i1)) < 6 Then
x_(i1) = String(6 - Len(x_(i1)), "0") & x_(i1)
ElseIf Len(y_(i1)) < 6 Then
y_(i1) = String(6 - Len(y_(i1)), "0") & y_(i1)
Else
x_(i1) = x_(i1): y_(i1) = y_(i1)
End If
Next
ReDim xr(0 To (Len(a) - 6) \ 6): ReDim yr(0 To (Len(B) - 6) \ 6): ReDim zr(0 To (Len(B) - 6) \ 6)
If Len(a) = 6 Then
xr(0) = a: yr(0) = B
Else
Dim I As Long, J As Long, mn As Long, lh As Long, k As Long
'位序倒置
n = sb '求数组大小，其值必须是2的幂
lh = n / 2
J = n / 2
For I = 1 To n - 2
Debug.Print I, J
k = lh '下面是向右进位算法
Do
If k > J Then Exit Do '高位是1吗
J = J - k '是的，高位置0
k = k / 2 '准备次高位的权
Loop Until k = 0 '次高位的权若非0，则检查新的次高位
J = J + k '非则若最高位是0，则置1
xr(I + 1) = x_(J + 1): yr(I + 1) = y_(J + 1)
Next
xr(0) = x_(1): xr(1) = x_(1 + sb / 2)
yr(0) = y_(1): yr(1) = y_(1 + sb / 2)
End If
Dim xi() As Currency: Dim yi() As Currency: Dim zi() As Currency
n = sb '求数组大小，其值必须是2的幂
m = 0
l = 2
pi = "3.141592653589793238"
Do
l = l + l
m = m + 1
Loop Until l > n
n = l / 2
ReDim xi(n - 1): ReDim yi(n - 1): ReDim zi(n - 1)
l = 1
Do
le = 2 ^ l
le1 = le / 2
wr = 1
wi = 0
If l = 1 Then
t = 0
Else
t = pi / le1
End If
w1r = Cos(t)
w1i = -Sin(t)
r = 0
Do
p = r
Do
q = p + le1
tr = xr(q) * wr - xi(q) * wi
ti = xr(q) * wi + xi(q) * wr
tr1 = yr(q) * wr - yi(q) * wi
ti1 = yr(q) * wi + yi(q) * wr
xr(q) = xr(p) - tr
xi(q) = xi(p) - ti
xr(p) = xr(p) + tr
xi(p) = xi(p) + ti
yr(q) = yr(p) - tr1
yi(q) = yi(p) - ti1
yr(p) = yr(p) + tr1
yi(p) = yi(p) + ti1
xr(p) = Format(Val(xr(p)), "0.000000"): xi(p) = Format(Val(xi(p)), "0.000000")
yr(p) = Format(Val(yr(p)), "0.000000"): yi(p) = Format(Val(yi(p)), "0.000000")
p = p + le
Loop Until p > n - 1
wr2 = wr * w1r - wi * w1i
wi2 = wr * w1i + wi * w1r
wr = wr2
wi = wi2
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m
For I = 0 To n - 1 '仅输出模
zr(I) = xr(I) * yr(I) - xi(I) * yi(I): zi(I) = xr(I) * yi(I) + xi(I) * yr(I)
zr(I) = Format(Val(zr(I)), "0.000000"): zi(I) = Format(Val(zi(I)), "0.000000")
's = s & "/" & zr(I)
's1 = s1 & "/" & zi(I)
Next
J = sb
ReDim x_(1 To sb): ReDim y_(1 To sb)
For k = 1 To J
n1 = n1 + 1
ReDim Preserve x_(1 To n1)
x_(n1) = zr(n1 - 1): y_(n1) = zi(n1 - 1)
x_(n1) = Format(Val(x_(n1)), "0.000000"): y_(n1) = Format(Val(y_(n1)), "0.000000")
Next
'位序倒置
n = sb '求数组大小，其值必须是2的幂
lh = n / 2
J = n / 2
For I = 1 To n - 2
Debug.Print I, J
k = lh '下面是向右进位算法
Do
If k > J Then Exit Do '高位是1吗
J = J - k '是的，高位置0
k = k / 2 '准备次高位的权
Loop Until k = 0 '次高位的权若非0，则检查新的次高位
J = J + k '非则若最高位是0，则置1
xr(I + 1) = x_(J + 1): yr(I + 1) = y_(J + 1)
'js = js & "/" & x_(J + 1)
'js1 = js1 & "/" & y_(J + 1)
Next
'sx1 = "/" & x_(1) & "/" & x_(1 + sb / 2) & js
'sy1 = "/" & y_(1) & "/" & y_(1 + sb / 2) & js1
xr(0) = x_(1): xr(1) = x_(1 + sb / 2)
yr(0) = y_(1): yr(1) = y_(1 + sb / 2)
ns = Len(a) \ 6: Jn = ns
ReDim zr(0 To ns - 1)
m = 0
l = 2
pi = "3.141592653589793238"
Do
l = l + l
m = m + 1
Loop Until l > ns
ns = l / 2
ReDim xi(ns - 1): ReDim yi(ns - 1): ReDim zi(ns - 1)
l = 1
Do
le = 2 ^ l
le1 = le / 2
wr = 1
wi = 0
If l = 1 Then
t = 0
Else
t = -1 * pi / le1
End If
w1r = Cos(t)
w1i = -Sin(t)
r = 0
Do
p = r
Do
q = p + le1
tr = xr(q) * wr - xi(q) * wi
ti = xr(q) * wi + xi(q) * wr
tr1 = yr(q) * wr - yi(q) * wi
ti1 = yr(q) * wi + yi(q) * wr
xr(q) = xr(p) - tr
xi(q) = xi(p) - ti
xr(p) = xr(p) + tr
xi(p) = xi(p) + ti
yr(q) = yr(p) - tr1
yi(q) = yi(p) - ti1
yr(p) = yr(p) + tr1
yi(p) = yi(p) + ti1
xr(p) = Format(Val(xr(p)), "0.000000"): xi(p) = Format(Val(xi(p)), "0.000000")
yr(p) = Format(Val(yr(p)), "0.000000"): yi(p) = Format(Val(yi(p)), "0.000000")
p = p + le
Loop Until p > ns - 1
wr2 = wr * w1r - wi * w1i
wi2 = wr * w1i + wi * w1r
wr = wr2
wi = wi2
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m
For I = 0 To ns - 1 '仅输出模
zr(I) = (xr(I) - yi(I)) / n
zr(I) = Format(Val(zr(I) + 0.5), "0.000000")
If InStr(zr(I), ".") = 0 Then
s121 = zr(I)
Else
s121 = Left(zr(I), InStr(zr(I), ".") - 1)
End If
's0 = "/" & s121 & s0
zr(I) = s121
Next
For i1 = 1 To Val(Jn - sb1 + 1)
zr(sb1 + i1 - 2) = 0
Next
For i1 = 0 To n - 1
If zr(i1) < 0 Then
zr(i1) = "000000"
ElseIf Len(zr(i1)) < 6 Then
zr(i1) = String(6 - Len(zr(i1)), "0") & zr(i1)
Else
zr(i1) = zr(i1)
End If
's5 = s5 & "/" & zr(i1)
If i1 = 0 Then
s6 = Val(Left(zr(i1), Len(zr(i1)) - 6))
If Len(s6) < 6 Then
s6 = String(6 - Len(s6), "0") & s6
Else
s6 = s6
End If
s8 = Right(zr(i1), 6)
ElseIf Val(zr(i1)) >= 0 Then
s7 = Val(zr(i1)) + Val(s6)
If Len(s7) = 6 Or Len(s7) = 12 Or Len(s7) = 18 Then
s7 = s7
Else
s7 = String(24 - Len(s7), "0") & s7
End If
s10 = Right(s7, 6)
s11 = s10 & s11
If Len(s7) < 6 Then
s7 = String(6 - Len(s7), "0") & s7
ElseIf Len(s7) = 6 Then
s6 = "000000"
Else
s7 = s7
s6 = Val(Left(s7, Len(s7) - 6))
End If
Else
s6 = s6
End If
Next
s9 = s6 & s11 & s8
s9 = qdqd0(Trim(s9))
's2 = nifft(dxcx1(Trim(s)), dxcx1(Trim(s1)), Trim(sb1))
's3 = nifft(Trim(sx1), Trim(sy1), Trim(sb1))
MbC4 = s9
End Function
Private Function qdqd0(sa As String) As String
a = sa
Do While Left(a, 1) = "0"
a = Mid(a, 2)
Loop
If a = "" Then
a = 0
Else
a = a
End If
qdqd0 = a
End Function
2022-02-16 21:20
