再谈 BigInteger - 使用快速傅里叶变换

我在“浅谈 BigInteger”的随笔中实现了一个 Skyiv.Numeric.BigInteger 类，那时乘法是使用常规的 O(N²) 的算法，所以比 .NET Framework 3.5 Base Class Library 中的 System.Numeric.BigInteger 类稍慢，后者的乘法是使用 Karatsuba 算法，其时间复杂度约为 O(N^1.585)。

现在，我已经实现了一个提供任意精度的算术运算的静态类: BigArithmetic，该类使用快速傅里叶变换计算大整数乘法，其时间复杂度降低到 O(N logN loglogN)。所以，Skyiv.Numeric.BigInteger 类也重新改写调用 BigArithmetic 类的静态方法来进行算术运算。

下面就是改写后的 Skyiv.Numeric.BigInteger 类的源程序代码：

  1 using System;
  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 为基的字节数组来存储。原来是使用 10⁹ 为基的 List<int> 来存储的。

我在改写过程中曾经使用过以 10⁸ 为基的 List<int> 来存储，在需要调用 BigArithmetic 类的静态方法进行算术运算时将参数转换为以 100 为基的字节数组，运算完成后再将运算结果转换回来。这样做的好处是加法和减法运算比较快(还使用原来的方法，不调用 BigArithmetic 类的静态方法)，并且除了 operator *、DivRem 和 Sqrt 方法以外，其他所有的方法都不用改写。不过经过综合考虑，还是采用现在的方案。

 

 14     BigInteger()
 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 

私有的构造函数，和原来的大同小异。

 

 39     public static BigInteger Parse(string s)
 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 方法以及单目 +、-、++ 和 -- 运算符，以及减法(-)运算符，和原来的相同。

 

 96     public static BigInteger operator +(BigInteger x, BigInteger y)
 97     {
 98       if (x == null || y == null) return null;
 99       if (x.AbsCompareTo(y) < 0) Utility.Swap(ref x, ref y);
       BigInteger z = new BigInteger();
       z.sign = x.sign;
       byte[] bs = Utility.Expand(y.data, x.data.Length);
       bool isAdd = x.sign * y.sign == 1;
      z.data = new byte[x.data.Length + (isAdd ? 1 : 0)];
      if (isAdd) BigArithmetic.Add(z.data, x.data, bs, bs.Length);
      else BigArithmetic.Subtract(z.data, x.data, bs, bs.Length);
      z.Shrink();
      return z;
    }

加法(+)运算符，调用 BigArithmetic 类的 Add 和 Subtract 方法。

     public static BigInteger operator *(BigInteger x, BigInteger y)
     {
       if (x == null || y == null) return null;
       if (x.sign * y.sign == 0) return 0;
       BigInteger z = new BigInteger();
       z.sign = (sbyte)(x.sign * y.sign);
       z.data = new byte[x.data.Length + y.data.Length];
       BigArithmetic.Multiply(z.data, x.data, x.data.Length, y.data, y.data.Length);
       z.Shrink();
       return z;
     }
 

乘法(*)运算，直接调用 BigArithmetic 类的 Multiply 方法。

程序的第 129 到 141 行是 / 和 % 运算符，和原来的相同。

 

     public static BigInteger DivRem(BigInteger dividend, BigInteger divisor, out BigInteger remainder)
     {
       remainder = null;
       if (dividend == null || divisor == null) return null;
       if (divisor.sign == 0) throw new DivideByZeroException();
       if (dividend.AbsCompareTo(divisor) < 0)
       {
         remainder = new BigInteger(dividend);
         return 0;
       }
       BigInteger quotient = new BigInteger();
       remainder = new BigInteger();
       quotient.data = new byte[dividend.data.Length - divisor.data.Length + 1];
       remainder.data = new byte[divisor.data.Length];
       BigArithmetic.DivRem(quotient.data, remainder.data, dividend.data,
         dividend.data.Length, divisor.data, divisor.data.Length);
       quotient.sign = (sbyte)(dividend.sign * divisor.sign);
       remainder.sign = dividend.sign;
       quotient.Shrink();
       remainder.Shrink();
       return quotient;
     }
 

