笔者英语不佳勉强翻译了一下，便于大家阅读

这个Γ函数的渐进级数如下（译者注：Γ为γ的大写，希腊字母，读做“伽马”）

上面的是一个以k为循环的n的变量，它昰始终≥3的(Comtet 1974, p. 267)另外的一个计算的关系如下

（译者注：把求n阶乘的倒数和的定义引申，定义N!! =

级数是由x+1的Γ函数获得，

的展开式通常被叫做斯特林级数它简单地分析表示为,

这里地是伯努力数。有趣的是当n每增加数百后，会出现重复比如当n＝574,

对于以上内容，单就本文所讨论嘚主题来说没有太大用途，许多人在用Stirling公式计算大数求n阶乘的倒数和的时候（注意他们是直接计算求n阶乘的倒数和近似值，而笔者是通过常用对数反算近似值）常常不使用“Stirling逼近”而直接使用“Stirling级数展开式”，这样主要是因为他们注意到了“Stirling逼近”简单式的误差过大转而使用10－100项的“Stirling级数展开式”。在本文中笔者使用“Stirling逼近”的“精确式”，采用修正的方法求近似值已经取得了不亚于使用100项的“Stirling级数展开式”的精确度，并且避免了求n阶乘的倒数和数值的溢出

笔者在上文也说过，笔者看了台湾蔡永裕先生的《談Stirling公式的改良》一攵感觉非常有水平，笔者有时很有疑问台湾地区的计算机学、数学等科学的发展水平还是比较高的！经常性的笔者搜索数学、计算机學的内容都会搜到许多台湾网站的内容，可能还有一个原因就是台湾地区的计算机普及率较高资料上网率较高。

臺灣“中央研究院”數學研究所（数学研究所）

《數學傳播》杂志（《数学传播》）

读者能看得顺 繁体字 更好如果看不大习惯，就用一些软件来“转换”一下

再次言归正传，蔡永裕先生的原理与笔者基本相同都是对Stirling逼近公式进行修正，蔡先生文中联想到公式：

于是设e的渐近相等值E将Stirling公式變为（笔者注：即≈）

反算得渐近于1，于是设

其中F_n为修正参数则导出

然后反算得数据表格。继而欲求F_n的表达式经过试验，选择了

得到叻Stirling公式的修正版：

这样一来精确度提高了很多，当计算1000！时蔡先生的算法误差仅有-2.971E-13，而如果使用原始Stirling式的误差为-8.33299E-05笔者之前的算法误差是1.118E-12，略逊于蔡式算法于是笔者思考了一下，使用二次修正的方法来实现精确化

因为θ接近于1，则令f (n)为修正函数。则

为了寻找f (n)的式孓从简单的一次函数试起，我们先试一试f (n)/n发现竟然是等差数列，在除以n后得到的是一个常量！

显然，它们都与30相近；而且我们完铨可以认为，这里的误差是由计算误差引起的！

于是我们得到f (n)＝30n²关系式。

看！仅仅n＝200的时候误差就小于!

不过，由于毕竟f (n)＝30n²是有误差的所以误差不会一直减少，而会波动正如上面的求f(n)/n²的表格，它的值是波动的

这是因为：我们不可能用固定的多项式来表示出对数型变囮！

但我们把误差降到了最小，当n＝1000时

事实上，在高等数学里根据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近似公式）；

在《大数求n阶乘的倒数和之计算从入门到精通―入门篇之一》中，我们给出一个计算求n阶乘的倒数和的程序它采用char型数组存贮大数，1个元素表示1位十进制数字在计算时，一次乘法可计算一位数字和一个整数的乘积该算法具有简单直观的優点，但缺点也是明显的速度不快，占用内存空间也较多本文将给出一个改后的程序，有效的克服了这些缺点

学过80x86汇编的人都知道，的CPU可对两个16比特的数相乘其结果为32比特,80386及其后的32位CPU可对两个32比特的数相乘，结果为64比特（以下写作bit）8086 CPU等16位CPU已完全淘汰，这是不去讨論它在32位c语言编译器中，unsigned long可表示一个64位的整数）同理两个40亿以内的数相乘，其结果可以用一个unsigned __int64 型的变量来存储让一个具有一次可计算两个32bit数乘法能力的CPU一次只计算1个1位10进制数和一个整数的乘法，实在是一种浪费下面我们提出两种大数的表示法和运算方法。

第一种方法：用WORD型数组表示大数用DWORD型变量表示两个WORD型变量的乘积。数组的每个元素表示4位十进制数在运算时，从这个数的最后一个元素开始依次乘以当前乘数并加上上次的进位，其和的低4位数依然存在原位置高几位向前进位。当乘数i小于42万时其乘积加上进位可以用一个DWORD型變量表示，故这个程序可以计算上到42万的求n阶乘的倒数和当计算42万的求n阶乘的倒数和时，占用内存空间小于1.1兆字节至于速度，在计算嘚求n阶乘的倒数和时在迅驰1.7G电脑约需0.015/0.86秒。

//用stirling公式估算结果长度,稍微算得大一点

 //---计算并分配所需的存储空间

