数学中国

 找回密码
 注册
搜索
热搜: 活动 交友 discuz
楼主: ysr

大数乘法与快速傅里叶变换

[复制链接]
 楼主| 发表于 2021-3-19 20:36 | 显示全部楼层
速度基本无法再提高了,改一下算法,改为8位一组试试,如何呢?如下是优化后的代码:
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 Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
  Dim pi As Double, t As Double, tr1 As Double

Private Sub Command1_Click()
  Dim xr() As Double, a As String
  a = Trim(Text1)
  b = Trim(Text3)
  ts = Timer
  sb1 = Len(a) + Len(b)
  sb2 = Log(sb1) / Log(2)
  If InStr(sb2, ".") = 0 Then
  sb2 = sb2
  Else
  sb2 = Int(sb2) + 1
  End If
  sb = 2 ^ sb2
  Print sb
  If Len(a) = Len(b) And 2 ^ (Int(Log(Len(a)) / Log(2))) = Len(a) Then
    a = String(Val(sb) - Len(a), "0") & a
  b = String(Val(sb) - Len(b), "0") & b
  a = dxcx0(Trim(a), Val(sb)): b = dxcx0(Trim(b), Val(sb))
  Else

  a = String(Val(sb) - Len(a), "0") & a
  b = String(Val(sb) - Len(b), "0") & b
  a = dxcx0(Trim(a), Val(sb)): b = dxcx0(Trim(b), Val(sb))
  End If
  ReDim xr(0 To Len(a) - 1): ReDim yr(0 To Len(b) - 1): ReDim zr(0 To Len(b) - 1)
  For i1 = 0 To Len(a) - 1
  xr(i1) = Mid(a, i1 + 1, 1)
  yr(i1) = Mid(b, i1 + 1, 1)

     Next
  
  Dim xi(): Dim yi(): Dim zi()
  n = Len(a) '求数组大小,其值必须是2的幂
m = 0
  l = 2
  pi = 3.14159265358979
  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
     
      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)
      

      s = s & "/" & zr(i)
      s1 = s1 & "/" & zi(i)
      Next
     s2 = nifft(dxcx1(Trim(s)), dxcx1(Trim(s1)), Trim(sb1))
     
      Text2 = s2 & "有" & Len(s2) & "位,用时" & Timer - ts & "秒"
  End Sub

  Private Sub Command2_Click()
  Text1 = ""
  Text2 = ""
  Text3 = ""
  Form1.Cls
  End Sub

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


  Private Function nifft(sa As String, sb As String, sb1 As String) As String
  
  Dim xi(): Dim yi(): Dim zi()
  Dim xr(), yr(), zr()
  s2 = Split(sa, "/")
  s3 = Split(sb, "/")
      j = UBound(s2)
      n = j
     For k = 1 To j
         n1 = n1 + 1
          ReDim Preserve xr(0 To n1 - 1)
          ReDim Preserve yr(0 To n1 - 1)
         xr(n1 - 1) = s2(n1): yr(n1 - 1) = s3(n1)
       Next
     

  ReDim zr(0 To j - 1)

  m = 0
  l = 2
  pi = 3.14159265358979
  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 = -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
     
      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) - yi(i)) / n
      
     
      s1 = Int(Val(zr(i) + 0.5))
      s = "/" & s1 & s
      zr(i) = s1
      Next
      For i1 = 1 To Val(j - sb1 + 1)
      zr(sb1 + i1 - 2) = 0
      Next
      
     
     
      For i1 = 0 To n - 1
      If zr(i1) < 0 Then
      zr(i1) = 0
      Else
      zr(i1) = zr(i1)
      End If
      
      
      If i1 = 0 Then
      s6 = Int(zr(i1)) \ 10
      s8 = Int(zr(i1)) Mod 10
      ElseIf Val(zr(i1)) >= 0 Then
      s7 = Int(zr(i1)) + Val(s6)
      s10 = Val(s7) Mod 10
      s11 = s10 & s11
      s6 = Val(s7) \ 10
      Else
      s6 = Val(s6)
      End If
     
      Next
      s9 = s6 & s11 & s8
     
  nifft = qdqd0(Trim(s9))

  End Function

  Private Function dxcx0(sa As String, sb As String) As String

  Dim x_() As Double, a As String
    a = Trim(sa)
    ReDim x_(1 To sb)
    For i1 = 1 To sb
    x_(i1) = Mid(a, sb - i1 + 1, 1)
      Next
    Dim n As Integer, i As Long, j As Long, mn As Long, lh As Long, t As Double, 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
    s = s & x_(j + 1)
    Next
    dxcx0 = x_(1) & x_(1 + sb / 2) & s
   

  End Function

  Private Function dxcx1(sa As String) As String

  Dim x_() As Double, a As String
    a = Trim(sa)
     

  s2 = Split(sa, "/")
  s3 = Split(sb, "/")
      j = UBound(s2)
      sb = j
     
       ReDim x_(1 To sb)
     For k = 1 To j
         n1 = n1 + 1
          ReDim Preserve x_(1 To n1)
        
         x_(n1) = s2(n1)
       Next
    Dim n As Integer, i As Long, mn As Long, lh As Long, t As Double
    '位序倒置
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
    s = s & "/" & x_(j + 1)
    Next
    dxcx1 = "/" & x_(1) & "/" & x_(1 + sb / 2) & s
   
    End Function
