数学中国

 找回密码
 注册
搜索
热搜: 活动 交友 discuz
查看: 10095|回复: 45

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

[复制链接]
发表于 2020-8-30 23:27 | 显示全部楼层 |阅读模式
貌似一种新的乘法快速计算方法已经提交论文,理论上可以达到大数乘法的效率极限。
文章链接:多项式乘法到快速傅里叶变换;此文介绍的非常详细,极力推荐。
文章链接:使用快速傅里叶变换计算大整数乘法;快速傅里叶变换,使用算法设计思想中的分治法,降低傅里叶变换的时间复杂度到 O(N logN)。
傅里叶变换,算法的时间复杂度还是 O(N2)。关键在于:直接进行离散傅里叶变换的计算复杂度是 O(N2)。快速傅里叶变换可以计算出与直接计算相同的结果,但只需要 O(N logN) 的计算复杂度。 N logN 和 N2 之间的差别是巨大的。例如,当 N = 106 时,在一个每秒运算百万次的计算机上,粗略地说,它们之间就是占用 30 秒 CPU 时间和两星期 CPU 时间的差别。
快速傅里叶变换的要点如下:一个界长为 N 的离散傅里叶变换可以重新写成两个界长各为 N/2 的离散傅里叶变换之和。其中一个变换由原来 N 个点中的偶数点构成,另一个变换由奇数点构成。这个过程可以递归地进行下去,直到我们将全部数据细分为界长为 1 的变换。
什么是界长为 1 的傅里叶变换呢?它正是把一个输入值复制成它的一个输出值的恒等运算。要实现以上算法,最容易的情况是原始的 N 为 2 的整幂次项,如果数据集的界长不是 2 的幂次时,则可添上一些零值,直到 2 的下一幂次。在这个算法中,每递归一次需 N 阶运算,共需要 log N 次递归,所以快速傅里叶变换算法的时间复杂度是 O(N logN)。

关键是把各位数字看成了多项式系数?如何把多项式表示法转换为点值表示法?
系数的数字都不是相等的,不规则变化的,与单位圆上均匀分布的点值不同,半径相等?
如何利用对称性,周期性?

不懂,大大的不懂!
 楼主| 发表于 2020-9-5 22:12 | 显示全部楼层
好像明白了一点:
我们以两个 3 ( N = 3 ) 位数 a = 678 和 b = 432 的乘积来说明以上过程吧。
a = (678)'10 = 6x10^2 + 7x10^1 + 8x10^0
b = (432)'10 = 4x10^2 + 3x10^1 + 2x10^0
由此:
c = a x b = 10^4xc4 + 10^3xc3 + 10^2xc2 + 10^1xc1 + 10^0xc0
   = 10000x24 + 1000x46 + 100x65 + 10x38 + 1x16
   = 292896
如果按以上方法计算大整数的乘法,时间复杂度是 O(N^2)。
在我们的例子中,乘积 c = 292896,共 6 位数字,N 需要扩展到 2^3 = 8。那么,向量 {ai} 和向量 {bj} 如下所示:
{ a7, a6, a5, a4, a3, a2, a1, a0 } = { 0, 0, 0, 0, 0, 6, 7, 8 }
{ b7, b6, b5, b4, b3, b2, b1, b0 } = { 0, 0, 0, 0, 0, 4, 3, 2 }

ω^0~ω^7分别是1,0.7-0.7i,-i,-0.7-0.7i,-1,-0.7+0.7i,i,0.7+0.7i.
可见ω^1和ω^7是共轭的,ω^2和ω^6是共轭的,等,这就是对称性。
ω^0和ω^8是相等的,ω^4和ω^12是相等的,等,这就是周期性。

向量A0就是8,7,6依次乘以ω^0,再加起来得到21.
以此类推得到:
{ A7, A6, A5, A4, A3, A2, A1, A0 } = { 12.9+10.9i, 2+7i, 3.1-1.1i, 7, 3.1+1.1i, 2-7i, 12.9-10.9i, 21 }
{ B7, B6, B5, B4, B3, B2, B1, B0 } = { 4.1+6.1i, -2+3i, -0.1-1.9i, 3, -0.1+1.9i, -2-3i, 4.1-6.1i, 9 }

现在,将她们逐项相乘得到向量 {Ck},即 { C7, C6, C5, C4, C3, C2, C1, C0 }
= { -13.6+123.4i, -25-8i, -2.4-5.8i, 21, -2.4+5.8i, -25+8i, -13.6-123.4i, 189 }

这样,我们就使用离散傅里叶变换和逆变换计算出了向量 {ai} 和向量 {bj} 的卷积向量 {ck},如下所示:
{ c7, c6, c5, c4, c3, c2, c1, c0 } = { 0, 0, 0, 0, 24, 46, 65, 38, 16 }

