在上一篇随笔“使用快速傅里叶变换计算大整数乘法”中，已经讲述了使用快速傅里叶变换计算大整数乘法的原理。在这一篇随笔中，我们就使用快速傅里叶变换来实现一个提供任意精度的算术运算的静态类：BigArithmetic。
下面就是 BigArithmetic 类源程序代码：
1 using System;
2 using System.Diagnostics;
3
4 namespace Skyiv.Numeric
5 {
6 /// <summary>
7 /// 提供任意精度的算术运算的静态类。使用快速傅里叶变换。
8 /// 本类对字节数组进行算术运算，字节数组以 100 为基。
9 /// 字节数组中第一个元素存储的数字是最高有效位。
10 /// </summary>
11 static class BigArithmetic
12 {
13 //= C语言数值算法程序大全(第二版)，ISBN 7-5053-2931-6 / TP 993
14 //= Numerical Recipes in C, The Art of Scientific Computing, Second Edition
15 //= Cambridge University Press 1988, 1992
16 //= [美] W.H.Press, S.A.Teukolsky, W.T.Vetterling, B.P.Flannery 著
17 //= 傅祖芸 赵梅娜 丁岩 等译，傅祖芸 校，电子工业出版社，1995年10月第一版
18
19 static readonly byte Len = 2; // 字节数组的元素包含的十进制数字的个数
20 static readonly byte Base = (byte)Math.Pow(10, Len); // 字节数组的基
21 static readonly byte MaxValue = (byte)(Base - 1); // 字节数组的元素的最大值
22
23 //= pp.431-432, four1, 12.2 快速傅里叶变换(FFT)
24 /// <summary>
25 /// 复函数的快速傅里叶变换
26 /// 变量 nn 是复数据点的个数，实型数组 data[1..2*nn] 的实际界长是两倍 nn，
27 /// 而每个复数值占据了两个相继的存储单元。nn 必须是 2 的整数幂
28 /// </summary>
29 /// <param name="data">实型数组 data[1..2*nn]。注意，下标从 1 开始</param>
30 /// <param name="isInverse">是否逆变换。注意: 逆变换未乘上归一化因子 1/nn</param>
31 public static void ComplexFFT(double[] data, bool isInverse)
32 {
33 int n = data.Length - 1; // n 必须是 2 的正整数幂
34 int nn = n >> 1; // 变量 nn 是复数据点的个数
35 for (int i = 1, j = 1; i < n; i += 2) // 这个循环实现位序颠倒
36 {
37 if (j > i)
38 {
39 Utility.Swap(ref data[j], ref data[i]);
40 Utility.Swap(ref data[j + 1], ref data[i + 1]);
41 }
42 int m = nn;
43 for (; m >= 2 && j > m; m >>= 1) j -= m;
44 j += m;
45 }
46 for (int mmax = 2, istep = 4; n > mmax; mmax = istep) // 执行 log2(nn) 次外循环
47 {
48 istep = mmax << 1; // 下面是关于三角递归的初始赋值
49 double theta = (isInverse ? -2 : 2) * Math.PI / mmax;
50 double wtemp = Math.Sin(0.5 * theta);
51 double wpr = -2 * wtemp * wtemp;
52 double wpi = Math.Sin(theta);
53 double wr = 1;
54 double wi = 0;
55 for (int m = 1; m < mmax; m += 2)
56 {
57 for (int i = m; i <= n; i += istep)
58 {
59 int j = i + mmax; // 下面是 Danielson-Lanczos 公式
60 double tempr = wr * data[j] - wi * data[j + 1];
61 double tempi = wr * data[j + 1] + wi * data[j];
62 data[j] = data[i] - tempr;
63 data[j + 1] = data[i + 1] - tempi;
64 data[i] += tempr;
65 data[i + 1] += tempi;
66 }
67 wr = (wtemp = wr) * wpr - wi * wpi + wr; // 三角递归
68 wi = wi * wpr + wtemp * wpi + wi;
69 }
70 }
71 }
72
2 using System.Diagnostics;
3
4 namespace Skyiv.Numeric
5 {
6 /// <summary>
7 /// 提供任意精度的算术运算的静态类。使用快速傅里叶变换。
8 /// 本类对字节数组进行算术运算，字节数组以 100 为基。
9 /// 字节数组中第一个元素存储的数字是最高有效位。
10 /// </summary>
11 static class BigArithmetic
12 {
13 //= C语言数值算法程序大全(第二版)，ISBN 7-5053-2931-6 / TP 993
14 //= Numerical Recipes in C, The Art of Scientific Computing, Second Edition
15 //= Cambridge University Press 1988, 1992
16 //= [美] W.H.Press, S.A.Teukolsky, W.T.Vetterling, B.P.Flannery 著
17 //= 傅祖芸 赵梅娜 丁岩 等译，傅祖芸 校，电子工业出版社，1995年10月第一版
18
19 static readonly byte Len = 2; // 字节数组的元素包含的十进制数字的个数
20 static readonly byte Base = (byte)Math.Pow(10, Len); // 字节数组的基
21 static readonly byte MaxValue = (byte)(Base - 1); // 字节数组的元素的最大值
22
23 //= pp.431-432, four1, 12.2 快速傅里叶变换(FFT)
24 /// <summary>
25 /// 复函数的快速傅里叶变换
26 /// 变量 nn 是复数据点的个数，实型数组 data[1..2*nn] 的实际界长是两倍 nn，
27 /// 而每个复数值占据了两个相继的存储单元。nn 必须是 2 的整数幂
28 /// </summary>
29 /// <param name="data">实型数组 data[1..2*nn]。注意，下标从 1 开始</param>
30 /// <param name="isInverse">是否逆变换。注意: 逆变换未乘上归一化因子 1/nn</param>
31 public static void ComplexFFT(double[] data, bool isInverse)
32 {
33 int n = data.Length - 1; // n 必须是 2 的正整数幂
34 int nn = n >> 1; // 变量 nn 是复数据点的个数
35 for (int i = 1, j = 1; i < n; i += 2) // 这个循环实现位序颠倒
36 {
37 if (j > i)
38 {
39 Utility.Swap(ref data[j], ref data[i]);
40 Utility.Swap(ref data[j + 1], ref data[i + 1]);
41 }
42 int m = nn;
43 for (; m >= 2 && j > m; m >>= 1) j -= m;
44 j += m;
45 }
46 for (int mmax = 2, istep = 4; n > mmax; mmax = istep) // 执行 log2(nn) 次外循环
47 {
48 istep = mmax << 1; // 下面是关于三角递归的初始赋值
49 double theta = (isInverse ? -2 : 2) * Math.PI / mmax;
50 double wtemp = Math.Sin(0.5 * theta);
51 double wpr = -2 * wtemp * wtemp;
52 double wpi = Math.Sin(theta);
53 double wr = 1;
54 double wi = 0;
55 for (int m = 1; m < mmax; m += 2)
56 {
57 for (int i = m; i <= n; i += istep)
58 {
59 int j = i + mmax; // 下面是 Danielson-Lanczos 公式
60 double tempr = wr * data[j] - wi * data[j + 1];
61 double tempi = wr * data[j + 1] + wi * data[j];
62 data[j] = data[i] - tempr;
63 data[j + 1] = data[i + 1] - tempi;
64 data[i] += tempr;
65 data[i + 1] += tempi;
66 }
67 wr = (wtemp = wr) * wpr - wi * wpi + wr; // 三角递归
68 wi = wi * wpr + wtemp * wpi + wi;
69 }
70 }
71 }
72
该类的第一个方法是 ComplexFFT，计算复函数的快速傅里叶变换。注意，ComplexFFT 并没有使用复数(不象 C++，C# 也没有提供复数)，而是让每个复数值占据实型数组的两个相继的存储单元。还有，要求输入的复数据点的个数必须是 2 的整数幂。该方法也能够计算复函数的快速傅里叶逆变换。
该程序的算法是使用 1942 年 Danielson 和 Lanczos 证明的引理：一个界长为 N 的离散傅里叶变换可以重新写成两个界长各为 N/2 的离散傅里叶变换之和。在算法的第一部分，将数据整理成位序颠倒的次序。而在第二部分，有一个执行 log2N 次的外循环。它依次计算界长为 2, 4, 8, ..., N 的变换。对于这一过程的每一步来说，为了履行 Danielson-Lanczos 引理，有两个嵌套的内循环，其涉及到已计算的子变换和每个变换的元素。通过限制外部调用正弦和余弦到外层循环，可以使运算更有效，在外层循环中只要调用它们 log2N 次。倍角的正弦和余弦的计算是通过内循环中简单的递归关系进行的，如下所示：
cos(θ + δ) = cosθ - [ α cosθ + βsinθ ]
sin(θ + δ) = sinθ - [ α sinθ - βcosθ ]
其中 α, β 是预先计算的系数：α = 2 sin2(δ/2)， β = sinδ 。
73 //= pp.436, realft, 12.3.2 单个实函数的 FFT
74 /// <summary>
75 /// 单个实函数的快速傅里叶变换
76 /// 计算一组 n 个实值数据点的傅里叶变换。用复傅里叶变换的正半频率替换这些数据，
77 /// 它存储在数组 data[1..n] 中。复变换的第一个和最后一个分量的实数值分别返回
78 /// 单元 data[1] 和 data[2] 中。n 必须是 2 的幂次。这个程序也能计算复数据数组
79 /// 的逆变换，只要该数组是实值数据的变换(在这种情况下，其结果必须乘以 1/n)即可。
80 /// </summary>
81 /// <param name="data">实型数组 data[1..n]。注意，下标从 1 开始</param>
82 /// <param name="isInverse">是否逆变换。注意: 逆变换未乘上归一化因子 1/n</param>
83 public static void RealFFT(double[] data, bool isInverse)
84 {
85 int n = data.Length - 1; // n 必须是 2 的整数幂
86 if (!isInverse) ComplexFFT(data, isInverse); // 此处是正向变换
87 double theta = (isInverse ? -2 : 2) * Math.PI / n; // 递归的初始赋值
88 double wtemp = Math.Sin(0.5 * theta);
89 double wpr = -2 * wtemp * wtemp;
90 double wpi = Math.Sin(theta);
91 double wr = 1 + wpr;
92 double wi = wpi;
93 double c1 = 0.5;
94 double c2 = isInverse ? 0.5 : -0.5;
95 int n3 = n + 3;
96 int n4 = n >> 2;
97 for (int i = 2; i <= n4; i++)
98 {
99 int i1 = i + i - 1, i2 = i1 + 1, i3 = n3 - i2, i4 = i3 + 1;
100 double h1r = c1 * (data[i1] + data[i3]); // 两个分离变换是
101 double h1i = c1 * (data[i2] - data[i4]); // 从 data 中分离出来
102 double h2r = -c2 * (data[i2] + data[i4]);
103 double h2i = c2 * (data[i1] - data[i3]);
104 data[i1] = h1r + wr * h2r - wi * h2i; // 此处重新组合以形成
105 data[i2] = h1i + wr * h2i + wi * h2r; // 原始实型数据的真实变换
106 data[i3] = h1r - wr * h2r + wi * h2i;
107 data[i4] = -h1i + wr * h2i + wi * h2r;
108 wr = (wtemp = wr) * wpr - wi * wpi + wr; // 递归式
109 wi = wi * wpr + wtemp * wpi + wi;
110 }
111 double tmp = data[1];
112 if (!isInverse)
113 {
114 data[1] = tmp + data[2]; // 同时挤压第一个和最后一个数据
115 data[2] = tmp - data[2]; // 使它们都在原始数组中
116 }
117 else
118 {
119 data[1] = c1 * (tmp + data[2]);
120 data[2] = c1 * (tmp - data[2]);
121 ComplexFFT(data, isInverse); // 此处是逆变换
122 }
123 }
124
74 /// <summary>
75 /// 单个实函数的快速傅里叶变换
76 /// 计算一组 n 个实值数据点的傅里叶变换。用复傅里叶变换的正半频率替换这些数据，
77 /// 它存储在数组 data[1..n] 中。复变换的第一个和最后一个分量的实数值分别返回
78 /// 单元 data[1] 和 data[2] 中。n 必须是 2 的幂次。这个程序也能计算复数据数组
79 /// 的逆变换，只要该数组是实值数据的变换(在这种情况下，其结果必须乘以 1/n)即可。
80 /// </summary>
81 /// <param name="data">实型数组 data[1..n]。注意，下标从 1 开始</param>
82 /// <param name="isInverse">是否逆变换。注意: 逆变换未乘上归一化因子 1/n</param>
83 public static void RealFFT(double[] data, bool isInverse)
84 {
85 int n = data.Length - 1; // n 必须是 2 的整数幂
86 if (!isInverse) ComplexFFT(data, isInverse); // 此处是正向变换
87 double theta = (isInverse ? -2 : 2) * Math.PI / n; // 递归的初始赋值
88 double wtemp = Math.Sin(0.5 * theta);
89 double wpr = -2 * wtemp * wtemp;
90 double wpi = Math.Sin(theta);
91 double wr = 1 + wpr;
92 double wi = wpi;
93 double c1 = 0.5;
94 double c2 = isInverse ? 0.5 : -0.5;
95 int n3 = n + 3;
96 int n4 = n >> 2;
97 for (int i = 2; i <= n4; i++)
98 {
99 int i1 = i + i - 1, i2 = i1 + 1, i3 = n3 - i2, i4 = i3 + 1;
100 double h1r = c1 * (data[i1] + data[i3]); // 两个分离变换是
101 double h1i = c1 * (data[i2] - data[i4]); // 从 data 中分离出来
102 double h2r = -c2 * (data[i2] + data[i4]);
103 double h2i = c2 * (data[i1] - data[i3]);
104 data[i1] = h1r + wr * h2r - wi * h2i; // 此处重新组合以形成
105 data[i2] = h1i + wr * h2i + wi * h2r; // 原始实型数据的真实变换
106 data[i3] = h1r - wr * h2r + wi * h2i;
107 data[i4] = -h1i + wr * h2i + wi * h2r;
108 wr = (wtemp = wr) * wpr - wi * wpi + wr; // 递归式
109 wi = wi * wpr + wtemp * wpi + wi;
110 }
111 double tmp = data[1];
112 if (!isInverse)
113 {
114 data[1] = tmp + data[2]; // 同时挤压第一个和最后一个数据
115 data[2] = tmp - data[2]; // 使它们都在原始数组中
116 }
117 else
118 {
119 data[1] = c1 * (tmp + data[2]);
120 data[2] = c1 * (tmp - data[2]);
121 ComplexFFT(data, isInverse); // 此处是逆变换
122 }
123 }
124
第二个方法是 RealFFT，计算单个实函数的快速傅里叶变换。因为在很多情况下期望求快速傅里叶变换的数据由实值样本组成。如果将这些样本放入一个复型数组中，并令其所有虚部为零的话，从执行时间和对存储需求二者来看，其效率是很低的。所以，我们将原始数据放入一个界长只有一半的复型数组中，其偶数项放入该数组的实部，奇数项放入它的虚部。然后调用 ComplexFFT 来计算快速傅里叶变换。
125 /// <summary>
126 /// 比较 x[0..n-1] 和 y[0..n-1]
127 /// </summary>
128 /// <param name="x">第一操作数 x[0..n-1]</param>
129 /// <param name="y">第二操作数 y[0..n-1]</param>
130 /// <param name="n">两个操作数 x 和 y 的精度</param>
131 /// <returns>比较结果：-1:小于 1:大于 0:等于</returns>
132 public static int Compare(byte[] x, byte[] y, int n)
133 {
134 Debug.Assert(x.Length >= n && y.Length >= n);
135 for (int i = 0; i < n; i++)
136 if (x[i] != y[i])
137 return (x[i] < y[i]) ? -1 : 1;
138 return 0;
139 }
140
141 //= pp.775, mpneg, 20.6 任意精度的运算
142 /// <summary>
143 /// 求补码。注意，被操作数被修改。
144 /// </summary>
145 /// <param name="data">被操作数 data[0..n-1]</param>
146 /// <param name="n">被操作数 data 的精度</param>
147 /// <returns>被操作数的补码 data[0..n-1]</returns>
148 public static byte[] Negative(byte[] data, int n)
149 {
150 Debug.Assert(data.Length >= n);
151 for (int k = Base, i = n - 1; i >= 0; i--)
152 data[i] = (byte)((k = MaxValue + k / Base - data[i]) % Base);
153 return data;
154 }
155
156 //= pp.774, mpsub, 20.6 任意精度的运算
157 /// <summary>
158 /// 减法。从 minuend[0..n-1] 中减去 subtrahend[0..n-1]，得到 difference[0..n-1]
159 /// </summary>
160 /// <param name="difference">差 difference[0..n-1]</param>
161 /// <param name="minuend">被减数 minuend[0..n-1]</param>
162 /// <param name="subtrahend">减数 subtrahend[0..n-1]</param>
163 /// <param name="n">被减数 minuend 和减数 subtrahend 的精度</param>
164 /// <returns>差 difference[0..n-1]</returns>
165 public static byte[] Subtract(byte[] difference, byte[] minuend, byte[] subtrahend, int n)
166 {
167 Debug.Assert(minuend.Length >= n && subtrahend.Length >= n && difference.Length >= n);
168 for (int k = Base, i = n - 1; i >= 0; i--)
169 difference[i] = (byte)((k = MaxValue + k / Base + minuend[i] - subtrahend[i]) % Base);
170 return difference;
171 }
172
173 //= pp.774, mpadd, 20.6 任意精度的运算
174 /// <summary>
175 /// 加法。augend[0..n-1] 与 addend[0..n-1] 相加，得到 sum[0..n]
176 /// </summary>
177 /// <param name="sum">和 sum[0..n]</param>
178 /// <param name="augend">被加数 augend[0..n-1]</param&g
126 /// 比较 x[0..n-1] 和 y[0..n-1]
127 /// </summary>
128 /// <param name="x">第一操作数 x[0..n-1]</param>
129 /// <param name="y">第二操作数 y[0..n-1]</param>
130 /// <param name="n">两个操作数 x 和 y 的精度</param>
131 /// <returns>比较结果：-1:小于 1:大于 0:等于</returns>
132 public static int Compare(byte[] x, byte[] y, int n)
133 {
134 Debug.Assert(x.Length >= n && y.Length >= n);
135 for (int i = 0; i < n; i++)
136 if (x[i] != y[i])
137 return (x[i] < y[i]) ? -1 : 1;
138 return 0;
139 }
140
141 //= pp.775, mpneg, 20.6 任意精度的运算
142 /// <summary>
143 /// 求补码。注意，被操作数被修改。
144 /// </summary>
145 /// <param name="data">被操作数 data[0..n-1]</param>
146 /// <param name="n">被操作数 data 的精度</param>
147 /// <returns>被操作数的补码 data[0..n-1]</returns>
148 public static byte[] Negative(byte[] data, int n)
149 {
150 Debug.Assert(data.Length >= n);
151 for (int k = Base, i = n - 1; i >= 0; i--)
152 data[i] = (byte)((k = MaxValue + k / Base - data[i]) % Base);
153 return data;
154 }
155
156 //= pp.774, mpsub, 20.6 任意精度的运算
157 /// <summary>
158 /// 减法。从 minuend[0..n-1] 中减去 subtrahend[0..n-1]，得到 difference[0..n-1]
159 /// </summary>
160 /// <param name="difference">差 difference[0..n-1]</param>
161 /// <param name="minuend">被减数 minuend[0..n-1]</param>
162 /// <param name="subtrahend">减数 subtrahend[0..n-1]</param>
163 /// <param name="n">被减数 minuend 和减数 subtrahend 的精度</param>
164 /// <returns>差 difference[0..n-1]</returns>
165 public static byte[] Subtract(byte[] difference, byte[] minuend, byte[] subtrahend, int n)
166 {
167 Debug.Assert(minuend.Length >= n && subtrahend.Length >= n && difference.Length >= n);
168 for (int k = Base, i = n - 1; i >= 0; i--)
169 difference[i] = (byte)((k = MaxValue + k / Base + minuend[i] - subtrahend[i]) % Base);
170 return difference;
171 }
172
173 //= pp.774, mpadd, 20.6 任意精度的运算
174 /// <summary>
175 /// 加法。augend[0..n-1] 与 addend[0..n-1] 相加，得到 sum[0..n]
176 /// </summary>
177 /// <param name="sum">和 sum[0..n]</param>
178 /// <param name="augend">被加数 augend[0..n-1]</param&g