[转]浅谈 BigInteger

http://www.cnblogs.com/skyivben/archive/2008/07/13/1241681.html

最近一段时间，我在 Timus Online Judge 网站做 ACM 题。发现其中不少题目需要用到 BigInteger，例如：

• 1531. Zones on a plane
• 1013. K-based numbers. Version 3
• 1012. K-based numbers. Version 2
• 1153. Supercomputer
• 1108. Heritage

于是，我就开始寻找可以用于 C# 的 BigInteger 类。结果如下表所示：

项目	速度	自动增长	进制	ModPow	Pow	Sqrt	源文件大小	备注
Fcl35Biginteger	最快	是	232	是	是	否	1,755行，47.1 KB	是 struct，不是 class
SkyivBigInteger	快	是	109	否	是	是	335行，9.16 KB	源文件短小精悍，功能全
ChewBigInteger	慢	否	232	是	否	是	3,335行，102 KB	源文件中有大量的注释
JavaBigInteger	极慢	是	28?	是	是	否	(vjslib.dll)	需要 J# run-time library

在上表中，JavaBigInteger 是指 java.math.BigInteger，她在 vjslib.dll 中，需要 J# run-time library (Visual Studio 2005/2008 中已经包括 J# run-time library)。用 .NET Reflector 打开 C:\Windows\Microsoft.NET\Framework\v2.0.50727\vjslib.dll 文件，找到 java.math.BigInteger，查看其源程序代码如下：

1 [Serializable, JavaInterfaces("2;System/IConvertible;java/lang/Comparable;")]
2 public class BigInteger : Number, IConvertible, Comparable
3 {
4   // Fields
5   private BigInteger __absOfthis;
6   [JavaFlags(0x1002)]
7   internal sbyte[] __bits;
8   // 
9 }

根据第 7 行的 internal sbyte[] __bits; 猜想 java.math.BigInteger 的内部实现是基于 2⁸ 进制的。

ChewBigInteger 就是 Chew Keong TAN 在 CodeProject 中的 C# BigInteger Class 。这是个很有名的 C# Biginteger 类，它写于2002年8月到9月，已经有很久的历史了。但是她有个致命的缺点，就是不能自动增长，需要事先指定最大允许有多大的数，如果在运算过程中产生超过你事先指定大小的数，就会引发异常。

她的部分源程序代码如下所示：

 1 namespace ChewKeongTAN
 2 {
 3   /*public */class BigInteger
 4   {
 5     // maximum length of the BigInteger in uint (4 bytes)
 6     // change this to suit the required level of precision.
 7 
 8     private const int maxLength = 26591; // == Math.Ceiling ( digits / log10(2^32) )
 9 
10     private uint[] data = null;          // stores bytes from the Big Integer
11     public int dataLength;               // number of actual chars used
12     
13     // 
14     
15   }
16 }

她内部使用一个 uint[] 来存储，是基于 2³² ≈ 4.3 x 10⁹ 进制的。需要事先设定 maxLength 的值(程序中第 8 行)，用以指定最多允许多少个 2³² 进制的数字，换算用的公式是：

maxLength = Math.Ceiling ( digits / log₁₀(2³²) ) ≈ digits / 9.6

上式中 digits 表示十进制数字的个数。

Fcl35Biginteger 就是 .NET Framework 3.5 Base Class Library 中的 System.Numeric.BigInteger，她在 System.Core.dll 中。注意，她是一个 internal struct。首先，她是一个 struct，不是一个 class，也就是说她是值类型而不是引用类型。其次，致命的是，她是 internal 的，而不是 public 的，也就是说，她就在那里，但是你就是不能使用她。实际上，在 .NET Framework 3.5 的 Beta 和 latest CTP 版本中，她还是 public 的，只不过是在 .NET Framework 3.5 正式发布时变成了 internal 的。这里有一篇 BCL Team 的 Melitta Andersen 于2008年1月4日发表的 Blog: Where did BigInteger go? 谈到这件事，要点是：

