学习笔记 - 快速傅里叶变换 / 大数A * B的另一种解法
阅读原文时间:2023年07月09日阅读:2

转:

学习笔记 - 快速傅里叶变换 / 大数A * B的另一种解法

###

文章目录

  • 前言
  • * Fast Fast TLE
  • 一、FFT是什么?
  • 二、FFT可以干什么?
  • * 1.多项式乘法
    • 2.大数乘法
  • 三、FFT怎么做?
  • * 1. 系数表示法和点值表示法
    • 3.如何巧妙地进行DFT? - FFT的高明之处
    • 4. *单位复根* - FFT的高明之处(二)
    • 5. Inverse Fast Fourier Transform - 从点值表示法回到系数表示法
  • 四、FFT的递归实现
  • 五、大数乘法的实现(FFT版)
  • 后话(废话)

前言

原本做大数乘法的时候是想偷懒的, 就百度了下大数 A * B 的代码, 无意中发现有使用FFT的做法, 于是便开始了学习 (受苦) , 本文用以记录算法应用层面上我对FFT的理解和一个手搓的大数乘法模版.

需要一定的复数知识, 但不多.

笔者对FFT的理解较为粗浅, 因此在下文可能大量略过某些部分的详细证明, 请见谅!

一、FFT是什么?

快速傅里叶变换( Fast Fourier Transform ), 是一种对离散傅里叶变换( Discrete Fourier Transform )的优化, 可令离散傅里叶变换的运算量大大降低, 在取样点数量较多的情况下可显著提升效率

二、FFT可以干什么?

实际上, FFT可以大大提升卷积1的计算效率, 而当卷积的两个函数为多项式的时候, FFT则能在多项式乘法上大显神威.

关于乘法, 朴素的算法则是模拟我们自小学会的竖式乘法, 不难证明这样的算法的复杂度为O(n2), 然而这样的算法不仅容易精度溢出, 效率也不够高, 应用FFT则可以使算法的复杂度降低至O(nlogn).

观察一个数字, 如123456, 换种角度思考, 可以将其转化为以下的式子:

123456

=

1

1

0

5

+

2

1

0

4

+

3

1

0

3

+

4

1

0

2

+

5

1

0

1

+

6

1

0

0

123456 = 1 *10^5+2*10^4+3*10^3+4*10^2+5*10^1+6*10^0

123456=1∗105+2∗104+3∗103+4∗102+5∗101+6∗100再转化一次, 令x = 10, 则可以得到:

1

x

5

+

2

x

4

+

3

x

3

+

4

x

2

+

5

x

1

+

6

x

0

1 *x^5+2*x^4+3*x^3+4*x^2+5*x^1+6*x^0

1∗x5+2∗x4+3∗x3+4∗x2+5∗x1+6∗x0
从这个角度思考, 大数乘法不正是x = 10的情况下的多项式乘法吗?

三、FFT怎么做?

PS:强烈建议您通过这个视频进行学习并在需要的情况下结合下文理解, 个人认为这个视频讲得相当易懂

任取一个n-1阶2 多项式, 将其表示为

a

0

+

a

1

x

+

a

2

x

2

+

.

.

.

+

a

n

1

x

n

1

a_0+a_1x+a_2x^2+…+a_{n-1}x^{n-1}

a0​+a1​x+a2​x2+…+an−1​xn−1实际上我们已经得到了它的系数表示, 即

[

a

0

a

1

a

2

.

.

.

a

n

1

]

left[ begin{array}{l} a_0&a_1&a_2&…&a_{n-1} end{array} right]

[a0​​a1​​a2​​…​an−1​​]
那么什么是它的点值表示法呢?

对一个多项式

f

(

x

)

=

a

0

+

a

1

x

+

a

2

x

2

+

.

.

.

+

a

n

1

x

n

1

f(x)= a_0+a_1x+a_2x^2+…+a_{n-1}x^{n-1}

f(x)=a0​+a1​x+a2​x2+…+an−1​xn−1任取不小于n个xi带入并计算得到f(xi)即为其点值表示法:

