【NA】误差与有效数字

误差与有效数字"/>

【NA】误差与有效数字

文章目录

• 误差简述.
• • 绝对误差.
  • 相对误差.
  • 有效数字.
• 减少误差的原则.
• • 减少误差的实例.

误差简述.

绝对误差.

• 【绝对误差】设 x ∗ x^* x∗ 为准确值， x x x 是 x ∗ x^* x∗ 的近似值，则称 e = x ∗ − x e=x^*-x e=x∗−x 为近似值 x x x 的绝对误差，简称误差。
• 由于准确值 x ∗ x^* x∗ 未知，实际工作中一般根据相关领域的知识、经验以及测量工具的精度，事先估计绝对误差的绝对值，即 ∣ x ∗ − x ∣ |x^*-x| ∣x∗−x∣ 不超过某个正数 ϵ \epsilon ϵ，所以有 ∣ e ∣ = ∣ x ∗ − x ∣ ≤ ϵ |e|=|x^*-x|≤\epsilon ∣e∣=∣x∗−x∣≤ϵ，该正数 ϵ \epsilon ϵ 称为近似值 x x x 的绝对误差限，简称为误差限或精度。
• 依据上面关于绝对误差限的定义，能够得出下式： x − ϵ ≤ x ∗ ≤ x + ϵ x-\epsilon≤x^*≤x+\epsilon x−ϵ≤x∗≤x+ϵ，因此有时可以将准确值 x ∗ x^* x∗ 写为 x ± ϵ x±\epsilon x±ϵ.
• 但近似值 x x x 的优劣通过绝对误差并不能完全反应出来，一个很直观的例子就是如果你年收入1000000元，财务少算了100元；和你年收入1000元少算了100元，基数大小的不同导致相同数量的计算失误的结果不同。

相对误差.

• 【相对误差】所以除了考虑绝对误差的大小，还需要考虑准确值本身的大小，定义近似值 x x x 的相对误差为 e r = e x ∗ = x ∗ − x x ∗ e_r=\frac{e}{x^*}=\frac{x^*-x}{x^*} er​=x∗e​=x∗x∗−x​，实际计算中由于 x ∗ x^* x∗ 未知，所以将定义改为 e r = e x = x ∗ − x x e_r=\frac{e}{x}=\frac{x^*-x}{x} er​=xe​=xx∗−x​，称为近似值 x x x 的相对误差。
• 上述例子中，前者的相对误差为 1 0 − 4 10^{-4} 10−4，而后者的相对误差为 0.1 0.1 0.1，一般来说相对误差越小，表明 x x x 对 x ∗ x^* x∗ 的近似程度越好。
• 仿照绝对误差限 ∣ e ∣ |e| ∣e∣ 的定义，可以定义相对误差限 ∣ e r ∣ = ∣ x ∗ − x x ∣ |e_r|=|\frac{x^*-x}{x}| ∣er​∣=∣xx∗−x​∣.

有效数字.

• 【定义】设近似值 x = ± 0. a 1 a 2 . . . a n × 1 0 m x=±0.a_1a_2...a_n×10^m x=±0.a1​a2​...an​×10m，其中 a i a_i ai​ 都是 0 0 0 到 9 9 9 之间的自然数，并且 a 1 ≠ 0 a_1≠0 a1​​=0， m m m 为整数，如果 ∣ x ∗ − x ∣ ≤ 1 2 × 1 0 m − n |x^*-x|≤\frac{1}{2}×10^{m-n} ∣x∗−x∣≤21​×10m−n，则称近似值 x x x 具有 n n n 位有效数字，也称 x x x 是有 n n n 位有效数字的近似值。
• 【定理一】设近似值 x = ± 0. a 1 a 2 . . . a n × 1 0 m x=±0.a_1a_2...a_n×10^m x=±0.a1​a2​...an​×10m 有 n n n 位有效数字，则其相对误差限 ∣ e r ∣ = 1 2 a 1 × 1 0 − n + 1 |e_r|=\frac{1}{2a_1}×10^{-n+1} ∣er​∣=2a1​1​×10−n+1.
• 【定理二】设近似值 x = ± 0. a 1 a 2 . . . a n × 1 0 m x=±0.a_1a_2...a_n×10^m x=±0.a1​a2​...an​×10m 的相对误差限 ϵ r = 1 2 ( a 1 + 1 ) × 1 0 − n + 1 \epsilon_r=\frac{1}{2(a_1+1)}×10^{-n+1} ϵr​=2(a1​+1)1​×10−n+1，则它至少具有 n n n 位有效数字。
• 从定理一和定理二的表述中不难看出，有效数字的位数越多，其相对误差限就越小。

• 【例题】设 x ∗ = 3.200169 x^*=3.200169 x∗=3.200169，则它的近似值 x 1 = 3.2001 x_1=3.2001 x1​=3.2001， x 2 = 3.2002 x_2=3.2002 x2​=3.2002 分别具有几位有效数字？
• 【解答】根据定义， x 1 = 0.32001 × 1 0 1 x_1=0.32001×10^1 x1​=0.32001×101，所以 m 1 = 1 m_1=1 m1​=1， ∣ x ∗ − x 1 ∣ = 0.000069 = 0.069 × 1 0 − 3 ≤ 0.5 × 1 0 − 3 |x^*-x_1|=0.000069=0.069×10^{-3}≤0.5×10^{-3} ∣x∗−x1​∣=0.000069=0.069×10−3≤0.5×10−3，所以 m 1 − n 1 = − 3 m_1-n_1=-3 m1​−n1​=−3，即 n 1 = 4 n_1=4 n1​=4，所以近似值 x 1 x_1 x1​ 的有效位数为4，最后一位的数字1不是有效数字；同理 x 2 = 0.32002 × 1 0 1 x_2=0.32002×10^1 x2​=0.32002×101， m 2 = 1 m_2=1 m2​=1， ∣ x ∗ − x 2 ∣ = 0.000031 = 0.31 × 1 0 − 4 ≤ 0.5 × 1 0 − 4 |x^*-x_2|=0.000031=0.31×10^{-4}≤0.5×10^{-4} ∣x∗−x2​∣=0.000031=0.31×10−4≤0.5×10−4，所以 m 2 − n 2 = − 4 m_2-n_2=-4 m2​−n2​=−4， n 2 = 5 n_2=5 n2​=5，所以 x 2 x_2 x2​ 的有效位数为5位。
• 【定理一证明】近似值 x x x 有 n n n 位有效数字，所以由定义可知 ∣ x ∗ − x ∣ ≤ 1 2 × 1 0 m − n |x^*-x|≤\frac{1}{2}×10^{m-n} ∣x∗−x∣≤21​×10m−n，因为 x = ± 0. a 1 a 2 . . . a n × 1 0 m = ± a 1 . a 2 . . . a n × 1 0 m − 1 x=±0.a_1a_2...a_n×10^m=±a_1.a_2...a_n×10^{m-1} x=±0.a1​a2​...an​×10m=±a1​.a2​...an​×10m−1，易得 ∣ x ∣ ≥ a 1 × 1 0 m − 1 |x|≥a_1×10^{m-1} ∣x∣≥a1​×10m−1，所以相对误差限 ∣ e r ∣ = ∣ x ∗ − x ∣ ∣ x ∣ ≤ 1 2 × 1 0 m − n a 1 × 1 0 m − 1 = 1 2 a 1 × 1 0 − n + 1 = ϵ r |e_r|=\frac{|x^*-x|}{|x|}≤\frac{\frac{1}{2}×10^{m-n}}{a_1×10^{m-1}}=\frac{1}{2a_1}×10^{-n+1}=\epsilon_r ∣er​∣=∣x∣∣x∗−x∣​≤a1​×10m−121​×10m−n​=2a1​1​×10−n+1=ϵr​.
• 【定理二证明】因为 x = ± 0. a 1 a 2 . . . a n × 1 0 m = ± a 1 . a 2 . . . a n × 1 0 m − 1 x=±0.a_1a_2...a_n×10^m=±a_1.a_2...a_n×10^{m-1} x=±0.a1​a2​...an​×10m=±a1​.a2​...an​×10m−1，易得 ∣ x ∣ ≤ ( a 1 + 1 ) × 1 0 m − 1 |x|≤(a_1+1)×10^{m-1} ∣x∣≤(a1​+1)×10m−1，所以 ∣ x ∗ − x ∣ = ∣ x ∗ − x ∣ ∣ x ∣ ⋅ ∣ x ∣ ≤ 1 2 ( a 1 + 1 ) × 1 0 − n + 1 ⋅ ( a 1 + 1 ) × 1 0 m − 1 = 0.5 × 1 0 m − n |x^*-x|=\frac{|x^*-x|}{|x|}·|x|≤\frac{1}{2(a_1+1)}×10^{-n+1}·(a_1+1)×10^{m-1}=0.5×10^{m-n} ∣x∗−x∣=∣x∣∣x∗−x∣​⋅∣x∣≤2(a1​+1)1​×10−n+1⋅(a1​+1)×10m−1=0.5×10m−n，根据定义可知近似值 x x x 至少具有 n n n 位有效数字。

减少误差的原则.

• 避免两个相近的数相减；
• 避免重要的小数被大数吞没；
• 避免除数的绝对值远小于被除数的绝对值；
• 注意算法的稳定性；

如果输入数据有误差，而算法在计算过程中误差不增张，那么称此算法是数值稳定的，否则为数值不稳定的。

• 简化计算步骤，减少计算次数。

减少误差的实例.

• ①当 x = 10003 x=10003 x=10003 时，计算 x + 1 − x \sqrt{x+1}-\sqrt{x} x+1 ​−x ​ 的近似值。(避免两个相近的数相减)
• 【分析】采用6位十进制浮点运算，如果直接相减，由于 x > > 1 x>>1 x>>1，所以减数与被减数相差极小， x + 1 − x = 100.020 − 100.015 = 0.005 \sqrt{x+1}-\sqrt{x}=100.020-100.015=0.005 x+1 ​−x ​=100.020−100.015=0.005，该式的精确值为 0.004999125231... 0.004999125231... 0.004999125231...，按照有效数字位数的定义可知，上述直接相减的结果只有 1 1 1 位有效数字，损失了 5 5 5 位有效数字。如果对算式进行等价变形，即 x + 1 − x = 1 x + 1 + x = 1 100.020 + 100.015 = 0.00499913 \sqrt{x+1}-\sqrt{x}=\frac{1}{\sqrt{x+1}+\sqrt{x}}=\frac{1}{100.020+100.015}=0.00499913 x+1 ​−x ​=x+1 ​+x ​1​=100.020+100.0151​=0.00499913，有 6 6 6 位有效数字。

• 求二次方程 x 2 − ( 1 0 9 + 1 ) x + 1 0 9 = 0 x^2-(10^9+1)x+10^9=0 x2−(109+1)x+109=0 的根。(避免重要的小数被大数吞没)
• 【分析】使用因式分解法可知该方程两个解为 1 1 1 和 1 0 9 10^9 109，但如果编制程序，使用二次方程求根公式 x = − b ± b 2 − 4 a c 2 a x=\frac{-b±\sqrt{b^2-4ac}}{2a} x=2a−b±b2−4ac ​​，如果在一台只能将数表示到小数点后 8 8 8 位的计算机上运算，首先进行对阶操作， − b = 1 0 9 + 1 = 0.1000000 × 1 0 10 + 0.0000000001 × 1 0 10 -b=10^9+1=0.1000000×10^{10}+0.0000000001×10^{10} −b=109+1=0.1000000×1010+0.0000000001×1010，由于只能表示到小数点后 8 8 8 位，加号之后的数值就被舍去，所以 − b = 1 0 9 -b=10^9 −b=109；类似地 4 a c = 0 4ac=0 4ac=0，所以最终得到的两个根为 1 0 9 10^9 109 和 0 0 0，产生这两个解的原因是大数在运算中吞没了小数。产生解 0 0 0 的计算式为 0 = − b − b 2 − 4 a c 2 a 0=\frac{-b-\sqrt{b^2-4ac}}{2a} 0=2a−b−b2−4ac ​​，如果对其进行等价变形，得到 x 2 = 2 c − b + b 2 − 4 a c = 2 × 1 0 9 1 0 9 + 1 0 9 = 1 x_2=\frac{2c}{-b+\sqrt{b^2-4ac}}=\frac{2×10^9}{10^9+10^9}=1 x2​=−b+b2−4ac ​2c​=109+1092×109​=1，在计算机精度不足的情况下也能够得到正确结果。
• 类似的情况，若 x = 1 0 12 x=10^{12} x=1012， y = 3 y=3 y=3， z = − 1 0 12 z=-10^{12} z=−1012，求 x + y + z x+y+z x+y+z. 如果按照从左向右的加法次序来编制程序，同样在一台只能将数表示到小数点后 8 8 8 位的计算机上运算， x + y x+y x+y时的对阶操作就会导致大数 x x x 将小数 y y y 吞没，从而最后的结果为 0 0 0；而如果改变次序，按照 ( x + z ) + y (x+z)+y (x+z)+y 的次序进行运算，就能够得到正确的结果 3 3 3.
• 【Digression】出自Anany Levitin《Introduction to The Design and Analysis of Algorithms》，其中公式(11.14)即 x = − b ± b 2 − 4 a c 2 a x=\frac{-b±\sqrt{b^2-4ac}}{2a} x=2a−b±b2−4ac ​​.
  
  
  

• 在 4 4 4 位浮点十进制数下，用消去法解线性方程组 { 0.00003 x 1 − 3 x 2 = 0.6 x 1 + 2 x 2 = 1 \begin{cases} 0.00003x_1-3x_2=0.6\\ x_1+2x_2=1\end{cases} {0.00003x1​−3x2​=0.6x1​+2x2​=1​(避免除数的绝对值远小于被除数的绝对值)
• 【分析】在计算机的 4 4 4 位浮点十进制数表示中，上述方程组写为： { 0.3000 × 1 0 − 4 x 1 − 0.3000 × 1 0 1 x 2 = 0.6000 × 1 0 0 − − − ① 0.1000 × 1 0 1 x 1 + 0.2000 × 1 0 1 x 2 = 0.1000 × 1 0 1 − − − ② \begin{cases} 0.3000×10^{-4}x_1-0.3000×10^1x_2=0.6000×10^0---①\\ 0.1000×10^1x_1+0.2000×10^1x_2=0.1000×10^1---②\end{cases} {0.3000×10−4x1​−0.3000×101x2​=0.6000×100−−−①0.1000×101x1​+0.2000×101x2​=0.1000×101−−−②​
• 第一种消去法思路将 ① 式除以 0.00003 0.00003 0.00003 后减去 ② 式，从而消去 x 1 x_1 x1​，得到 − 0.1000 × 1 0 6 x 2 = 0.2000 × 1 0 5 -0.1000×10^6x_2=0.2000×10^5 −0.1000×106x2​=0.2000×105，解得 x 2 = − 0.2 x_2=-0.2 x2​=−0.2，回代入 ① 式得 x 1 = 0 x_1=0 x1​=0，显然这组解偏差极大。原方程组的准确解为 x 1 = 1.399... x_1=1.399... x1​=1.399...， x 2 = − 0.199... x_2=-0.199... x2​=−0.199...，在 x 1 x_1 x1​ 处发生了严重失真。
• 第二种消去思路保留 ② 式，消去 ① 式中的 x 1 x_1 x1​，即将 ② 式乘以 0.00003 0.00003 0.00003 后减去 ① 式，得到 − 0.3000 × 1 0 1 x 2 = 0.6000 × 1 0 1 -0.3000×10^1x_2=0.6000×10^1 −0.3000×101x2​=0.6000×101，所以解得 x 2 = − 0.2 x_2=-0.2 x2​=−0.2，回代入 ② 式得 x 1 = 1.4 x_1=1.4 x1​=1.4，相较于准确解来说是一个很好的近似解。

• 当 n = 0 , 1 , 2 , . . . , 11 n=0,1,2,...,11 n=0,1,2,...,11 时，计算积分 I n = ∫ 0 1 x n x + 9 d x I_n=\int_{0}^{1}\frac{x^n}{x+9}dx In​=∫01​x+9xn​dx 的近似值。(注意算法的稳定性)
• 【分析】 I n + 9 I n − 1 = ∫ 0 1 ( x n x + 9 + 9 x n − 1 x + 9 ) d x = ∫ 0 1 x n − 1 d x = x n n ∣ 0 1 = 1 n I_n+9I_{n-1}=\int_{0}^{1}(\frac{x^n}{x+9}+\frac{9x^{n-1}}{x+9})dx=\int_{0}^{1}x^{n-1}dx=\frac{x^n}{n}|_0^1=\frac{1}{n} In​+9In−1​=∫01​(x+9xn​+x+99xn−1​)dx=∫01​xn−1dx=nxn​∣01​=n1​，所以有 I n = 1 n − 9 I n − 1 I_n=\frac{1}{n}-9I_{n-1} In​=n1​−9In−1​. 当 n = 0 n=0 n=0 时， I 0 = ∫ 0 1 1 x + 9 d x = l n ( x + 9 ) ∣ 0 1 ≈ 0.105361 = I 0 ‾ I_0=\int^1_0\frac{1}{x+9}dx=ln(x+9)|_0^1≈0.105361=\overline{I_0} I0​=∫01​x+91​dx=ln(x+9)∣01​≈0.105361=I0​​，从而有 { I 0 ‾ = 0.105361 I n ‾ = 1 n − 9 I n − 1 ‾ ， n = 1 , 2 , . . . , 11 \begin{cases} \overline{I_0}=0.105361\\ \overline{I_n}=\frac{1}{n}-9\overline{I_{n-1}}，n=1,2,...,11\end{cases} {I0​​=0.105361In​​=n1​−9In−1​​，n=1,2,...,11​
• 根据上式计算出 I 1 ‾ = 0.051751 \overline{I_1}=0.051751 I1​​=0.051751，直到 I 5 ‾ = − 0.011689 \overline{I_5}=-0.011689 I5​​=−0.011689，问题出现了，上述积分在 [ 0 , 1 ] [0,1] [0,1] 内必然是正数，显然误差已经过大了。我们记初始误差为 ϵ 0 = I 0 − I 0 ‾ \epsilon_0=I_0-\overline{I_0} ϵ0​=I0​−I0​​，根据迭代式，得到 ϵ 1 = I 1 − I 1 ‾ = ( 1 − 9 I 0 ) − ( 1 − 9 I 0 ‾ ) = − 9 ϵ 0 \epsilon_1=I_1-\overline{I_1}=(1-9I_0)-(1-9\overline{I_0})=-9\epsilon_0 ϵ1​=I1​−I1​​=(1−9I0​)−(1−9I0​​)=−9ϵ0​，所以最终 ∣ ϵ n ∣ = 9 n ∣ ϵ 0 ∣ |\epsilon_n|=9^n|\epsilon_0| ∣ϵn​∣=9n∣ϵ0​∣，其误差绝对值会随着计算过程的进行不断增大，显然该算法不是数值稳定的。
• 考虑另一种算法，由 I n = 1 n − 9 I n − 1 I_n=\frac{1}{n}-9I_{n-1} In​=n1​−9In−1​ 可得 I n − 1 = 1 9 ( 1 n − I n ) I_{n-1}=\frac{1}{9}(\frac{1}{n}-I_n) In−1​=91​(n1​−In​)，直接估计出 I 12 I_{12} I12​ 的近似值，而后将 n n n 递减。因为 I n = ∫ 0 1 x n x + 9 d x I_n=\int_{0}^{1}\frac{x^n}{x+9}dx In​=∫01​x+9xn​dx，所以 ∫ 0 1 x n 10 d x ≤ I n ≤ ∫ 0 1 x n 9 d x \int^1_0\frac{x^n}{10}dx≤I_n≤\int^1_0\frac{x^n}{9}dx ∫01​10xn​dx≤In​≤∫01​9xn​dx，所以 1 10 ( n + 1 ) ≤ I n ≤ 1 9 ( n + 1 ) \frac{1}{10(n+1)}≤I_n≤\frac{1}{9(n+1)} 10(n+1)1​≤In​≤9(n+1)1​，令 n = 12 n=12 n=12 可以得到 I 12 ≈ 1 2 ( 1 10 × 13 + 1 9 × 13 ) ≈ 0.008120 = I 12 ‾ I_{12}≈\frac{1}{2}(\frac{1}{10×13}+\frac{1}{9×13})≈0.008120=\overline{I_{12}} I12​≈21​(10×131​+9×131​)≈0.008120=I12​​，可得 { I 12 ‾ = 0.008120 I n − 1 ‾ = 1 9 ( 1 n − I n ‾ ) ， n = 12 , 11 , . . . , 1 \begin{cases} \overline{I_{12}}=0.008120\\ \overline{I_{n-1}}=\frac{1}{9}(\frac{1}{n}-\overline{I_{n}})，n=12,11,...,1\end{cases} {I12​​=0.008120In−1​​=91​(n1​−In​​)，n=12,11,...,1​
• 采用相同的误差分析方法，记初始误差 ϵ 12 = I 12 − I 12 ‾ \epsilon_{12}=I_{12}-\overline{I_{12}} ϵ12​=I12​−I12​​，则 ∣ ϵ n ∣ = ∣ ϵ 12 ∣ 9 n |\epsilon_{n}|=\frac{|\epsilon_{12}|}{9^n} ∣ϵn​∣=9n∣ϵ12​∣​，在计算过程中，其误差绝对值是不断减小的，因而该算法是数值稳定的。

• 离散Fourier变换与快速Fourier变换(简化计算步骤)
• 【离散Fourier变换】已知一离散信号 ( f 0 , f 1 , . . . , f N − 1 ) T (f_0,f_1,...,f_{N-1})^T (f0​,f1​,...,fN−1​)T，则其对应的Fourier系数为 a j = 1 N ∑ k = 0 N − 1 f k ⋅ e − i j k 2 π N , j = 0 , 1 , 2 , . . . , N − 1 a_j=\frac{1}{N}\sum^{N-1}_{k=0}f_k·e^{-ijk\frac{2\pi}{N}},j=0,1,2,...,N-1 aj​=N1​k=0∑N−1​fk​⋅e−ijkN2π​,j=0,1,2,...,N−1
• 相应地，如果已知Fourier系数 ( a 0 , a 1 , . . . , a N − 1 ) T (a_0,a_1,...,a_{N-1})^T (a0​,a1​,...,aN−1​)T，那么可以反推出 f j = ∑ k = 0 N − 1 a k ⋅ e i j k 2 π N , j = 0 , 1 , 2 , . . . , N − 1 f_j=\sum^{N-1}_{k=0}a_k·e^{ijk\frac{2\pi}{N}},j=0,1,2,...,N-1 fj​=k=0∑N−1​ak​⋅eijkN2π​,j=0,1,2,...,N−1
• 上述两式分别是离散Fourier变换和离散Fourier 逆变换，观察后发现主要是计算自然指数部分，若令 ω = e − i 2 π N \omega=e^{-i\frac{2\pi}{N}} ω=e−iN2π​(着眼于正向变换)，即我们想得到 c j = ∑ k = 0 N − 1 x k ω k j , j = 0 , 1 , . . . , N − 1 c_j=\sum^{N-1}_{k=0}x_k\omega^{kj},j=0,1,...,N-1 cj​=k=0∑N−1​xk​ωkj,j=0,1,...,N−1
• 上式 c j c_j cj​ 计算需要进行 N N N 次复数乘法，完成整个变换共需要 N 2 N^2 N2 次复数乘法，一般取 N = 2 p N=2^p N=2p.
• 【快速Fourier变换】FFT的基本思想是利用 a b + a c = a ( b + c ) ab+ac=a(b+c) ab+ac=a(b+c) 减少乘法运算的次数，将 N 2 N^2 N2 量级的乘法次数降低到 N ⋅ l o g N N·logN N⋅logN. 在上面介绍的计算 c j c_j cj​ 式中，将相同的 ω k j \omega^{kj} ωkj 提取，对应的 x k x_k xk​ 相加，就能够减少乘法运算的次数。在 ω k j = e − i k j 2 π N \omega^{kj}=e^{-ikj\frac{2\pi}{N}} ωkj=e−ikjN2π​ 中，设 k j = q N + r kj=qN+r kj=qN+r，则有 ω k j = e − i ( q N + r ) 2 π N = e − i q 2 π ⋅ e − i 2 π r N = ω r \omega^{kj}=e^{-i(qN+r)\frac{2\pi}{N}}=e^{-iq2\pi}·e^{^{-i\frac{2\pi r}{N}}}=\omega^r ωkj=e−i(qN+r)N2π​=e−iq2π⋅e−iN2πr​=ωr
• 以 N = 2 3 N=2^3 N=23 为例， c j = ∑ k = 0 N − 1 x k ω k j , j = 0 , 1 , . . . , 7 c_j=\sum^{N-1}_{k=0}x_k\omega^{kj},j=0,1,...,7 cj​=∑k=0N−1​xk​ωkj,j=0,1,...,7，将 k , j k,j k,j 以二进制表示，即 k = k 0 2 0 + k 1 2 1 + k 2 2 2 = k 2 k 1 k 0 k=k_02^0+k_12^1+k_22^2=k_2k_1k_0 k=k0​20+k1​21+k2​22=k2​k1​k0​，同理 j = j 2 j 1 j 0 j=j_2j_1j_0 j=j2​j1​j0​，记 c j = c ( j 2 , j 1 , j 0 ) , x k = x ( k 2 , k 1 , k 0 ) c_j=c(j_2,j_1,j_0),x_k=x(k_2,k_1,k_0) cj​=c(j2​,j1​,j0​),xk​=x(k2​,k1​,k0​)，则有 c ( j 2 , j 1 , j 0 ) = ∑ k 0 = 0 1 ∑ k 1 = 0 1 ∑ k 2 = 0 1 x ( k 2 , k 1 , k 0 ) ω ( k 0 2 0 + k 1 2 1 + k 2 2 2 ) ( j 0 2 0 + j 1 2 1 + j 2 2 2 ) = ∑ k 0 = 0 1 ∑ k 1 = 0 1 ∑ k 2 = 0 1 x ( k 2 , k 1 , k 0 ) ω j 0 2 0 ( k 0 2 0 + k 1 2 1 + k 2 2 2 ) ω j 1 2 1 ( k 0 2 0 + k 1 2 1 + k 2 2 2 ) ω j 2 2 2 ( k 0 2 0 + k 1 2 1 + k 2 2 2 ) = ∑ k 0 = 0 1 ∑ k 1 = 0 1 ∑ k 2 = 0 1 x ( k 2 , k 1 , k 0 ) ω j 0 ( k 0 2 0 + k 1 2 1 + k 2 2 2 ) ω j 1 ( k 0 2 1 + k 1 2 2 ) ω j 2 ( k 0 2 2 ) = ∑ k 0 = 0 1 [ ∑ k 1 = 0 1 ( ∑ k 2 = 0 1 x ( k 2 , k 1 , k 0 ) ω j 0 ( k 0 2 0 + k 1 2 1 + k 2 2 2 ) ) ω j 1 ( k 0 2 1 + k 1 2 2 ) ] ω j 2 ( k 0 2 2 ) = ∑ k 0 = 0 1 [ ∑ k 1 = 0 1 ( ∑ k 2 = 0 1 x ( k 2 , k 1 , k 0 ) ω j 0 ( k 2 , k 1 , k 0 ) ω j 1 ( k 1 , k 0 , 0 ) ] ω j 2 ( k 2 , 0 , 0 ) c(j_2,j_1,j_0)=\sum^{1}_{k_0=0}\sum^1_{k_1=0}\sum^1_{k_2=0}x(k_2,k_1,k_0)\omega^{(k_02^0+k_12^1+k_22^2)(j_02^0+j_12^1+j_22^2)}\\=\sum^{1}_{k_0=0}\sum^1_{k_1=0}\sum^1_{k_2=0}x(k_2,k_1,k_0)\omega^{j_02^0(k_02^0+k_12^1+k_22^2)}\omega^{j_12^1(k_02^0+k_12^1+k_22^2)}\omega^{j_22^2(k_02^0+k_12^1+k_22^2)}\\=\sum^{1}_{k_0=0}\sum^1_{k_1=0}\sum^1_{k_2=0}x(k_2,k_1,k_0)\omega^{j_0(k_02^0+k_12^1+k_22^2)}\omega^{j_1(k_02^1+k_12^2)}\omega^{j_2(k_02^2)}\\=\sum^{1}_{k_0=0}[\sum^1_{k_1=0}(\sum^1_{k_2=0}x(k_2,k_1,k_0)\omega^{j_0(k_02^0+k_12^1+k_22^2)})\omega^{j_1(k_02^1+k_12^2)}]\omega^{j_2(k_02^2)}\\=\sum^{1}_{k_0=0}[\sum^1_{k_1=0}(\sum^1_{k_2=0}x(k_2,k_1,k_0)\omega^{j_0(k_2,k_1,k_0})\omega^{j_1(k_1,k_0,0)}]\omega^{j_2(k_2,0,0)} c(j2​,j1​,j0​)=k0​=0∑1​k1​=0∑1​k2​=0∑1​x(k2​,k1​,k0​)ω(k0​20+k1​21+k2​22)(j0​20+j1​21+j2​22)=k0​=0∑1​k1​=0∑1​k2​=0∑1​x(k2​,k1​,k0​)ωj0​20(k0​20+k1​21+k2​22)ωj1​21(k0​20+k1​21+k2​22)ωj2​22(k0​20+k1​21+k2​22)=k0​=0∑1​k1​=0∑1​k2​=0∑1​x(k2​,k1​,k0​)ωj0​(k0​20+k1​21+k2​22)ωj1​(k0​21+k1​22)ωj2​(k0​22)=k0​=0∑1​[k1​=0∑1​(k2​=0∑1​x(k2​,k1​,k0​)ωj0​(k0​20+k1​21+k2​22))ωj1​(k0​21+k1​22)]ωj2​(k0​22)=k0​=0∑1​[k1​=0∑1​(k2​=0∑1​x(k2​,k1​,k0​)ωj0​(k2​,k1​,k0​)ωj1​(k1​,k0​,0)]ωj2​(k2​,0,0)
• 将FFT表示如下： A 0 ( k 2 , k 1 , k 0 ) = x ( k 2 , k 1 , k 0 ) A 1 ( k 1 , k 0 , j 0 ) = ∑ k 2 = 0 1 A 0 ( k 2 , k 1 , k 0 ) ω j 0 ( k 2 , k 1 , k 0 ) A 2 ( k 0 , j 1 , j 0 ) = ∑ k 1 = 0 1 A 1 ( k 1 , k 0 , j 0 ) ω j 1 ( k 1 , k 0 , 0 ) A 3 ( j 2 , j 1 , j 0 ) = ∑ k 0 = 0 1 A 2 ( k 0 , j 1 , j 0 ) ω j 2 ( k 0 , 0 , 0 ) c ( j 2 , j 1 , j 0 ) = A 3 ( j 2 , j 1 , j 0 ) A_0(k_2,k_1,k_0)=x(k_2,k_1,k_0)\\A_1(k_1,k_0,j_0)=\sum^1_{k_2=0}A_0(k_2,k_1,k_0)\omega^{j_0(k_2,k_1,k_0)}\\A_2(k_0,j_1,j_0)=\sum^1_{k_1=0}A_1(k_1,k_0,j_0)\omega^{j_1(k_1,k_0,0)}\\A_3(j_2,j_1,j_0)=\sum^1_{k_0=0}A_2(k_0,j_1,j_0)\omega^{j_2(k_0,0,0)}\\c(j_2,j_1,j_0)=A_3(j_2,j_1,j_0) A0​(k2​,k1​,k0​)=x(k2​,k1​,k0​)A1​(k1​,k0​,j0​)=k2​=0∑1​A0​(k2​,k1​,k0​)ωj0​(k2​,k1​,k0​)A2​(k0​,j1​,j0​)=k1​=0∑1​A1​(k1​,k0​,j0​)ωj1​(k1​,k0​,0)A3​(j2​,j1​,j0​)=k0​=0∑1​A2​(k0​,j1​,j0​)ωj2​(k0​,0,0)c(j2​,j1​,j0​)=A3​(j2​,j1​,j0​)
• 当 N = 2 p N=2^p N=2p 时，单个系数的FFT分为 p p p 步，每步进行2次复数乘法运算，共需要 N ⋅ 2 p = 2 N ⋅ l o g N N·2p=2N·logN N⋅2p=2N⋅logN 次复数乘法运算。