So why was BigInteger cut?  The basic rationale behind making BigInteger internal was that it just wasn't ready to ship.  We thought our implementation met the needs for a BigInteger type.  But then we had some other teams take a look at it and they pointed out some performance and compatibility issues that we just didn't have time to fix before we shipped.

于是， .NET Reflector 又上场了，用她打开 C:\Windows\assembly\GAC_MSIL\System.Core\3.5.0.0__b77a5c561934e089\System.Core.dll 文件，找到 System.Numeric.BigInteger，在 Disassembler 窗口中，点击最下方的 Expand Methods，然后把源程序代码全部复制到 Visual Studio 2008 的 IDE 窗口(注意，目标 Framework 设定为 .NET Framework 2.0)中。按下编译按钮，天啊，出现了好多错误。于是，就开始修改源程序代码：

• 实现两个空的 sealed class: ImmutableAttribute 和 PureAttribute。(前一个好理解，后一个不知是什么东东。当然，实际上她们不应该是空的)
• 实现一个 static class: Contract，其中包括两个重载的 public static void Requires 方法。(当然，实际上的 Contract.Requires 应该更复杂)
• 实现一个 static class: Res，其中包括若干个 public static readonly string 字段。(当然，实际上 Res 应该从资源文件中获得本地化的数据)
• 将程序中出现的某些常数值还原为 private const 字段。注意，这不是必须的，只是为了使程序更容易看懂。如果不做这一步，就可以将程序中所有的 private const 字段全部删除。
  
• 修改程序中某些变量的名称。注意，这也不是必须的，也只是为了使程序更容易看懂。
  

然后，再次编译，耶，通过了。但是，还是高兴得太早了一点，运行一些测试，发现答案根本不对。

经检查，发现是 .NET Reflector 在反汇编为 C# 代码时会丢失一些强制类型转换，这在我的上一篇随笔 浅谈 Math.BigMul 方法 中已经提到过了：她把 Math.BigMul 方法中的唯一的语句 return (long)a * b; 错误地反汇编为 return (a * b); 了，但是反汇编为 IL 代码是正确的。

于是，只好逐条语句检查，如果需要进行强制类型转换，就自己加上，同时还去掉一些不必要的强制类型转换和括号，以便使源程序代码更易读。真是一个很烦人和困难和工作，热切盼望 .NET Reflector 能够尽快修复这个 BUG。:)

好了，经过调试，程序终于运行正常了。下面就是修改后的部分源程序代码：

 1 namespace Fcl35.Numeric
 2 {
 3   [Serializable, StructLayout(LayoutKind.Sequential), Immutable, ComVisible(false)]
 4   public /*internal*/ struct BigInteger : IFormattable,
 5     IEquatable<BigInteger>, IComparable<BigInteger>, IComparable
 6   {
 7     private const ulong Base = uint.MaxValue + 1UL; // 0x100000000L;
 8     private readonly short _sign;
 9     private readonly uint[] _data;
10     private int _length;
11     // 
12   }
13 }

和 ChewKeongTAN.BigInteger 类一样，Fcl35.Numberic.BigInteger 结构内部也使用一个 uint[] 来存储，同样也是基于 2³² ≈ 4.3 x 10⁹ 进制的。但是，最重要的一点，Fcl35.Numberic.BigInteger 结构是动态增长的，不需要事先指定最大允许有多大的数。

现在来讨论我最近写的一个 Skyiv.Numeric.BigInteger 类，她是从我在做 ACM 题时写的极简单的如下所示的 BigInteger 类发展而来的：

1 sealed class BigInteger
2 {
3   int[] digits = new int[1800];
4   // 
5 }

这个极简单的 BigInteger 类的全部源程序代码可以在本随笔开头给出的 URL 中找到，只有五十多行。她是基于 10 进制的，内部使用一个 int[] 来存储，需要事先指定该数组的大小，不能动态增长，而且只能表示非负整数。 

