在上一篇隨筆“使用快速傅裡葉變換計算大整數乘法”中,已經講述了使用快速傅裡葉變換計算大整數乘法的原理。在這一篇隨筆中,我們就使用快速傅裡葉變換來實現一個提供任意精度的算術運算的靜態類: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
復制代碼
該類的第一個方法是 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
復制代碼
第二個方法是 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>
179 /// <param name="addend">加數 addend[0..n-1]</param>
180 /// <param name="n">被加數 augend 和加數 addend 的精度</param>
181 /// <returns>和 sum[0..n]</returns>
182 public static byte[] Add(byte[] sum, byte[] augend, byte[] addend, int n)
183 {
184 Debug.Assert(augend.Length >= n && addend.Length >= n && sum.Length >= n + 1);
185 int k = 0;
186 for (int i = n - 1; i >= 0; i--)
187 sum[i + 1] = (byte)((k = k / Base + augend[i] + addend[i]) % Base);
188 sum[0] += (byte)(k / Base);
189 return sum;
190 }
191
192 //= pp.774, mpadd, 20.6 任意精度的運算
193 /// <summary>
194 /// 捷加法。augend[0..n-1] 與整數 addend 相加,得到 sum[0..n]
195 /// </summary>
196 /// <param name="sum">和 sum[0..n]</param>
197 /// <param name="augend">被加數 augend[0..n-1]</param>
198 /// <param name="n">被加數 augend 的精度</param>
199 /// <param name="addend">加數 addend</param>
200 /// <returns>和 sum[0..n]</returns>
201 public static byte[] Add(byte[] sum, byte[] augend, int n, byte addend)
202 {
203 Debug.Assert(augend.Length >= n && sum.Length >= n + 1);
204 int k = Base * addend;
205 for (int i = n - 1; i >= 0; i--)
206 sum[i + 1] = (byte)((k = k / Base + augend[i]) % Base);
207 sum[0] += (byte)(k / Base);
208 return sum;
209 }
210
211 //= pp.775, mpsdv, 20.6 任意精度的運算
212 /// <summary>
213 /// 捷除法。dividend[0..n-1] 除以整數 divisor,得到 quotient[0..n-1]
214 /// </summary>
215 /// <param name="quotient">商 quotient[0..n-1]</param>
216 /// <param name="dividend">被除數 dividend[0..n-1]</param>
217 /// <param name="n">被除數 dividend 的精度</param>
218 /// <param name="divisor">除數 divisor</param>
219 /// <returns>商 quotient[0..n-1]</returns>
220 public static byte[] Divide(byte[] quotient, byte[] dividend, int n, byte divisor)
221 {
222 Debug.Assert(quotient.Length >= n && dividend.Length >= n);
223 for (int r = 0, k = 0, i = 0; i < n; i++, r = k % divisor)
224 quotient[i] = (byte)((k = Base * r + dividend[i]) / divisor);
225 return quotient;
226 }
227
復制代碼
接下來是比較、求補碼、減法、加法、捷加法、捷除法,都相當簡單,程序中已經有詳細的注釋了。
復制代碼
228 //= pp.776-777, mpmul, 20.6 任意精度的運算
229 /// <summary>
230 /// 乘法。multiplicand[0..n-1] 與 multiplier[0..m-1] 相乘,得到 product[0..n+m-1]
231 /// </summary>
232 /// <param name="product">積 product[0..n+m-1]</param>
233 /// <param name="multiplicand">被乘數 multiplicand[0..n-1]</param>
234 /// <param name="n">被乘數 multiplicand 的精度</param>
235 /// <param name="multiplier">乘數 multiplier[0..m-1]</param>
236 /// <param name="m">乘數 multiplier 的精度</param>
237 /// <returns>積 product[0..n+m-1]</returns>
238 public static byte[] Multiply(byte[] product, byte[] multiplicand, int n, byte[] multiplier, int m)
239 {
240 int mn = m + n, nn = 1;
241 Debug.Assert(product.Length >= mn && multiplicand.Length >= n && multiplier.Length >= m);
242 while (nn < mn) nn <<= 1; // 為變換找出最小可用的 2 的冪次
243 double[] a = new double[nn + 1], b = new double[nn + 1];
244 for (int i = 0; i < n; i++) a[i + 1] = multiplicand[i];
245 for (int i = 0; i < m; i++) b[i + 1] = multiplier[i];
246 RealFFT(a, false); // 執行卷積,首先求出二個傅裡葉變換
247 RealFFT(b, false);
248 b[1] *= a[1]; // 復數相乘的結果(實部和虛部)
249 b[2] *= a[2];
250 for (int i = 3; i <= nn; i += 2)
251 {
252 double t = b[i];
253 b[i] = t * a[i] - b[i + 1] * a[i + 1];
254 b[i + 1] = t * a[i + 1] + b[i + 1] * a[i];
255 }
256 RealFFT(b, true); // 進行傅裡葉逆變換
257 byte[] bs = new byte[nn + 1];
258 long cy = 0; // 執行最後完成所有進位的過程
259 for (int i = nn, n2 = nn / 2; i >= 1; i--)
260 {
261 long t = (long)(b[i] / n2 + cy + 0.5);
262 bs[i] = (byte)(t % Base); // 原書中這句使用循環,有嚴重的性能問題
263 cy = t / Base;
264 }
265 if (cy >= Base) throw new OverflowException("FFT Multiply");
266 bs[0] = (byte)cy;
267 Array.Copy(bs, product, n + m);
268 return product;
269 }
270
復制代碼
接下來的方法是 Multiply,乘法。其算法是使用 RealFFT 求被乘數和乘數的快速傅裡葉變換,將結果相乘,然後進行傅裡葉逆變換得到卷積,最後執行適當的進位。其原理已經在上一篇隨筆“使用快速傅裡葉變換計算大整數乘法”中講述得很清楚了。
復制代碼
271 //= pp.778, mpdiv, 20.6 任意精度的運算
272 /// <summary>
273 /// 除法。dividend[0..n-1] 除以 divisor[0..m-1],m ≤ n,
274 /// 得到:商 quotient[0..n-m],余數 remainder[0..m-1]。
275 /// </summary>
276 /// <param name="quotient">商 quotient[0..n-m]</param>
277 /// <param name="remainder">余數 remainder[0..m-1]</param>
278 /// <param name="dividend">被除數 dividend[0..n-1]</param>
279 /// <param name="n">被除數 dividend 的精度</param>
280 /// <param name="divisor">除數 divisor[0..m-1]</param>
281 /// <param name="m">除數 divisor 的精度</param>
282 /// <returns>商 quotient[0..n-m]</returns>
283 public static byte[] DivRem(byte[] quotient, byte[] remainder, byte[] dividend, int n, byte[] divisor, int m)
264 {
285 Debug.Assert(m <= n && dividend.Length >= n && divisor.Length >= m && quotient.Length >= n - m + 1 && remainder.Length >= m);
286 int MACC = 3;
287 byte[] s = new byte[n + MACC], t = new byte[n - m + MACC + n];
288 Inverse(s, n - m + MACC, divisor, m); // s = 1 / divisor
289 Array.Copy(Multiply(t, s, n - m + MACC, dividend, n), 1, quotient, 0, n - m + 1); // quotient = dividend / divisor
290 Array.Copy(Multiply(t, quotient, n - m + 1, divisor, m), 1, s, 0, n); // s = quotient * divisor
291 Subtract(s, dividend, s, n); // s = dividend - quotient * divisor
292 Array.Copy(s, n - m, remainder, 0, m);
293 if (Compare(remainder, divisor, m) >= 0) // 調整商和余數
294 {
295 Subtract(remainder, remainder, divisor, m);
296 Add(s, quotient, n - m + 1, 1);
297 Array.Copy(s, 1, quotient, 0, n - m + 1);
298 }
299 return quotient;
300 }
301
302 //= pp.777 - 778, mpinv, 20.6 任意精度的運算
303 /// <summary>
304 /// 求倒數。
305 /// </summary>
306 /// <param name="inverse">倒數 inverse[0..n-1],在 inverse[0] 後有基數的小數點</param>
307 /// <param name="n">倒數 inverse 的精度</param>
308 /// <param name="data">被操作數 data[0..m-1],data[0] > 0,在 data[0] 後有基數的小數點</param>
309 /// <param name="m">被操作數 data 的精度</param>
310 /// <returns>倒數 inverse[0..n-1],在 inverse[0] 後有基數的小數點</returns>
311 public static byte[] Inverse(byte[] inverse, int n, byte[] data, int m)
312 {
313 Debug.Assert(inverse.Length >= n && data.Length >= m);
314 InitialValue(inverse, n, data, m, false);
315 if (n == 1) return inverse;
316 byte[] s = new byte[n], t = new byte[n + n];
317 for (; ; ) // 牛頓迭代法: inverse = inverse * ( 2 - data * inverse ) => inverse = 1 / data
318 {
319 Array.Copy(Multiply(t, inverse, n, data, m), 1, s, 0, n); // s = data * inverse
320 Negative(s, n); // s = -(data * inverse)
321 s[0] -= (byte)(Base - 2); // s = 2 - data * inverse
322 Array.Copy(Multiply(t, s, n, inverse, n), 1, inverse, 0, n); // inverse = inverse * s
323 int i = 1;
324 for (; i < n - 1 && s[i] == 0; i++) ; // 判斷 s 的小數部分是否為零
325 if (i == n - 1) return inverse; // 若 s 收斂到 1 則返回 inverse = 1 / data
326 }
327 }
328
復制代碼
接下來的方法是 DivRem,除法,以及 Inverse,求倒數。
除法用除數的倒數乘以被除數來計算,倒數值用牛頓迭代法:
Ui+1 = Ui (2 - VUi)
來計算,這導致 U∞二次收斂於 1/V。實際上,許多超級計算機和 RISC 機都是使用這種迭代法來實現除法的。
要注意在 DivRem 中,求得的余數可能大於等於除數,此時商比實際的值要小一。所以在程序的 293 到 298 行調整商和余數。
復制代碼
329 //= pp.778-779, mpsqrt, 20.6 任意精度的運算
330 /// <summary>
331 /// 求平方根 sqrt,以及平方根的倒數 invSqrt。invSqrt 也可設為 sqrt,此時,invSqrt 也是平方根。
332 /// </summary>
333 /// <param name="sqrt">平方根 sqrt[0..n-1],在 sqrt[0] 後有基數的小數點</param>
334 /// <param name="invSqrt">平方根的倒數 invSqrt[0..n-1],在 invSqrt[0] 後有基數的小數點</param>
335 /// <param name="n">平方根的精度</param>
336 /// <param name="data">被操作數 data[0..m-1],data[0] > 0,在 data[0] 後有基數的小數點</param>
337 /// <param name="m">被操作數 data 的精度</param>
338 /// <returns>平方根 sqrt[0..n-1],在 sqrt[0] 後有基數的小數點</returns>
339 public static byte[] Sqrt(byte[] sqrt, byte[] invSqrt, int n, byte[] data, int m)
340 {
341 Debug.Assert(sqrt.Length >= n && invSqrt.Length >= n && data.Length >= m);
342 if (n <= 1) throw new ArgumentOutOfRangeException("n", "must greater than 1");
343 InitialValue(invSqrt, n, data, m, true);
344 byte[] s = new byte[n], t = new byte[n + Math.Max(m, n)];
345 for (; ; ) // invSqrt = invSqrt * (3 - data * invSqrt * invSqrt) / 2 => invSqrt = 1 / sqrt(data)
346 {
347 Array.Copy(Multiply(t, invSqrt, n, invSqrt, n), 1, s, 0, n); // s = invSqrt * invSqrt
348 Array.Copy(Multiply(t, s, n, data, m), 1, s, 0, n); // s = data * invSqrt * invSqrt
349 Negative(s, n); // s = -(data * invSqrt * invSqrt)
350 s[0] -= (byte)(Base - 3); // s = 3 - data * invSqrt * invSqrt
351 Divide(s, s, n, 2); // s = (3 - data * invSqrt * invSqrt) / 2
352 Array.Copy(Multiply(t, s, n, invSqrt, n), 1, invSqrt, 0, n); // invSqrt = invSqrt * s
353 int i = 1;
354 for (; i < n - 1 && s[i] == 0; i++) ; // 判斷 s 的小數部分是否為零
355 if (i < n - 1) continue; // 若 s 沒有收斂到 1 則繼續迭代
356 Array.Copy(Multiply(t, invSqrt, n, data, m), 1, sqrt, 0, n); // sqrt = invSqrt * data = sqrt(data)
357 return sqrt;
358 }
359 }
360
361 /// <summary>
362 /// 采用浮點運算以得到一個初始近似值 u[0..n-1]: u = 1 / data 或者 u = 1 / sqrt(data)
363 /// </summary>
364 /// <param name="u">初始近似值 u[0..n-1]</param>
365 /// <param name="n">所需的精度</param>
366 /// <param name="data">被操作數 data[0..m-1]</param>
367 /// <param name="m">被操作數 data 的精度</param>
368 /// <param name="isSqrt">是否求平方根</param>
369 /// <returns>初始近似值 u[0..n-1]</returns>
370 static byte[] InitialValue(byte[] u, int n, byte[] data, int m, bool isSqrt)
371 {
372 Debug.Assert(u.Length >= n && data.Length >= m);
373 int scale = 16 / Len; // double 可達到 16 位有效數字
374 double fu = 0;
375 for (int i = Math.Min(scale, m) - 1; i >= 0; i--) fu = fu / Base + data[i];
376 fu = 1 / (isSqrt ? Math.Sqrt(fu) : fu);
377 for (int i = 0; i < Math.Min(scale + 1, n); i++)
378 {
379 int k = (int)fu;
380 u[i] = (byte)k;
381 fu = Base * (fu - k);
382 }
383 return u;
384 }
385 }
386 }
復制代碼
最後的方法是 Sqrt,求平方根。用牛頓迭代法計算平方根和除法類似。若:
Ui+1 = Ui (3 - VUi2) / 2
則 U∞二次收斂於 1/√V,最後乘以 V 就得到√V。
在上一篇隨筆“使用快速傅裡葉變換計算大整數乘法”中提到:
由於快速傅裡葉變換是采用了浮點運算,因此我們需要足夠的精度,以使在出現捨入誤差時,結果中每個組成部分的准確整數值仍是可辨認的。長度為 N 的 B 進制數可產生大到 B2N 階的卷積分量。我們知道,雙精度浮點數的尾數是 53 個二進位,所以:
2 x log2B + log2N + 幾個 x log2log2N < 53
上式中左邊最後一項是為了快速傅裡葉變換的捨入誤差。
我們假設上式中的“幾個”為“三個”,那麼,經過簡單的計算,得到以下結果:
注意,基數 B 選取得越小,計算速度就越慢。
在下一篇隨筆中,我將使用 BigArithmetic 類重新實現 Skyiv.Numeric.BigInteger 類。