{

(

x

0

,

f

(

x

0

)

,

(

x

1

,

f

(

x

1

)

,

,

(

x

n

1

,

f

(

x

n

1

)

}

left{ (x_0,f(x_0),(x_1,f(x_1), dots,(x_{n-1},f(x_{n-1}) right}

{(x0​,f(x0​),(x1​,f(x1​),…,(xn−1​,f(xn−1​)}

将其化为矩阵形式, 有

[

f

(

x

1

)

f

(

x

2

)

f

(

x

1

)

]

=

[

1

x

0

x

0

2

x

0

n

1

1

x

1

x

1

2

x

1

n

1

1

x

n

1

x

n

1

2

x

n

1

n

1

]

[

a

0

a

1

a

n

1

]

left[ begin{array}{l} f(x_1)\ f(x_2)\ vdots \ f(x-1) end{array} right] = left[ begin{array}{l} 1&x_0&x_0^2&dots&x_{0}^{n-1}\ 1&x_1&x_1^2&dots&x_{1}^{n-1}\ dots&dots&dots&vdots&dots\ 1&x_{n-1}&x_{n-1}^2&dots&x_{n-1}^{n-1} end{array} right] left[ begin{array}{l} a_0\ a_1\ vdots \ a_{n-1} end{array} right]

⎣⎢⎢⎢⎡​f(x1​)f(x2​)⋮f(x−1)​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​11…1​x0​x1​…xn−1​​x02​x12​…xn−12​​……⋮…​x0n−1​x1n−1​…xn−1n−1​​⎦⎥⎥⎥⎤​⎣⎢⎢⎢⎡​a0​a1​⋮an−1​​⎦⎥⎥⎥⎤​
PS:从 系数表示法 得到 点值表示法 的这一步即为DFT(离散傅里叶变换)

由范德蒙德行列式易证得当x0 ~ xn-1互不相同时矩阵可逆, 即当所取的xi不小于矩阵的阶加一时, 点值表示法和系数表示法一一对应.

h

(

x

)

=

f

(

x

)

g

(

x

)

h(x)=f(x)g(x)

h(x)=f(x)g(x)
在点值表示法的基础上, 要计算h(x)的表达式, 只需要取相同的xi得出f(xi)和g(xi)并相乘, 再从(xi ,f(xi)g(xi))的点值表示法转换成系数表示法即可得到h(x)的值.

那么如何从点值表示法转换成系数表示法呢? 仔细观察上面的矩阵表达式不难得出:

{

(

x

0

,

f

(

x

0

)

,

(

x

1

,

f

(

x

1

)

,

,

(

x

n

1

,

f

(

x

n

1

)

}

left{ (x_0,f(x_0),(x_1,f(x_1), dots,(x_{n-1},f(x_{n-1}) right}

{(x0​,f(x0​),(x1​,f(x1​),…,(xn−1​,f(xn−1​)}

将其化为矩阵形式, 有

[

a

0

a

1

a

n

1

]

=

[

1

x

0

x

0

2

x

0

n

1

1

x

1

x

1

2

x

1

n

1

1

x

n

1

x

n

1

2

x

n

1

n

1

]

1

[

f

(

x

1

)

f

(

x

2

)

f

(

x

1

)

]

left[ begin{array}{l} a_0\ a_1\ vdots \ a_{n-1} end{array} right] = left[ begin{array}{l} 1&x_0&x_0^2&dots&x_{0}^{n-1}\ 1&x_1&x_1^2&dots&x_{1}^{n-1}\ dots&dots&dots&vdots&dots\ 1&x_{n-1}&x_{n-1}^2&dots&x_{n-1}^{n-1} end{array} right]^{-1} left[ begin{array}{l} f(x_1)\ f(x_2)\ vdots \ f(x-1) end{array} right]

⎣⎢⎢⎢⎡​a0​a1​⋮an−1​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​11…1​x0​x1​…xn−1​​x02​x12​…xn−12​​……⋮…​x0n−1​x1n−1​…xn−1n−1​​⎦⎥⎥⎥⎤​−1⎣⎢⎢⎢⎡​f(x1​)f(x2​)⋮f(x−1)​⎦⎥⎥⎥⎤​
PS:从 点值表示法 得到 系数表示法 的这一步即为Inverse DFT(逆离散傅里叶变换)

对上述的DFT过程, 一个通常的想法是在x轴上随便挑几个点代入计算, 比如说x0= 0, x1 = 1, …然而这样做的话又重新掉入了O(n2)的窠臼之中, 我们依旧需要对n个甚至以上的点计算n次, 而要探讨FFT的做法, 让我们先着眼于两个简单的多项式

P

1

=

x

2

P_1=x^2

P1​=x2
因为这是一个偶函数, 当我们确定f(xi)时通过偶函数的性质可以立刻确定f(-xi) = f(xi), 这样我们任取一个点, 立刻就可以确定其相反数的函数值了.

而对于

P

2

=

x

3

P_2= x^3

P2​=x3
和P1不同, 这是个奇函数, 因此f(-xi) = -f(xi).

据此, 我们不妨将一个一般的多项式(设n为偶数)分解为奇函数的部分和偶函数的部分, 再从奇函数部分提取一个x:

P

(

x

)

=

a

0

+

a

1

x

+

a

2

x

2

+

.

.

.

+

a

n

1

x

n

1

P

e

(

x

)

=

a

0

+

a

2

x

2

+

.

.

.

+

a

n

2

x

n

2

,

P

o

(

x

)

=

a

1

+

a

3

x

2

+

.

.

.

+

a

n

1

x

n

2

P

(

x

i

)

=

P

e

(

x

i

)

+

x

i

P

o

(

x

i

)

,

P

(

x

i

)

=

P

e

(

x

i

)

x

i

P

o

(

x

i

)

P(x) = a_0+a_1x+a_2x^2+…+a_{n-1}x^{n-1}\ 令P_e(x) = a_0+ a_2x^2+…+a_{n-2}x^{n-2} , P_o(x) = a_1+a_3x^2+…+a_{n-1}x^{n-2}\ 则P(x_i)=P_e(x_i)+x_iP_o(x_i), \P(-x_i)=P_e(x_i)-x_iP_o(x_i)

P(x)=a0​+a1​x+a2​x2+…+an−1​xn−1令Pe​(x)=a0​+a2​x2+…+an−2​xn−2,Po​(x)=a1​+a3​x2+…+an−1​xn−2则P(xi​)=Pe​(xi​)+xi​Po​(xi​),P(−xi​)=Pe​(xi​)−xi​Po​(xi​)
若令 p = x2, 则Pe(x)和Po(x)则简化为以下式子:

P

e

(

x

)

=

a

0

+

a

2

x

2

+

.

.

.

+

a

n

2

x

n

2

P

o

(

x

)

=

a

1

+

a

3

x

2

+

.

.

.

+

a

n

1

x

n

2

p

=

x

2

P

e

(

p

)

=

a

0

+

a

2

p

+

.

.

.

+

a

n

2

p

n

2

2

P

o

(

p

)

=

a

1

+

a

3

p

+

.

.

.

+

a

n

1

p

n

2

2

P_e(x) = a_0+ a_2x^2+…+a_{n-2}x^{n-2}\ P_o(x) = a_1+a_3x^2+…+a_{n-1}x^{n-2}\ Updownarrow 令p=x^2\ P_e(p) = a_0 + a_2p+…+a_{n-2}p^{frac{n-2}{2}}\ P_o(p) = a_1 + a_3p+…+a_{n-1}p^{frac{n-2}{2}}\

Pe​(x)=a0​+a2​x2+…+an−2​xn−2Po​(x)=a1​+a3​x2+…+an−1​xn−2⇕令p=x2Pe​(p)=a0​+a2​p+…+an−2​p2n−2​Po​(p)=a1​+a3​p+…+an−1​p2n−2​
如果此时将Pe和Po视为新的式子, 重复上述操作则可将Pe和Po继续降阶直到Pe和Po的式子的阶降至0, 那么此时Pe和Po的值便为常数a0和a1.

以多项式P(x) = a0 + a1x + a2x2 + a3x2为例, DFT过程需要构造4个根x1, x2, x3, x4并带入求值, 以上述方法降阶则有:

第一步:

P

0

(

x

)

=

a

0

+

a

1

x

+

a

2

x

2

+

a

3

x

3

t

=

x

2

P

e

0

(

t

)

=

a

0

+

a

2

t

P

o

0

(

t

)

=

a

1

+

a

3

t

P

0

(

x

i

)

=

P

e

0

(

x

i

2

)

+

x

i

P

o

0

(

x

i

2

)

,

P

0

(

x

i

)

=

P

e

0

(

x

i

2

)

x

i

P

o

0

(

x

i

2

)

P_0(x) = a_0+a_1x+a_2x^2+a_3x^3\ Updownarrow令t = x^2\ P_{e0}(t) = a_0+a_2t\ P_{o0}(t)= a_1 + a_3t\\ P_0(x_i)=P_{e0}(x_i^2)+x_iP_{o0}(x_i^2), \P_0(-x_i)=P_{e0}(x_i^2)-x_iP_{o0}(x_i^2)

P0​(x)=a0​+a1​x+a2​x2+a3​x3⇕令t=x2Pe0​(t)=a0​+a2​tPo0​(t)=a1​+a3​tP0​(xi​)=Pe0​(xi2​)+xi​Po0​(xi2​),P0​(−xi​)=Pe0​(xi2​)−xi​Po0​(xi2​)


第二步:

P

1

=

P

e

0

(

t

)

=

a

0

+

a

2

t

,

P

2

=

P

o

0

(

t

)

=

a

1

+

a

3

t

,

P

e

1

(

t

)

=

a

0

,

P

e

2

(

t

)

=

a

1

,

P

o

1

(

t

)

=

a

2

P

e

2

(

t

)

=

a

3

P

1

(

x

i

2

)

=

P

e

1

(

x

i

2

)

+

x

i

2

P

o

1

(

x

i

2

)

=

a

0

+

x

i

2

a

2

,

P

2

(

x

i

2

)

=

P

e

2

(

x

i

2

)

+

x

i

2

P

o

2

(

x

i

2

)

=

a

1

+

x

i

2

a

3

,

P

1

(

x

i

2

)

=

P

e

1

(

x

i

2

)

x

i

2

P

o

1

(

x

i

2

)

=

a

0

x

i

2

a

2

P

2

(

x

i

2

)

=

P

e

2

(

x

i

2

)

x

i

2

P

o

2

(

x

i

2

)

=

a

1

x

i

2

a

3

begin{array}{c|c} 令P_1=P_{e0}(t) = a_0+a_2t, & 令P_2 = P_{o0}(t)= a_1 + a_3t,\ 则P_{e1}(t) = a_0, & 则P_{e2}(t) = a_1,\ P_{o1}(t) = a_2 & P_{e2}(t) = a_3\ P_1(x_i^2) = P_{e1}(x_i^2)+x_i^2P_{o1}(x_i^2) = a_0 + x_i^2cdot a_2, & P_2(x_i^2) = P_{e2}(x_i^2)+x_i^2P_{o2}(x_i^2) = a_1 + x_i^2cdot a_3,\ P_1(-x_i^2) = P_{e1}(-x_i^2)-x_i^2P_{o1}(x_i^2) = a_0 - x_i^2cdot a_2 & P_2(-x_i^2) = P_{e2}(-x_i^2)-x_i^2P_{o2}(x_i^2)=a_1 - x_i^2cdot a_3\ end{array}

令P1​=Pe0​(t)=a0​+a2​t,则Pe1​(t)=a0​,Po1​(t)=a2​P1​(xi2​)=Pe1​(xi2​)+xi2​Po1​(xi2​)=a0​+xi2​⋅a2​,P1​(−xi2​)=Pe1​(−xi2​)−xi2​Po1​(xi2​)=a0​−xi2​⋅a2​​令P2​=Po0​(t)=a1​+a3​t,则Pe2​(t)=a1​,Pe2​(t)=a3​P2​(xi2​)=Pe2​(xi2​)+xi2​Po2​(xi2​)=a1​+xi2​⋅a3​,P2​(−xi2​)=Pe2​(−xi2​)−xi2​Po2​(xi2​)=a1​−xi2​⋅a3​​
如此, 对P0(x)的快速傅里叶变换过程已经呼之欲出了, 然而这里却有一个严重的问题亟待解决……

也许你已经注意到了其中的问题, 让我们回到上述第二步, 该如何取xi的值才能使得x02与x12互为相反数呢, 显然在实数域中这个问题是无法被解决的.

这就是FFT的第二个高明之处, 使用n次单位复根来解决问题.

让我们继续使用上面那个多项式P0(x), 我们需要4个数作为x代入多项式并求值, 从而得出P0(x)的点值表达式, 那么我们不妨使用4次单位复根(1, -1, i, -i)作为 xi :

P

0

(

x

i

)

,

P

0

(

x

)

=

a

0

+

a

1

x

+

a

2

x

2

+

a

3

x

3

,

x

i

=

[

1

,

1

,

i

,

i

]

P

0

(

1

)

=

P

e

0

(

1

)

+

1

P

o

0

(

1

)

,

P

0

(

1

)

=

P

e

0

(

1

)

1

P

o

0

(

1

)

P

0

(

i

)

=

P

e

0

(

1

)

+

i

P

o

0

(

1

)

,

P

0

(

1

)

=

P

e

0

(

1

)

i

P

o

0

(

1

)

求P_0(x_i)的值, 其中P_0(x) = a_0+a_1x+a_2x^2+a_3x^3 ,x_i = left[ 1,-1,i,-i right]\ P_0(1)=P_{e0}(1)+1 cdot P_{o0}(1), \P_0(-1)=P_{e0}(1)-1 cdot P_{o0}(1)\ P_0(i)=P_{e0}(-1)+i cdot P_{o0}(-1), \P_0(-1)=P_{e0}(-1)-i cdot P_{o0}(-1)\

求P0​(xi​)的值,其中P0​(x)=a0​+a1​x+a2​x2+a3​x3,xi​=[1,−1,i,−i]P0​(1)=Pe0​(1)+1⋅Po0​(1),P0​(−1)=Pe0​(1)−1⋅Po0​(1)P0​(i)=Pe0​(−1)+i⋅Po0​(−1),P0​(−1)=Pe0​(−1)−i⋅Po0​(−1)


P

e

0

(

1

)

=

P

1

(

1

)

=

P

e

1

(

1

)

+

1

P

o

1

(

1

)

=

a

0

+

1

a

2

,

P

e

0

(

1

)

=

P

1

(

1

)

=

P

e

1

(

1

)

1

P

o

1

(

1

)

=

a

0

1

a

2

P

o

0

(

1

)

=

P

1

(

1

)

=

P

e

2

(

1

)

+

1

P

o

2

(

1

)

=

a

1

+

1

a

3

,

P

o

0

(

1

)

=

P

1

(

1

)

=

P

e

2

(

1

)

1

P

o

2

(

1

)

=

a

1

1

a

3

P_{e0}(1) = P_1(1) = P_{e1}(1)+1cdot P_{o1}(1) = a_0 + 1cdot a_2, \ P_{e0}(-1) = P_1(-1) = P_{e1}(-1)-1cdot P_{o1}(1) = a_0 - 1cdot a_2\ P_{o0}(1) = P_1(1) = P_{e2}(1)+1cdot P_{o2}(1) = a_1 + 1cdot a_3, \ P_{o0}(-1) = P_1(-1) = P_{e2}(-1)-1cdot P_{o2}(1) = a_1 - 1cdot a_3

Pe0​(1)=P1​(1)=Pe1​(1)+1⋅Po1​(1)=a0​+1⋅a2​,Pe0​(−1)=P1​(−1)=Pe1​(−1)−1⋅Po1​(1)=a0​−1⋅a2​Po0​(1)=P1​(1)=Pe2​(1)+1⋅Po2​(1)=a1​+1⋅a3​,Po0​(−1)=P1​(−1)=Pe2​(−1)−1⋅Po2​(1)=a1​−1⋅a3​
再将值代入上式, 整个FFT过程便结束了.

当然单位复根的其他(我并不熟知)的性质(如可约引理, 等分引理, 求和引理等), 造就了其十分适合用于快速傅里叶变换中.

简单来说, 对于n次单位复根(其中n是2的正整数次幂), 他们之间两两正负配对(第k个根与第k + n/2个根)且平方后同样两两正负配对, 这令其十分适合用于递归当中.

当然, 这些性质在n是2的正整数次幂时才生效, 当n不是2的正整数次幂时, 可自行将多项式的次数补足, 此时该项系数为0即可.

当然由于我十分的菜, 具体的证明就此略过了, 各位可以看看我推荐的那个视频或者自行寻找资料进行证明(

这一步骤之简单你绝对无法想象, 观察上文的IDFT的实现:

[

a

0

a

1

a

n

1

]

=

[

1

x

0

x

0

2

x

0

n

1

1

x

1

x

1

2

x

1

n

1

1

x

n

1

x

n

1

2

x

n

1

n

1

]

1

[

f

(

x

1

)

f

(

x

2

)

f

(

x

1

)

]

left[ begin{array}{l} a_0\ a_1\ vdots \ a_{n-1} end{array} right] = left[ begin{array}{l} 1&x_0&x_0^2&dots&x_{0}^{n-1}\ 1&x_1&x_1^2&dots&x_{1}^{n-1}\ dots&dots&dots&vdots&dots\ 1&x_{n-1}&x_{n-1}^2&dots&x_{n-1}^{n-1} end{array} right]^{-1} left[ begin{array}{l} f(x_1)\ f(x_2)\ vdots \ f(x-1) end{array} right]

⎣⎢⎢⎢⎡​a0​a1​⋮an−1​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​11…1​x0​x1​…xn−1​​x02​x12​…xn−12​​……⋮…​x0n−1​x1n−1​…xn−1n−1​​⎦⎥⎥⎥⎤​−1⎣⎢⎢⎢⎡​f(x1​)f(x2​)⋮f(x−1)​⎦⎥⎥⎥⎤​
让我们对第二个矩阵进行求逆:

[

1

x

0

x

0

2

x

0

n

1

1

x

1

x

1

2

x

1

n

1

1

x

n

1

x

n

1

2

x

n

1

n

1

]

1

=

[

1

1

1

1

1

ω

ω

2

ω

n

1

1

ω

n

1

ω

2

(

n

1

)

ω

(

n

1

)

(

n

1

)

]

1

=

1

n

[

1

1

1

1

1

ω

1

ω

2

ω

(

n

1

)

1

ω

(

n

1

)

ω

2

(

n

1

)

ω

(

n

1

)

(

n

1

)

]

left[ begin{array}{l} 1&x_0&x_0^2&dots&x_{0}^{n-1}\ 1&x_1&x_1^2&dots&x_{1}^{n-1}\ dots&dots&dots&vdots&dots\ 1&x_{n-1}&x_{n-1}^2&dots&x_{n-1}^{n-1} end{array} right]^{-1} = left[ begin{array}{l} 1&1 &1&dots&1\ 1&omega &omega ^2 &dots&omega ^{n-1}\ dots&dots&dots&vdots&dots\ 1&omega ^{n-1}&omega ^{2(n-1)}&dots&omega^{(n-1)(n-1)} end{array} right]^{-1}=frac{1}{n}left[ begin{array}{l} 1&1 &1&dots&1\ 1&omega^{-1} &omega ^{-2} &dots&omega ^{-(n-1)}\ dots&dots&dots&vdots&dots\ 1&omega ^{-(n-1)}&omega ^{-2(n-1)}&dots&omega^{-(n-1)(n-1)} end{array} right]

⎣⎢⎢⎢⎡​11…1​x0​x1​…xn−1​​x02​x12​…xn−12​​……⋮…​x0n−1​x1n−1​…xn−1n−1​​⎦⎥⎥⎥⎤​−1=⎣⎢⎢⎢⎡​11…1​1ω…ωn−1​1ω2…ω2(n−1)​……⋮…​1ωn−1…ω(n−1)(n−1)​⎦⎥⎥⎥⎤​−1=n1​⎣⎢⎢⎢⎡​11…1​1ω−1…ω−(n−1)​1ω−2…ω−2(n−1)​……⋮…​1ω−(n−1)…ω−(n−1)(n−1)​⎦⎥⎥⎥⎤​其中

x

k

=

ω

k

,

ω

=

e

2

π

i

n

x_k = omega_k, omega = e^{frac{2pi i}{n}}

xk​=ωk​,ω=en2πi​
也就是说:

[

p

0

p

1

p

n

1

]

=

1

n

[

1

1

1

1

1

ω

ω

2

ω

n

1

1

ω

n

1

ω

2

(

n

1

)

ω

(

n

1

)

(

n

1

)

]

[

P

(

x

1

)

P

(

x

2

)

P

(

x

1

)

]

left[ begin{array}{l} p_0\ p_1\ vdots \ p_{n-1} end{array} right] = frac{1}{n}left[ begin{array}{l} 1&1 &1&dots&1\ 1&omega &omega ^2 &dots&omega ^{n-1}\ dots&dots&dots&vdots&dots\ 1&omega ^{n-1}&omega ^{2(n-1)}&dots&omega^{(n-1)(n-1)} end{array} right] left[ begin{array}{l} P(x_1)\ P(x_2)\ vdots \ P(x-1) end{array} right]

⎣⎢⎢⎢⎡​p0​p1​⋮pn−1​​⎦⎥⎥⎥⎤​=n1​⎣⎢⎢⎢⎡​11…1​1ω…ωn−1​1ω2…ω2(n−1)​……⋮…​1ωn−1…ω(n−1)(n−1)​⎦⎥⎥⎥⎤​⎣⎢⎢⎢⎡​P(x1​)P(x2​)⋮P(x−1)​⎦⎥⎥⎥⎤​
不难发现, 要实现IFFT, 只需在FFT的基础上将ω取倒数, 并在最终结果乘以n分之一即可.

PS:我觉得自己写的FFT解释很烂, 在此建议您如果 (肯定) 有看不懂的部分结合刚才推荐的视频学习

四、FFT的递归实现

据此, 我们终于可以写出FFT的代码了, 我手写的递归实现如下:

#include 

using namespace std;

typedef unsigned long long ull;
const double Pi = acos(-1.0);    //将常数Pi的值设为arccos(-1)

vector<complex<double> > FFT(vector<complex<double> > &a, int flag) {
    //第一个参数为一个多项式的系数, 以次数从小到大的顺序, 向量中每一项的实部为该项系数
    //flag == 1->FFT, flag == -1->IFFT
    int n = a.size();
    if(n == 1)        //如果当前多项式仅有常数项时直接返回多项式的值
        return a;
    vector<complex<double> > Pe, Po;        //即原文中的Pe与Po的系数表示法
    complex<double> omega(cos(2.0 * Pi / n), sin(flag * 2.0 * Pi / n)), cur(1, 0);        //omega为第一个n次复根, cur为第零个n次复根即1
    for (unsigned int i = 0; i < n; i++) {
        if(i & 1)
            Po.push_back(a[i]);
        else
            Pe.push_back(a[i]);
    }
    vector<complex<double> > ye = FFT(Pe, flag), yo = FFT(Po, flag), y(n);        //递归求得ye=Pe(xi), yo=Po(xi)
    for (int i = 0; i < n / 2; i++) {
        y[i] = ye[i] + cur * yo[i];                //求出P(xi)
        y[i + n / 2] = ye[i] - cur * yo[i];        //由单位复根的性质可知第k个根与第k + n/2个根互为相反数
        cur *= omega;                            //cur * omega得到下一个复根
    }
    return y;            //返回最终的系数
}

当然需要注意的是, 如果当前进行的是IFFT时, 需要将最终结果除以n, 其中n为不小于原多项式阶+1的最小2的整数次幂.

因此你需要在上述代码的基础上增加一个代码, 向原系数向量补0直到向量中元素的个数为二的整数幂为止.

五、大数乘法的实现(FFT版)

了解了FFT的实现方法, 大数乘法的实现可以说是再简单不过了(毕竟前面的鬼东西都看懂了), 这里便直接贴上源代码了,当然有一个非常需要注意的地方, 最终获取的多项式系数是有可能大于等于10的, 此时需要做进位处理!

#include 

#define sync ios::sync_with_stdio(false)

using namespace std;

typedef unsigned long long ull;
const double Pi = acos(-1.0);    //将常数Pi的值设为arccos(-1)

vector<complex<double> > FFT(vector<complex<double> > &a, int flag) {
    //第一个参数为一个多项式的系数, 以次数从小到大的顺序, 向量中每一项的实部为该项系数
    //flag == 1->FFT, flag == -1->IFFT
    int n = a.size();
    if(n == 1)        //如果当前多项式仅有常数项时直接返回多项式的值
        return a;
    vector<complex<double> > Pe, Po;        //即原文中的Pe与Po的系数表示法
    complex<double> omega(cos(2.0 * Pi / n), sin(flag * 2.0 * Pi / n)), cur(1, 0);        //omega为第一个n次复根, cur为第零个n次复根即1
    for (unsigned int i = 0; i < n; i++) {
        if(i & 1)
            Po.push_back(a[i]);
        else
            Pe.push_back(a[i]);
    }
    vector<complex<double> > ye = FFT(Pe, flag), yo = FFT(Po, flag), y(n);        //递归求得ye=Pe(xi), yo=Po(xi)
    for (int i = 0; i < n / 2; i++) {
        y[i] = ye[i] + cur * yo[i];                //求出P(xi)
        y[i + n / 2] = ye[i] - cur * yo[i];        //由单位复根的性质可知第k个根与第k + n/2个根互为相反数
        cur *= omega;                            //cur * omega得到下一个复根
    }
    return y;            //返回最终的系数
}

vector<complex<double> > read() {    //获取系数, 注意需要从右向左获取(从0次幂的系数开始)
    string num;
    cin >> num;
    vector<complex<double> > result;
    for (int i = num.size() - 1; i >= 0; i--) {
        complex<double> tmp(num[i] - '0', 0);
        result.push_back(tmp);
    }
    return result;
}

void solve(vector<complex<double> > &a, vector<complex<double> > &b) {        //多项式系数表达式的修饰
    complex<double> tmp(0, 0);
    int sum = a.size() + b.size();
    while (a.size() < sum)
        a.push_back(tmp);
    while (b.size() < sum)
        b.push_back(tmp);
    //如果两式的阶不同, 先补齐
    int temp = 1;
    while (temp < a.size())
        temp <<= 1;
    //获取不小于n的最小2的整数次幂
    while (a.size() < temp) {
        a.push_back(tmp);
        b.push_back(tmp);
    }
    //补齐
}

int main() {
    sync;
    vector<complex<double> > num1 = read(), num2 = read(), tmp1, tmp2, mid, ans;
    solve(num1, num2);
    tmp1 = FFT(num1, 1), tmp2 = FFT(num2, 1);
    num1.clear();
    num2.clear();
    for (int i = 0; i < tmp1.size(); i++)
        mid.push_back(tmp1[i] * tmp2[i]);
    tmp1.clear();
    tmp2.clear();
    ans = FFT(mid, -1);
    bool Ans = false;
    int add = 0;
    string final;
    //进位处理
    for (int i = 0; i < ans.size(); i++) {
        int op = round(round(ans[i].real()) / ans.size()) + add;
        add = 0;
        if(op >= 10)
            add = op / 10;
        final += op % 10 + '0';
    }
    if(add > 0)
        final += add % 10 + '0';
    for (int i = final.size() - 1; i >= 0; i--) {
        if(final[i] != '0')
            Ans = true;
        else if(!Ans)
            continue;
        cout << final[i];
    }
    cout << 'n';
}

后话(废话)

  1. 历时一个寒假! 终于把这篇万字长文(x)写完了, 多谢某csome同学的鼓励(催更) !
  2. 理解FFT足足花了2天时间, 好难……
  3. 如果博客有错漏, 敬请指出并不吝赐教, 感谢!

转:

学习笔记 - 快速傅里叶变换 / 大数A * B的另一种解法

手机扫一扫

移动阅读更方便

阿里云服务器
腾讯云服务器
七牛云服务器

你可能感兴趣的文章