这些两位数错位相加,就得到292896。

问题是如何利用对称性和周期性?使这个算法的时间复杂度由 O(N^2),变成O(N logN)?

各位老师,请给与指点!谢谢!

注意到(可以观察到)向量Ai和向量Bj,都有共轭复数,对称出现的。

由向量Ai和向量Bj,得到向量Ck仍然需要相乘?虽然向量Ck也有共轭复数,且是对称出现的。

各位老师,请给与指点!非常感谢!
回复 支持 反对

使用道具 举报

 楼主| 发表于 2020-9-5 22:24 | 显示全部楼层
注意到:3*7=21,9*21=189.
所以,可能是对应项相乘的。
回复 支持 反对

使用道具 举报

 楼主| 发表于 2020-10-9 18:00 | 显示全部楼层
本帖最后由 ysr 于 2020-10-9 14:56 编辑

我们以两个 3 ( N = 3 ) 位数 a = 678 和 b = 432 的乘积来说明以上过程吧。
a = (678)'10 = 6x10^2 + 7x10^1 + 8x10^0
b = (432)'10 = 4x10^2 + 3x10^1 + 2x10^0
由此:
c = a x b = 10^4xc4 + 10^3xc3 + 10^2xc2 + 10^1xc1 + 10^0xc0
   = 10000x24 + 1000x46 + 100x65 + 10x38 + 1x16
   = 292896
如果按以上方法计算大整数的乘法,时间复杂度是 O(N^2)。
在我们的例子中,乘积 c = 292896,共 6 位数字,N 需要扩展到 2^3 = 8。那么,向量 {ai} 和向量 {bj} 如下所示:
{ a7, a6, a5, a4, a3, a2, a1, a0 } = { 0, 0, 0, 0, 0, 6, 7, 8 }
{ b7, b6, b5, b4, b3, b2, b1, b0 } = { 0, 0, 0, 0, 0, 4, 3, 2 }

ω^0~ω^7分别是1,0.7-0.7i,-i,-0.7-0.7i,-1,-0.7+0.7i,i,0.7+0.7i.
可见ω^1和ω^7是共轭的,ω^2和ω^6是共轭的,等,这就是对称性。
ω^0和ω^8是相等的,ω^4和ω^12是相等的,等,这就是周期性。

向量A0就是8,7,6依次乘以ω^0,再加起来得到21.
以此类推得到:
{ A7, A6, A5, A4, A3, A2, A1, A0 } = { 12.9+10.9i, 2+7i, 3.1-1.1i, 7, 3.1+1.1i, 2-7i, 12.9-10.9i, 21 }
{ B7, B6, B5, B4, B3, B2, B1, B0 } = { 4.1+6.1i, -2+3i, -0.1-1.9i, 3, -0.1+1.9i, -2-3i, 4.1-6.1i, 9 }

现在,将她们逐项相乘得到向量 {Ck},即 { C7, C6, C5, C4, C3, C2, C1, C0 }
= { -13.6+123.4i, -25-8i, -2.4-5.8i, 21, -2.4+5.8i, -25+8i, -13.6-123.4i, 189 }
这一步好像就是把前两个向量的对应项相乘得到的。

这样,我们就使用离散傅里叶变换和逆变换计算出了向量 {ai} 和向量 {bj} 的卷积向量 {ck},如下所示:
{ c7, c6, c5, c4, c3, c2, c1, c0 } = { 0, 0, 0, 0, 24, 46, 65, 38, 16 }

这些两位数错位相加,就得到292896。

问题是如何利用对称性和周期性?使这个算法的时间复杂度由 O(N^2),变成O(N logN)?

ω是如何计算的呢?
为了计算离散傅里叶变换,如当N=8时,我们令:
ω = e^(-2πi/N )= e^(-2πi/8) = e^(-πi/4) = cos(-π/4) + i sin(-π/4) = √2 / 2 - i √2 / 2 ≈ 0.7-0.7i
而ω^0~ω^8都以此为基础,以此类推,如ω^0=e^0=1,
ω^1=e^(-πi/4) = cos(-π/4) + i sin(-π/4) = √2 / 2 - i √2 / 2 ≈ 0.7-0.7i,
ω^2=e^(-πi/2)=0+i*(-1)=-i,
……。


这就是快速离散傅里叶变换在大数乘法中的应用。

回复 支持 反对

使用道具 举报

 楼主| 发表于 2020-10-9 18:10 | 显示全部楼层
总计进行了两次快速傅里叶变换(被乘数和乘数各一次),得到:
{ A7, A6, A5, A4, A3, A2, A1, A0 } = { 12.9+10.9i, 2+7i, 3.1-1.1i, 7, 3.1+1.1i, 2-7i, 12.9-10.9i, 21 }
{ B7, B6, B5, B4, B3, B2, B1, B0 } = { 4.1+6.1i, -2+3i, -0.1-1.9i, 3, -0.1+1.9i, -2-3i, 4.1-6.1i, 9 }