回复 支持 反对

使用道具 举报

 楼主| 发表于 2021-3-20 10:31 | 显示全部楼层
修改了一下程序,这回可能对了,程序速度和模仿手工计算的相似,还能提高,希望老师能优化一下程序,前面的一位一组是对的,速度太慢能否优化?这个待优化的快速程序代码如下:

Private Sub Command1_Click()
  Dim xr() As String, a As String
  a = Trim(Text1)
  b = Trim(Text3)
  ts = Timer
  x = Len(a) \ 4: Y = Len(b) \ 4
  If Val(4 * x) = Len(a) Then
  a = a
  ElseIf Val(4 * Y) = Len(b) Then
  b = b
  Else
  a = String(Val(x * 4 + 4 - Len(a)), "0") & a
  b = String(Val(Y * 4 + 4 - Len(b)), "0") & b
  x = x + 1: Y = Y + 1
  End If
  
  
  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
  Print sb
  

  a = String(Val(sb) * 4 - Len(a), "0") & a
  b = String(Val(sb) * 4 - Len(b), "0") & b
  a = dxcx0(Trim(a), Val(sb)): b = dxcx0(Trim(b), Val(sb))
  Print a
  ReDim xr(0 To (Len(a) - 4) \ 4): ReDim yr(0 To (Len(b) - 4) \ 4): ReDim zr(0 To (Len(b) - 4) \ 4)
  If Len(a) = 4 Then
  xr(0) = a: yr(0) = b
  Else
  For i1 = 0 To (Len(a) - 4) \ 4
  xr(i1) = Mid(a, (i1 + 1) * 4 - 3, 4)
  yr(i1) = Mid(b, (i1 + 1) * 4 - 3, 4)

     Next
     End If
  
  Dim xi(): Dim yi(): Dim zi()
  n = Len(a) \ 4 '求数组大小,其值必须是2的幂
