我的高精度浮点数类，我用它计算出了1亿位的PI(3.1415926...)值

housisong 2007-12-24 11:22:41
转自blog文章: http://blog.csdn.net/housisong/archive/2007/12/24/1965652.aspx

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

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

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

超高精度浮点数类TLargeFloat的完整源代码参见: http://blog.csdn.net/housisong/archive/2005/11/08/525215.aspx ；在导读文章中列出的代码(仅作示例)很可能和实际的代码有区别；导读文章会更关注于算法和原理，而更多的细节需要读者去阅读源代码；

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长度为2，m_Digits[0]==1234,m_Digits[1]=5678; 则这个TLargeFloat代表的是浮点数：-0.12345678E-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的数也不会产生任何误差；

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进制数时长度不能超过约17百万位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/4
重复计算:y=a a=(a+b)/2 b=sqrt(by) c=c-x(a-y)^2 x=2x
最后:pi=(a+b)^2/(4c)

实现它的函数是: TLargeFloat GetAGMPI(unsigned long uiDigitsLength);
我也实现了Borwein四次收敛迭代式：TLargeFloat GetBorweinPI(unsigned long uiDigitsLength);



...全文
2969 46 打赏 收藏 转发到动态 举报
写回复
用AI写文章
46 条回复
切换为时间正序
请发表友善的回复…
发表回复
yaos 2008-02-14
  • 打赏
  • 举报
 回复
晕，是41.5M
推翻上帖的文字
这么做看合适不
少开内存缓存，多开内存映射文件
理由同上
yaos 2008-02-14
  • 打赏
  • 举报
 回复
觉得内存占用忒大

应该有办法压缩的

一个1亿位数字刚计算了，不过18M多点而已

你开10个数字缓存足够了啊
再不行，可以考虑内存映射文件的
后期一次计算，足够读多少次内存映射文件的

计算过，读20M内存映射文件绝对是不超过1秒的时间
housisong 2008-01-13
  • 打赏
  • 举报
 回复
结贴 我的计算PI算抛砖引玉 等待gxqcn的实现:)
goodmrning 2008-01-03
  • 打赏
  • 举报
 回复
lz 好厉害!!!!!!!!!!!!!!!!!!!!!!!!
up!!!!!!!!!!!!1
jf!!!!!!!!!!!!
gxqcn 2008-01-03
  • 打赏
  • 举报
 回复
可以预见的是，如果在四核上对比测试，HugeCalc 将领先 GMP 更多；可惜我没有这样的测试环境，哪位热心网友来帮这个忙？:)

更详细深入的探讨见上的帖子：http://emath.5d6d.com/thread-27-1-1.html

有了大数乘法计算PI确实很容易:) 不过近段时间有点忙，编写以HugeCalc内核的超高精度PI暂未进行任何实质性的工作。
housisong 2008-01-03
  • 打赏
  • 举报
 回复
恭喜 有了大数乘法 计算PI很容易的
gxqcn 2008-01-03
  • 打赏
  • 举报
 回复
好消息，在双核上测试，我的 HugeCalc 核心算法效率全面超出 GMP！
[code=INIFile]
GMP HugeCalc-HI HugeCalc-HX
Fib[100000000]( 69424191 Bits ) 4.349451 s 3.702678 s 3.479291 s
Fib[200000000]( 138848382 Bits ) 9.424257 s 7.230688 s 6.871886 s
Fib[400000000]( 277696765 Bits ) 20.659469 s 15.330399 s 14.545750 s
[/code]

测试代码如下：


// DESCRIPTION for Fib.cpp :

// Defines the entry point for the console application.

// 目的：测试 GMP 与 HugeCalc 两大高精度算法库核心算法效率

//       因为计算 Fibonacci 本身可优化的地方不多，关键是大数乘法及大数平方，

//       故，通过测试 Fibonacci 效率而检测算法库最核心的大数乘法（平方）效率。

//

// 设计：郭先强 ( gxqcn@163.com; HugeCalc@Gmail.com )

// 日期：2008-01-03

//

// Web: http://www.emath.ac.cn/

// BBS: http://bbs.emath.ac.cn/

//

//////////////////////////////////////////////////////////////////////////



// Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:

//      Win32 Debug:    Debug Multithreaded DLL

//      Win32 Release:  Multithreaded DLL



#include <iostream.h>



#include <gmp.h>



#include <HugeCalc.h>   // 公共接口

#include <HugeInt.h>    // 10进制系统

#include <HugeIntX.h>   // 16进制系统



#pragma message( "automatic link to .../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )

#pragma comment( lib, "gmp.lib" )

#pragma comment( lib, "HugeCalc.lib" )



#if 0

#   define FIRST_INDEX  10000000L   /*1e7*/

#else

#   define FIRST_INDEX  100000000L  /*1e8*/

#endif



int main(int argc, char* argv[])

{

    cout << "Call " << HugeCalc::GetVer() << endl;



    if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel() )

    {

        cout << "\r\n警告：您未通过 HugeCalc.dll 的许可认证！\r\n" \

            << "\r\n解决方案可选下列方案之一：" \

            << "\r\n    一、请将本程序移动到“/CopyrightByGuoXianqiang/[../]”目录下运行；" \

            << "\r\n    二、或请在 HugeCalc.chm 中进行注册（一劳永逸）。" \

            << endl;

    }

    else

    {

        for ( UINT32 u32Index = FIRST_INDEX; u32Index <= FIRST_INDEX * 4; u32Index <<= 1 )

        {

            cout << "\r\n--------------------------------\r\n" << endl;

            cout << "Calculate Fib[" << u32Index << "]:\r\n" << endl;



#if 1

            // 测试 GMP

            {

                mpz_t hugeFIB;

                mpz_init( hugeFIB );



                HugeCalc::EnableTimer( TRUE );

                HugeCalc::ResetTimer();

                mpz_fib_ui( hugeFIB, u32Index );

                HugeCalc::EnableTimer( FALSE );



                cout << "GMP used time: " << HugeCalc::GetTimerStr( FT_DOT06SEC_s ) << " ( " \

                    << mpz_sizeinbase( hugeFIB, 2 ) << " Bits )" << endl;



                mpz_clear( hugeFIB );

            }

#endif



            // 测试 HI

            {

                CHugeInt hugeFIB;



                HugeCalc::EnableTimer( TRUE );

                HugeCalc::ResetTimer();

                hugeFIB.Fibonacci( u32Index );

                HugeCalc::EnableTimer( FALSE );



                cout << "HI used time: " << HugeCalc::GetTimerStr( FT_DOT06SEC_s ) << " ( " \

                    << hugeFIB.GetDigits() << " Digits )" << endl;

            }



            // 测试 HX

            {

                CHugeIntX hugeFIB;



                HugeCalc::EnableTimer( TRUE );

                HugeCalc::ResetTimer();

                hugeFIB.Fibonacci( u32Index );

                HugeCalc::EnableTimer( FALSE );



                cout << "HX used time: " << HugeCalc::GetTimerStr( FT_DOT06SEC_s ) << " ( " \

                    << hugeFIB.GetBits() << " Bits )" << endl;

            }

        }

    }



    cout << endl;

    system( "pause" );



    return 0;

}


huangxw000 2008-01-03
  • 打赏
  • 举报
 回复
很好，很强大！
Eddie005 2008-01-03
  • 打赏
  • 举报
 回复
再来接分~~
Eddie005 2008-01-02
  • 打赏
  • 举报
 回复
@_@
grellen 2008-01-02
  • 打赏
  • 举报
 回复
顶
liangbch 2007-12-29
  • 打赏
  • 举报
 回复
不好意思，31楼 写错了
“
位数 ooura FFT用时 super-pI 用时
104万 17 sec 65 sec
419万 321 sec 90 sec
”
应该是
"位数 ooura FFT用时 super-pI 用时
104万 17 sec 65 sec
419万 90 sec 321 sec"


ooura 的计算圆周率的代码 晦涩难懂还是其次，倒是它的FFT代码极为难懂，简直就是天书。
housisong 2007-12-29
  • 打赏
  • 举报
 回复