下面是 Skyiv.Numeric.BigInteger 的源程序代码，非常短小精悍，只有区区 335 行，该有的功能基本上都有了。实现了基本的算术运算和逻辑运算，还实现了 Pow 和 Sqrt 方法。

  1 using System;
  2 using System.Text;
  3 using System.Collections.Generic;
  4 
  5 namespace Skyiv.Numeric
  6 {
  7   sealed class BigInteger : IEquatable<BigInteger>, IComparable<BigInteger>
  8   {
  9     static readonly int Len = 9; // int.MaxValue = 2,147,483,647
 10     static readonly int Base = (int)Math.Pow(10, Len);
 11 
 12     int sign;
 13     List<int> data;
 14 

她内部使用一个 List<int> 来存储(程序中第 13 行)，是动态增长的，基 于 10⁹ 进制(第 9 和 10 行)。

 15     BigInteger(long x)
 16     {
 17       sign = (x == 0) ? 0 : ((x > 0) ? 1 : -1);
 18       data = new List<int>();
 19       for (ulong z = Abs(x); z != 0; z /= (ulong)Base) data.Add((int)(z % (ulong)Base));
 20     }
 21 
 22     BigInteger(BigInteger x)
 23     {
 24       sign = x.sign;  // x != null
 25       data = new List<int>(x.data);
 26     }
 27 
 28     public static implicit operator BigInteger(long x)
 29     {
 30       return new BigInteger(x);
 31     }
 32 

她只有两个构造函数，而且都是 private 的。要获得该类的实例，主要途径之一是通过从 long、int 等类型隐式转换(第 28 到 31 行)而来，例如 BigInteger a = 2; 语句。

 33     public static BigInteger Parse(string s)
 34     {
 35       if (s == null) return null;
 36       s = s.Trim().Replace(",", null);
 37       if (s.Length == 0) return 0;
 38       BigInteger z = 0;
 39       z.sign = (s[0] == '-') ? -1 : 1;
 40       if (s[0] == '-' || s[0] == '+') s = s.Substring(1);
 41       int r = s.Length % Len;
 42       z.data = new List<int>(new int[s.Length / Len + ((r != 0) ? 1 : 0)]);
 43       int i = z.data.Count - 1;
 44       if (r != 0) z.data[i--] = int.Parse(s.Substring(0, r));
 45       for (; i >= 0; i--, r += Len) z.data[i] = int.Parse(s.Substring(r, Len));
 46       z.Shrink();
 47       return z;
 48     }
 49 

获得该类的实例的主要途径之二是调用静态的 Parse 方法，例如：var a = BigInteger.Parse("-123,456,789"); 语句。

 50     public static void Swap(ref BigInteger x, ref BigInteger y)
 51     {
 52       BigInteger z = x;
 53       x = y;
 54       y = z;
 55     }
 56 
 57     public static ulong Abs(long x)
 58     {
 59       return (x < 0) ? (ulong)-x : (ulong)x;
 60     }
 61 
 62     public static BigInteger Abs(BigInteger x)
 63     {
 64       if (x == null) return null;
 65       BigInteger z = new BigInteger(x);
 66       z.sign = Math.Abs(x.sign);
 67       return z;
 68     }
 69 

上面三个静态的公共方法都很简单，主要是为了方便，如果不提供这些方法也没关系。要注意的是第二个 Abs(long) 方法返回 ulong 类型，而且从不抛出异常。而 Math.Abs(long) 返回 long 类型，而且当参数为 long.MinValue 时会抛出异常。

 70     public static BigInteger Pow(BigInteger x, int y)
 71     {
 72       if (x == null) return null;
 73       BigInteger z = 1, n = x;
 74       for (; y > 0; y >>= 1, n *= n) if ((y & 1) != 0) z *= n;
 75       return z;
 76     }
 77 
 78     public static BigInteger operator +(BigInteger x)
 79     {
 80       if (x == null) return null;
 81       return new BigInteger(x);
 82     }
 83 
 84     public static BigInteger operator -(BigInteger x)
 85     {
 86       if (x == null) return null;
 87       BigInteger z = new BigInteger(x);
 88       z.sign = -x.sign;
 89       return z;
 90     }
 91 

注意，重载的单目 + 操作符不仅是为了和重载的单目 - 操作符对应，而且可以用于克隆。 BigInteger x = +y; 和 BigInteger x = y; 语句是完全不同的，前者为克隆一份新的对象，后者只是直接引用同一对象。之所以不实现 IClone 接口，是因为该接口的 Clone 方法的返回类型是 object，使用起来很不方便：var x = y.Clone() as BigInteger; 需要一次强制类型转换。如果 FCL 中有泛型版本的 IClone<T> 接口就好了。

 92     public static BigInteger operator ++(BigInteger x)
 93     {
 94       return x + 1;
 95     }
 96 
 97     public static BigInteger operator --(BigInteger x)
 98     {
 99       return x - 1;
100     }
101 

在 C# 语言中重载 ++ 和 -- 运算符比在 C++ 语言中容易，C# 编译器会处理好前缀和后缀的情况。

102     public static BigInteger operator +(BigInteger x, BigInteger y)
103     {
104       if (x == null || y == null) return null;
105       if (x.AbsCompareTo(y) < 0) Swap(ref x, ref y);
106       BigInteger z = new BigInteger(x);
107       if (x.sign * y.sign == -1) z.AbsSubtract(y);
108       else z.AbsAdd(y);
109       return z;
110     }
111 
112     public static BigInteger operator -(BigInteger x, BigInteger y)
113     {
114       if (x == null || y == null) return null;
115       return x + (-y);
116     }
117 
118     public static BigInteger operator *(BigInteger x, BigInteger y)
119     {
120       if (x == null || y == null) return null;
121       BigInteger z = 0;
122       z.sign = x.sign * y.sign;
123       z.data = new List<int>(new int[x.data.Count + y.data.Count]);
124       for (int i = x.data.Count - 1; i >= 0; i--)
125         for (int j = y.data.Count - 1; j >= 0; j--)
126         {
127           long n = Math.BigMul(x.data[i], y.data[j]);
128           z.data[i + j] += (int)(n % Base);
129           z.CarryUp(i + j);
130           z.data[i + j + 1] += (int)(n / Base);
131           z.CarryUp(i + j + 1);
132         }
133       z.Shrink();
134       return z;
135     }
136 
137     public static BigInteger operator /(BigInteger dividend, BigInteger divisor)
138     {
139       BigInteger remainder;
140       return DivRem(dividend, divisor, out remainder);
141     }
142 
143     public static BigInteger operator %(BigInteger dividend, BigInteger divisor)
144     {
145       BigInteger remainder;
146       DivRem(dividend, divisor, out remainder);
147       return remainder;
148     }
149 

在 C# 中重载加减乘除等运算符，C# 编译器也会自动添上相应的赋值运算符(+=, -=, ...)。

150     public static BigInteger DivRem(BigInteger dividend, BigInteger divisor, out BigInteger remainder)
151     {
152       remainder = null;
153       if (dividend == null || divisor == null) return null;
154       if (divisor.sign == 0) throw new Exception("division by zero");
155       if (dividend.AbsCompareTo(divisor) < 0)
156       {
157         remainder = new BigInteger(dividend);
158         return 0;
159       }
160       remainder = 0;
161       BigInteger quotient = 0;
162       quotient.sign = dividend.sign * divisor.sign;
163       quotient.data = new List<int>(new int[dividend.data.Count]);
164       divisor = Abs(divisor); // NOT: divisor.sign = Math.Abs(divisor.sign);
165       for (int i = dividend.data.Count - 1; i >= 0; i--)
166       {
167         remainder = remainder * Base + dividend.data[i];
168         int iQuotient = remainder.BinarySearch(divisor, -1);
169         quotient.data[i] = iQuotient;
170         remainder -= divisor * iQuotient;
171       }
172       quotient.Shrink();
173       if (remainder.sign != 0) remainder.sign = dividend.sign;
174       return quotient;
175     }
176 
177     public static BigInteger Sqrt(BigInteger x)
178     {
179       if (x == null || x.sign < 0) return null;
180       BigInteger root = 0;
181       root.sign = 1;
182       root.data = new List<int>(new int[x.data.Count / 2 + 1]);
183       for (int i = root.data.Count - 1; i >= 0; i--) root.data[i] = x.BinarySearch(root, i);
184       root.Shrink();
185       return root;
186     }
187 
188     int BinarySearch(BigInteger divisor, int i)
189     {
190       int low = 0, high = Base - 1, mid = 0, cmp = 0;
191       while (low <= high)
192       {
193         mid = (low + high) / 2;
194         cmp = CompareTo(divisor, mid, i);
195         if (cmp > 0) low = mid + 1;
196         else if (cmp < 0) high = mid - 1;
197         else return mid;
198       }
199       return (cmp < 0) ? (mid - 1) : mid;
200     }
201 
202     int CompareTo(BigInteger divisor, int mid, int i)
203     {
204       if (i >= 0) divisor.data[i] = mid;
205       return AbsCompareTo(divisor * ((i >= 0) ? divisor : mid));
206     }
207 
208     void AbsAdd(BigInteger x)
209     {
210       for (int i = 0; i < x.data.Count; i++)
211       {
212         data[i] += x.data[i];
213         CarryUp(i);
214       }
215     }
216 
217     void AbsSubtract(BigInteger x)
218     {
219       for (int i = 0; i < x.data.Count; i++)
220       {
221         data[i] -= x.data[i];
222         CarryDown(i);
223       }
224       Shrink();
225     }
226 
227     void CarryUp(int n)
228     {
229       for (; data[n] >= Base; n++)
230       {
231         if (n == data.Count - 1) data.Add(data[n] / Base);
232         else data[n + 1] += data[n] / Base;
233         data[n] %= Base;
234       }
235     }
236 
237     void CarryDown(int n)
238     {
239       for (; data[n] < 0; n++)
240       {
241         data[n + 1]--;
242         data[n] += Base;
243       }
244       Shrink();
245     }
246 
247     void Shrink()
248     {
249       for (int i = data.Count - 1; i >= 0 && data[i] == 0; i--) data.RemoveAt(i);
250       if (data.Count == 0) sign = 0;
251     }
252 

除法、取模和求平方根都是通过二分搜索来进行的，每次搜索时都有进行一次乘法和比较运算，以决定下次要搜索的区间。

253     public static bool operator ==(BigInteger x, BigInteger y)
254     {
255       if (object.ReferenceEquals(x, null)) return object.ReferenceEquals(y, null);
256       return x.Equals(y);
257     }
258 
259     public static bool operator !=(BigInteger x, BigInteger y)
260     {
261       if (object.ReferenceEquals(x, null)) return !object.ReferenceEquals(y, null);
262       return !x.Equals(y);
263     }
264 
265     public static bool operator <(BigInteger x, BigInteger y)
266     {
267       if (object.ReferenceEquals(x, null)) return !object.ReferenceEquals(y, null);
268       return x.CompareTo(y) < 0;
269     }
270 
271     public static bool operator >(BigInteger x, BigInteger y)
272     {
273       if (object.ReferenceEquals(x, null)) return false;
274       return x.CompareTo(y) > 0;
275     }
276 
277     public static bool operator <=(BigInteger x, BigInteger y)
278     {
279       if (object.ReferenceEquals(x, null)) return true;
280       return x.CompareTo(y) <= 0;
281     }
282 
283     public static bool operator >=(BigInteger x, BigInteger y)
284     {
285       if (object.ReferenceEquals(x, null)) return object.ReferenceEquals(y, null);
286       return x.CompareTo(y) >= 0;
287     }
288 
289     public override string ToString()
290     {
291       StringBuilder sb = new StringBuilder();
292       if (sign < 0) sb.Append('-');
293       sb.Append((data.Count == 0) ? 0 : data[data.Count - 1]);
294       for (int i = data.Count - 2; i >= 0; i--) sb.Append(data[i].ToString("D" + Len));
295       return sb.ToString();
296     }
297 
298     public override int GetHashCode()
299     {
300       int hash = sign;
301       foreach (int n in data) hash ^= n;
302       return hash;
303     }
304 
305     public override bool Equals(object other)
306     {
307       if (other == null || GetType() != other.GetType()) return false;
308       return Equals(other as BigInteger);
309     }
310 
311     public bool Equals(BigInteger other)
312     {
313       return CompareTo(other) == 0;
314     }
315 
316     public int CompareTo(BigInteger other)
317     {
318       if (object.ReferenceEquals(other, null)) return 1;
319       if (sign < other.sign) return -1;
320       if (sign > other.sign) return 1;
321       if (sign == 0) return 0;
322       return sign * AbsCompareTo(other);
323     }
324 
325     int AbsCompareTo(BigInteger other)
326     {
327       if (data.Count < other.data.Count) return -1;
328       if (data.Count > other.data.Count) return 1;
329       for (int i = data.Count - 1; i >= 0; i--)
330         if (data[i] != other.data[i])
331           return (data[i] < other.data[i]) ? -1 : 1;
332       return 0;
333     }
334   }
335 }