m = 0
  l = 2
  pi = 3.14159265358979
  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")
     Print xr(p), xr(q)
      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
     s2 = nifft(dxcx1(Trim(s)), dxcx1(Trim(s1)), Trim(sb1))
     
      Text2 = s2 & "有" & Len(s2) & "位,用时" & Timer - ts & "秒"
  End Sub

  Private Sub Command2_Click()
  Text1 = ""
  Text2 = ""
  Text3 = ""
  Form1.Cls
  End Sub

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


  Private Function nifft(sa As String, sb As String, sb1 As String) As String
  
  Dim xi(): Dim yi(): Dim zi()
  Dim xr(), yr()
  Dim zr() As String
  
  s2 = Split(sa, "/")
  s3 = Split(sb, "/")
      J = UBound(s2)
      n = J
     For k = 1 To J
         n1 = n1 + 1
          ReDim Preserve xr(0 To n1 - 1)
          ReDim Preserve yr(0 To n1 - 1)
         xr(n1 - 1) = s2(n1): yr(n1 - 1) = s3(n1)
         xr(n1 - 1) = Format(Val(xr(n1 - 1)), "0.000000"): yr(n1 - 1) = Format(Val(yr(n1 - 1)), "0.000000")

       Next
     

  ReDim zr(0 To J - 1)

  m = 0
  l = 2
  pi = 3.14159265358979
  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 = -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 > 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) - yi(I)) / n
      zr(I) = Format(Val(zr(I) + 0.5), "0.000000")
     If InStr(zr(I), ".") = 0 Then
     s1 = zr(I)
     Else
     s1 = Left(zr(I), InStr(zr(I), ".") - 1)
      End If
      s = "/" & s1 & s
      zr(I) = s1
      Next
      For i1 = 1 To Val(J - sb1 + 1)
      zr(sb1 + i1 - 2) = 0
      Next
      
     
     
      For i1 = 0 To n - 1
      If zr(i1) < 0 Then
      zr(i1) = "00000000"
      Else
      zr(i1) = zr(i1)
      End If
      s5 = s5 & "/" & zr(i1)
      
      If i1 = 0 Then
      If Len(zr(i1)) < 4 Then
      zr(i1) = String(4 - Len(zr(i1)), "0") & zr(i1)
      Else
      zr(i1) = zr(i1)
      End If
      s6 = Val(Left(zr(i1), Len(zr(i1)) - 4))
      If Len(s6) < 4 Then
      s6 = String(4 - Len(s6), "0") & s6
      Else
      s6 = s6
      End If
      s8 = Right(zr(i1), 4)
      ElseIf Val(zr(i1)) >= 0 Then
      s7 = MPC1(Trim(zr(i1)), Trim(s6))
      s10 = Right(s7, 4)
      s11 = s10 & s11
      If Len(s7) < 4 Then
      s7 = String(4 - Len(s7), "0") & s7
      ElseIf Len(s7) = 4 Then
      s6 = "00000000"
      Else
      s7 = s7
      s6 = Val(Left(s7, Len(s7) - 4))
      End If
      Else
      s6 = s6
      End If
     
      Next
      s9 = s6 & s11 & s8
     
  nifft = qdqd0(Trim(s9))

  End Function

  Private Function dxcx0(sa As String, sb As String) As String

  Dim x_() As String, a As String
    a = Trim(sa)
    ReDim x_(1 To sb)
    For i1 = 1 To sb
    x_(i1) = Mid(a, (sb - i1 + 1) * 4 - 3, 4)
    If Len(x_(i1)) < 4 Then
    x_(i1) = String(4 - Len(x_(i1)), "0") & x_(i1)
    Else
    x_(i1) = x_(i1)
    End If
   
      Next
    Dim n As Integer, I As Long, J As Long, mn As Long, lh As Long, t As Double, 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
    s = s & x_(J + 1)
    Next
    dxcx0 = x_(1) & x_(1 + sb / 2) & s
   
   

  End Function

  Private Function dxcx1(sa As String) As String

  Dim x_() As Double, a As String
    a = Trim(sa)
     

  s2 = Split(sa, "/")
  s3 = Split(sb, "/")
      J = UBound(s2)
      sb = J
     
       ReDim x_(1 To sb)
     For k = 1 To J
         n1 = n1 + 1
          ReDim Preserve x_(1 To n1)
        
         x_(n1) = Format(Val(s2(n1)), "0.000000")
       Next
    Dim n As Integer, I As Long, mn As Long, lh As Long, t As Double
    '位序倒置
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
x_(J + 1) = Format(Val(x_(J + 1)), "0.000000")
    s = s & "/" & x_(J + 1)
    Next
    x_(1) = Format(Val(x_(1)), "0.000000"): x_(1 + sb / 2) = Format(Val(x_(1 + sb / 2)), "0.000000")
    dxcx1 = "/" & x_(1) & "/" & x_(1 + sb / 2) & s
   
    End Function

Public Function MPC1(D1 As String, D2 As String) As String 'jiafa
Dim x, Y '两数长度

If Len(D1) >= Len(D2) Then
D4 = String(Len(D1) - Len(D2), "0") & D2
d3 = D1
Else
D4 = D2
d3 = String(Len(D2) - Len(D1), "0") & D1
End If
x = Len(d3): Y = Len(D4)
Dim a() As Integer, B1() As Integer, C1() As Integer, E1() As Integer
ReDim a(1 To x)
ReDim B1(1 To Y)
ReDim C1(1 To x)
ReDim E1(1 To x)
Dim I, J, C2, CJ, JW
For J = Y To 1 Step -1 'D2
JW = 0 '进位清0
B1(J) = Mid$(D4, J, 1) '每位数
For I = x To 1 Step -1  'D1
   a(I) = Mid$(d3, I, 1) '每位数
   C1(I) = a(I) + B1(I) + JW '计算jia
   JW = C1(I) \ 10
   E1(I) = C1(I) Mod 10
  Next
  Next
  For r = 1 To x
  If JW = 0 Then
  MPC1 = MPC1 & E1(r)
  Else
  jc = jc & E1(r)
  MPC1 = JW & jc
  End If
  Next
  
End Function

999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999(270位)*9=8999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999991有271位,用时0.8242188秒(哈哈哈这回可能对了)
回复 支持 反对

使用道具 举报

 楼主| 发表于 2021-3-21 00:00 | 显示全部楼层
4位一组的程序重发如下,是修改好的,还有提升空间,希望老师指点,谢谢!代码如下:

Private Sub Command1_Click()
  Dim xr() As String, a As String
  a = Trim(Text1)
  b = Trim(Text3)
  ts = Timer
  x = Len(a) \ 4: Y = Len(b) \ 4
  If Val(4 * x) = Len(a) Then
  a = a
  ElseIf Val(4 * Y) = Len(b) Then
  b = b
  Else
  a = String(Val(x * 4 + 4 - Len(a)), "0") & a
  b = String(Val(Y * 4 + 4 - Len(b)), "0") & b
  x = x + 1: Y = Y + 1
  End If
  
  
  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
  Print sb
  

  a = String(Val(sb) * 4 - Len(a), "0") & a
  b = String(Val(sb) * 4 - Len(b), "0") & b
  a = dxcx0(Trim(a), Val(sb)): b = dxcx0(Trim(b), Val(sb))
  Print a
  ReDim xr(0 To (Len(a) - 4) \ 4): ReDim yr(0 To (Len(b) - 4) \ 4): ReDim zr(0 To (Len(b) - 4) \ 4)
  If Len(a) = 4 Then
  xr(0) = a: yr(0) = b
  Else
  For i1 = 0 To (Len(a) - 4) \ 4
  xr(i1) = Mid(a, (i1 + 1) * 4 - 3, 4)
  yr(i1) = Mid(b, (i1 + 1) * 4 - 3, 4)

     Next
     End If
  
  Dim xi(): Dim yi(): Dim zi()
  n = Len(a) \ 4 '求数组大小,其值必须是2的幂
m = 0
  l = 2
  pi = 3.14159265358979
  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")
     Print xr(p), xr(q)
      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
     s2 = nifft(dxcx1(Trim(s)), dxcx1(Trim(s1)), Trim(sb1))
     
      Text2 = s2 & "有" & Len(s2) & "位,用时" & Timer - ts & "秒"
  End Sub

  Private Sub Command2_Click()
  Text1 = ""
  Text2 = ""
  Text3 = ""
  Form1.Cls
  End Sub

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


  Private Function nifft(sa As String, sb As String, sb1 As String) As String
  
  Dim xi(): Dim yi(): Dim zi()
  Dim xr(), yr()
  Dim zr() As String
  
  s2 = Split(sa, "/")
  s3 = Split(sb, "/")
      J = UBound(s2)
      n = J
     For k = 1 To J
         n1 = n1 + 1
          ReDim Preserve xr(0 To n1 - 1)
          ReDim Preserve yr(0 To n1 - 1)
         xr(n1 - 1) = s2(n1): yr(n1 - 1) = s3(n1)
         xr(n1 - 1) = Format(Val(xr(n1 - 1)), "0.000000"): yr(n1 - 1) = Format(Val(yr(n1 - 1)), "0.000000")

       Next
     

  ReDim zr(0 To J - 1)

  m = 0
  l = 2
  pi = 3.14159265358979
  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 = -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 > 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) - yi(I)) / n
      zr(I) = Format(Val(zr(I) + 0.5), "0.000000")
     If InStr(zr(I), ".") = 0 Then
     s1 = zr(I)
     Else
     s1 = Left(zr(I), InStr(zr(I), ".") - 1)
      End If
      s = "/" & s1 & s
      zr(I) = s1
      Next
      For i1 = 1 To Val(J - sb1 + 1)
      zr(sb1 + i1 - 2) = 0
      Next
      
     
     
      For i1 = 0 To n - 1
      If zr(i1) < 0 Then
      zr(i1) = "0000"
      ElseIf Len(zr(i1)) < 4 Then
      zr(i1) = String(4 - 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)) - 4))
      If Len(s6) < 4 Then
      s6 = String(4 - Len(s6), "0") & s6
      Else
      s6 = s6
      End If
      s8 = Right(zr(i1), 4)
      ElseIf Val(zr(i1)) >= 0 Then
      s7 = MPC1(Trim(zr(i1)), Trim(s6))
      s10 = Right(s7, 4)
      s11 = s10 & s11
      If Len(s7) < 4 Then
      s7 = String(4 - Len(s7), "0") & s7
      ElseIf Len(s7) = 4 Then
      s6 = "0000"
      Else
      s7 = s7
      s6 = Val(Left(s7, Len(s7) - 4))
      End If
      Else
      s6 = s6
      End If
     
      Next
      s9 = s6 & s11 & s8
     
  nifft = qdqd0(Trim(s9))

  End Function

  Private Function dxcx0(sa As String, sb As String) As String

  Dim x_() As String, a As String
    a = Trim(sa)
    ReDim x_(1 To sb)
    For i1 = 1 To sb
    x_(i1) = Mid(a, (sb - i1 + 1) * 4 - 3, 4)
    If Len(x_(i1)) < 4 Then
    x_(i1) = String(4 - Len(x_(i1)), "0") & x_(i1)
    Else
    x_(i1) = x_(i1)
    End If
   
      Next
    Dim n As Integer, I As Long, J As Long, mn As Long, lh As Long, t As Double, 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
    s = s & x_(J + 1)
    Next
    dxcx0 = x_(1) & x_(1 + sb / 2) & s
   
   

  End Function

  Private Function dxcx1(sa As String) As String

  Dim x_() As Double, a As String
    a = Trim(sa)
     

  s2 = Split(sa, "/")
  s3 = Split(sb, "/")
      J = UBound(s2)
      sb = J
     
       ReDim x_(1 To sb)
     For k = 1 To J
         n1 = n1 + 1
          ReDim Preserve x_(1 To n1)
        
         x_(n1) = Format(Val(s2(n1)), "0.000000")
       Next
    Dim n As Integer, I As Long, mn As Long, lh As Long, t As Double
    '位序倒置
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
x_(J + 1) = Format(Val(x_(J + 1)), "0.000000")
    s = s & "/" & x_(J + 1)
    Next
    x_(1) = Format(Val(x_(1)), "0.000000"): x_(1 + sb / 2) = Format(Val(x_(1 + sb / 2)), "0.000000")
    dxcx1 = "/" & x_(1) & "/" & x_(1 + sb / 2) & s
   
    End Function