to: liangbch
谢谢你的测试
你提到的ooura中减小误差的方法很有价值，这两天试试，看看效果
什么时候集成ooura的FFT试试；
他们的代码写的都有点晦涩;计算PI的代码写得巨复杂; (我的计算PI的代码可以用很少量的代码很清晰的写出来，就如写公式一样) 还是看着自己的代码舒服:D

FFTW以前也知道，但没有用过，不过intel一直宣称它的FFT优化库的代码速度超过FFTW很多(而且接口基本兼容)；

一生有爱1980 2007-12-29
  • 打赏
  • 举报
 回复
学习学习
liangbch 2007-12-28
  • 打赏
  • 举报
 回复
关于FFTW和其它FFT程序的速度对比见：http://www.fftw.org/speed/Pentium4-2.4GHz-icc/
liangbch 2007-12-28
  • 打赏
  • 举报
 回复
刚刚对ooura FFT的例子程序（计算PI）测试了一下.下面给出 它 和SUPER PI 的速度对比。（fft的代码使用fftsg.c）

测试环境：PIV 2.6G, 768M RAM, windows XP

位数 ooura FFT用时 super-pI 用时
104万 17 sec 65 sec
419万 321 sec 90 sec

从上面的数据可以看出 ooura FFT的速度为super PI 的3-4倍左右。

下面介绍一点儿 ooura FFT 计算圆周率些程序的一些信息：
1. 该程序使用的是 实数 FFT 变换，占用内存较复数FFT少了一半，速度也快了一倍。
2. 为了减小误差， ooura FFT 并不是将4位数直接转化为double型的，而是做了一些处理，即：本位 减去 0.5 × RADIX, 高位加上0.5,这个做法可以将每个浮点数的绝对值拉小，从而减小误差。
示例如下：
整数序列： 0 1234 9999 1111
直接转化为double型序列为： 0 1234.00 9999.00 1111.00
ooura FFT 的做法： 0.5 -3765.5 4999.5 -3999.00

3.用FFT算法计算大数乘法，累计误差不能超过0.5, ooura FFT 提前评估误差，如果误差较大，采用每个浮点数表示更少的10进制数位，当计算104万，419万，3355万时， ooura FFT 均使用每个浮点数 表示 4位10进制数的方法。楼主提到“采用4个10进制数时长度不能超过约17百万位10进制位”，但是我们在ooura FFT 中看到，即使计算到 3355万时，仍然采用每个浮点数表示4位10进制的做法，这说明ooura FFT 误差控制平衡得很好。这个程序显示，当计算104万，419万，3355万时， 误差上界分别为： 0.005，0.021，0.16，并没有超过0.5.

4. 西方最快的FFT程序是FFTW,我刚才查了一下，实数FFT变换，当序列长度为262144时，FFTW的速度指数是564,而ooura FFT的3种算法（4g,8g,sg)的速度指数分别为266,314,438.即FFTW比sg快了29%，这意味着，使用同样的计算圆周率算法，如果FFT换用fftw的代码，速度可以提高近三分之一。

5. 对照使用GMP计算圆周率的程序，如果GMP使用的是同样的计算圆周率的算法，它的大数乘法的实现应该比最快的FFT( FFTW )的大数乘法更快。故我以为超过GMP算法是一个非常困难的事情。退一步，我们想要写一个超过ooura FFT的大数乘法也很困难。
mathe 2007-12-28
  • 打赏
  • 举报
 回复
呵呵，正因为模幂算法是如此重要，所以更加需要好好衡量一下，我相信gmp对于它也会做不少的优化，gmp也是有直接的模幂函数可以调用的。你可以两者都评估一下看看。
gxqcn 2007-12-28
  • 打赏
  • 举报
 回复
to mathe:
模幂算法是大数素性判定和RSA算法的核心，非常重要，研究者众，提交的论文也非常多。我的模幂算法参考了不少论文，并糅合了自己原创的部分算法，效率也非常高，大家可以从我的网站下载压缩包，用里面的HugeCalc\testDLL\bin\NumberTheoretic.exe进行测试；相对计算 Fibonacci 来说，模幂算法牵涉到自身的算法更多，很难直接评判到核心算法的效率。
我已下载了gmp库文件4.2.1，将抽空去做一下对比测试，我估计HugeCalc不会输于它。:)
mathe 2007-12-28
  • 打赏
  • 举报
 回复
