BigArithmetic - 提供任意精度的算术运算的静态类

最新推荐文章于 2024-06-29 17:47:30 发布

weixin_33940102

最新推荐文章于 2024-06-29 17:47:30 发布

阅读量209

点赞数 1

文章标签： ui c/c++ c#

在上一篇随笔“使用快速傅里叶变换计算大整数乘法”中，已经讲述了使用快速傅里叶变换计算大整数乘法的原理。在这一篇随笔中，我们就使用快速傅里叶变换来实现一个提供任意精度的算术运算的静态类：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 的离散傅里叶变换之和。在算法的第一部分，将数据整理成位序颠倒的次序。而在第二部分，有一个执行 log₂N 次的外循环。它依次计算界长为 2, 4, 8, ..., N 的变换。对于这一过程的每一步来说，为了履行 Danielson-Lanczos 引理，有两个嵌套的内循环，其涉及到已计算的子变换和每个变换的元素。通过限制外部调用正弦和余弦到外层循环，可以使运算更有效，在外层循环中只要调用它们 log₂N 次。倍角的正弦和余弦的计算是通过内循环中简单的递归关系进行的，如下所示：

cos(θ + δ) = cosθ - [ α cosθ + βsinθ ]
sin(θ + δ) = sinθ - [ α sinθ - βcosθ ]

其中 α, β 是预先计算的系数：α = 2 sin²(δ/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，求倒数。

除法用除数的倒数乘以被除数来计算，倒数值用牛顿迭代法：

U_i+1 = U_i (2 - VU_i)

来计算，这导致 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，求平方根。用牛顿迭代法计算平方根和除法类似。若：

U_i+1 = U_i (3 - VU_i²) / 2

则 U_∞二次收敛于 1/√V，最后乘以 V 就得到√V。

在上一篇随笔“使用快速傅里叶变换计算大整数乘法”中提到：

由于快速傅里叶变换是采用了浮点运算，因此我们需要足够的精度，以使在出现舍入误差时，结果中每个组成部分的准确整数值仍是可辨认的。长度为 N 的 B 进制数可产生大到 B²N 阶的卷积分量。我们知道，双精度浮点数的尾数是 53 个二进位，所以：

2 x log₂B + log₂N + 几个 x log₂log₂N < 53

上式中左边最后一项是为了快速傅里叶变换的舍入误差。

我们假设上式中的“几个”为“三个”，那么，经过简单的计算，得到以下结果：