再谈 BigInteger - 使用快速傅里叶变换难点.doc

再谈 BigInteger - 使用快速傅里叶变换我在“浅谈 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 类的源程序代码：复制代码 1 using System; 2 using System.Text; 3 4 namespace Skyiv.Numeric 5 6 sealed class BigInteger : IEquatable, IComparable 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 来存储的。我在改写过程中曾经使用过以 108 为基的 List 来存储，在需要调用 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 byte10; / 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) datai = (byte)(z % Base); 24 Shrink(); 25 26 27 BigInteger(BigInteger x) 28 29 sign = x.sign; / x != null 30 data = new bytex.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)(s0 = -) ? -1 : 1); 46 if (s0 = - | s0 = +) s = s.Substring(1); 47 int r = s.Length % Len; 48 z.data = new bytes.Length / Len + (r != 0) ? 1 : 0); 49 int i = 0; 50 if (r != 0) z.datai+ = byte.Parse(s.Substring(0, r); 51 for (; i z.data.Length; i+, r += Len) z.datai = 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);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 bytex.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 方法。复制代码117 public static BigInteger operator *(BigInteger x, BigInteger y)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 bytex.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 行是 / 和 % 运算符，和原来的相同。 复制代码142 public static BigInteger DivRem(BigInteger dividend, BigInteger divisor, out BigInteger remainder)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 bytedividend.data.Length - divisor.data.Length + 1;155 remainder.data = new 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 方法。 复制代码165 public static BigInteger Sqrt(BigInteger x)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.data0);170 BigInteger z = new BigInteger();171 z.sign = 1;172 z.data = new bytex.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 = 10) throw new OverflowException(sqrt adjust);182 byte nbs = new byte(digits + 1) / 2;183 if (digits % 2 = 0)184 for (int k = bs0, i = 0; i nbs.Length; i+, k = bsi % 10)185 nbsi = (byte)(k * 10 + bsi + 1 / 10);186 else Array.Copy(bs, nbs, nbs.Length);187 return nbs;188 189 复制代码求平方根(Sqrt)，调用 BigArithmetic 类的 Sqrt 方法。但是要注意，BigArithmetic 类的 Sqrt 方法求得的平方根是小数，需要使用 Adjust 方法调整为整数。并且平方根有可能比实际的小 1，也需要判断和调整。 复制代码190 void Shrink()191 192 int i;193 for (i = 0; i data.Length; i+) if (datai != 0) break;194 if (i != 0)195 196 byte bs = new bytedata.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 行是各种逻辑运算符，和原来的完全一样。复制代码239 public override string ToString()240 241 StringBuilder sb = new StringBuilder();242 if (sign 0) Append(-);243 sb.Append(data.Length = 0) ? 0 : (int)data0);244 for (int i = 1; i data.Length; i+) sb.Append(datai.ToString(D + Len);245 return sb.ToString();246 247 复制代码从基类继承的 ToString 方法，和原来的大同小异。程序的第 248 到 274 行是从基类继承的 GetHashCode、Equals 方法，以及用于实现接口的 Equals 和 CompareTo 方法，和原来的相同。 复制代码275 int AbsCompareTo(BigInteger other)276 277 if (data.Length other.data.Length) return 1;279