摘要

本文分为四个部分:
第1部分介绍大数阶乘的一般计算方法;
第2部分讨论阶乘位数的精确计算和估计;
第3部分介绍如何计算阶乘末尾的0的个数;
第4部分讨论大数阶乘计算的优化方法。

阶乘是什么我想就不用我说了,既然你对这篇文章感兴趣,那么我假设你已经会了简单的阶乘计算,并且想进一步研究一下大数阶乘的计算。

研究阶乘有什么用呢?说实话,对于大数阶乘精确计算的应用,笔者还未曾有所了解,或许它本来就没什么实际用处。那为什么还要研究它呢?可以这样问一下你自己,“为什么要学习?”就像有人地日复一日、乐此不疲地计算着圆周率$\pi$的第$n$位小数一样,研究大数阶乘的计算无非有两个目的:

  1. 为着自己的乐趣并想看下能否从中获得一些新发现。
  2. 在实践中学习,在学习中实践,检验一下自己对于某些算法的理解程度

大数阶乘的简单计算

C语言的整数类型可表示数据范围有限,对于64位long long类型,超过$2^{63}-1$就会发生溢出,那么如何用C语言来计算10000的阶乘呢?下面提供一个笔者用C语言写的打印10000以内的数字的阶乘的程序。在下面的代码中对于大数运算的进位的处理,笔者采用了10亿进制($10^9$),10000的阶乘的十进制位数为35660,除以9稍微小于4000,是故笔者使用了一个长度为4000的int数组。

打印10000的阶乘

#include <stdio.h>

// 采用十亿进制,每个int表示存储9位十进制
#define scale 1000000000
#define print_cell(cell) printf("%09d", cell)

void print_factorial(const int N)
{
    int i, j, len;
    int fac[4000] = {0};
    long long carry, factor;
    j = len = 0; fac[0] = 1;

    for (factor = 2; factor <= N; ++factor) {
        for (carry = 0, i = j; i <= len; ++i) {
            carry += factor * fac[i];
            fac[i] = carry % scale;  // 取余数
            carry /= scale;          // 进位
        }
        if (fac[j] == 0) ++j;
        if (carry > 0) {
            fac[i] = carry;
            ++len;
        }
    }
    printf("%d", fac[len]);
    while (--len >= 0)print_cell(fac[len]);
    printf("\n");
}

int main()
{
    // for (int i = 0; i < 40; ++i)
        // print_factorial(i);
    print_factorial(10000);

    return 0;
}

在上面给的代码中,我们把一个int作为一个单元来处理,在现代的计算机上,int的一般表示范围为$-2147483648(-2^{31})\sim2147483647(2^{31}-1)$,最大可以表示$10^9$量级的整数,采用$10^9$进制的目的是减少循环的次数,提高效率。

代码中采用一个数组来表示一个大整数,按$10^9$进制倒着存储,向后进位,高位向后存储。也即是说,数组的第一个元素是所表示的大整数的末位(十进制表示的最后几位)。第一层循环的作用是不断地将每一个因子乘到数组表示的大整数中去,第二层循环是代码的核心部分,处理一个因子与数组表示的大整数的乘法,其中的关键部分是取余和进位。第二层循环后面的两个if语句的作用分别是

  1. 跳过阶乘末尾的0(将任何数乘到0中还是0)。
  2. 计算数组表示的大整数的位数。

阶乘的位数

随着N的增大,$N!$的增长速度比指数函数还快。定义函数$D(N)$表示$N!$的十进制位数,那么
$$
\begin{align*}
D(N) &= \lceil\log_{10}(1×2× \dots ×N-1×N)\rceil \\
&= \lceil\sum_{n=1}^{N}{\log_{10}(n)}\rceil
\end{align*}
$$
我们以$D(10^7)$为例,如果你会一点简单的编程并且此刻你身边有一台可编程的计算机,计算这个数字不是什么难事。 经过笔者编程计算,$D(10^7)$的精确值为65657060。 这意味着,一千万的阶乘是一个六千多万位的数字!事实上,当N很大时,我们可以用积分代替求和近似计算
$$
\begin{align*}
D(N) &= \lceil\sum_{n=1}^{N}{\log_{10}(n)}\rceil\\
&≈ \int_1^{N}\log_{10}(x)dx\\
&= {\frac{x(\ln(x) - 1)}{\ln(10)}}{\Bigg|}_1^N \\
&= {\frac{N(\ln(N) - 1)}{\ln(10)}} + 0.4343
\end{align*}
$$
那么,一个简单的计算器就可以得出结果
$$
\begin{align*}
D(10^7) &≈ {\frac{10^7(\ln(10^7) - 1)}{\ln(10)}} + 0.4343\\
&≈ 65657055.6153
\end{align*}
$$
与精确值相差无几。

