使用快速傅里叶变换计算大整数乘法-代码

在上一篇随笔“使用快速傅里叶变换计算大整数乘法”中，已经讲述了使用快速傅里叶变换计算大整数乘法的原理。在这一篇随笔中，我们就使用快速傅里叶变换来实现一个提供任意精度的算术运算的静态类：BigArithmetic。

下面就是 BigArithmetic 类源程序代码： 

  1 using System;
  2 using System.Diagnostics;
  3 
  4 namespace Skyiv.Numeric
  5 {
  6   /// <summary>
  7   /// 提供任意精度的算术运算的静态类。使用快速傅里叶变换。
  8   /// 本类对字节数组进行算术运算，字节数组以 100 为基。
  9   /// 字节数组中第一个元素存储的数字是最高有效位。
 10   /// </summary>
 11   static class BigArithmetic
 12   {
 13     //= C语言数值算法程序大全(第二版)，ISBN 7-5053-2931-6 / TP 993
 14     //= Numerical Recipes in C, The Art of Scientific Computing, Second Edition
 15     //= Cambridge University Press 1988, 1992
 16     //= [美] W.H.Press, S.A.Teukolsky, W.T.Vetterling, B.P.Flannery 著
 17     //= 傅祖芸 赵梅娜 丁岩 等译，傅祖芸 校，电子工业出版社，1995年10月第一版
 18 
 19     static readonly byte Len = 2; // 字节数组的元素包含的十进制数字的个数
 20     static readonly byte Base = (byte)Math.Pow(10, Len); // 字节数组的基
 21     static readonly byte MaxValue = (byte)(Base - 1);    // 字节数组的元素的最大值
 22 
 23     //= pp.431-432, four1, 12.2 快速傅里叶变换(FFT)
 24     /// <summary>
 25     /// 复函数的快速傅里叶变换
 26     /// 变量 nn 是复数据点的个数，实型数组 data[1..2*nn] 的实际界长是两倍 nn，
 27     /// 而每个复数值占据了两个相继的存储单元。nn 必须是 2 的整数幂
 28     /// </summary>
 29     /// <param name="data">实型数组 data[1..2*nn]。注意，下标从 1 开始</param>
 30     /// <param name="isInverse">是否逆变换。注意: 逆变换未乘上归一化因子 1/nn</param>
 31     public static void ComplexFFT(double[] data, bool isInverse)
 32     {
 33       int n = data.Length - 1; // n 必须是 2 的正整数幂
 34       int nn = n >> 1;         // 变量 nn 是复数据点的个数
 35       for (int i = 1, j = 1; i < n; i += 2) // 这个循环实现位序颠倒
 36       {
 37         if (j > i)
 38         {
 39           Utility.Swap(ref data[j], ref data[i]);
 40           Utility.Swap(ref data[j + 1], ref data[i + 1]);
 41         }
 42         int m = nn;
 43         for (; m >= 2 && j > m; m >>= 1) j -= m;
 44         j += m;
 45       }
 46       for (int mmax = 2, istep = 4; n > mmax; mmax = istep) // 执行 log2(nn) 次外循环
 47       {
 48         istep = mmax << 1; // 下面是关于三角递归的初始赋值
 49         double theta = (isInverse ? -2 : 2) * Math.PI / mmax;
 50         double wtemp = Math.Sin(0.5 * theta);
 51         double wpr = -2 * wtemp * wtemp;
 52         double wpi = Math.Sin(theta);
 53         double wr = 1;
 54         double wi = 0;
 55         for (int m = 1; m < mmax; m += 2)
 56         {
 57           for (int i = m; i <= n; i += istep)
 58           {
 59             int j = i + mmax; // 下面是 Danielson-Lanczos 公式
 60             double tempr = wr * data[j] - wi * data[j + 1];
 61             double tempi = wr * data[j + 1] + wi * data[j];
 62             data[j] = data[i] - tempr;
 63             data[j + 1] = data[i + 1] - tempi;
 64             data[i] += tempr;
 65             data[i + 1] += tempi;
 66           }
 67           wr = (wtemp = wr) * wpr - wi * wpi + wr; // 三角递归
 68           wi = wi * wpr + wtemp * wpi + wi;
 69         }
 70       }
 71     }
 72 

该类的第一个方法是 ComplexFFT，计算复函数的快速傅里叶变换。注意，ComplexFFT 并没有使用复数(不象 C++，C# 也没有提供复数)，而是让每个复数值占据实型数组的两个相继的存储单元。还有，要求输入的复数据点的个数必须是 2 的整数幂。该方法也能够计算复函数的快速傅里叶逆变换。

该程序的算法是使用 1942 年 Danielson 和 Lanczos 证明的引理：一个界长为 N 的离散傅里叶变换可以重新写成两个界长各为 N/2 的离散傅里叶变换之和。在算法的第一部分，将数据整理成位序颠倒的次序。而在第二部分，有一个执行 log₂N 次的外循环。它依次计算界长为 2, 4, 8, ..., N 的变换。对于这一过程的每一步来说，为了履行 Danielson-Lanczos 引理，有两个嵌套的内循环，其涉及到已计算的子变换和每个变换的元素。通过限制外部调用正弦和余弦到外层循环，可以使运算更有效，在外层循环中只要调用它们 log₂N 次。倍角的正弦和余弦的计算是通过内循环中简单的递归关系进行的，如下所示：

cos(θ + δ) = cosθ - [ α cosθ  +  βsinθ ]
sin(θ + δ) = sinθ - [ α sinθ  -  βcosθ ]

其中 α, β 是预先计算的系数：α = 2 sin²(δ/2)，  β = sinδ 。

 

 73     //= pp.436, realft, 12.3.2 单个实函数的 FFT
 74     /// <summary>
 75     /// 单个实函数的快速傅里叶变换
 76     /// 计算一组 n 个实值数据点的傅里叶变换。用复傅里叶变换的正半频率替换这些数据，
 77     /// 它存储在数组 data[1..n] 中。复变换的第一个和最后一个分量的实数值分别返回
 78     /// 单元 data[1] 和 data[2] 中。n 必须是 2 的幂次。这个程序也能计算复数据数组
 79     /// 的逆变换，只要该数组是实值数据的变换(在这种情况下，其结果必须乘以 1/n)即可。
 80     /// </summary>
 81     /// <param name="data">实型数组 data[1..n]。注意，下标从 1 开始</param>
 82     /// <param name="isInverse">是否逆变换。注意: 逆变换未乘上归一化因子 1/n</param>
 83     public static void RealFFT(double[] data, bool isInverse)
 84     {
 85       int n = data.Length - 1; // n 必须是 2 的整数幂
 86       if (!isInverse) ComplexFFT(data, isInverse); // 此处是正向变换
 87       double theta = (isInverse ? -2 : 2) * Math.PI / n; // 递归的初始赋值
 88       double wtemp = Math.Sin(0.5 * theta);
 89       double wpr = -2 * wtemp * wtemp;
 90       double wpi = Math.Sin(theta);
 91       double wr = 1 + wpr;
 92       double wi = wpi;
 93       double c1 = 0.5;
 94       double c2 = isInverse ? 0.5 : -0.5;
 95       int n3 = n + 3;
 96       int n4 = n >> 2;
 97       for (int i = 2; i <= n4; i++)
 98       {
 99         int i1 = i + i - 1, i2 = i1 + 1, i3 = n3 - i2, i4 = i3 + 1;
         double h1r = c1 * (data[i1] + data[i3]); // 两个分离变换是
         double h1i = c1 * (data[i2] - data[i4]); // 从 data 中分离出来
         double h2r = -c2 * (data[i2] + data[i4]);
         double h2i = c2 * (data[i1] - data[i3]);
         data[i1] = h1r + wr * h2r - wi * h2i; // 此处重新组合以形成
         data[i2] = h1i + wr * h2i + wi * h2r; // 原始实型数据的真实变换
         data[i3] = h1r - wr * h2r + wi * h2i;
         data[i4] = -h1i + wr * h2i + wi * h2r;
         wr = (wtemp = wr) * wpr - wi * wpi + wr; // 递归式
         wi = wi * wpr + wtemp * wpi + wi;
       }
       double tmp = data[1];
       if (!isInverse)
       {
         data[1] = tmp + data[2]; // 同时挤压第一个和最后一个数据
         data[2] = tmp - data[2]; // 使它们都在原始数组中
       }
       else
       {
         data[1] = c1 * (tmp + data[2]);
         data[2] = c1 * (tmp - data[2]);
         ComplexFFT(data, isInverse); // 此处是逆变换
       }
     }
 

第二个方法是 RealFFT，计算单个实函数的快速傅里叶变换。因为在很多情况下期望求快速傅里叶变换的数据由实值样本组成。如果将这些样本放入一个复型数组中，并令其所有虚部为零的话，从执行时间和对存储需求二者来看，其效率是很低的。所以，我们将原始数据放入一个界长只有一半的复型数组中，其偶数项放入该数组的实部，奇数项放入它的虚部。然后调用 ComplexFFT 来计算快速傅里叶变换。

     /// <summary>
     /// 比较 x[0..n-1] 和 y[0..n-1]
     /// </summary>
     /// <param name="x">第一操作数 x[0..n-1]</param>
     /// <param name="y">第二操作数 y[0..n-1]</param>
     /// <param name="n">两个操作数 x 和 y 的精度</param>
     /// <returns>比较结果：-1:小于 1:大于 0:等于</returns>
     public static int Compare(byte[] x, byte[] y, int n)
     {
       Debug.Assert(x.Length >= n && y.Length >= n);
       for (int i = 0; i < n; i++)
         if (x[i] != y[i])
           return (x[i] < y[i]) ? -1 : 1;
       return 0;
     }
 
     //= pp.775, mpneg, 20.6 任意精度的运算
     /// <summary>
     /// 求补码。注意，被操作数被修改。
     /// </summary>
     /// <param name="data">被操作数 data[0..n-1]</param>
     /// <param name="n">被操作数 data 的精度</param>
     /// <returns>被操作数的补码 data[0..n-1]</returns>
     public static byte[] Negative(byte[] data, int n)
     {
       Debug.Assert(data.Length >= n);
       for (int k = Base, i = n - 1; i >= 0; i--)
         data[i] = (byte)((k = MaxValue + k / Base - data[i]) % Base);
       return data;
     }
 
     //= pp.774, mpsub, 20.6 任意精度的运算
     /// <summary>
     /// 减法。从 minuend[0..n-1] 中减去 subtrahend[0..n-1]，得到 difference[0..n-1]
     /// </summary>
     /// <param name="difference">差 difference[0..n-1]</param>
     /// <param name="minuend">被减数 minuend[0..n-1]</param>
     /// <param name="subtrahend">减数 subtrahend[0..n-1]</param>
     /// <param name="n">被减数 minuend 和减数 subtrahend 的精度</param>
     /// <returns>差 difference[0..n-1]</returns>
     public static byte[] Subtract(byte[] difference, byte[] minuend, byte[] subtrahend, int n)
     {
       Debug.Assert(minuend.Length >= n && subtrahend.Length >= n && difference.Length >= n);
       for (int k = Base, i = n - 1; i >= 0; i--)
         difference[i] = (byte)((k = MaxValue + k / Base + minuend[i] - subtrahend[i]) % Base);
       return difference;
     }
 
     //= pp.774, mpadd, 20.6 任意精度的运算
     /// <summary>
     /// 加法。augend[0..n-1] 与 addend[0..n-1] 相加，得到 sum[0..n]
     /// </summary>
     /// <param name="sum">和 sum[0..n]</param>
     /// <param name="augend">被加数 augend[0..n-1]</param>
     /// <param name="addend">加数 addend[0..n-1]</param>
     /// <param name="n">被加数 augend 和加数 addend 的精度</param>
     /// <returns>和 sum[0..n]</returns>
     public static byte[] Add(byte[] sum, byte[] augend, byte[] addend, int n)
     {
       Debug.Assert(augend.Length >= n && addend.Length >= n && sum.Length >= n + 1);
       int k = 0;
       for (int i = n - 1; i >= 0; i--)
         sum[i + 1] = (byte)((k = k / Base + augend[i] + addend[i]) % Base);
       sum[0] += (byte)(k / Base);
       return sum;
     }
 
     //= pp.774, mpadd, 20.6 任意精度的运算
     /// <summary>
     /// 捷加法。augend[0..n-1] 与整数 addend 相加，得到 sum[0..n]
     /// </summary>
     /// <param name="sum">和 sum[0..n]</param>
     /// <param name="augend">被加数 augend[0..n-1]</param>
     /// <param name="n">被加数 augend 的精度</param>
     /// <param name="addend">加数 addend</param>
     /// <returns>和 sum[0..n]</returns>
     public static byte[] Add(byte[] sum, byte[] augend, int n, byte addend)
     {
       Debug.Assert(augend.Length >= n && sum.Length >= n + 1);
       int k = Base * addend;
       for (int i = n - 1; i >= 0; i--)
         sum[i + 1] = (byte)((k = k / Base + augend[i]) % Base);
       sum[0] += (byte)(k / Base);
       return sum;
     }
 
     //= pp.775, mpsdv, 20.6 任意精度的运算
     /// <summary>
     /// 捷除法。dividend[0..n-1] 除以整数 divisor，得到 quotient[0..n-1]
     /// </summary>
     /// <param name="quotient">商 quotient[0..n-1]</param>
     /// <param name="dividend">被除数 dividend[0..n-1]</param>
     /// <param name="n">被除数 dividend 的精度</param>
     /// <param name="divisor">除数 divisor</param>
     /// <returns>商 quotient[0..n-1]</returns>
     public static byte[] Divide(byte[] quotient, byte[] dividend, int n, byte divisor)
     {
       Debug.Assert(quotient.Length >= n && dividend.Length >= n);
       for (int r = 0, k = 0, i = 0; i < n; i++, r = k % divisor)
         quotient[i] = (byte)((k = Base * r + dividend[i]) / divisor);
       return quotient;
     }
 

接下来是比较、求补码、减法、加法、捷加法、捷除法，都相当简单，程序中已经有详细的注释了。

     //= pp.776-777, mpmul, 20.6 任意精度的运算
     /// <summary>
     /// 乘法。multiplicand[0..n-1] 与 multiplier[0..m-1] 相乘，得到 product[0..n+m-1]
     /// </summary>
     /// <param name="product">积 product[0..n+m-1]</param>
     /// <param name="multiplicand">被乘数 multiplicand[0..n-1]</param>
     /// <param name="n">被乘数 multiplicand 的精度</param>
     /// <param name="multiplier">乘数 multiplier[0..m-1]</param>
     /// <param name="m">乘数 multiplier 的精度</param>
     /// <returns>积 product[0..n+m-1]</returns>
     public static byte[] Multiply(byte[] product, byte[] multiplicand, int n, byte[] multiplier, int m)
     {
       int mn = m + n, nn = 1;
       Debug.Assert(product.Length >= mn && multiplicand.Length >= n && multiplier.Length >= m);
       while (nn < mn) nn <<= 1; // 为变换找出最小可用的 2 的幂次
       double[] a = new double[nn + 1], b = new double[nn + 1];
       for (int i = 0; i < n; i++) a[i + 1] = multiplicand[i];
       for (int i = 0; i < m; i++) b[i + 1] = multiplier[i];
       RealFFT(a, false); // 执行卷积，首先求出二个傅里叶变换
       RealFFT(b, false);
       b[1] *= a[1]; // 复数相乘的结果(实部和虚部)
       b[2] *= a[2];
       for (int i = 3; i <= nn; i += 2)
       {
         double t = b[i];
         b[i] = t * a[i] - b[i + 1] * a[i + 1];
         b[i + 1] = t * a[i + 1] + b[i + 1] * a[i];
       }
       RealFFT(b, true); // 进行傅里叶逆变换
       byte[] bs = new byte[nn + 1];
       long cy = 0; // 执行最后完成所有进位的过程
       for (int i = nn, n2 = nn / 2; i >= 1; i--)
       {
         long t = (long)(b[i] / n2 + cy + 0.5);
         bs[i] = (byte)(t % Base); // 原书中这句使用循环，有严重的性能问题
         cy = t / Base;
       }
       if (cy >= Base) throw new OverflowException("FFT Multiply");
       bs[0] = (byte)cy;
       Array.Copy(bs, product, n + m);
       return product;
     }
 

接下来的方法是 Multiply，乘法。其算法是使用 RealFFT 求被乘数和乘数的快速傅里叶变换，将结果相乘，然后进行傅里叶逆变换得到卷积，最后执行适当的进位。其原理已经在上一篇随笔“使用快速傅里叶变换计算大整数乘法”中讲述得很清楚了。

     //= pp.778, mpdiv, 20.6 任意精度的运算
     /// <summary>
     /// 除法。dividend[0..n-1] 除以 divisor[0..m-1]，m ≤ n，
     /// 得到：商 quotient[0..n-m]，余数 remainder[0..m-1]。
     /// </summary>
     /// <param name="quotient">商 quotient[0..n-m]</param>
     /// <param name="remainder">余数 remainder[0..m-1]</param>
     /// <param name="dividend">被除数 dividend[0..n-1]</param>
     /// <param name="n">被除数 dividend 的精度</param>
     /// <param name="divisor">除数 divisor[0..m-1]</param>
     /// <param name="m">除数 divisor 的精度</param>
     /// <returns>商 quotient[0..n-m]</returns>
     public static byte[] DivRem(byte[] quotient, byte[] remainder, byte[] dividend, int n, byte[] divisor, int m)
     {
       Debug.Assert(m <= n && dividend.Length >= n && divisor.Length >= m && quotient.Length >= n - m + 1 && remainder.Length >= m);
       int MACC = 3;
       byte[] s = new byte[n + MACC], t = new byte[n - m + MACC + n];
       Inverse(s, n - m + MACC, divisor, m); // s = 1 / divisor
       Array.Copy(Multiply(t, s, n - m + MACC, dividend, n), 1, quotient, 0, n - m + 1); // quotient = dividend / divisor
       Array.Copy(Multiply(t, quotient, n - m + 1, divisor, m), 1, s, 0, n); //  s = quotient * divisor
       Subtract(s, dividend, s, n); // s = dividend - quotient * divisor
       Array.Copy(s, n - m, remainder, 0, m);
       if (Compare(remainder, divisor, m) >= 0) // 调整商和余数
       {
         Subtract(remainder, remainder, divisor, m);
         Add(s, quotient, n - m + 1, 1);
         Array.Copy(s, 1, quotient, 0, n - m + 1);
       }
       return quotient;
     }
 
     //= pp.777 - 778, mpinv, 20.6 任意精度的运算
     /// <summary>
     /// 求倒数。
     /// </summary>
     /// <param name="inverse">倒数 inverse[0..n-1]，在 inverse[0] 后有基数的小数点</param>
     /// <param name="n">倒数 inverse 的精度</param>
     /// <param name="data">被操作数 data[0..m-1]，data[0] > 0，在 data[0] 后有基数的小数点</param>
     /// <param name="m">被操作数 data 的精度</param>
     /// <returns>倒数 inverse[0..n-1]，在 inverse[0] 后有基数的小数点</returns>
     public static byte[] Inverse(byte[] inverse, int n, byte[] data, int m)
     {
       Debug.Assert(inverse.Length >= n && data.Length >= m);
       InitialValue(inverse, n, data, m, false);
       if (n == 1) return inverse;
       byte[] s = new byte[n],  t = new byte[n + n];
       for (; ; ) // 牛顿迭代法: inverse = inverse * ( 2 - data * inverse )  =>  inverse = 1 / data
       {
         Array.Copy(Multiply(t, inverse, n, data, m), 1, s, 0, n); // s = data * inverse
         Negative(s, n);                                         // s = -(data * inverse)
         s[0] -= (byte)(Base - 2);                             // s = 2 - data * inverse
         Array.Copy(Multiply(t, s, n, inverse, n), 1, inverse, 0, n); // inverse = inverse * s
         int i = 1;
         for (; i < n - 1 && s[i] == 0; i++) ; // 判断 s 的小数部分是否为零
         if (i == n - 1) return inverse; // 若 s 收敛到 1 则返回 inverse = 1 / data
       }
     }
 

 

接下来的方法是 DivRem，除法，以及 Inverse，求倒数。

除法用除数的倒数乘以被除数来计算，倒数值用牛顿迭代法：

U_i+1 = U_i (2 - VU_i)

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

要注意在 DivRem 中，求得的余数可能大于等于除数，此时商比实际的值要小一。所以在程序的 293 到 298 行调整商和余数。

 

     //= pp.778-779, mpsqrt, 20.6 任意精度的运算
     /// <summary>
     /// 求平方根 sqrt，以及平方根的倒数 invSqrt。invSqrt 也可设为 sqrt，此时，invSqrt 也是平方根。
     /// </summary>
     /// <param name="sqrt">平方根 sqrt[0..n-1]，在 sqrt[0] 后有基数的小数点</param>
     /// <param name="invSqrt">平方根的倒数 invSqrt[0..n-1]，在 invSqrt[0] 后有基数的小数点</param>
     /// <param name="n">平方根的精度</param>
     /// <param name="data">被操作数 data[0..m-1]，data[0] > 0，在 data[0] 后有基数的小数点</param>
     /// <param name="m">被操作数 data 的精度</param>
     /// <returns>平方根 sqrt[0..n-1]，在 sqrt[0] 后有基数的小数点</returns>
     public static byte[] Sqrt(byte[] sqrt, byte[] invSqrt, int n, byte[] data, int m)
     {
       Debug.Assert(sqrt.Length >= n && invSqrt.Length >= n && data.Length >= m);
       if (n <= 1) throw new ArgumentOutOfRangeException("n", "must greater than 1");
       InitialValue(invSqrt, n, data, m, true);
       byte[] s = new byte[n], t = new byte[n + Math.Max(m, n)];
       for (; ; ) // invSqrt = invSqrt * (3 - data * invSqrt * invSqrt) / 2 => invSqrt = 1 / sqrt(data)
       {
         Array.Copy(Multiply(t, invSqrt, n, invSqrt, n), 1, s, 0, n); // s = invSqrt * invSqrt
         Array.Copy(Multiply(t, s, n, data, m), 1, s, 0, n);   // s = data * invSqrt * invSqrt
         Negative(s, n);                                     // s = -(data * invSqrt * invSqrt)
         s[0] -= (byte)(Base - 3);                         // s = 3 - data * invSqrt * invSqrt
         Divide(s, s, n, 2);                              // s = (3 - data * invSqrt * invSqrt) / 2
         Array.Copy(Multiply(t, s, n, invSqrt, n), 1, invSqrt, 0, n);   // invSqrt = invSqrt * s
         int i = 1;
         for (; i < n - 1 && s[i] == 0; i++) ; // 判断 s 的小数部分是否为零
         if (i < n - 1) continue; // 若 s 没有收敛到 1 则继续迭代
        Array.Copy(Multiply(t, invSqrt, n, data, m), 1, sqrt, 0, n); // sqrt = invSqrt * data = sqrt(data)
         return sqrt;
       }
     }
 
     /// <summary>
     /// 采用浮点运算以得到一个初始近似值 u[0..n-1]: u = 1 / data 或者 u = 1 / sqrt(data)
     /// </summary>
     /// <param name="u">初始近似值 u[0..n-1]</param>
     /// <param name="n">所需的精度</param>
     /// <param name="data">被操作数 data[0..m-1]</param>
     /// <param name="m">被操作数 data 的精度</param>
     /// <param name="isSqrt">是否求平方根</param>
     /// <returns>初始近似值 u[0..n-1]</returns>
     static byte[] InitialValue(byte[] u, int n, byte[] data, int m, bool isSqrt)
     {
       Debug.Assert(u.Length >= n && data.Length >= m);
       int scale = 16 / Len; // double 可达到 16 位有效数字
       double fu = 0;
       for (int i = Math.Min(scale, m) - 1; i >= 0; i--) fu = fu / Base + data[i];
       fu = 1 / (isSqrt ? Math.Sqrt(fu) : fu);
       for (int i = 0; i < Math.Min(scale + 1, n); i++)
       {
         int k = (int)fu;
         u[i] = (byte)k;
         fu = Base * (fu - k);
       }
       return u;
     }
   }
 }

最后的方法是 Sqrt，求平方根。用牛顿迭代法计算平方根和除法类似。若：

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

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

 

在上一篇随笔“使用快速傅里叶变换计算大整数乘法”中提到：

由于快速傅里叶变换是采用了浮点运算，因此我们需要足够的精度，以使在出现舍入误差时，结果中每个组成部分的准确整数值仍是可辨认的。长度为 N 的 B 进制数可产生大到 B²N 阶的卷积分量。我们知道，双精度浮点数的尾数是 53 个二进位，所以：

2 x log₂B + log₂N + 几个 x log₂log₂N < 53

上式中左边最后一项是为了快速傅里叶变换的舍入误差。

我们假设上式中的“几个”为“三个”，那么，经过简单的计算，得到以下结果：

基数 B	长度 N	十进制数字个数
256	约 107	约 2.4 x 107
100	约 108	约 2 x 108
10	约 109	约 109

注意，基数 B 选取得越小，计算速度就越慢。

转自：http://www.cnblogs.com/skyivben/archive/2008/07/25/1250891.html