渐近估阶的威力：欧拉-麦克劳林公式的应用

摘要: 用欧拉-麦克劳林公式对调和级数估阶

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

---

各位好，本文我们来讨论一下算法分析中的问题。

由于递归和迭代是最常见的两种程序结构。因此时间复杂度 $T(n)$ 最终往往归结到求解和式或递归式，例如：

• 素数筛的最坏时间复杂度最终归结为和式：

$$T(n) = n\sum\limits_{p}\frac{1}{p}$$

• 快速排序的平均时间复杂度最终归结为递归式：

$$T(n) = n + 1 + \frac{1}{n}\sum\limits_{j=1}\limits^{N}(T(j-1)+T(n-j))$$

一般来讲我们总是希望能够解出精确的结果。但很多时候得出精确结果很难，或者没必要，我们的关注点就会由求出精确解改为对解的渐近估计。

对于递归式来讲，某些递归式用生成函数表达之后，经过一些变换可以得到生成函数的微分方程或函数方程，这个方程很可能很难解出精确解，即使有精确解 $C(z)$，也很可能给不出 $C(z)$ 的展开中 $z^{n}$ 的系数的精确解。

对于和式来说，有些和式可以给出精确解，例如 $\sum\limits_{i=1}\limits^{n}i = \frac{1}{2}n(n+1)$，更多的给不出精确解，例如前面提到的 $\sum\limits_{p}\frac{1}{p}$。

因此，在对实际算法的分析中，最终很可能会归结到对某个量的渐近估阶，在文章 函数增长与渐近分析入门 中我们介绍了渐近分析中的一些基本概念。

本文我们通过一个例子，即 $\sum\limits_{n=1}\limits^{\infty}\frac{1}{n}$，来看一下渐近估阶中的一个重要方法也就是欧拉-麦克劳林公式的应用。推导出这个量的渐近估计，我们就可以进一步得到快速排序平均情况下的时间复杂度，也就是我们在八股文中背过的 $O(n\log n)$。

首先给出结论，当 $n\rightarrow \infty$ 时，$\sum\limits_{k=1}\limits^{n}\frac{1}{k}$ 的渐近估阶如下，其中 $\gamma$ 为欧拉常数：

$$\sum\limits_{k=1}\limits^{n}\frac{1}{k} = \ln n + \gamma + o(1)$$

下面我们通过渐近估阶的方法推导出这个结论。

本文涉及到一些数学知识，大致有生成函数、微分方程、泰勒级数、用积分逼近求和、伯努利多项式、分部积分、欧拉-麦克劳林公式。

---

问题背景：快速排序的时间复杂度

在前面的快速排序的递归式中，整理后得到 $nT(n) = n(n + 1) + 2\sum\limits_{j=1}\limits^{N}T(j-1)$。记 $T(n)$ 的生成函数为 $C(z) = \sum\limits_{n=0}\limits^{\infty}T(n)z^{n}$，两边同时乘以 $z^{n}$，然后对 $n$ 求和，得到：

$$\sum\limits_{n=1}\limits^{\infty}nT(n)z^{n} = \sum\limits_{n=1}\limits^{\infty}n(n+1)z^{n} + 2 \sum\limits_{n=1}\limits^{\infty} \sum\limits_{k=1}\limits^{n}T(k-1)z^{n}$$

由生成函数的性质，上式最终转化为了关于生成函数 $C(z)$ 的微分方程：

$$C'(z) = \frac{2}{(1-z)^{3}} + \frac{2C(z)}{1-z}$$

求解该微分方程，得到生成函数的精确解 $C(z) = \frac{2}{(1-z)^{2}}\ln \frac{1}{1-z}$。

由生成函数的性质，以及泰勒级数，可以得到生成函数 $C(z)$ 的展开，其 $z^{n}$ 项的系数为：$2(n+1)(H_{n+1} - 1)$。

