使用快速傅裡葉變換計算大整數乘法

 

我們知道，兩個 N 位數字的整數的乘法，如果使用常規的算法，時間復雜度是 O(N²)。然而，使用快速傅裡葉變換，時間復雜度可以降低到 O(N logN loglogN)。

 

假設我們要計算以下兩個 N 位數字的乘積：

a = (a_N-1a_N-2...a₁a₀)₁₀ = a_N-1x10^N-1 + a_N-2x10^N-2 + ... + a₁x10¹ + a₀x10⁰

b = (b_N-1b_N-2...b₁b₀)₁₀ = b_N-1x10^N-1 + b_N-2x10^N-2 + ... + b₁x10¹ + b₀x10⁰

將上面兩個式子相乘，得到以下公式 (共 2N - 1 項)：

c = a x b = c_2N-2x10^2N-2 + c_2N-3x10^2N-3 + ... + c₁x10¹ + c₀x10⁰

非常容易驗證，上式中的 c_k ( 0 ≤ k ≤ 2N-2 ) 滿足以下公式：

c_k = a₀xb_k + a₁xb_k-1 + ... + a_k-2xb₂ + a_k-1xb₁
      + a_kxb₀ + a_k+1xb_-1 + ... + a_N-2xb_-(N-2-k) + a_N-1xb_-(N-1-k)

上式共有 N 項，a_i 和 b_j 的下標 i 和 j 滿足 i + j = k。若不滿足 0 ≤ i, j ≤ N-1 時，則令 a_i = b_j = 0。

 

 

我們以兩個 3 ( N = 3 ) 位數 a = 678 和 b = 432 的乘積來說明以上過程吧。

a = (678)₁₀ = 6x10² + 7x10¹ + 8x10⁰

b = (432)₁₀ = 4x10² + 3x10¹ + 2x10⁰

由此：
 

c₀ = a₀xb₀ + a₁xb_-1 + a₂xb_-2 = 8x2 + 7x0 + 6x0 = 16 + 0 + 0 = 16
 

c₁ = a₀xb₁ + a₁xb₀ + a₂xb_-1 = 8x3 + 7x2 + 6x0 = 24 + 14 + 0 = 38
 

c₂ = a₀xb₂ + a₁xb₁ + a₂xb₀ = 8x4 + 7x3 +6x2 = 32 + 21 + 12 = 65
 

c₃ = a₀xb₃ + a₁xb₂ + a₂xb₁ = 8x0 + 7x4 + 6x3 = 0 + 28 + 18 = 46
 

c₄ = a₀xb₄ + a₁xb₃ + a₂xb₂ = 8x0 + 7x0 + 6x4 = 0 + 0 + 24 = 24

最後：

c = a x b = 10⁴xc₄ + 10³xc₃ + 10²xc₂ + 10¹xc₁ + 10⁰xc₀
   = 10000x24 + 1000x46 + 100x65 + 10x38 + 1x16
   = 292896

如果按以上方法計算大整數的乘法，時間復雜度是 O(N²)。

 

但是，我們注意到，向量 {c_k} 是向量 {a_i} 和向量 {b_j} 的卷積。根據卷積定理，向量卷積的離散傅裡葉變換是向量離散傅裡葉變換的乘積。於是，我們可以按照以下步驟來計算大整數乘法：

1. 分別求出向量 {a_i} 和向量 {b_j} 的離散傅裡葉變換 {A_i} 和 {B_j}。
2. 將 {A_i} 和 {B_j} 逐項相乘得到向量 {C_k}。
3. 對 {C_k} 求離散傅裡葉逆變換，得到的向量 {c_k} 就是向量 {a_i} 和向量 {b_j} 的卷積。
4. 對的向量 {c_k} 進行適當的進位就得到了大整數 a 和 b 的乘積 c。

對於復數向量 { x_N-1, ..., x₁, x₀ }，離散傅裡葉變換公式為：

離散傅裡葉逆變換公式為：