已经相当不错了:) 不过看复杂度增长关系，104万倍只有慢一倍多一些，但是1677万位慢的多很多，我估计可能内存方面优化的没那么好
housisong 2007-12-28
  • 打赏
  • 举报
 回复
to: mathe
时间复杂度方面的改进很小了;用agm公式计算PI的算法里起决定因素的是大数乘法，我的是接近N*log(N)(略大于) 现在理论界应该还没有超越这个复杂度的可用算法；
但实际上我估计GMP的代码速度是我的代码的几十倍:D

作为对比,在我的酷睿4400(2.0G)上，计算PI (位数不完全相同)
104万位 Surper PI 是 31秒 我的代码是 72秒
1677万位 Surper PI 是 758秒 我的代码是 2771秒
加载更多回复（26）
28.精讲JavaScript数字类型，你不知道的必考点。
文章目录JavaScript数字类型详解前言数字的写法分隔符省略的0十六、八、二进制toString(base)舍入精度偏差偏差现象偏差原因解决方法Infinity、NaNObject.isparseInt、parseFloatMath对象总结 JavaScript数字类型详解 你真的完全了解JavaScript数字类型吗，看完这篇在回答我！ 前言 JavaScript中存在两种表示数字的类型： Number类型，也就是常规意义上的数字类型，以64位的IEEE-754格式存储，属于“双精度浮点数”，截至.
java学习笔记-基础知识-2023.3.29
学习网站：https://www.sxt.cn/Java_jQuery_in_action/History_Direction.html 第一章 Java入门 计算机语言经历了三代：第一代是机器语言，第二代是汇编语言，第三代是高级语言 1.3.2 Java的核心优势–跨平台（JAVA虚拟机(JVM)是JAVA实现跨平台的核心） 1.3.3 Java各版本的含义 JavaSE(Java Standard Edition)：标准版，定位在个人计算机上的应用 JavaEE(Java Enterprise Edit
面向对象之14：开发中的常用类之一： Math、Random、System、Arrays、BigInteger、BigDecimal、Scanner 的使用总结：
数学相关的类：Math ，Math在 java.lang包下，直接使用，无需导包。 Math 的 UML类图： ，Math被final所修饰，不可被继承。 ，构造器私有化，不能创建Math的对象。（只能是 Math 类 自己在自己的类中使用自己的对象。） Math类内部所有的 属性、方法都是 静态的，被static 修饰，直接用类名.静态成员 访问即可。 ，如果本类中的静态成员有和静态导入的类的静态成员重名时，在运行程序调用静态成员的时候，优先调用的是本类中的静态成员！没有，则是使用静态导...
今天是圆周率节，我们一起来用C++来求π吧！
前言 3.1415926535897932384626433832795… 不好意思，小编我只能被到这里了，当然我觉得做题只要3.14就够了，没必要背那么长的啊。。。当时我上小学时总是和同学们比背π，后来发现也没什么必要哈~当时我们班有能背到100位的，还被我们膜拜了一番（虽然背到30几位我们也不知道他是不是在乱说）。当然了，如果哪位想锻炼以下记忆力的，你尽管背好了。（Haha）因为圆周率是一个...
深入浅出大数阶乘
深入浅出大数阶乘 大数阶乘的计算是一个有趣的话题，从中学生到大学教授，许多人都投入到这个问题的探索和研究之中，并发表了他们自己的研究成果。如果你用阶乘作关键字在google上搜索，会找到许多此类文章，另外，如果你使用google学术搜索，也能找到一些计算大数阶乘的学术论文。但这些文章和论文的深度有限，并没有给出一个高速的算法和程序 菜鸟篇 程序1，一个最直接的计算阶乘的程序 ...
数据结构与算法

33,008

社区成员

35,326

社区内容

发帖
与我相关
我的任务
社区描述
数据结构与算法相关内容讨论专区
社区管理员
  • 数据结构与算法社区
加入社区
  • 近7日
  • 近30日
  • 至今

加载中

查看更多榜单
社区公告
暂无公告

试试用AI创作助手写篇文章吧

+ 用AI写文章