用Stirling逼近近似计算阶乘的探讨与应用
江苏省赣榆高级中学 仲晨
【关键词】：Stirling逼近，阶乘，极限论，微积分，数学实验，计算机算法
“阶乘”（factorial）在信息学竞赛中具有重要角色，更广泛的说，“阶乘”在数学领域也是占有重要地位。在许多人刚刚学习计算机语言的时候，大多会被要求写一个算阶乘的程序，而在学习高精度算法的时候，也会写一个计算较大数字阶乘的程序。不过，在实际的运用之中，可能遇到更大数字的阶乘计算和不同要求的阶乘结果，例如：TOJ（同济大学ACM网络题库，http://acm.tongji.edu.cn/problem.php）的1016题——“求N！左边第二位的数字”，这就需要一定的精度思考了。
可是我们通常对于较大数字阶乘的要求是求结果位数或前几位数字，这怎么办呢？
在刘汝佳的《算法艺术与信息学竞赛》一书中，（Page241）介绍了Stirling公式：
其中的～符号是指“同阶”或“相当”，即两者随n增加的大致速度相同，在n较大时，两者极其相近。
这是一个极限的概念（现行教材高二下学期数学内容），属于微分学内容，准确写法为：
但遗憾的是在《算法艺术与信息学竞赛》书中只提供了这个算式，并无他物！
本人近日看到一本数学科普读物——《好玩的数学——不可思议的e》（陈任政著，科学出版社），其中5.12节简介了Stirling逼近近似算阶乘，本人感到好奇，于是对这种算法的具体步骤进行了分析，并研究了它的精确度，故为本文。在2005年8月7日完工之日，笔者上网搜索了一下，找到了一些关于Stirling逼近的文章，偶然地在臺灣亞洲聚合公司蔡永裕《談Stirling公式的改良》（刊自台湾《數學傳播》20卷4期，民国85年12月）一文中找到同感，蔡先生的做法于笔者方向不同，作出来的结果比笔者的算法精确一个数量级左右，惭愧，于是，笔者又再次研究，寻找更好算法，写于本文后部。
在1730年，棣莫弗(棣,音Dì)（法国数学家，Abraham De Moiver，1667～1754）发表的《分析杂论》中首先对n!地一个无穷级数展开式给出了近似公式：
但是，现在我们叫这个式子为“Stirling逼近”，中文叫做“斯特林逼近”，这是为什么呢？
因为棣莫弗的朋友苏格兰数学家斯特林(James Stirling,1696~1770)在同年的《微分法或无穷级数的简述》中也给出了等价的级数。
事实上，棣莫弗首先得到的式子是
，
但是，他没有把C求出来。而斯特林则利用棣莫弗的发现做了探讨，求出了
。
这些式子的来源是一个无穷级数展开式：
其中B2=1/6，B4=-1/30，B6=1/42
… B2k是雅格布·伯努力数。（具体内容请参见后文介绍）
这里介绍一下，还没上高中的同学还没有学到，“乘方”的逆运算有两种：开方和对数。
对于一个幂：
，其中a成为底数，n成为指数，b成为幂。已知a和n求b，就是乘方运算；已知b和n，求a，就是开方运算；而已知a和b求n，就是对数运算，写做：
，这里n就称为以a为底b的对数（logarithm）。
当底数为10的时候，叫做常用底数，简写做
；当底数为e的时候，叫做自然对数，简写做
。
至于e的含义：e是重要性仅次于π的数，是极限论的主要内容，具体的说，即：
意思是当n趋向于正的无限大的时候，
趋向于e 。
e是无理数，即无限不循环小数，也是超越数，即不能满足某个整数系数代数方程的数（不能满足某个整数系数代数方程的数叫做代数数）。
目前e只算到了几千位。
e＝2.718281828459045235360287471352662497757247093...
特别说明的是，在Pascal语言中，exp(n)函数就是e的n次方。
另外，有个著名的公式被成为“整个数学中最卓越的公式”：
其中的i为虚数的单位，
。来自算术的0、1，来自代数的i，来自几何的π，来自分析学的e，奇妙的组成了一个公式！
这是欧拉（瑞士数学家，Leonhard Euler，1707~1783）发现的！所以称作“欧拉公式”。
不过，真正的欧拉公式是：
那个“最卓越的公式”只是欧拉公式的一个推倒公式。
言归正传，由公式
两边同时取e的幂，得
即：
再经过近似等处理（这些处理比较麻烦，而且牵扯到微积分内容）,我们得到了Stirling公式：
注：在本文后部有Stirling公式的推倒过程。
当然，我们可以得到它得更具体形式：
其中的θ是不定的，在（0,1）区间之内。
讲到这里，我们有了公式，可以开始计算了。
但是，我们计算什么呢？难道我们要计算N！的值吗？
虽然这个式子比阶乘原始计算式简单，但是实际计算的时候仍然得到的将是上百上千位的数字，而且这个式子毕竟是近似公式，无法真正得到确切得数！
难道我们辛苦的到的式子无用了吗？
答案当然是否定的！我们可以求N！的位数！
求位数的方法是取常用对数（以10为底的对数）。那，这是为什么呢？看看下表：
|
n |
|
n位数 |
|
1 |
0.000000 |
1 |
|
10 |
1.000000 |
2 |
|
100 |
2.000000 |
3 |
|
1000 |
3.000000 |
4 |
|
10000 |
4.000000 |
5 |
|
100000 |
5.000000 |
6 |
|
1000000 |
6.000000 |
7 |
|
10000000 |
7.000000 |
8 |
|
100000000 |
8.000000 |
9 |
好了！我们知道了n的位数等于
，对吗？
不对！
不一定是整数，比如
，所以，我们得到的公式是：
n的位数＝
其中
是取整符号，指不大于n的最大整数，在Pascal语言中为trunc(n)函数。
如果我们用Stirling逼近算出n！再算它的位数就太不值得了，这里我们就要运用对数的一些特性来简化公式。
我们两边取常用对数：
进而化简为:
现在我们可以来求值了！
以求1000！的位数为例，
所以1000！的位数为2568位！
我们可以看到，这里的计算复杂度很低，都是O(n) 级别的！对于其计算机编程运用来说，计算过程中数字不大不会溢出，取对数后的精确度也可以保证。
这样我们就解决了计算阶乘位数问题！而且可以绝对放心，这个公式在位数方面的准确率是99.9999999999999％以上。(而且这个仅有的可能性也只是从正态性来推测的)
用Pascal语言来书写：
Trunc(0.5*Ln(2*n*3.1415926) /Ln(10) +n*Ln(n/Exp(1)) /Ln(10))+1
建议使用分步计算，这样能避免计算过程中的溢出现象。
这样，笔者使用批处理计算了几个值：
|
n |
n!位数 |
|
1 |
1 |
|
10 |
7 |
|
100 |
158 |
|
1000 |
2568 |
|
10000 |
35660 |
|
100000 |
456574 |
|
1000000 |
5565709 |
|
10000000 |
65657060 |
|
100000000 |
756570557 |
|
1000000000 |
8565705523 |
由于Longint范围的限制,无法计算再大的数值，但是可以通过使用int64或者extended来解决。
不过，我们有点不甘心，只计算位数难免用途不广，其实，我们可以用来计算n！的前几位！
什么原理呢？
例如计算1000！时，我们得到的数为2567.6046080272230971441871361093945……
而我们从下表中可以看出点东西来:
|
n |
|
|
5 |
0.698970004336019…… |
|
50 |
1.698970004336019…… |
|
500 |
2.698970004336019…… |
|
5000 |
3.698970004336019…… |
我们可以得知：一个数增大10倍，它的常用对数整数部分将增加1，而小数部分不变。
我们得到的更重要的结果是：一个数的常用对数的小数部分所对应的幂的各位数字（不考虑小数点）与原数各位数字相同。
所以，对于1000！来说:
2567.6046080272230971441871361093945……的小数部分为
0.6046080272230971441871361093945……，
它的10为底的幂为4.0235372577197005393836315765635……，
也就是说1000！这个2568位数的前几位为40235372577197005393836315765635……。
Well！我们可以求出一个阶乘的前几位数了！
但是，不要高兴太早！这里的精度无法保证！
下面让我们来看看精确度方面：
特别说明的是，这里说的精确度是指n!与Stirling逼近公式的精确度，而不是位数计算的精确度，因为位数计算基本上是精确的！
|
n |
逼近值 |
n！准确值 |
相对误差 |
|
10 |
3598695.619 |
3628800 |
0.008365359 |
|
20 |
2.42279E+18 |
2.4329E+18 |
0.004175011 |
|
30 |
2.64517E+32 |
2.65253E+32 |
0.002781536 |
|
40 |
8.14E+47 |
8.16E+47 |