注意到離散傅裡葉逆變換除了指數的符號相反以及結果需要乘以歸一化因子 1/N 外，與離散傅裡葉變換是相同的。所以計算離散傅裡葉變換的程序稍做修改也可以用於計算逆變換。
 

 

在我們的例子中，乘積 c = 292896，共 6 位數字，N 需要擴展到 2³ = 8。那麼，向量 {a_i} 和向量 {b_j} 如下所示：

{ a₇, a₆, a₅, a₄, a₃, a₂, a₁, a₀ } = { 0, 0, 0, 0, 0, 6, 7, 8 }

{ b₇, b₆, b₅, b₄, b₃, b₂, b₁, b₀ } = { 0, 0, 0, 0, 0, 4, 3, 2 }

為了求出以上向量的離散傅裡葉變換，我們令

ω = 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

 

為了方便計算，我們預先求出 ω 的各次方，如下：

ω⁸ = ω⁰ = e⁰ = 1

ω⁹ = ω¹ = e^-πi/4 = cos(-π/4) + i sin(-π/4) ≈ 0.7-0.7i

ω¹⁰ = ω² = e^-πi/2 = cos(-π/2) + i sin(-π/2) = -i

ω¹¹ = ω³ = e^-3πi/4 = cos(-3π/4) + i sin(-3π/4) ≈ -0.7-0.7i

ω¹² = ω⁴ = e^-πi = cos(-π) + i sin(-π) = -1

ω¹³ = ω⁵ = e^-5πi/4 = cos(-5π/4) + i sin(-5π/4) ≈ -0.7+0.7i

ω¹⁴ = ω⁶ = e^-3πi/2 = cos(-3π/2) + i sin(-3π/2) = i

ω¹⁵ = ω⁷ = e^-7πi/4 = cos(-7π/4) + i sin(-7π/4) ≈ 0.7+0.7i

注意到當 n > 2 時，a_n = 0，於是：
 

A₀ = a₀xω^0x0 + a₁xω^1x0 + a₂xω^2x0 = 8xω⁰ + 7xω⁰ + 6xω⁰ = 8x1 + 7x1 + 6x1 = 21
 

A₁ = a₀xω^0x1 + a₁xω^1x1 + a₂xω^2x1 = 8xω⁰ + 7xω¹ + 6xω² ≈ 8x1 + 7x(0.7 - 0.7i) + 6x(-i) = 12.9-10.9i
 

A₂ = a₀xω^0x2 + a₁xω^1x2 + a₂xω^2x2 = 8xω⁰ + 7xω² + 6xω⁴ = 8x1 + 7x(-i) + 6x(-1) = 2-7i
 

A₃ = a₀xω^0x3 + a₁xω^1x3 + a₂xω^2x3 = 8xω⁰ + 7xω³ + 6xω⁶ ≈ 8x1 + 7x(-0.7 - 0.7i) + 6xi = 3.1+1.1i

A₄ = a₀xω^0x4 + a₁xω^1x4 + a₂xω^2x4 = 8xω⁰ + 7xω⁴ + 6xω⁸ = 8x1 + 7x(-1) + 6x1 = 7
 

A₅ = a₀xω^0x5 + a₁xω^1x5 + a₂xω^2x5 = 8xω⁰ + 7xω⁵ + 6xω¹⁰ ≈ 8x1 + 7x(-0.7 + 0.7i) + 6x(-i) = 3.1-1.1i

A₆ = a₀xω^0x6 + a₁xω^1x6 + a₂xω^2x6 = 8xω⁰ + 7xω⁶ + 6xω¹² = 8x1 + 7xi + 6x(-1) = 2+7i
 

A₇ = a₀xω^0x7 + a₁xω^1x7 + a₂xω^2x7 = 8xω⁰ + 7xω⁷ + 6xω¹⁴ ≈ 8x1 + 7x(0.7 + 0.7i) + 6xi = 12.9+10.9i
 

同樣，當 n > 2 時，b_n = 0，於是：
 

B₀ = b₀xω^0x0 + b₁xω^1x0 + b₂xω^2x0 = 2xω⁰ + 3xω⁰ + 4xω⁰ = 2x1 + 3x1 + 4x1 = 9
 

B₁ = b₀xω^0x1 + b₁xω^1x1 + b₂xω^2x1 = 2xω⁰ + 3xω¹ + 4xω² ≈ 2x1 + 3x(0.7 - 0.7i) + 4x(-i) = 4.1-6.1i
 

B₂ = b₀xω^0x2 + b₁xω^1x2 + b₂xω^2x2 = 2xω⁰ + 3xω² + 4xω⁴ = 2x1 + 3x(-i) + 4x(-1) = -2-3i
 

B₃ = b₀xω^0x3 + b₁xω^1x3 + b₂xω^2x3 = 2xω⁰ + 3xω³ + 4xω⁶ ≈ 2x1 + 3x(-0.7 - 0.7i) + 4xi = -0.1+1.9i

B₄ = b₀xω^0x4 + b₁xω^1x4 + b₂xω^2x4 = 2xω⁰ + 3xω⁴ + 4xω⁸ = 2x1 + 3x(-1) + 4x1 = 3
 

B₅ = b₀xω^0x5 + b₁xω^1x5 + b₂xω^2x5 = 2xω⁰ + 3xω⁵ + 4xω¹⁰ ≈ 2x1 + 3x(-0.7 + 0.7i) + 4x(-i) = -0.1-1.9i

B₆ = b₀xω^0x6 + b₁xω^1x6 + b₂xω^2x6 = 2xω⁰ + 3xω⁶ + 4xω¹² = 2x1 + 3xi + 4x(-1) = -2+3i
 

B₇ = b₀xω^0x7 + b₁xω^1x7 + b₂xω^2x7 = 2xω⁰ + 3xω⁷ + 4xω¹⁴ ≈ 2x1 + 3x(0.7 + 0.7i) + 4xi = 4.1+6.1i

 

這樣，向量 {a_i} 和向量 {b_j} 的離散傅裡葉變換 {A_i} 和 {B_j} 如下所示：

 

{ A₇, A₆, A₅, A₄, A₃, A₂, A₁, A₀ } = { 12.9+10.9i, 2+7i, 3.1-1.1i, 7, 3.1+1.1i, 2-7i, 12.9-10.9i, 21 }

{ B₇, B₆, B₅, B₄, B₃, B₂, B₁, B₀ } = { 4.1+6.1i, -2+3i, -0.1-1.9i, 3, -0.1+1.9i, -2-3i, 4.1-6.1i, 9 }

 

現在，將她們逐項相乘得到向量 {C_k}，即 { C₇, C₆, C₅, C₄, C₃, C₂, C₁, C₀ }

= { -13.6+123.4i, -25-8i, -2.4-5.8i, 21, -2.4+5.8i, -25+8i, -13.6-123.4i, 189 }

 

為了求出向量 {C_k} 的離散傅裡葉逆變換，我們令

ω = 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

 

為了方便計算，我們預先求出 ω 的各次方(注意 ω^k+8 = ω^k)，如下：

ω⁰ = e⁰ = 1

ω¹ = e^πi/4 = cos(π/4) + i sin(π/4) ≈ 0.7+0.7i

ω² = e^πi/2 = cos(π/2) + i sin(π/2) = i

ω³ = e^3πi/4 = cos(3π/4) + i sin(3π/4) ≈ -0.7+0.7i

ω⁴ = e^πi = cos(π) + i sin(π) = -1

ω⁵ = e^5πi/4 = cos(5π/4) + i sin(5π/4) ≈ -0.7-0.7i

ω⁶ = e^3πi/2 = cos(3π/2) + i sin(3π/2) = -i

ω⁷ = e^7πi/4 = cos(7π/4) + i sin(7π/4) ≈ 0.7-0.7i

 

於是：
 

c₀ = (1/N) x ( C₀xω^0x0 + C₁xω^1x0 + C₂xω^2x0 + C₃xω^3x0
                  + C₄xω^4x0 + C₅xω^5x0 + C₆xω^6x0 + C₇xω^7x0 )
    = (1/8) x ( 189xω⁰ + (-13.6-123.4i)xω⁰ + (-25+8i)xω⁰ + (-2.4+5.8i)xω⁰
                  + 21xω⁰ + (-2.4-5.8i)xω⁰ + (-25-8i)xω⁰ + (-13.6+123.4i)xω⁰ )
    = 0.125 x ( 189x1 + (-13.6-123.4i)x1 + (-25+8i)x1 + (-2.4+5.8i)x1
                  + 21x1 + (-2.4-5.8i)x1 + (-25-8i)x1 + (-13.6+123.4i)x1 )
    = 0.125 x 128 = 16
 

c₁ = (1/N) x ( 8xc₁ = C₀xω^0x1 + C₁xω^1x1 + C₂xω^2x1 + C₃xω^3x1
                  + C₄xω^4x1 + C₅xω^5x1 + C₆xω^6x1 + C₇xω^7x1 )
    = (1/8) x ( 189xω⁰ + ( -13.6-123.4i)xω¹ + (-25+8i)xω² + (-2.4+5.8i)xω³
                  + 21xω⁴ + (-2.4-5.8i)xω⁵ + (-25-8i)xω⁶ + (-13.6+123.4i)xω⁷ )
    ≈ 0.125 x ( 189x1 + (-13.6-123.4i)x(0.7+0.7i) + (-25+8i)x(i) + (-2.4+5.8i)x(-0.7+0.7i)
                  + 21x(-1) + (-2.4-5.8i)x(-0.7-0.7i) + (-25-8i)x(-i) + (-13.6+123.4i)x(0.7-0.7i) )
    = 0.125 x 300.96 = 37.62 ≈ 38

c₂ = (1/N) x ( C₀xω^0x2 + C₁xω^1x2 + C₂xω^2x2 + C₃xω^3x2
                  + C₄xω^4x2 + C₅xω^5x2 + C₆xω^6x2 + C₇xω^7x2 )
    = (1/8) x ( 189xω⁰ + (-13.6-123.4i)xω² + (-25+8i)xω⁴ + (-2.4+5.8i)xω⁶
                  + 21xω⁸ + (-2.4-5.8i)xω¹⁰ + (-25-8i)xω¹² + (-13.6+123.4i)xω¹⁴ )
    = 0.125 x ( 189x1 + (-13.6-123.4i)x(i) + (-25+8i)x(-1) + (-2.4+5.8i)x(-i)
                  + 21x1 + (-2.4-5.8i)x(i) + (-25-8i)x(-1) + (-13.6+123.4i)x(-i) )
    ≈ 0.125 x 518.4 = 64.8 ≈ 65
 

c₃ = (1/N) x ( C₀xω^0x3 + C₁xω^1x3 + C₂xω^2x3 + C₃xω^3x3
                  + C₄xω^4x3 + C₅xω^5x3 + C₆xω^6x3 + C₇xω^7x3 )
    = (1/8) x ( 189xω⁰ + (-13.6-123.4i)xω³ + (-25+8i)xω⁶ + (-2.4+5.8i)xω⁹
                  + 21xω¹² + (-2.4-5.8i)xω¹⁵ + (-25-8i)xω¹⁸ + (-13.6+123.4i)xω²¹ )
    ≈ 0.125 x ( 189x1 + (-13.6-123.4i)x(-0.7+0.7i) + (-25+8i)x(-i) + (-2.4+5.8i)x(0.7+0.7i)
                  + 21x(-1) + (-2.4-5.8i)x(0.7-0.7i) + (-25-8i)x(i) + (-13.6+123.4i)x(-0.7-0.7i) )
    = 0.125 x 364.32 = 45.54 ≈ 46

c₄ = (1/N) x ( C₀xω^0x4 + C₁xω^1x4 + C₂xω^2x4 + C₃xω^3x4
                  + C₄xω^4x4 + C₅xω^5x4 + C₆xω^6x4 + C₇xω^7x4 )
    = (1/8) x ( 189xω⁰ + (-13.6-123.4i)xω⁴ + (-25+8i)xω⁸ + (-2.4+5.8i)xω¹²
                  + 21xω¹⁶ + (-2.4-5.8i)xω²⁰ + (-25-8i)xω²⁴ + (-13.6+123.4i)xω²⁸ )
    = 0.125 x ( 189x1 + (-13.6-123.4i)x(-1) + (-25+8i)x1 + (-2.4+5.8i)x(-1)
                  + 21x1 + (-2.4-5.8i)x(-1) + (-25-8i)x1 + (-13.6+123.4i)x(-1) )
    = 0.125 x 192 = 24

c₅ = (1/N) x ( C₀xω^0x5 + C₁xω^1x5 + C₂xω^2x5 + C₃xω^3x5
                  + C₄xω^4x5 + C₅xω^5x5 + C₆xω^6x5 + C₇xω^7x5 )
    = (1/8) x ( 189xω⁰ + (-13.6-123.4i)xω⁵ + (-25+8i)xω¹⁰ + (-2.4+5.8i)xω¹⁵
                  + 21xω²⁰ + (-2.4-5.8i)xω²⁵ + (-25-8i)xω³⁰ + (-13.6+123.4i)xω³⁵ )
    ≈ 0.125 x ( 189x1 + (-13.6-123.4i)x(-0.7-0.7i) + (-25+8i)x(i) + (-2.4+5.8i)x(0.7-0.7i)
                  + 21x(-1) + (-2.4-5.8i)x(0.7+0.7i) + (-25-8i)x(-i) + (-13.6+123.4i)x(-0.7+0.7i) )
    = 0.125 x 3.04 = 0.38 ≈ 0
 

c₆ = (1/N) x ( C₀xω^0x6 + C₁xω^1x6 + C₂xω^2x6 + C₃xω^3x6
                  + C₄xω^4x6 + C₅xω^5x6 + C₆xω^6x6 + C₇xω^7x6 )
    = (1/8) x ( 189xω⁰ + (-13.6-123.4i)xω⁶ + (-25+8i)xω¹² + (-2.4+5.8i)xω¹⁸
                  + 21xω²⁴ + (-2.4-5.8i)xω³⁰ + (-25-8i)xω³⁶ + (-13.6+123.4i)xω⁴² )
    = 0.125 x ( 189x1 + (-13.6-123.4i)x(-i) + (-25+8i)x(-1) + (-2.4+5.8i)x(i)
                  + 21x1 + (-2.4-5.8i)x(-i) + (-25-8i)x(-1) + (-13.6+123.4i)x(i) )
    = 0.125 x 1.6 = 0.2 ≈ 0
 

 

c₇ = (1/N) x ( C₀xω^0x7 + C₁xω^1x7 + C₂xω^2x7 + C₃xω^3x7
                  + C₄xω^4x7 + C₅xω^5x7 + C₆xω^6x7 + C₇xω^7x7 )
    = (1/8) x ( 189xω⁰ + (-13.6-123.4i)xω⁷ + (-25+8i)xω¹⁴ + (-2.4+5.8i)xω²¹
                  + 21xω²⁸ + (-2.4-5.8i)xω³⁵ + (-25-8i)xω⁴² + (-13.6+123.4i)xω⁴⁹ )
    ≈ 0.125 x ( 189x1 + (-13.6-123.4i)x(0.7-0.7i) + (-25+8i)x(-i) + (-2.4+5.8i)x(-0.7-0.7i)
                  + 21x(-1) + (-2.4-5.8i)x(-0.7+0.7i) + (-25-8i)x(i) + (-13.6+123.4i)x(0.7+0.7i) )
    = 0.125 x 3.68 = 0.46 ≈ 0
 

這樣，我們就使用離散傅裡葉變換和逆變換計算出了向量 {a_i} 和向量 {b_j} 的卷積向量 {c_k}，如下所示：
 

{ c₇, c₆, c₅, c₄, c₃, c₂, c₁, c₀ } = { 0, 0, 0, 0, 24, 46, 65, 38, 16 }
 

這和我們在前面直接使用向量 {a_i} 和向量 {b_j} 來計算卷積的結果是一樣的。

但是，這個算法的時間復雜度還是 O(N²)。我們繞了這麼一大圈，不是白費勁了嗎？

 

現在就到了關鍵時刻，關鍵在於：直接進行離散傅裡葉變換的計算復雜度是 O(N²)。快速傅裡葉變換可以計算出與直接計算相同的結果，但只需要 O(N logN) 的計算復雜度。 N logN 和 N² 之間的差別是巨大的。例如，當 N = 10⁶ 時，在一個每秒運算百萬次的計算機上，粗略地說，它們之間就是占用 30 秒 CPU 時間和兩星期 CPU 時間的差別。
 

快速傅裡葉變換的要點如下：一個界長為 N 的離散傅裡葉變換可以重新寫成兩個界長各為 N/2 的離散傅裡葉變換之和。其中一個變換由原來 N 個點中的偶數點構成，另一個變換由奇數點構成。這個過程可以遞歸地進行下去，直到我們將全部數據細分為界長為 1 的變換。什麼是界長為 1 的傅裡葉變換呢？它正是把一個輸入值復制成它的一個輸出值的恆等運算。要實現以上算法，最容易的情況是原始的 N 為 2 的整冪次項，如果數據集的界長不是 2 的冪次時，則可添上一些零值，直到 2 的下一冪次。在這個算法中，每遞歸一次需 N 階運算，共需要 log N 次遞歸，所以快速傅裡葉變換算法的時間復雜度是 O(N logN)。

 

 
 

由於快速傅裡葉變換是采用了浮點運算，因此我們需要足夠的精度，以使在出現捨入誤差時，結果中每個組成部分的准確整數值仍是可辨認的。長度為 N 的 B 進制數可產生大到 B²N 階的卷積分量。我們知道，雙精度浮點數的尾數是 53 個二進位，所以：

2 x log₂B + log₂N + 幾個 x log₂log₂N < 53

上式中左邊最後一項是為了快速傅裡葉變換的捨入誤差。

所以，為了能夠計算盡量大的整數，一般 B 不會取得太大。在計算機程序中經常使用 256 進制進行運算。但是如果經常需要將計算結果和十進制互相轉換，則往往使用 100 進制進行運算。

 

關於快速傅裡葉變換以及卷積定理的更深入的知識，請參閱文末的參考文獻。這一篇隨筆主要是講述相關的原理，在下一篇隨筆中，我將給出一個使用快速傅裡葉變換進行任意精度的算術運算的 C# 程序。

 

順便說一句，我在准備正文的例題的時候，是使用 google 計算器來進行復雜的復數運算的。發現她非常好用。以計算 c₂ 為例， 只要將要計算的表達式復制到 goole 搜索欄，然後按回車，就能得到計算結果：
 

(189 x 1) + (((-13.6) - (123.4 * i)) x i) + (((-25) + (8 * i)) x (-1)) + (((-2.4) + (5.8 * i)) x (-i)) + (21 x 1) + (((-2.4) - (5.8 * i)) x i) + (((-25) - (8 * i)) x (-1)) + (((-13.6) + (123.4 * i)) x (-i)) = 518.4 - 1.77635684 × 10^-15 i
Google 計算器詳情

找不到和您的查詢 "189x1 + (-13.6-123.4i)x(i) + (-25+8i)x(-1) + (-2.4+5.8i)x(-i) + 21x1 + (-2.4-5.8i)x(i) + (-25-8i)x(-1) + (-13.6+123.4i)x(-i)" 相符的網頁。