Public Function MPC1(D1 As String, D2 As String) As String 'jiafa
Dim x, Y '两数长度

If Len(D1) >= Len(D2) Then
D4 = String(Len(D1) - Len(D2), "0") & D2
d3 = D1
Else
D4 = D2
d3 = String(Len(D2) - Len(D1), "0") & D1
End If
x = Len(d3): Y = Len(D4)
Dim a() As Integer, B1() As Integer, C1() As Integer, E1() As Integer
ReDim a(1 To x)
ReDim B1(1 To Y)
ReDim C1(1 To x)
ReDim E1(1 To x)
Dim I, J, C2, CJ, JW
For J = Y To 1 Step -1 'D2
JW = 0 '进位清0
B1(J) = Mid$(D4, J, 1) '每位数
For I = x To 1 Step -1  'D1
   a(I) = Mid$(d3, I, 1) '每位数
   C1(I) = a(I) + B1(I) + JW '计算jia
   JW = C1(I) \ 10
   E1(I) = C1(I) Mod 10
  Next
  Next
  For r = 1 To x
  If JW = 0 Then
  MPC1 = MPC1 & E1(r)
  Else
  jc = jc & E1(r)
  MPC1 = JW & jc
  End If
  Next
  
End Function
回复 支持 反对

使用道具 举报

 楼主| 发表于 2021-3-21 00:01 | 显示全部楼层