最后就是各种逻辑运算符以及从基类继承的 ToString、GetHashCode、Equals 方法，以及用于实现接口的 Equals 和 CompareTo 方法。

下面，让我们来比较这四种 BigInteger 的效率吧。测试程序如下所示：

TestMain.cs：

 1 using System;
 2 
 3 namespace Skyiv
 4 {
 5   class TestMain
 6   {
 7     static void Main()
 8     {
 9       string s1 = new TestFcl35BigInteger().Run(20);
10       string s2 = new TestSkyivBigInteger().Run(20);
11       string s3 = new TestChewBigInteger().Run(20);
12       string s4 = new TestJavaBigInteger().Run(14);
13       Console.WriteLine(s1 == s2 && s1 == s3);
14     }
15   }
16 }

TestBase.cs:

 1 using System;
 2 using System.Diagnostics;
 3 
 4 namespace Skyiv
 5 {
 6   class TestBase
 7   {
 8     public string Run(int n)
 9     {
10       var stopWatch = Stopwatch.StartNew();
11       var s = RunTest(n);
12       stopWatch.Stop();
13       Console.WriteLine("{0,-25} {1,11:F7} {2,7:N0}",
14         GetType(), stopWatch.Elapsed.TotalSeconds, s.Length);
15       return s;
16     }
17 
18     protected virtual string RunTest(int n)
19     {
20       throw new InvalidOperationException("TestBase.RunTest()");
21     }
22   }
23 }