一次对应项相乘,得到向量 {Ck},即 { C7, C6, C5, C4, C3, C2, C1, C0 }
= { -13.6+123.4i, -25-8i, -2.4-5.8i, 21, -2.4+5.8i, -25+8i, -13.6-123.4i, 189 }

一次逆变换计算出了向量 {ai} 和向量 {bj} 的卷积向量 {ck},如下所示:
{ c7, c6, c5, c4, c3, c2, c1, c0 } = { 0, 0, 0, 0, 24, 46, 65, 38, 16 }

最后错位相加得到结果。
回复 支持 反对

使用道具 举报

 楼主| 发表于 2020-10-23 18:42 | 显示全部楼层
例如如下向量:
{ A7, A6, A5, A4, A3, A2, A1, A0 } = { 12.9+10.9i, 2+7i, 3.1-1.1i, 7, 3.1+1.1i, 2-7i, 12.9-10.9i, 21 }
显然以7为对称中心,对称项都是共轭复数,只有7和21是不对称的,如何快速得到?
先得到7,还是先得到 3.1+1.1i?然后由对称性和其它算法(而不再进行乘法)推演出来其他项?
回复 支持 反对

使用道具 举报

 楼主| 发表于 2020-12-25 00:38 | 显示全部楼层
改造后输出复数的蝶形运算程序代码及结果:
Private Sub Command1_Click()
Dim xr() As Double, a As String
a = Trim(Text1)
ReDim xr(0 To Len(a) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
  Next
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
Dim xi()
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)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  wr = 1
  wi = 0
  t = pi / le1
  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
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   p = p + le
Loop Until p > n - 1

wr = wr * w1r - wi * w1i
wi = wr * w1i - wi * w1r
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
   
   Text2 = Text2 & "  " & xr(i) & "+" & xi(i) & "i"
   Next
End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
End Sub
计算结果:
Text1=00000678,Text2= 21+0i  -4.24264068711929+3i  -1.31801948466055+-2.25i  -1.68198051533947+2.25i  -21+0i  4.24264068711929+-3i  1.31801948466055+2.25i  1.68198051533947+-2.25i
回复 支持 反对

使用道具 举报

 楼主| 发表于 2020-12-25 00:39 | 显示全部楼层
跟人家的例子(正确值)比较如下:
21+0i  -4.24264068711929+3i  -1.31801948466055+-2.25i  -1.68198051533947+2.25i  -21+0i  4.24264068711929+-3i  1.31801948466055+2.25i  1.68198051533947+-2.25i

12.9+10.9i, 2+7i, 3.1-1.1i, 7, 3.1+1.1i, 2-7i, 12.9-10.9i, 21

不一样,咋回事呢?希望老师指点!谢谢!
Text1=00000432,Text2=
9+0i  -2.82842712474619+2i  -0.146446609406727+-0.25i  -1.12132034355964+1.5i  -9+0i  2.82842712474619+-2i  0.146446609406727+0.25i  1.12132034355964+-1.5i
和人家的例子比较:
4.1+6.1i, -2+3i, -0.1-1.9i, 3, -0.1+1.9i, -2-3i, 4.1-6.1i, 9

不一样,显然有错误,咋回事呢?
回复 支持 反对

使用道具 举报

 楼主| 发表于 2020-12-25 00:46 | 显示全部楼层
仅输出模的程序,如下为代码和计算结果(有问题,不知道哪错了):
Private Sub Command1_Click()
Dim xr() As Double, a As String
a = Trim(Text1)
ReDim xr(0 To Len(a) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
  Next
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
Dim xi()
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)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  wr = 1
  wi = 0
  t = pi / le1
  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
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   p = p + le
Loop Until p > n - 1

wr = wr * w1r - wi * w1i
wi = wr * w1i - wi * w1r
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
   xr(i) = Sqr(xr(i) ^ 2 + xi(i) ^ 2)
   Text2 = Text2 & "  " & xr(i)
   Next
End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
End Sub
计算结果:
Text1=00000678,Text2= 21  5.19615242270664  2.60761871483253  2.80919177949488  21  5.19615242270664  2.60761871483253  2.80919177949488.
回复 支持 反对

使用道具 举报

 楼主| 发表于 2020-12-25 01:02 | 显示全部楼层
当Text1=00000678时,各点值对应的角度分别为:Text2=  3.14159265358979  1.5707963267949  0.785398163397448,
缺少0度角了,实数值21和9又是如何来的?可能是有问题了。
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2024-3-29 16:56 , Processed in 0.079102 second(s), 16 queries .

Powered by Discuz! X3.4

Copyright © 2001-2020, Tencent Cloud.

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