[转]浅谈 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 类。结果如下表所示：

在上表中，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? 谈到这件事，要点是：

于是， .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 方法。

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

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

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

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

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

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

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

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

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

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

TestMain.cs：

TestBase.cs:

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

这个测试程序主要是测试 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)，如下所示：

其实如果使用快速傅里叶变换(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 - 优化。