其中 $H_{n+1} = \sum\limits_{k=1}\limits^{n+1}\frac{1}{k}$ 正是我们要研究的和式，称其为调和数。

通过积分逼近求和的方法，可以得到渐近估阶 $H_{n} = O(\ln n)$，由此已经可以知道快速排序平均情况下的时间复杂度 $T(n) = O(n\ln n)$ 。进一步，通过欧拉-麦克劳林公式，还可以得到对 $H_{N}$ 更精确的渐近估阶。

用积分逼近求和

一个自然的想法是用积分 $\int_{a}^{b}f(x)\mathrm{d}x$ 来逼近和式 $\sum\limits_{k=a}\limits^{b}f(k)$ ，但问题在于误差有多大，这个问题依赖于 $f(k)$ 有多光滑。

如果记 $\delta_{k} = \max\limits_{k\leq x < k+1}|f(x) - f(k)|$，可以得出总误差的大致估计：

$$\sum\limits_{k=a}\limits^{b}f(k) = \int_{a}^{b}f(x)\mathrm{d}x + \Delta, \quad |\Delta| \leq \sum\limits_{k=a}\limits^{b-1}\delta_{k}$$

若 $f(x)$ 在 $[a, b]$ 单调，那么各个小区间 $[k, k+1)$ 的误差项叠加到一起得到 $|\Delta| \leq |f(a) - f(b)|$，例如对于调和数，用积分逼近求和就可以给出以下估计：

$$H_{n} = \sum\limits_{k=1}\limits^{n}\frac{1}{k} = \int_{1}^{n}\frac{1}{x}\mathrm{d}x + \Delta = \ln n + \Delta$$

估计的误差 $|\Delta| \leq |f(a) - f(b)| = 1 - \frac{1}{N}$。这样我们实际上已经得到了 $H_{n}$ 的估阶，即 $H_{n} \sim \ln n$。

误差项如果想要估计的更精确，将取决于 $f(x)$ 的各阶导数，这就是欧拉-麦克劳林公式的动机，它也是渐近估阶中最强大的工具。

欧拉-麦克劳林公式

用积分逼近求和想要减少误差，有两种不同的路线，一个是区间固定，随着步长越来越小，和与积分的差值收敛到 0，类似于黎曼积分的思想。另一种步长固定，让积分区间越来越大，和与积分的差值将收敛到一个常数。

第一形式

考虑积分 $\int_{0}^{1}g(x)\mathrm{d}x$，在定积分的分部积分 $\int_{a}^{b}u\mathrm{d}v = uv\big|_{a}^{b} - \int_{a}^{b}v\mathrm{d}u$ 中，取 $u = g(x)$，$v = B_{1}(x) = x - \frac{1}{2}$，$a = 0, b = 1$。有：

$$\int_{0}^{1}g(x)\mathrm{d}x = (x - \frac{1}{2})g(x)\bigg|_{0}^{1} - \int_{0}^{1}(x - \frac{1}{2})g'(x)\mathrm{d}x$$

其中 $B_{1}(x) = 1 - \frac{1}{2}$ 为伯努利多项式。注意这里对 $v = B_{1}(x)$ 的取法，为什么要这么取，实际上是有线索的，感兴趣的可以参考《具体数学》Chap9。

在上式中取 $g(x) = f(x + k)$，得到：

$$\int_{k}^{k+1}f(x)\mathrm{d}x = \frac{1}{2}(f(x + k) + f(k)) - \int_{k}^{k+1}(x - k - \frac{1}{2})f'(x)\mathrm{d}x$$

记 $x - k$ 为 $\{x\} = x = [x]$ 表示 $x$ 的小数部分，取所有 $k \in [a, b)$，并相加，得到：

$$\int_{a}^{b}f(x)\mathrm{d}x = \sum\limits_{a \leq k \leq b}f(k) - \frac{1}{2}(f(a) + f(b)) - \int_{a}^{b}(\{x\} - \frac{1}{2})f'(x)\mathrm{d}x$$

整理之后即可得和式 $\sum\limits_{a \leq k \leq b}f(k)$ 与积分 $\int_{a}^{b}f(x)\mathrm{d}x$ 的一个精确关系：

$$\sum\limits_{a \leq k \leq b}f(k) = \int_{a}^{b}f(x)\mathrm{d}x + \frac{1}{2}(f(a) + f(b)) + \int_{a}^{b}(\{x\} - \frac{1}{2})f'(x)\mathrm{d}x$$

想要知道误差是多少，需要知道上式尾部积分的界。而前面的推导过程还可以继续迭代，使得误差项越来越小。沿着这个思路可以得到欧拉-麦克劳林公式的第一种形式。下面直接给出结论，证明过程参考潘承洞《阶的估计基础》或菲赫金哥尔茨《微积分学教程》。

定理(欧拉-麦克劳林公式，第一种形式)

设 $f(x)$ 为定义在区间 $[a, b]$ 上的函数，$a, b \in \mathbb{Z}$，假设当 $1 \leq i \leq 2m$ 时，导数 $f^{(i)}(x)$ 存在且连续，其中 $m$ 为常数，则：

$$\sum\limits_{a \leq k \leq b}f(k) = \int_{a}^{b}f(x)\mathrm{d}x + \frac{1}{2}(f(a) + f(b)) + \sum\limits_{1 \leq i \leq m}\frac{B_{2i}}{(2i)!}f^{(2i-1)}(x)\bigg|_{a}^{b} + R_{m}$$

其中 $B_{2i}$ 为伯努利数，$R_{m}$ 为余项，且满足：

$$|R_{m}| \leq \frac{|B_{2m}|}{(2m)!}\int_{a}^{b}|f^{(2m)}(x)|\mathrm{d}x < \frac{4}{(2\pi)^{2m}}\int_{a}^{b}|f^{(2m)}(x)|\mathrm{d}x$$

其中的伯努利数 $B_{n}$ 是指具有指数生成函数 $\sum\limits_{k}B_{k}\frac{z^{k}}{k!} = \frac{z}{e^{z}-1}$ 的序列。关于伯努利数和伯努利多项式，可以参考《具体数学》或菲赫金哥尔茨《微积分学教程》。

以上公式可以解决一些问题，比如自然数幂和问题 $\sum\limits_{i=1}^{n}i^{k}$ 用以上公式可以完整解决。

第二形式

当求和/积分的区间变大且步长固定时，以上定理就不能给出足够的准确的估计了，例如试图用积分 $\int_{a}^{b}f(x)\mathrm{d}x$ 来估计 $H_{N} = \sum\limits_{k=1}^{N}\frac{1}{k}$ 时，当 $N \rightarrow \infty$ 时，会发现和与积分的差将趋于一个未知常数。

在前面推导第一形式的欧拉-麦克劳林公式过程中，取 $a=1, b=N$，得到：

$$\int_{1}^{N}f(x)\mathrm{d}x = \sum\limits_{1 \leq k \leq N}f(k) - \frac{1}{2}(f(1) + f(N)) - \int_{1}^{N}(\{x\} - \frac{1}{2})f'(x)\mathrm{d}x$$

当 $N\rightarrow\infty$ 时，$f’(x)$ 趋于 0 的速度如果足够快，那么上式将把和式与积分关联到一个常数因子上，特别地，若以下积分存在：

$$\int_{1}^{\infty}(\{x\} - \frac{1}{2})f'(x)\mathrm{d}x$$

则可以定义 $f$ 的欧拉-麦克劳林常数：$C_{f} = \frac{1}{2}f(1) + \int_{1}^{\infty}({x} - \frac{1}{2})f’(x)\mathrm{d}x$，也就是：

$$C_{f} = \lim\limits_{N\rightarrow\infty}(\sum\limits_{1 \leq k \leq N}f(k) - \int_{1}^{N}f(x)\mathrm{d}x - \frac{1}{2}f(N))$$