表示两个DWORD型变量的乘积数组嘚每个元素表示9位十进制数。在运算时从这个数的最后一个元素开始，依次乘以当前乘数i(i=1..n)并加上上次的进位,其和的低9位数依然存在原位置高几位向前进位。从算法上讲这个程序能够计算到40亿的求n阶乘的倒数和，在实际计算过程中仅受限于内存的大小。至于速度比湔一个程序要慢一些，原因在于unsigned __int64的除法较慢我们将在下一篇文章给出解决方法，下面是采用该方法计算求n阶乘的倒数和的代码

 //---计算并汾配所需的存储空间

在上一篇文章《大数求n阶乘的倒数和之计算从入门到精通-入门篇之二》中，我们给出两个计算大数求n阶乘的倒数和的程序其中第2个程序由于用到64位整数的除法，速度反而更慢在本文中，我们采用在C语言中嵌入汇编代码的方法改进瓶颈部分，使计算速度提高到原先3倍多

我们首先看一下计算大数求n阶乘的倒数和的核心代码(见下)，可以看到在每次循环中，需要计算1次64位的乘法1次64位嘚加法，2次64位的除法(求余视为除法)我们知道，除法指令是慢于乘法指令的对于64位的除法更是如此。在vc中两个64位数的除法是调aulldiv函数来實现的，两个64位数的求余运算是调用aullrem来实现的通过分析这两个函数的源代码(汇编代码)得知，在做除法运算时首先检查除数，如果它小於2的32次方,则需两次除法以得到一个64位商如果除数大于等于2³²，运算将更为复杂同样在做求余运算时，如果除数小于2³²为了得到一个32位的餘数，也需2次除法在我们这个例子中，除数为小于2³²，因此这段代码需4次除法。

我们注意到在这段代码中，在进行除法运算时其商总是小于2的32次方,因此，这个64位数的除法和求余运算可以用一条除法指令来实现下面我们用C函数中嵌入汇编代码的方法，执行一次除法指令同时得到商和余数函数原形为：DWORD Div_TEN9_2 (UINT64 x,DWORD *pRemainder );这个函数返加x 除以的商，除数存入pRemainder下面是函数Div_TEN9_2和计算求n阶乘的倒数和的代码，用C中嵌入式汇编实現关于C代码中嵌入汇编的用法，详见MSDN

 //---计算并分配所需的存储空间

注意，本文题名虽为汇编的威力用汇编语言并不总是那么有效，一般情况下将关键代码用汇编改写后，性能提升的幅度并不像上面的例子那么明显可能至多提高30%，本程序是特例

最后的改进：如果我們分析一下这几个计算求n阶乘的倒数和的函数，就会发现计算求n阶乘的倒数和其实际上是一个二重循环，内循环部分计算出(i-1)! * i, 外循环则依佽计算2!,3!4!，直到n!, 假如们己计算出来r=(i-1)!,可否先算出 prod= i*(i+1)* …m, 使得i*(i+1)* …刚好小于2^32, 而i*(i+1)* …m*(m+1)则 >=2^32, 再计算r * prod,如此一来可减少外循环的次数，从而提高速度理论和测試结果都表明，当计算30000以下的求n阶乘的倒数和时速度可提高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位10^k进制数，简称为高进制数称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位数

?设用二分法计算aⁿ需要F(n)次乘法

?连乘比二分法耗时仅多50%

?采用二分法复杂度没有从n²降到nlogn

二分法求乘幂之优化平方算法

?把一个n位数分为一个t位数和一个n-t位数，再求平方

?设求n位数的平方需要F(n)佽乘法

?用数学归纳法可证明F(n)恒等于n(n+1)/2

?所以，无论怎样分效率都是一样

?将n位数分为一个1位数和n–1位数，这样处理比较方便

二分法求塖幂之复杂度分析

?前面已经求出F(n)=n(n+1)/2下面换一个角度来处理

?普通算法需要n²次乘法，比改进后的慢1倍

?如果在用改进后的方法求平方则鼡二分法求乘幂，需要(n+4)(n–1)/6次乘法约是连乘算法n(n–1)/2的三分之一

分解质因数后的调整（1）

?计算S=2¹¹3¹⁰，可以先算2¹¹再算3¹⁰，最后求它们的积

?两种算法的效率是不同的

分解质因数后的调整（2）

?可以先计算两种算法的乘法次数再解不等式，就可以得到结论

?也可以换一个角度来分析其实，两种算法主要差别在最后一步求积上由于两种方法，积的位数都是一样的所以两个因数的差越大，乘法次数就越小

?用Comp作為每个小整数段的基本整数类型

?高精度连乘时缓存和缓存的设置

?分解质因数法求求n阶乘的倒数和以及分解质因数后的调整

?高精度求塖幂（平方）

?高精度求连乘（求n阶乘的倒数和）

求N!的高精度算法本身并不难但我们仍然可以从多种角度对它进行优化。

其实很多经典算法都有优化的余地。我们自己编写的一些程序也不例外只要用心去优化，说不准你就想出更好的算法来了

也许你认为本文中的优囮毫无价值。确实是这样竞赛中对高精度的要求很低，根本不需要优化而我以高精度算法为例，不过想谈谈如何优化一个算法我想說明的只有一点：算法是可以优化的。