51750801837147361345408953922231823615475578427966187002956389087112242842559611794590895524485015222232340190036677951157401229518412775512617768569186939216558814320044037671525512763073762127357272238470370174050144130962253437215369727381909754581075278888137687950749577751967083233227932391898438489520107788418593104216764745143648165875043861634459264809523254076330115365651752264033578829280651927862062277153553168278640509846511729608378923480331705134467105387853440058108864733429085916392927954109944314725590758950098556991513109760871599045355054258610189920009222629391227784118536578933067532162200768112408111142115021110731900107445187734955403922470144187293970602024056652942406064101184615113779243566842311336047149767396006622923831778833173731302325162741660053390002572456069410383734906953419854919650114722834689335022576764661864630532595250136417800041415827379233503655173380741283758557001102498493237113089198106283487003030225058520047964863929279622015825870972146222540614667928697270131499340362072956956518213266361286372027834000838185920579641835952404678850816456467012890793914108776148483471455515695596694573035435841962871983614397882366460940068418989292218104650923374259244879331938925143014431488635252724366437874702141802712286424754174951797598162531093933772708024046029028652425912112434874797404878317364682580259997737809632860143160557722236459697738869857870982308547756781808568834445404196611833859447257611721058384661943492980180884154939581533269700568975940675513953240040220315659784214936947443784598154082110443748171557676558299300252057475022409541083454882202785077569281268346332901357829843950754331120530978456454396022570647840977328115888711408099894489014427119253061839770890965114853246067272044090282103668762164047374624856119676063310166566519499552342907590495696658359542582469176737197423550084705962160477441546453931948920780829433665262356431973422034925001341330921812806939990688094288773216182971976433480031438271163051341645350474631300965486818841553843112688792803367397276204832936049162335586140151875524537720425345701193713093955432559579831822084761702059241698184717469616732582600658006352069963017377661059372648121969448644394344799369193952197519795756182552461838014896498816263546256621347748100920475049252091842372001283217811864786225692344297842649982716736752004086642881662614773154743879410328420184057019873899911520017723788801*
51750801837147361345408953922231823615475578427966187002956389087112242842559611794590895524485015222232340190036677951157401229518412775512617768569186939216558814320044037671525512763073762127357272238470370174050144130962253437215369727381909754581075278888137687950749577751967083233227932391898438489520107788418593104216764745143648165875043861634459264809523254076330115365651752264033578829280651927862062277153553168278640509846511729608378923480331705134467105387853440058108864733429085916392927954109944314725590758950098556991513109760871599045355054258610189920009222629391227784118536578933067532162200768112408111142115021110731900107445187734955403922470144187293970602024056652942406064101184615113779243566842311336047149767396006622923831778833173731302325162741660053390002572456069410383734906953419854919650114722834689335022576764661864630532595250136417800041415827379233503655173380741283758557001102498493237113089198106283487003030225058520047964863929279622015825870972146222540614667928697270131499340362072956956518213266361286372027834000838185920579641835952404678850816456467012890793914108776148483471455515695596694573035435841962871983614397882366460940068418989292218104650923374259244879331938925143014431488635252724366437874702141802712286424754174951797598162531093933772708024046029028652425912112434874797404878317364682580259997737809632860143160557722236459697738869857870982308547756781808568834445404196611833859447257611721058384661943492980180884154939581533269700568975940675513953240040220315659784214936947443784598154082110443748171557676558299300252057475022409541083454882202785077569281268346332901357829843950754331120530978456454396022570647840977328115888711408099894489014427119253061839770890965114853246067272044090282103668762164047374624856119676063310166566519499552342907590495696658359542582469176737197423550084705962160477441546453931948920780829433665262356431973422034925001341330921812806939990688094288773216182971976433480031438271163051341645350474631300965486818841553843112688792803367397276204832936049162335586140151875524537720425345701193713093955432559579831822084761702059241698184717469616732582600658006352069963017377661059372648121969448644394344799369193952197519795756182552461838014896498816263546256621347748100920475049252091842372001283217811864786225692344297842649982716736752004086642881662614773154743879410328420184057019873899911520017723788801=
2678145490787694710138406683675886762424440548647327525594133873038347266950777818028529804877891353347244255046176220880989519574806598935727647635506939520211209794328712328126947740602255036264392628020291100698432900022333877728824270241063030019536831205264726059774673235680796770108334187808344470141093691232926021244583649747591446066921535333854157441330339217364515328375889937313870826987520723914810358707312468857074455945610038118763873992681688857200178452843740327143990247797893489008808544535536081089994800692174759199578596435124153437689121057209523619989593900355037400039597165241729693472001589594642879166693443994939359639600558001937547209387175085119128575833396037478944824322315197771750730768383267287661855880237549141571520374996029107650341071074081501632143441365461865431362023790463469688071768208320226561675630348708531708110373889241713511809029065048200924851629507265095151843687852114948736130517024578482059771284515915597919026074110811882700307150650858080695566672584075215926450360584738705106755945308199980660478183563586945902705730231409948938478910176292179410952310592864745025373902249844607078729230098512959498522851990040833671698873183076996740726479481352429756295966452788332789259316653805673821716832204368483794501474641036177256005810416507649770210480747056764697918800240201155395628561684853684747063434500664449114760179367086926249813006056735202795763742271472136500817313178050955623517148960973128433455798444839070215800699075256371341531892841721456473301085595261371522078416420399702683504171019191363372883472415693921509069484448792009627325245658854979159002549542769741504596589556474585723101403486046952706847309188488390899096649210853675228464611999630360164703837131620319658277117300208326944026347808645954842782309633205645471959842390075986572806483755223497552549948718545978420793425036666212191832327909389117021539426917904222971515424022606157676379417252939026404363841301595348953340657745965591468370831371449860774296757556463049237248472613472436666799615694672489600751763956185893123579218797539971253828850163437053164723845196089466298000486042246621634263064752906475198879616484926366053597071510886820878425907736453843197228282179943365206039514850477377542684361730162273246709577442085159313947668552462343596819299335060863266863059421286080919061657175757802930378698062290163092494909661982639303946271889619023417245433500782685534406291391527022657608608130680134939380126048353657648559735671562230714451258764762350451448379438849319710473631625533527602085050858633346936315115309333355986556546113134787213610042970405680162084678021281774169131725296977361816657904233850711837826724976507876061043779177845550646612327663283839673950405472039432309339110648301347672597370034346896334064135099613033528785924976781865447506814375627647103020761819104412269362847848884053177378762432543285984193091289854949394023586451139393716370633759317698908834326333339634557238526607966375935960004339535161203546369910853171524516945554862019134684985868547068709986010719120259421723270478586992153989283736650948010894951434857207445141886625290777550599700738909777829412358941157674866550887608199469723859806486384262574182811908513849905040209776132986875560077482527819093584150609347560543627976011009349905362597935257895287412087233810106693141683912683308201407373637490835862452911394718399632093441745957478079754723679119343225931665785380948807114301164623717144211285944154289093317610519159281033583963015810556464180819928345256832572962092994993950655913067768151113345566231088759030148164710287289335782462899525700500822021472132532771242198844498579102733028256523810973405350508253183399668290339844052874956694050028526577775764032598023622192221226501334606994447870277619118208162133562541678250491647671042612088902952091163094215645902682086024337098910110924115381788866213335441814530141911754446270296630426675020592640816713810936850961517475097707179713115856513919631998747192606018612535436009162735908155778165736923393260822217415131632296810207091918499370220533980666393617256227781299273913099896641458499949922475383775878771011361417556799406151773618626138749998464514570482965286184387350245762317271666329142667987131139265167425606863654107927682419958167046956962340419295394058625063862272314570914511355913547470055072878835349142049813368019960713974852351833144227544851625677753419141348585495260230041111323576432417502998832485027071744397105338790364666312325483218209150368759120749671419554589928540715813531144241705501211397506662448276056294585639532800634727428900861538486753011006379264618338782036690340859441153926551954112364946638768223012718907190897165613382555036872673785444132692292217184974818914501886689736449912131900762592363447945036085692999163815947845753538999001633767986447907729462453017601有4888位,用时3.351563秒(这是4位一组的算法,速度还可以,还有提升空间)
回复 支持 反对