于是有：

$$\sum\limits_{1 \leq k \leq N}f(k) = \int_{1}^{N}f(x)\mathrm{d}x + \frac{1}{2}f(N) + C_{f} - \int_{N}^{\infty}(\{x\} - \frac{1}{2})f'(x)\mathrm{d}x$$

想要知道误差是多少，需要知道上式尾部积分的界。而前面的推导过程还可以继续迭代，使得误差项越来越小。沿着这个思路可以得到欧拉-麦克劳林公式的第二种形式。下面直接给出结论，证明过程参考潘承洞《阶的估计基础》或菲赫金哥尔茨《微积分学教程》。

定理(欧拉-麦克劳林公式，第二种形式)

设 $f(x)$ 为定义在区间 $[1, \infty)$ 上的函数，假设当 $1 \leq i \leq 2m$ 时，导数 $f^{(i)}$ 存在且绝对可积，其中 $m$ 为常数，则：

$$\sum\limits_{1 \leq k \leq N}f(k) = \int_{1}^{N}f(x)\mathrm{d}x + \frac{1}{2}f(N) + C_{f} + \sum\limits_{1 \leq k \leq m}\frac{B_{2k}}{(2k)!}f^{(2k-1)}(N) + R_{m}$$

其中 $C_{f}$ 为与 $f$ 有关的常数，称为欧拉-麦克劳林常数，$R_{m}$ 为余项，且满足：

$$|R_{m}| = O(\int_{N}^{\infty}|f^{(2m)}(x)|\mathrm{d}x)$$

调和数的渐近估计

取 $f(x) = \frac{1}{x}$，其欧拉-麦克劳林常数为：

$$\begin{align} C_{f} &= \frac{1}{2}f(1) + \int_{1}^{\infty}(\{x\} - \frac{1}{2})f'(x)\mathrm{d}x \\ &= \frac{1}{2} - \int_{1}^{\infty}(\{x\} - \frac{1}{2})\frac{\mathrm{d}x}{x^{2}} \\ \end{align}$$

这里常数 $C_{f}$ 大约为 $0.57721$，正是我们熟知的欧拉常数，记为 $\gamma$。

在第二形式的欧拉-麦克劳林公式中，将 $f(x) = \frac{1}{x}$ 和 $C_{f} = \gamma$ 代入，得：

$$\begin{align} H_{N} &= \ln N + \gamma + \frac{1}{2N} - \sum\limits_{1\leq k < M}\frac{B_{2k}}{2kN^{2k}} + O(\frac{1}{N^{2m}}) \\ &= \ln N + \gamma + o(1) \end{align}$$

由此我们得到了比 $H_{N} = O(\ln N)$ 更精确的渐近估阶。

总结

本文讨论了调和级数的渐近估阶，这个例子非常经典，也是渐近估阶的入门必会题，在算法分析中非常常见，例如快速排序的平均时间复杂度，素数筛的时间复杂度等。

首先从快速排序平均情况时间复杂度的递归式出发，导出其生成函数满足的微分方程，求解之后得到生成函数的表达式。进一步结合生成函数的性质以及泰勒级数，将生成函数展开，即得到 $z^{n}$ 项的系数为一个含有调和数的公式。

在推导调和数的渐近估阶的过程中，我们应用了通过积分对级数进行估计的思想，结合所考察的函数的单调性，得到了比较粗略的渐近估阶结果。这也是我们在算法分析中熟知的结论。

进一步地，我们介绍了欧拉-麦克劳林公式的引入动机，以及简略的推导。基于欧拉-麦克劳林公式，我们得到了调和级数更加细致的估计，并且得到了欧拉常数。

这只是入门的例子，在算法分析中，还有很多更加精彩的例子，用到更多的数学工具，欧拉-麦克劳林公式依然是有力手段。

---