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

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的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}$ 项的系数为一个含有调和数的公式。

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

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

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

最后宣传一下：2024 届或 2025 届正在找实习或想走校招提前批的同学可以看一下拼多多的机会：拼多多2025届校园招聘研发岗实习生+提前批内推。

【算法分析】渐近估阶与欧拉-麦克劳林公式