使用道具 举报

 楼主| 发表于 2021-3-21 10:40 | 显示全部楼层
改成了3位一组的程序,省略了调用大数加法程序,但速度仍然不快:

Private Sub Command1_Click()
  Dim xr() As String, a As String
  a = Trim(Text1)
  b = Trim(Text3)
  ts = Timer
  x = Len(a) \ 3: Y = Len(b) \ 3
  If Val(3 * x) = Len(a) Then
  a = a
  ElseIf Val(3 * Y) = Len(b) Then
  b = b
  Else
  a = String(Val(x * 3 + 3 - Len(a)), "0") & a
  b = String(Val(Y * 3 + 3 - Len(b)), "0") & b
  x = x + 1: Y = Y + 1
  End If
  
  
  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
  Print sb
  

  a = String(Val(sb) * 3 - Len(a), "0") & a
  b = String(Val(sb) * 3 - Len(b), "0") & b
  a = dxcx0(Trim(a), Val(sb)): b = dxcx0(Trim(b), Val(sb))
  Print a
  ReDim xr(0 To (Len(a) - 3) \ 3): ReDim yr(0 To (Len(b) - 3) \ 3): ReDim zr(0 To (Len(b) - 3) \ 3)
  If Len(a) = 3 Then
  xr(0) = a: yr(0) = b
  Else
  For i1 = 0 To (Len(a) - 3) \ 3
  xr(i1) = Mid(a, (i1 + 1) * 3 - 2, 3)
  yr(i1) = Mid(b, (i1 + 1) * 3 - 2, 3)

     Next
     End If
  
  Dim xi(): Dim yi(): Dim zi()
  n = Len(a) \ 3 '求数组大小,其值必须是2的幂
