我在“浅谈 BigInteger”的随笔中实现了一个 Skyiv.Numeric.BigInteger 类，那时乘法是使用常规的 O(N2) 的算法，所以比 .NET Framework 3.5 Base Class Library 中的 System.Numeric.BigInteger 类稍慢，后者的乘法是使用 Karatsuba 算法，其时间复杂度约为 O(N1.585)。
现在，我已经实现了一个提供任意精度的算术运算的静态类: BigArithmetic，该类使用快速傅里叶变换计算大整数乘法，其时间复杂度降低到 O(N logN loglogN)。所以，Skyiv.Numeric.BigInteger 类也重新改写调用 BigArithmetic 类的静态方法来进行算术运算。
下面就是改写后的 Skyiv.Numeric.BigInteger 类的源程序代码：
2 using System.Text;
3
4 namespace Skyiv.Numeric
5 {
6 sealed class BigInteger : IEquatable<BigInteger>, IComparable<BigInteger>
7 {
8 static readonly byte Len = 2;
9 static readonly byte Base = (byte)Math.Pow(10, Len);
10
11 sbyte sign; // 符号，取值：-1, 0, 1。
12 byte[] data; // 字节数组以 100 为基，字节数组中第一个元素存储的数字是最高有效位。
13
因为 BigArithmetic 类是对字节数组进行算术运算，字节数组以 100 为基，为了方便调用 BigArithmetic 类的静态方法进行算术运算，BigInteger 也改为使用以 100 为基的字节数组来存储。原来是使用 109 为基的 List<int> 来存储的。
我在改写过程中曾经使用过以 108 为基的 List<int> 来存储，在需要调用 BigArithmetic 类的静态方法进行算术运算时将参数转换为以 100 为基的字节数组，运算完成后再将运算结果转换回来。这样做的好处是加法和减法运算比较快(还使用原来的方法，不调用 BigArithmetic 类的静态方法)，并且除了 operator *、DivRem 和 Sqrt 方法以外，其他所有的方法都不用改写。不过经过综合考虑，还是采用现在的方案。
15 {
16 }
17
18 BigInteger(long x)
19 {
20 sign = (sbyte)((x == 0) ? 0 : ((x > 0) ? 1 : -1));
21 data = new byte[10]; // long.MinValue = -9,223,372,036,854,775,808
22 ulong z = (x < 0) ? (ulong)-x : (ulong)x;
23 for (int i = data.Length - 1; z != 0; i--, z /= Base) data[i] = (byte)(z % Base);
24 Shrink();
25 }
26
27 BigInteger(BigInteger x)
28 {
29 sign = x.sign; // x != null
30 data = new byte[x.data.Length];
31 Array.Copy(x.data, data, data.Length);
32 }
33
34 public static implicit operator BigInteger(long x)
35 {
36 return new BigInteger(x);
37 }
38
私有的构造函数，和原来的大同小异。
40 {
41 if (s == null) return null;
42 s = s.Trim().Replace(",", null);
43 if (s.Length == 0) return 0;
44 BigInteger z = new BigInteger();
45 z.sign = (sbyte)((s[0] == '-') ? -1 : 1);
46 if (s[0] == '-' || s[0] == '+') s = s.Substring(1);
47 int r = s.Length % Len;
48 z.data = new byte[s.Length / Len + ((r != 0) ? 1 : 0)];
49 int i = 0;
50 if (r != 0) z.data[i++] = byte.Parse(s.Substring(0, r));
51 for (; i < z.data.Length; i++, r += Len) z.data[i] = byte.Parse(s.Substring(r, Len));
52 z.Shrink();
53 return z;
54 }
55
静态的 Parse 方法，也和原来的大同小异。
程序的第 56 到 95 行以及 111 到 116 行是 Abs、Pow 方法以及单目 +、-、++ 和 -- 运算符，以及减法(-)运算符，和原来的相同。
97 {
98 if (x == null || y == null) return null;
99 if (x.AbsCompareTo(y) < 0) Utility.Swap(ref x, ref y);
100 BigInteger z = new BigInteger();
101 z.sign = x.sign;
102 byte[] bs = Utility.Expand(y.data, x.data.Length);
103 bool isAdd = x.sign * y.sign == 1;
104 z.data = new byte[x.data.Length + (isAdd ? 1 : 0)];
105 if (isAdd) BigArithmetic.Add(z.data, x.data, bs, bs.Length);
106 else BigArithmetic.Subtract(z.data, x.data, bs, bs.Length);
107 z.Shrink();
108 return z;
109 }
110
加法(+)运算符，调用 BigArithmetic 类的 Add 和 Subtract 方法。
118 {
119 if (x == null || y == null) return null;
120 if (x.sign * y.sign == 0) return 0;
121 BigInteger z = new BigInteger();
122 z.sign = (sbyte)(x.sign * y.sign);
123 z.data = new byte[x.data.Length + y.data.Length];
124 BigArithmetic.Multiply(z.data, x.data, x.data.Length, y.data, y.data.Length);
125 z.Shrink();
126 return z;
127 }
128
乘法(*)运算，直接调用 BigArithmetic 类的 Multiply 方法。
程序的第 129 到 141 行是 / 和 % 运算符，和原来的相同。
143 {
144 remainder = null;
145 if (dividend == null || divisor == null) return null;
146 if (divisor.sign == 0) throw new DivideByZeroException();
147 if (dividend.AbsCompareTo(divisor) < 0)
148 {
149 remainder = new BigInteger(dividend);
150 return 0;
151 }
152 BigInteger quotient = new BigInteger();
153 remainder = new BigInteger();
154 quotient.data = new byte[dividend.data.Length - divisor.data.Length + 1];
155 remainder.data = new byte[divisor.data.Length];
156 BigArithmetic.DivRem(quotient.data, remainder.data, dividend.data,
157 dividend.data.Length, divisor.data, divisor.data.Length);
158 quotient.sign = (sbyte)(dividend.sign * divisor.sign);
159 remainder.sign = dividend.sign;
160 quotient.Shrink();
161 remainder.Shrink();
162 return quotient;
163 }
164
静态的 DevRem 方法，也是调用 BigArithmetic 类的 DivRem 方法。
166 {
167 if (x == null || x.sign < 0) return null;
168 if (x.sign == 0) return 0;
169 if (x.data.Length == 1) return new BigInteger((long)Math.Sqrt(x.data[0]));
170 BigInteger z = new BigInteger();
171 z.sign = 1;
172 z.data = new byte[x.data.Length / 2 + 3];
173 z.data = Adjust(BigArithmetic.Sqrt(z.data, z.data, z.data.Length, x.data, x.data.Length), x.data.Length);
174 z.Shrink();
175 BigInteger z1 = z + 1; // 平方根有可能比实际小 1。
176 return (z1 * z1 <= x) ? z1 : z;
177 }
178
179 static byte[] Adjust(byte[] bs, int digits)
180 {
181 if (bs[0] >= 10) throw new OverflowException("sqrt adjust");
182 byte[] nbs = new byte[(digits + 1) / 2];
183 if (digits % 2 == 0)
184 for (int k = bs[0], i = 0; i < nbs.Length; i++, k = bs[i] % 10)
185 nbs[i] = (byte)(k * 10 + bs[i + 1] / 10);
186 else Array.Copy(bs, nbs, nbs.Length);
187 return nbs;
188 }
189
求平方根(Sqrt)，调用 BigArithmetic 类的 Sqrt 方法。但是要注意，BigArithmetic 类的 Sqrt 方法求得的平方根是小数，需要使用 Adjust 方法调整为整数。并且平方根有可能比实际的小 1，也需要判断和调整。
191 {
192 int i;
193 for (i = 0; i < data.Length; i++) if (data[i] != 0) break;
194 if (i != 0)
195 {
196 byte[] bs = new byte[data.Length - i];
197 Array.Copy(data, i, bs, 0, bs.Length);
198 data = bs;
199 }
200 if (data.Length == 0) sign = 0;
201 }
202
私有的 Shrink 方法用于在需要时收缩存储数据的字节数组。
程序的第 203 到 238 行是各种逻辑运算符，和原来的完全一样。
240 {
241 StringBuilder sb = new StringBuilder();
242 if (sign < 0) sb.Append('-');
243 sb.Append((data.Length == 0) ? 0 : (int)data[0]);
244 for (int i = 1; i < data.Length; i++) sb.Append(data[i].ToString("D" + Len));
245 return sb.ToString();
246 }
247
从基类继承的 ToString 方法，和原来的大同小异。
程序的第 248 到 274 行是从基类继承的 GetHashCode、Equals 方法，以及用于实现接口的 Equals 和 CompareTo 方法，和原来的相同。
276 {
277 if (data.Length < other.data.Length) return -1;
278 if (data.Length > other.data.Length) return 1;
279 return BigArithmetic.Compare(data, other.data, data.Length);
280 }
281 }
282 }
最后是私有的 AbsCompareTo 方法，调用 BigArithmetic 类的 Compare 方法。
下面是程序中用到的辅助类：
2
3 namespace Skyiv
4 {
5 static class Utility
6 {
7 public static T[] Expand<T>(T[] x, int n)
8 {
9 T[] z = new T[n]; // assume n >= x.Length
10 Array.Copy(x, 0, z, n - x.Length, x.Length);
11 return z;
12 }
13
14 public static void Swap<T>(ref T x, ref T y)
15 {
16 T z = x;
17 x = y;
18 y = z;
19 }
20 }
21 }
经过运行测试程序，发现计算 256,142 位的乘积的运行时间从原来的 34.2497761 秒降低到 0.1749114 秒，速度提高了将近二百倍。如果计算的数字更大的话，速度将提高得更多。因为原来乘法的时间复杂度是 O(N2)。而现在乘法使用 FFT，时间复杂度是 O(N logN loglogN)。
我想，这个 BigInteger 类的运算速度已经是非常快的了。不知道有没有其他 C# 的 BigInteger 类的运算速度更快。
本文中的所有源代码可以到 https://bitbucket.org/ben.skyiv/biginteger 页面下载。
也可以使用 hg clone http://bitbucket.org/ben.skyiv/biginteger/ 命令下载。
关于 hg ，请参阅 Mercurial 备忘录。
2010-04-22 更新: 对 Skyiv.Numeric.BigInteger 的 GetHashCode、Parse 和 ToString 方法进行了优化。请参阅：再谈 BigInteger - 优化