静态的 DevRem 方法，也是调用 BigArithmetic 类的 DivRem 方法。

 

     public static BigInteger Sqrt(BigInteger x)
     {
       if (x == null || x.sign < 0) return null;
       if (x.sign == 0) return 0;
       if (x.data.Length == 1) return new BigInteger((long)Math.Sqrt(x.data[0]));
       BigInteger z = new BigInteger();
       z.sign = 1;
       z.data = new byte[x.data.Length / 2 + 3];
       z.data = Adjust(BigArithmetic.Sqrt(z.data, z.data, z.data.Length, x.data, x.data.Length), x.data.Length);
       z.Shrink();
       BigInteger z1 = z + 1; // 平方根有可能比实际小 1。
       return (z1 * z1 <= x) ? z1 : z;
     }
 
     static byte[] Adjust(byte[] bs, int digits)
     {
       if (bs[0] >= 10) throw new OverflowException("sqrt adjust");
       byte[] nbs = new byte[(digits + 1) / 2];
       if (digits % 2 == 0)
         for (int k = bs[0], i = 0; i < nbs.Length; i++, k = bs[i] % 10)
           nbs[i] = (byte)(k * 10 + bs[i + 1] / 10);
       else Array.Copy(bs, nbs, nbs.Length);
       return nbs;
     }
 

求平方根(Sqrt)，调用 BigArithmetic 类的 Sqrt 方法。但是要注意，BigArithmetic 类的 Sqrt 方法求得的平方根是小数，需要使用 Adjust 方法调整为整数。并且平方根有可能比实际的小 1，也需要判断和调整。

 

     void Shrink()
     {
       int i;
       for (i = 0; i < data.Length; i++) if (data[i] != 0) break;
       if (i != 0)
       {
         byte[] bs = new byte[data.Length - i];
         Array.Copy(data, i, bs, 0, bs.Length);
         data = bs;
       }
       if (data.Length == 0) sign = 0;
     }
 

私有的 Shrink 方法用于在需要时收缩存储数据的字节数组。

程序的第 203 到 238 行是各种逻辑运算符，和原来的完全一样。

     public override string ToString()
     {
       StringBuilder sb = new StringBuilder();
       if (sign < 0) sb.Append('-');
       sb.Append((data.Length == 0) ? 0 : (int)data[0]);
       for (int i = 1; i < data.Length; i++) sb.Append(data[i].ToString("D" + Len));
       return sb.ToString();
     }
 

从基类继承的 ToString 方法，和原来的大同小异。

程序的第 248 到 274 行是从基类继承的 GetHashCode、Equals 方法，以及用于实现接口的 Equals 和 CompareTo 方法，和原来的相同。

 

     int AbsCompareTo(BigInteger other)
     {
       if (data.Length < other.data.Length) return -1;
       if (data.Length > other.data.Length) return 1;
       return BigArithmetic.Compare(data, other.data, data.Length);
     }
   }
 }

最后是私有的 AbsCompareTo 方法，调用 BigArithmetic 类的 Compare 方法。

下面是程序中用到的辅助类：

 1 using System;
 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
       Array.Copy(x, 0, z, n - x.Length, x.Length);
       return z;
     }
 
     public static void Swap<T>(ref T x, ref T y)
     {
       T z = x;
       x = y;
       y = z;
     }
   }
 }

 

经过运行测试程序，发现计算 256,142 位的乘积的运行时间从原来的 34.2497761 秒降低到 0.1749114 秒，速度提高了将近二百倍。如果计算的数字更大的话，速度将提高得更多。因为原来乘法的时间复杂度是 O(N²)。而现在乘法使用 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 - 优化