m = 0
  l = 2
  pi = 3.14159265358979
  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")
     Print xr(p), xr(q)
      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
     s2 = nifft(dxcx1(Trim(s)), dxcx1(Trim(s1)), Trim(sb1))
     
      Text2 = s2 & "有" & Len(s2) & "位,用时" & Timer - ts & "秒"
  End Sub

  Private Sub Command2_Click()
  Text1 = ""
  Text2 = ""
  Text3 = ""
  Form1.Cls
  End Sub

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


  Private Function nifft(sa As String, sb As String, sb1 As String) As String
  
  Dim xi(): Dim yi(): Dim zi()
  Dim xr(), yr()
  Dim zr() As String
  
  s2 = Split(sa, "/")
  s3 = Split(sb, "/")
      J = UBound(s2)
      n = J
     For k = 1 To J
         n1 = n1 + 1
          ReDim Preserve xr(0 To n1 - 1)
          ReDim Preserve yr(0 To n1 - 1)
         xr(n1 - 1) = s2(n1): yr(n1 - 1) = s3(n1)
         xr(n1 - 1) = Format(Val(xr(n1 - 1)), "0.000000"): yr(n1 - 1) = Format(Val(yr(n1 - 1)), "0.000000")

       Next
     

  ReDim zr(0 To J - 1)

  m = 0
  l = 2
  pi = 3.14159265358979
  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 = -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 > 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) - yi(I)) / n
      zr(I) = Format(Val(zr(I) + 0.5), "0.000000")
     If InStr(zr(I), ".") = 0 Then
     s1 = zr(I)
     Else
     s1 = Left(zr(I), InStr(zr(I), ".") - 1)
      End If
      s = "/" & s1 & s
      zr(I) = s1
      Next
      For i1 = 1 To Val(J - sb1 + 1)
      zr(sb1 + i1 - 2) = 0
      Next
      
     
     
      For i1 = 0 To n - 1
      If zr(i1) < 0 Then
      zr(i1) = "000"
      ElseIf Len(zr(i1)) < 3 Then
      zr(i1) = String(3 - 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)) - 3))
      If Len(s6) < 3 Then
      s6 = String(3 - Len(s6), "0") & s6
      Else
      s6 = s6
      End If
      s8 = Right(zr(i1), 3)
      ElseIf Val(zr(i1)) >= 0 Then
      s7 = Val(zr(i1)) + Val(s6)
         If Len(s7) = 3 Or Len(s7) = 6 Or Len(s7) = 9 Then
          s7 = s7
          Else
          s7 = String(9 - Len(s7), "0") & s7
         End If
      s10 = Right(s7, 3)
      s11 = s10 & s11
      If Len(s7) < 3 Then
      s7 = String(3 - Len(s7), "0") & s7
      ElseIf Len(s7) = 3 Then
      s6 = "000"
      Else
      s7 = s7
      s6 = Val(Left(s7, Len(s7) - 3))
      End If
      Else
      s6 = s6
      End If
     
      Next
      s9 = s6 & s11 & s8
     
  nifft = qdqd0(Trim(s9))

  End Function

  Private Function dxcx0(sa As String, sb As String) As String

  Dim x_() As String, a As String
    a = Trim(sa)
    ReDim x_(1 To sb)
    For i1 = 1 To sb
    x_(i1) = Mid(a, (sb - i1 + 1) * 3 - 2, 3)
    If Len(x_(i1)) < 3 Then
    x_(i1) = String(3 - Len(x_(i1)), "0") & x_(i1)
    Else
    x_(i1) = x_(i1)
    End If
   
      Next
    Dim n As Integer, I As Long, J As Long, mn As Long, lh As Long, t As Double, 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
    s = s & x_(J + 1)
    Next
    dxcx0 = x_(1) & x_(1 + sb / 2) & s
   
   

  End Function

  Private Function dxcx1(sa As String) As String

  Dim x_() As Double, a As String
    a = Trim(sa)
     

  s2 = Split(sa, "/")
  s3 = Split(sb, "/")
      J = UBound(s2)
      sb = J
     
       ReDim x_(1 To sb)
     For k = 1 To J
         n1 = n1 + 1
          ReDim Preserve x_(1 To n1)
        
         x_(n1) = Format(Val(s2(n1)), "0.000000")
       Next
    Dim n As Integer, I As Long, mn As Long, lh As Long, t As Double
    '位序倒置
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
x_(J + 1) = Format(Val(x_(J + 1)), "0.000000")
    s = s & "/" & x_(J + 1)
    Next
    x_(1) = Format(Val(x_(1)), "0.000000"): x_(1 + sb / 2) = Format(Val(x_(1 + sb / 2)), "0.000000")
    dxcx1 = "/" & x_(1) & "/" & x_(1 + sb / 2) & s
   
    End Function

Public Function MPC1(D1 As String, D2 As String) As String 'jiafa
Dim x, Y '两数长度

If Len(D1) >= Len(D2) Then
D4 = String(Len(D1) - Len(D2), "0") & D2
d3 = D1
Else
D4 = D2
d3 = String(Len(D2) - Len(D1), "0") & D1
End If
x = Len(d3): Y = Len(D4)
Dim a() As Integer, B1() As Integer, C1() As Integer, E1() As Integer
ReDim a(1 To x)
ReDim B1(1 To Y)
ReDim C1(1 To x)
ReDim E1(1 To x)
Dim I, J, C2, CJ, JW
For J = Y To 1 Step -1 'D2
JW = 0 '进位清0
B1(J) = Mid$(D4, J, 1) '每位数
For I = x To 1 Step -1  'D1
   a(I) = Mid$(d3, I, 1) '每位数
   C1(I) = a(I) + B1(I) + JW '计算jia
   JW = C1(I) \ 10
   E1(I) = C1(I) Mod 10
  Next
  Next
  For r = 1 To x
  If JW = 0 Then
  MPC1 = MPC1 & E1(r)
  Else
  jc = jc & E1(r)
  MPC1 = JW & jc
  End If
  Next
  
End Function
回复 支持 反对

使用道具 举报

 楼主| 发表于 2021-3-23 12:52 | 显示全部楼层
乘法程序速度又提高了一点,2444位的数*2444位的数得到了4888位的数,用时1.546997秒(优化了一下程序,3位一组的,把3位当一位数计算,去掉了快速逆变换可调用程序,直接把逆变换加在主程序中,居然速度提高了。)除法慢一点4888位的数除以2444位的数得到2444位,用时548.1836秒(合9分钟多一点点,不比模仿手工除法快吗?)

可能是还有提升空间的!!
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|手机版|小黑屋|数学中国 ( 京ICP备05040119号 )

GMT+8, 2024-3-29 06:34 , Processed in 0.066406 second(s), 14 queries .

Powered by Discuz! X3.4

Copyright © 2001-2020, Tencent Cloud.

快速回复 返回顶部 返回列表