《自己动手打造“超高精度浮点数类”》源代码简要导读

                 《自己动手打造“超高精度浮点数类”》源代码简要导读
                           HouSisong@GMail.com

tag: PI,超高精度浮点数,TLargeFloat,FFT乘法,二分乘法,牛顿迭代法,borwein四次迭代,AGM二次迭代,源代码导读

摘要: 很多人可能都想自己写一个能够执行任意精度计算的浮点数；:D我写的第一个程序就是用qbasic计算自然数e到100万位(后来计算PI);  我的blog文章《自己动手打造“超高精度浮点数类”》里有一个C++类的实现TLargeFloat，它能够执行高精度的浮点数运算；演示代码里面有一个计算PI的Borwein四次迭代式和一个AGM二次迭代式（我用它计算出了上亿位的PI小数位:）  本文章是对其源代码的进一步解读；

  本系列文章的由来: 源于一次在abp论坛(现在的www.cpper.com)的帖子,那是2004年的时候,有初学者询问高精度计算的原理；我那时比较激动(哈哈)，因为这不就是我刚学编程的时候做的吗?! 就计划"重操旧业"用C++给初学者写一个高精度计算的Demo程序(决定实际动手写要归功于abp的codinggirl)；第一个版本实现前后花了10小时；但实际上完成的代码没有真正发布，也缺少源代码的分析和介绍；后来放在了csdn的blog上；由于实在太慢，摆在我的blog里“碍眼”:D , 最近对它做了一些速度改进；本文章是它的简要导读；
  超高精度浮点数类TLargeFloat的完整源代码参见: http://blog.csdn.net/housisong/archive/2005/11/08/525215.aspx ；在导读文章中列出的代码(仅作示例)很可能和实际的代码有区别；导读文章会更关注于算法和原理，而更多的细节需要读者去阅读源代码；

ps:题外话(因为是blog文章,说点题外话才是正常状态)，回想自己上大学的时候，不务正业的整天跑(泡)图书馆，对编写计算机程序产生了浓厚的兴趣；我的第一个程序是想计算出自然数e(2.718281828...)的上百万位;由于没有电脑，就把所有代码都写在纸上(用的qbasic语言)，觉得差不多了的时候就到学校的计算机房去把代码敲到电脑里；那时候对电脑上编写程序一点概念都没有(高中学习接触的根本不能算写程序)；折腾了一阵子后旁边的一个玩游戏中(好像)的学长实在看不下去了，转过头来告诉我怎么打开qbasic环境!  我笨拙的把纸上写好的代码敲入计算机，在修改了几处打字错误后运行成功，运行第一次就正确的输出了e的上千位小数，哈哈，当时很有成就感； 在纸上写程序和修改，并在脑袋中"运行"多遍后再让电脑来实际运行的那段记忆只能是美好的回忆了 (记得有一个程序我甚至写了几十页的) :D

a.高精度浮点数的表示方法(数据结构)

超高精度浮点数类TLargeFloat的数据结构定义：
class TLargeFloat 
{
    long               m_Sign;     //符号位
    long               m_Exponent; //保存10为底的指数
    std::vector<long>  m_Digits;   //小数部分
};
  其中:  m_Sign用来保存该浮点数的正负号(正1和负1)；我用0来代表浮点数0，这样对于某些判断比较方便；
    m_Exponent保存的是该浮点数的指数部分，其值代表的是10进制为底的指数 (指数部分可以使用int64整数)
    m_Digits数组保存的是该浮点数的小数尾数部分，每个元素保存4位10进制数，也就是取值[0--9999]; 对于m_Digits[0]，正常情况下的值取值范围为[1000,9999],否则该浮点数就是未规格化的(未规格化小数只存在于运算过程中的临时值)；
    比如当m_Sign=-1,m_Exponent=-100,m_Digits长度为3，m_Digits[0]==1234,m_Digits[1]=5678,m_Digits[2]=258; 则这个TLargeFloat代表的是浮点数：-0.123456780258E-100

  源代码中，为了捕获指数运算超出值域的异常情况我定义了一个TCatchIntError类用来包装m_Exponent所使用的整数类型；TCatchIntError实现了“+=” 、“-=” 、“*=”3个运算符，并处理可能的值域超界，如果发生值域超界将抛出TLargeFloatException类型的异常。

  小数部分数据结构的选择，将决定很多算法的具体实现细节；为了容易编码和阅读我舍弃了一些更节约内存的设计也舍弃了一些复杂的优化(比如汇编等)；选择4位10进制数将简化很多代码的实现(后面会看到)； 对于更专业的实现，建议使用8(或9)位10进制数或者以2进制为基础来实现；

b.其他数到超高精度浮点数类TLargeFloat的转换
  TLargeFloat实现了long double到TLargeFloat的多个转换和相互运算，这样那些能够自动隐式转换到long double的类型(比如int、float等)也能自动的参与TLargeFloat的运算；程序也提供了TLargeFloat与字符串类型之间的转换方法：AsString()和StrToLargeFloat(); 这对于超大和超高精度的数传递给TLargeFloat就比较有用，而且通过字符串的方式转换到TLargeFloat的数也不会产生任何误差！

  数值转换成TLargeFloat的构造函数，一般的声明是 TLargeFloat(数值,高精度浮点数的有效精度);但这样写容易混浠，比如 TLargeFloat(200,300); 所以我定义了一个类TDigits，要求这样写代码：TLargeFloat(200,TLargeFloat::TDigits(300)); 