事实上,借助斯特林公式,可以精确计算阶乘的位数。斯特林公式如下:
$$
N!=N^Ne^{-N}\sqrt{2\pi N}=(\frac{N}{e})^N\sqrt{2\pi N}
$$
运用斯特林公式,不难写出$N!$的位数为:
$$
\begin{align*}
D(N) &= \lceil\log_{10}{(N^Ne^{-N}\sqrt{2\pi N})}\rceil\\
&=\lceil\frac{\ln{(N^Ne^{-N}\sqrt{2\pi N})}}{\ln(10)}\rceil\\
&=\lceil\frac{N(\ln(N)-1) + \frac{1}{2}\ln(2\pi N)}{\ln(10)}\rceil
\end{align*}
$$
将$N=10^7$带入,可得$D(10^7)=65657060$,为1千万的阶乘的精确位数。

如果你嫌1千万这个数字太大,不妨试着将10,100等数字代入进行计算,亦可得到精确值(事实上,这一公式除了在$N=1$处,都能给出$N!$位数的准确值)。

这个公式有点类似于素数计数公式$\pi(x)$,通过实数计算来获取一个整数值。在《素数之恋·黎曼和数学中最大的未解之谜》一书中有$\pi(x)$精确的计算公式:
$$
\begin{align*}
\pi(x)=\sum_{n}{\frac{\mu(n)}{n}}J(\sqrt[n]{X})
\end{align*}
$$
$\pi(x)$的这一公式好像是黎曼给出的,这个公式比上面的阶乘位数公式强大深奥太多,需要复变函数、黎曼$\zeta$函数、莫比乌斯反演等知识基础。$\pi(N)$的计算结果刚好是一个整数,不需取整即可得到$N$以内素数的个数。鄙人愚笨,还未真正弄懂,所以暂不作深入探讨。

阶乘末尾0的个数

这一小节将介绍整数$N$的阶乘的末尾0的个数的计算方法。
$$
\lim_{N->\infty}\sum_{n=1}^{5^n\leqslant N} = \frac{1}{4}
$$
待絮!

阶乘计算的优化方法

乘法·作为基本的运算

乘法作为阶乘的基本运算,其算法必须高效,不然阶乘的算法设计得再怎么好,阶乘的速率也没有多大的提升空间。

乘法运算的普通算法时间复杂度为$O(M×N)$,空间复杂度为$O(M+N)$,其中$M$和$N$分别为两个乘数的位数。这样做乘法效率是极低的,尤其是遇到两个位数相近的整数。比如,拿两个位数$10^6$(100万)的整数相乘,你要花费$10^{12}$量级的计算次数才能得出一个只不过两百万位上下的整数!这在普通的个人计算机上至少要花费几十分钟的计算时间,是无法忍受的!

两个数相乘如此,何况是1千万个数相乘呢!或许你会说,我们可以预先分配好一段内存空间来表示作为阶乘计算结果的大整数,然后不断地将$1~10^7$中的每一个数字作为因子,乘到表示阶乘的大整数中去(上面打印1万的阶乘的代码正是这么做的)。不过,下面的计算结果将指出,即使采用了这种方法,用这种时间量级的算法来计算1千万的阶乘也不具备可行性,况且这种方法一点也不灵活,计算结果也不可复用!

根据上面的结论,$N!$的位数约为

$$
cN(\ln(N)-1)
$$

整数$n$的位数约为$c\ln(n)$,其中$c$为常数。

那么,将$n$乘到$(n-1)!$所需的计算次数约为
$$
c^2\ln(n)(\ln(n-1)-1)(n-1)
$$
于是计算$N!$所需的计算次数约为
$$
T(N) = c^2\sum_{n=2}^{N}\ln(n)(\ln(n-1)-1)(n-1)
$$

将$n$近似代替$n-1$并用积分代替求和,得
$$
\begin{align*}
T(N) &= c^2\sum_{n=2}^N[{n\ln^2(n)}-n\ln(n)]\\
&\approx c^2\int_2^N[x\ln^2(x)-x\ln(x)]dx\\
&\approx c^2N^2[\frac{1}{2}\ln^2(N)-\ln(N)+\frac{1}{2}]\\
&\approx c’N^2[\ln(N)-1]^2
\end{align*}
$$
取主项得
$$
T(N)=O(N^2\ln^2(N))
$$
如果采用该方法计算1千万的阶乘,那么$T(10^7)$的量级将在$10^{15}$以上。对于频率为几个$G$的个人电脑来说,这么大的计算量不花上几天的时间是搞不定的!

想要提高我们的阶乘算法的效率,作为基本运算的乘法的效率必须得提高。这时候就该快速傅里叶变换了(FFT,Fast Fourier Transformation)出场了。这里我们只说明,快速傅里叶变换算法作为大整数运算库用来加速乘法运算的必备手段,其时间复杂度为$O(n\ln(n))$,并不提供快速傅里叶变换的教程。关于快速傅里叶变换的教程,至少可以写一篇长长的文章。以下假设我们已经有了快速傅里叶变换加速的大数运算环境,并以此为基础来探讨和实现大数阶乘的高效算法。

素数筛选·阶乘加速之道

快速幂·减少乘法次数的关键

快速幂乘·进一步提升