《自己动手打造“超高精度浮点数类”》源代码简要导读
                           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);