c.高精度浮点数的加减法的实现
  加减法可以通过绝对值加和绝对值减来实现(Abs_Add(),Abs_Sub_Abs()函数)，先考虑参与运算的TLargeFloat浮点数的正负号然后调用绝对值加(符号相同)或绝对值减(符号不同)就可以了；
  要实现两个高精度数相加减，需要首先调整到指数相同(就如手工计算时的小数点对齐)；
  比如: (-0.1234E2) + (-0.5678E-1) == -(0.1234E2 + 0.5678E-1)== -( 0.1234 +0.0005678 )*1E2 == -0.1239678E2
   (-0.1234E2) + (0.5678E-1) == -(0.1234E2 - 0.5678E-1)== -( 0.1234 -0.0005678 )*1E2 == -0.1228322E2

  核心实现函数：
  小数部分的多位数对齐加法: void TLargeFloat_ArrayAdd(TInt32bit* result,long rsize,const TInt32bit* y,long ysize);
      它实现将y加到result上；
  小数部分的多位数对齐减法: void TLargeFloat_ArraySub(TInt32bit* result,long rsize,const TInt32bit* y,long ysize);
      它实现从result中减去y；
 

d.高精度浮点数的乘法的简单实现
  回想一下我们是怎么手工计算多位数乘法的，那高精度浮点数的乘法也可以这样实现；

  核心实现函数：
  小数部分的多位数乘法: TLargeFloat_ArrayMUL_ExE(TInt32bit* result,long rminsize,const TInt32bit* x,long xsize,const TInt32bit* y,long ysize)
  简要过程如下:
    for (long i=0;i<xsize;++i)
    {
        for (long j=0;j<ysize;++j)
        {
            result[j+i+1]+=(x[i]*y[j]);
            ...//进位
        }
    }
  实际的代码做了一些优化；为了减少进位的次数,当result快超出值域的时候,才会执行进位;
  该函数的时间复杂度为O(N*N);在进行高精度运算过程中(N很大的时候)，绝大部分时间都在做乘法运算；后面会对该实现进行优化；

e.高精度浮点数的除法和开方的实现 牛顿迭代法
  高精度浮点数的除法和开方可以用牛顿迭代法来实现，也就是利用乘法来实现；原理和推导(源代码中也有)可以参见我的blog文章《代码优化－之－优化除法》：http://blog.csdn.net/housisong/archive/2006/08/25/1116423.aspx  (包括牛顿迭代法的原理、图示 、 除法和开方的牛顿迭代公式的推导和公式的收敛速度证明)
  由于有效精度随迭代次数成倍增加，所以可以控制当前参与运算的高精度数的精度，从而优化速度；完成该自动精度调整功能的类是TDigitsLengthAutoCtrl；(开始用较低的运算精度，随着迭代次数的增加，增加运算的精度)这样实现后，除法和开方的时间复杂度和乘法一样(想当于几次高精度乘法)；

f.用二分法来优化乘法实现
  两个高精度的数相乘x*y ; 将x,y从数组中间分成两部分表示: x=a*E+b; y=c*E+d; 
  x*y=(a*E+b)*(c*E+d)=(E*E-E)*ac + E*(a+b)*(c+d) -(E-1)*bd；
  可见，长度为N的数相乘，可以分解为3次长度为N/2的数相乘(递推)； 时间上有递推式 f(N)=3f(N/2); 求解该方程可得该算法时间复杂度为(N^log2(3)); 比前面的O(N^2)复杂度低了很多；
 
  二分法乘法核心实现函数：void ArrayMUL_Dichotomy(TInt32bit* result,long rminsize,const TInt32bit* x,const TInt32bit* y,long MulSize);

g.用FFT来优化乘法实现
  高精度数的乘法其实可以看作一种卷积运算，可以转化为傅立叶变换运算 (进阶实现:转换成数论变换!) ，而它可以有O(N*log(N))复杂度的算法--快速傅立叶变换(FFT)；
  FFT核心函数: void FFT(TComplex* a,const TComplex* w,const long n);
  FFT实现的乘法核心函数：void MulFFT(const TInt32bit* ACoef,const long ASize,const TInt32bit* BCoef,const long BSize,TInt32bit* CCoef,const long CMinSize);
  这里的实现使用的是double实现的复数算法；当乘数的长度超过一定范围的时候，计算中产生的误差会放大到超过0.5从而造成取整错误；
所以MulFFT支持的运算长度是有限制的；采用4个10进制数为基时长度不能超过约16百万位多10进制位；当然可以使用更短的基从而提高允许的最大的位数；但实际上这样的内存消耗太恐怖了，所以不适合采用，而是在超过设定长度的时候转调用ArrayMUL_Dichotomy的实现；

  这里我们有了3个乘法的实现，它们在不同的情况下各有各的优势，当乘法位数较短的时候，TLargeFloat_ArrayMUL_ExE更快；再大一些时使用ArrayMUL_Dichotomy更快，再大一些的话使用MulFFT，当超过MulFFT允许的最大长度时调用ArrayMUL_Dichotomy来拆开乘法运算，使乘法长度适合MulFFT;
  这个综合函数就是：void ArrayMUL(TInt32bit* result,long rminsize,const TInt32bit* x,long xsize,const TInt32bit* y,long ysize);  它是整个高精度运算库的核心;

h.高精度PI值计算采用的公式
   我用超高精度浮点数类TLargeFloat实现的AGM算法计算出了PI的1亿位小数；
   AGM二次收敛迭代公式:
      初值:a=x=1; b=1/sqrt(2); c=1.0/4;
      重复计算:y=a; a=(a+b)/2; b=sqrt(b*y); c=c-x*sqr(a-y); x=2*x;
      最后:pi=sqr(a+b)/(4*c);

   实现它的函数是TLargeFloat GetAGMPI(unsigned long uiDigitsLength);  (参数是需要计算的有效10进制位数): 
   我也实现了Borwein四次收敛迭代式：TLargeFloat GetBorweinPI(unsigned long uiDigitsLength);