Test?????BigInteger.cs (以 TestSkyivBigInteger 为例，其它的大同小异):

 1 using System;
 2 using BigInteger = Skyiv.Numeric.BigInteger;
 3 
 4 namespace Skyiv
 5 {
 6   sealed class TestSkyivBigInteger : TestBase
 7   {
 8     protected sealed override string RunTest(int n)
 9     {
10       BigInteger a = 2, b;
11       for (int i = 2; i <= n; i++, a *= b) b = a + i;
12       return a.ToString();
13     }
14   }
15 }

这个测试程序主要是测试 BigInteger 的加法和乘法运算，其运行结果如下：

Skyiv.TestFcl35BigInteger  25.8100647 256,142
Skyiv.TestSkyivBigInteger  34.2497761 256,142
Skyiv.TestChewBigInteger  396.8586529 256,142
Skyiv.TestJavaBigInteger   61.9083458   4,003
True

上表中，第二栏是测试程序的运行时间，单位是秒。第三栏是最终运算结果的十进制数字的个数。前三个 BigInteger 都使用了同一测试用例，运算结果达 256,142 个十进制数字，并且经过比较这三个 BigInteger 的运算结果是相同的。最后的 java.math.BigInteger 由于运算速度极慢(估计是微软 J# 的问题，我想在真正的 java 语言中应该是比较快的)，测试用例的规模更小，运算结果只有 4,003 个十进制数字。从测试结果可以看出：

• Fcl35.Numeric.BigInteger 最快，推荐使用。
• Skyiv.Numeric.BigInteger 也很快，并且源程序代码短小精悍，也推荐使用。
• ChewKeongTAN.BigInteger 很慢，并且不能动态增长，不推荐使用。
• java.math.BigInteger 极慢，并且需要 J# run-time library，也不推荐使用。

Skyiv.Numeric.BigInteger 和 ChewKeongTAN.BigInteger 的乘法都是使用简单的乘法，时间复杂度是 O(N²)。Fcl35.Numeric.BigInteger 的乘法是使用 Karatsuba 算法，其时间复杂度约为 O(N^1.585)，如下所示：

 1 public static BigInteger Multiply(BigInteger x, BigInteger y)
 2 {
 3   int xLength = x.Length;
 4   int yLength = y.Length;
 5   if ((((xLength + yLength) >= UpperBoundForSchoolBookMultiplicationDigits /* 0x40 */)
 6     && (xLength >= ForceSchoolBookMultiplicationThresholdDigits /* 8 */))
 7     && (yLength >= ForceSchoolBookMultiplicationThresholdDigits /* 8 */))
 8   {
 9     return MultiplyKaratsuba(x, y);
10   }
11   return MultiplySchoolBook(x, y);
12 }
13 
14 private static BigInteger MultiplyKaratsuba(BigInteger x, BigInteger y)
15 {
16   int length = Math.Max(x.Length, y.Length) / 2;
17   if (length <= 2 * ForceSchoolBookMultiplicationThresholdDigits /* 0x10 */
18     || x.Length < 2 * ForceSchoolBookMultiplicationThresholdDigits /* 0x10 */
19     || y.Length < 2 * ForceSchoolBookMultiplicationThresholdDigits /* 0x10 */)
20     return MultiplySchoolBook(x, y);
21   int shift = BitsPerDigit /* 0x20 */ * length;
22   BigInteger i1 = RightShift(x, shift);
23   BigInteger i2 = x.RestrictTo(length);
24   BigInteger i3 = RightShift(y, shift);
25   BigInteger i4 = y.RestrictTo(length);
26   BigInteger i5 = Multiply(i1, i3);
27   BigInteger i6 = Multiply(i2, i4);
28   BigInteger i7 = Multiply(i1 + i2, i3 + i4) - (i5 + i6);
29   return (i6 + LeftShift(i7 + LeftShift(i5, shift), shift));
30 }

其实如果使用快速傅里叶变换(FFT)来进行乘法运算，时间复杂度可降低到 O(N logN loglogN)。而且除法可以用除数的倒数乘以被除数来计算，倒数值用牛顿迭代法：

U_i+1 = U_i (2 - VU_i)

来计算，这导致 U_∞二次收敛于 1/V。实际上，许多超级计算机和 RISC 机都是使用这种迭代法来实现除法的。用牛顿迭代法计算平方根和除法类似。若：

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

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

我可能会使用快速傅里叶变换来改写 Skyiv.Numeric.BigInteger 类。微软声称在下一版本的 .NET Framework 中会重新发布 System.Numeric.BigInteger，希望也使用快速傅里叶变换。:)

本文中的所有源代码可以到 http://bitbucket.org/ben.skyiv/biginteger/downloads 页面下载。
也可以使用 hg clone http://bitbucket.org/ben.skyiv/biginteger/ 命令下载。
关于 hg ，请参阅 Mercurial 备忘录。

---

1. 2008-07-25 更新: 请参阅：再谈 BigInteger - 使用快速傅里叶变换。
2. 2010-04-22 更新: 对 Skyiv.Numeric.BigInteger 的 GetHashCode、Parse 和 ToString 方法进行了优化。请参阅：再谈 BigInteger - 优化。