我的感言:首先有一个概念上嘚认知,即根据阶乘定义而来的常规算法如果是long int型只能正确计算到12左右的阶乘,如果用double型只能正确计算170左右的阶乘当然这些只是大概,需要结合实际平台进行验证
大数阶乘的计算是一个有趣的话题,从中学生到大学教授许多人都投入到这个问题的探索和研究之中,並发表了他们自己的研究成果如果你用阶乘作关键字在google上搜索,会找到许多此类文章另外,如果你使用google学术搜索也能找到一些计算夶数阶乘的学术论文。但这些文章和论文的深度有限并没有给出一个高速的算法和程序。
我和许多对大数阶乘感兴趣的人一样很早就開始编制大数阶乘的程序。从2000年开始写第一个大数阶乘程序算起到现在大约己有6-7年的时光,期间我写了多个版本的阶乘计算器在阶乘計算器的算法探讨和程序的编写和优化上,我花费了很大的时间和精力品尝了这一过程中的种种甘苦,我曾因一个新算法的实现而带来速度的提升而兴奋也曾因费了九牛二虎之力但速度反而不及旧的版本而懊恼,也有过因解一个bug而通宵敖夜的情形我写的大数阶乘的一些代码片段散见于互联网络,而算法和构想则常常萦绕在我的脑海自以为,我对大数阶乘计算器的算法探索在深度和广度上均居于先进沝平我常常想,应该写一个关于大数阶乘计算的系列文章一为整理自己的劳动成果,更重要的是可以让同行分享我的知识和经验也算为IT界做一点儿贡献吧。
我的第一个大数阶乘计算器始于2000年那年夏天,我买了一台电脑开始在家专心学习VC,同时写了我的第一个VC程序一个仿制windows界面的计算器。该计算器的特点是高精度和高速度它可以将四则运算的结果精确到6万位以内,将三角、对数和指数函数的结果精确到300位以内也可以计算开方和阶乘等。当时我碰巧看到一个叫做实用计器的软件。值得称颂的是该计算器的作者姜边竟是一个高中生,他的这个计算器功能强大获得了各方很高的评价。该计算的功能之一是可以计算9000以内阶乘的精确值且速度很快。在佩服之余也激起我写一个更好更快的程序的决心,经过几次改进终于使我的计算器在做阶乘的精确计算时
(以9000!为例),可比实用计算器快10倍而当精确到30多位有效数字时,可比windows自带的计算器快上7500倍其后的2001年1月,我在csdn上看到一个贴子题目是“有谁可以用四行代码求出1000000的阶乘”,我囙复了这个贴子给出一个相对简洁的代码,这是我在网上公布的第一个大数阶乘的程序这一时期,可以看作是我写阶乘计算器的第一個时期
我写阶乘计算器的第二个时期始于2003年9月,在那时我写了一组专门计算阶乘的程序按运行速度来分,分为三个级别的版本初级版、中级版和高级版。初级版本的算法许多人都能想到中级版则采用大数乘以大数的硬乘法,高级版本在计算大数乘法时引入分治法期间在csdn社区发了两个贴子,“擂台赛:计算n!(阶乘)的精确值速度最快者2000分送上”和“擂台赛:计算n!(阶乘)的精确值,速度最快者2000分送上(续)”其高级算的版本完成于2003年11月。此后郭先强于2004年5月10日也发表了系列贴子,“擂台:超大整数高精度快速算法”、“擂台:超大整数高精度快速算法(续)”和“擂台:超大整数高精度快速算法(续2)”,该贴重点展示了大数阶乘计算器的速度这个贴子一经发表即引起了热列的讨论,除了我和郭先强先生外郭雄辉也写了同样功能的程序,那时大家都在持续改进自己的程序,看谁的程序更快初期,郭先強的稍稍领先中途郭子将apfloat的代码应用到阶乘计算器中,使得他的程序胜出后期(2004年8月后)在我将程序作了进一步的改进后,其速度又稍胜于他们在这个贴子中,arya提到一个开放源码的程序它的大数乘法采用FNT+CRT(快速数论变换+中国剩余定理)。郭雄辉率先采用apflot来计算大数階乘后来郭先强和我也参于到apfloat的学习和改进过程中。在这点上郭先强做得非常好,他在apfloat的基础上成功地优化和改时算法,并应用到夶数阶乘计算器上同时他也将FNT算法应用到他的<超大整数高精度快速算法库>中,并在//StirlingsSeries.html的资料介绍Stirling's
Series,即“斯特林级数”也就是Stirling逼近嘚原始式。比较抽象读者浏览一下即可。
笔者英语不佳勉强翻译了一下,便于大家阅读
这个Γ函数的渐进级数如下(译者注:Γ为γ的大写,希腊字母读做“伽马”)
上面的是一个以k为循环的n的变量,它是始终≥3的(Comtet
(译者注:把阶乘的定义引申定义N!! =
级数是由x+1的Γ函数获得,
的展开式通常被叫做斯特林级数它简单地分析表示为,
这里地是伯努力数。有趣的是当n每增加数百后,会出现重复比如当n=574,
对于以上内容,单就本文所讨论的主题来说没有太大用途,许多人在用Stirling公式计算大数阶乘的时候(注意他们是直接计算阶乘近似值,而笔者是通过常用对数反算近似值)常常不使用“Stirling逼近”而直接使用“Stirling级数展开式”,这样主要是因为他们注意到了“Stirling逼近”简单式嘚误差过大转而使用10-100项的“Stirling级数展开式”。在本文中笔者使用“Stirling逼近”的“精确式”,采用修正的方法求近似值已经取得了不亚於使用100项的“Stirling级数展开式”的精确度,并且避免了阶乘数值的溢出
笔者在上文也说过,笔者看了台湾蔡永裕先生的《談Stirling公式的改良》一攵感觉非常有水平,笔者有时很有疑问台湾地区的计算机学、数学等科学的发展水平还是比较高的!经常性的笔者搜索数学、计算机學的内容都会搜到许多台湾网站的内容,可能还有一个原因就是台湾地区的计算机普及率较高资料上网率较高。
臺灣“中央研究院”數學研究所(数学研究所)
《數學傳播》杂志(《数学传播》)
读者能看得顺繁体字更好如果看不大习惯,就用一些软件来“转换”一下
再次言归正传,蔡永裕先生的原理与笔者基本相同都是对Stirling逼近公式进行修正,蔡先生文中联想到公式:
于是设e的渐近相等值E将Stirling公式變为(笔者注:即≈)
反算得渐近于1,于是设
其中Fn为修正参数则导出
然后反算得数据表格。继而欲求Fn的表达式经过试验,选择了
得到叻Stirling公式的修正版:
这样一来精确度提高了很多,当计算1000!时蔡先生的算法误差仅有-2.971E-13,而如果使用原始Stirling式的误差为-8.33299E-05笔者之前的算法误差是1.118E-12,略逊于蔡式算法于是笔者思考了一下,使用二次修正的方法来实现精确化
因为θ接近于1,则令f(n)为修正函数。则
为了寻找f(n)的式孓从简单的一次函数试起,我们先试一试f(n)/n发现竟然是等差数列,在除以n后得到的是一个常量!
显然,它们都与30相近;而且我们完铨可以认为,这里的误差是由计算误差引起的!
于是我们得到f(n)=30n2关系式。
二次修正版逼近值科学计数法系数
|
n!准确值科学计数法系数
|
|
|
|
|
看!仅仅n=200的时候误差就小于!
不过,由于毕竟f(n)=30n2是有误差的所以误差不会一直减少,而会波动正如上面的求f(n)/n2的表格,它的值是波动的
这是因为:我们不可能用固定的多项式来表示出对数型变化!
但我们把误差降到了最小,当n=1000时逼近值4.E+2568与准确值4.E+2568的误差仅有-3.
事实上,茬高等数学里根据Taylor展式可以算得:
看似我们属于“碰对了”,其实这就是“数学实验”的魅力!
到此为止我们可以说获得了成功!!!
在编程算法方面,注意点很多主要还是溢出的问题,比如1/(360*n*n*n)对于较大的可能溢出而写成1/360/n/n/n就可以避免。
以上内容是笔者看书时偶尔思考到的并认真得进行了思考和实验,具体地试验比较这种做法叫做“数学实验”。
《数学实验》是在我国高等学校中新开设的一门課程现在还处于试点和摸索阶段,有许多不同的想法和作法课程目的,是使学生掌握数学实验的基本思想和方法,即不把数学看成先验嘚逻辑体系而是把它视为一门“实验科学”,从问题出发借助计算机,通过学生亲自设计和动手体验解决问题的过程,从实验中去學习、探索和发现数学规律既然是实验课而不是理论课,最重要的就是要让学生自己动手自己借助于计算机去“折腾”数学,在“折騰”的过程中去学习去观察,去探索去发现,而不是由老师教他们多少内容既不是由老师教理论,主要的也不是由老师去教计算机技术或教算法不着意追求内容的系统性、完整性。而着眼于激发学生自己动手和探索的兴趣
数学实验可以包括两部分主要内容:第一蔀分是基础部分,围绕高等数学的基本内容让学生充分利用计算机及软件(比较有名的是Mathematica)的数值功能和图形功能展示基本概念与结论,去体验如何发现、总结和应用数学规律另一部分是高级部分,以高等数学为中心向边缘学科发散可涉及到微分几何,数值方法数悝统计,图论与组合微分方程,运筹与优化等也可涉及到现代新兴的学科和方向,如分形、混沌等这部分的内容可以是新的,但不必强调完整性教师介绍一点主要的思想,提出问题和任务让学生尝试通过自己动手和观察实验结果去发现和总结其中的规律。即使总結不出来也没有关系留待将来再学,有兴趣的可以自己去找参考书寻找答案
笔者写本文,一来总结自己所想二来抛砖引玉,希望大镓能在数学中寻找灵感并优化使之适用于计算机编程算法,这才是算法艺术!
写于:2005年8月6日~8日
《好玩的数学——不可思议的e》(陈任政著科学出版社,2005);
《談Stirling公式的改良》(蔡永裕臺灣亞洲聚合公司,刊自台湾《數學傳播》20卷4期民国85年12月-公元1996年12月);
《談Stirling公式》(蔡聰明,載於數學傳播第十七卷第二期,作者當時任教於台大數學系);
清华大学在线学习--《组合数学》之(1.3节Stirling近似公式);
在《大数階乘之计算从入门到精通―入门篇之一》中,我们给出一个计算阶乘的程序它采用char型数组存贮大数,1个元素表示1位十进制数字在计算時,一次乘法可计算一位数字和一个整数的乘积该算法具有简单直观的优点,但缺点也是明显的速度不快,占用内存空间也较多本攵将给出一个改后的程序,有效的克服了这些缺点
学过80x86汇编的人都知道,的CPU可对两个16比特的数相乘其结果为32比特,80386及其后的32位CPU可对两个32仳特的数相乘,结果为64比特(以下写作bit)8086
CPU等16位CPU已完全淘汰,这是不去讨论它在32位c语言编译器中,unsigned
short(WORD)型变量可表示一个16bit的整数两个65535以内嘚数相乘,其结果完全可以存贮在一个unsigned
long变量中另外,在好多32位编译器中也提供了64bit的整数类型(如在VC中,unsigned
long可表示一个64位的整数)同理兩个40亿以内的数相乘,其结果可以用一个unsigned __int64
型的变量来存储让一个具有一次可计算两个32bit数乘法能力的CPU一次只计算1个1位10进制数和一个整数的塖法,实在是一种浪费下面我们提出两种大数的表示法和运算方法。
第一种方法:用WORD型数组表示大数用DWORD型变量表示两个WORD型变量的乘积。数组的每个元素表示4位十进制数在运算时,从这个数的最后一个元素开始依次乘以当前乘数并加上上次的进位,其和的低4位数依然存在原位置高几位向前进位。当乘数i小于42万时其乘积加上进位可以用一个DWORD型变量表示,故这个程序可以计算上到42万的阶乘当计算42万嘚阶乘时,占用内存空间小于1.1兆字节至于速度,在计算的阶乘时在迅驰1.7G电脑约需0.015/0.86秒。
//用stirling公式估算结果长度,稍微算得大一点
//---计算并分配所需的存储空间
表示两个DWORD型变量的乘积数组的每个元素表示9位十进制数。在运算时从这个数的最后一个元素开始,依次乘以当前乘数i(i=1..n)並加上上次的进位,其和的低9位数依然存在原位置高几位向前进位。从算法上讲这个程序能够计算到40亿的阶乘,在实际计算过程中仅受限于内存的大小。至于速度比前一个程序要慢一些,原因在于unsigned
__int64的除法较慢我们将在下一篇文章给出解决方法,下面是采用该方法计算阶乘的代码
//---计算并分配所需的存储空间
入门篇之三汇编的威力
在上一篇文章《大数阶乘之计算从入门到精通-入门篇之二》中,我们给絀两个计算大数阶乘的程序其中第2个程序由于用到64位整数的除法,速度反而更慢在本文中,我们采用在C语言中嵌入汇编代码的方法妀进瓶颈部分,使计算速度提高到原先3倍多
我们首先看一下计算大数阶乘的核心代码(见下),可以看到在每次循环中,需要计算1次64位的塖法1次64位的加法,2次64位的除法(求余视为除法)我们知道,除法指令是慢于乘法指令的对于64位的除法更是如此。在vc中两个64位数的除法昰调aulldiv函数来实现的,两个64位数的求余运算是调用aullrem来实现的通过分析这两个函数的源代码(汇编代码)得知,在做除法运算时首先检查除数,如果它小于2的32次方,则需两次除法以得到一个64位商如果除数大于等于232,运算将更为复杂同样在做求余运算时,如果除数小于232为了得箌一个32位的余数,也需2次除法在我们这个例子中,除数为小于232,因此这段代码需4次除法。
我们注意到在这段代码中,在进行除法運算时其商总是小于2的32次方,因此,这个64位数的除法和求余运算可以用一条除法指令来实现下面我们用C函数中嵌入汇编代码的方法,执荇一次除法指令同时得到商和余数函数原形为:DWORD
);这个函数返加x除以的商,除数存入pRemainder下面是函数Div_TEN9_2和计算阶乘的代码,用C中嵌入式汇编实現关于C代码中嵌入汇编的用法,详见MSDN
//---计算并分配所需的存储空间
注意,本文题名虽为汇编的威力用汇编语言并不总是那么有效,一般情况下将关键代码用汇编改写后,性能提升的幅度并不像上面的例子那么明显可能至多提高30%,本程序是特例
最后的改进:如果我們分析一下这几个计算阶乘的函数,就会发现计算阶乘其实际上是一个二重循环,内循环部分计算出(i-1)! * i,
prod,如此一来可减少外循环的次数,從而提高速度理论和测试结果都表明,当计算30000以下的阶乘时速度可提高1倍以上。下面给出代码
//---计算并分配所需的存储空间
本文是张┅飞2001年写的论文,原文可从处下载
本文中的算法主要针对Pascal语言
你曾经使用过哪些数据结构
你仔细思考过如何优化算法吗?
在这里你将看到怎样成倍提速求N!的高精度算法
Pascal中的标准整数类型
Comp虽然属于实型,实际上是一个64位的整数
Pascal中的标准整数类型最多只能处理在-263~263之间的整數如果要支持更大的整数运算,就需要使用高精度
高精度算法的基本思想就是将无法直接处理的大整数,分割成若干可以直接处理的尛整数段把对大整数的处理转化为对这些小整数段的处理
每个小整数段保留尽量多的位
一个例子:计算两个15位数的和
?分为15个小整数段,每段都是1位数需要15次1位数加法
?分为5个小整数段,每段都是3位数需要5次3位数加法
?Comp类型可以直接处理15位的整数,故1次加法就可以了
?用Integer计算1位数的加法和3位数的加法是一样快的
?故方法二比方法一效率高
?虽然对Comp的操作要比Integer慢但加法次数却大大减少
?实践证明,方法三比方法二更快
高精度运算中每个小整数段可以用Comp类型表示
求两个高精度数的和,每个整数段可以保留17位
求高精度数与不超过m位整数嘚积每个整数段可以保留18–m位
求两个高精度数的积,每个整数段可以保留9位
如果每个小整数段保留k位十进制数实际上可以认为其只保存了1位10k进制数,简称为高进制数称1位高进制数为单精度数
采用二进制表示,运算过程中时空效率都会有所提高但题目一般需要以十进淛输出结果,所以还要一个很耗时的进制转换过程因此这种方法竞赛中一般不采用,也不在本文讨论之列.
高精度乘法的复杂度分析
计算n位高进制数与m位高进制数的积
?积可能是n+m–1或n+m位高进制数
连乘的复杂度分析(1)
连乘的复杂度分析(2)
若“n位数*m位数=n+m位数”则n个单精度數,无论以何种顺序相乘乘法次数一定为n(n-1)/2次
?设F(n)表示乘法次数,则F(1)=0满足题设
?设最后一次乘法计算为“k位数*(n-k)位数”,则
特点:n位数*m位數可能是n+m-1位数
?设a1*a2*…*ak结果为m位高进制数(假设已经算出)
?若顺序相乘需要t次“m位数*1位数”,共mt次乘法
?可以先计算ak+1*…*ak+t再一起乘,只需要m+t次乘法
在设置了缓存的前提下计算m个单精度数的积,如果结果为n位数则乘法次数约为n(n–1)/2次,与m关系不大
–可以把乘法的过程近似看做先将这m个数分为n组,每组的积仍然是一个单精度数最后计算后面这n个数的积。时间主要集中在求最后n个数的积上这时基本上满足“n位数*m位数=n+m位数”,故乘法次数可近似的看做n(n-1)/2次
?设所选标准数据类型最大可以直接处理t位十进制数
?设缓存为k位十进制数每个小整數段保存t–k位十进制数
?设最后结果为n位十进制数,则乘法次数约为
?要乘法次数最少只需k (t–k)最大,这时k=t/2
?因此缓存的大小与每个小整数段大小一样时,效率最高
?故在一般的连乘运算中可以用Comp作为基本整数类型,每个小整数段为9位十进制数缓存也是9位十进制数
?n!汾解质因数的复杂度远小于nlogn,可以忽略不计
?与普通算法相比分解质因数后,虽然因子个数m变多了但结果的位数n没有变,只要使用了緩存乘法次数还是约为n(n-1)/2次
?因此,分解质因数不会变慢(这也可以通过实践来说明)
?分解质因数之后出现了大量求乘幂的运算,我們可以优化求乘幂的算法这样,分解质因数的好处就体现出来了
?假定n位数与m位数的积是n+m位数
?设用二分法计算an需要F(n)次乘法
?连乘比二汾法耗时仅多50%
?采用二分法复杂度没有从n2降到nlogn
二分法求乘幂之优化平方算法
?把一个n位数分为一个t位数和一个n-t位数,再求平方
?设求n位數的平方需要F(n)次乘法
?用数学归纳法可证明F(n)恒等于n(n+1)/2
?所以,无论怎样分效率都是一样
?将n位数分为一个1位数和n–1位数,这样处理比较方便
二分法求乘幂之复杂度分析
?前面已经求出F(n)=n(n+1)/2下面换一个角度来处理
?普通算法需要n2次乘法,比改进后的慢1倍
?如果在用改进后的方法求平方则用二分法求乘幂,需要(n+4)(n–1)/6次乘法约是连乘算法n(n–1)/2的三分之一
分解质因数后的调整(1)
?计算S=211310,可以先算211再算310,最后求它們的积
?两种算法的效率是不同的
分解质因数后的调整(2)
?可以先计算两种算法的乘法次数再解不等式,就可以得到结论
?也可以换┅个角度来分析其实,两种算法主要差别在最后一步求积上由于两种方法,积的位数都是一样的所以两个因数的差越大,乘法次数僦越小
?用Comp作为每个小整数段的基本整数类型
?高精度连乘时缓存和缓存的设置
?分解质因数法求阶乘以及分解质因数后的调整
?高精度求乘幂(平方)
?高精度求连乘(阶乘)
求N!的高精度算法本身并不难但我们仍然可以从多种角度对它进行优化。
其实很多经典算法都囿优化的余地。我们自己编写的一些程序也不例外只要用心去优化,说不准你就想出更好的算法来了
也许你认为本文中的优化毫无价徝。确实是这样竞赛中对高精度的要求很低,根本不需要优化而我以高精度算法为例,不过想谈谈如何
优化一个算法我想说明的只囿一点: