[注：题图是 J. C. Schulze 1778 年刊行的书Neue und erweiterte Sammlung logarithmischer...中的一张对数表，由荷兰人 Wolfram 花费六年时间制成。表中包含从 1 到 10000 的所有整数的自然对数 (准确到小数点后 48 位)。Gauss 1791 年时从 Braunschweig 公爵Charles William Ferdinand, Duke of Brunswick-Wolfenbüttel处得到此书。]

虽说 Klein 把 Gauss 放在十九世纪数学家的行列当中，但是 Gauss 仍然具有十八世纪数学家的特质——也就是说，偏向用数值实验来归纳出各种定理。这种特质在十八世纪的数学家中普遍存在。具备这个特质的任何一个数学家都是计算狂魔，Euler 以及 Gauss 更是这些人中的佼佼者。根据高木贞治《近世数学史谈》第八章的内容，我们找到了 1852 年 10 月 11 日 Encke 写给 Gauss 的信件手稿。这一天 Eisenstein 去世，Encke 立即向 Gauss 通报了 Eisenstein 的死讯。信纸最后一页左下角的计算属于 Gauss 本人。他在 Encke 寄来的信件上计算Gotthold Eisenstein(1823.4.16 [这里《近世数学史谈》原文是错误的]—1852.10.11) 一共在世多少天。这一年 Gauss 75 岁。

[吐槽：Gauss 的计算瘾又犯了-_-||||]

从 1798 年开始，Gauss 的各种数学记录逐渐转移到所谓的Schedae(活页) 上。根据 Schlesinger 与 Klein 在 Gauss 全集第十卷，536-538 页的记载，Gauss 本人至少可能有两种数值计算 $P (u),Q (u)$Fourier 级数展开的方法：

• 用有限乘积近似无穷乘积，并将有限乘积展开为 $\sum x_n\sin^n\frac {\pi}{\varpi} u$ , 数值计算出各项系数，之后将 $\sin$ 的 $3,5,7,\cdots$ 次幂转化为 $\sin \frac {\pi}{\varpi} u,\sin \frac {3\pi}{\varpi} u,\sin \frac {5\pi}{\varpi} u,\sin \frac {7\pi}{\varpi} u,\cdots$ 的线性组合。(见 Gauss 全集第十卷第一册，537 页，这是Schedae上就有的记录) 这属于比较笨重的方法；

• 先计算出 $\log P (u)$ (或 $\log Q$) 的三角级数展开，然后借此可以通过数值方法算出 $P (u),Q (u)$ 的值。这是我们本篇的重点。

Gauss 在Scheda Aa(1798 年 7 月开始) 中 [Gauss 全集第三卷，413-432 页] 明确写下了 $\log P,\log Q$ Fourier 展开的解析表达式以及前若干项系数的数值计算结果。他的三角级数的展开来自这一节开头的公式 [Gauss 漏抄了下面公式等号左侧第二项，Gauss 全集的编辑Schering以及Schlesinger等人也没有把这一项补上]：

$$\log (1+\mu\cos\varphi)+\log (\frac {2}{1+\sqrt {1-\mu^2}})=2\sum_{n=1}^{\infty}\frac {(-1)^{n-1}}{n}\left (\frac {\mu}{1+\sqrt {1-\mu^2}}\right)^n\cos n\varphi$$

这个公式来自 $\log\left\vert 1+\frac {\mu}{1+\sqrt {1-\mu^2}} e^{i\varphi}\right\vert$ 的级数展开。那么我们有

$$\log (1-\mu^2\cos^2\varphi)+2\log (\frac {2}{1+\sqrt {1-\mu^2}})=-4\sum_{n=1}^{\infty}\frac {1}{2n}\left (\frac {\mu}{1+\sqrt {1-\mu^2}}\right)^{2n}\cos 2n\varphi$$

以及

$$\log (1+\mu^2\sin^2\varphi)+2\log (\frac {2}{1+\sqrt {1+\mu^2}})=-4\sum_{n=1}^{\infty}\frac {1}{2n}\left (\frac {\mu}{1+\sqrt {1+\mu^2}}\right)^{2n}\cos 2n\varphi$$

回忆

$$P (u)=\frac {\varpi}{\pi}\sin {\frac {\pi}{\varpi} u}\prod_{n=1}^{\infty}\left (1+\frac {\sin^2 {\frac {\pi}{\varpi} u}}{\sinh^2 {\pi n}}\right)$$

我们有

$$\begin {align}&\mu=\frac {1}{\sinh \pi n},\\&\frac {2}{1+\sqrt {1+\mu^2}}=(1-e^{-2\pi n}),\\&\frac {\mu}{1+\sqrt {1+\mu^2}}=e^{-2\pi n}.\end {align}$$

所以

$$\begin {align}&\log P (u)+2\sum_{n=1}^{\infty}\log (1-e^{-2\pi n})\\=&\log\frac {\varpi}{\pi}+\log\sin {\frac {\pi}{\varpi} u}-2\sum_{n=1}^{\infty}\frac {1}{n}\cos 2n\frac {\pi}{\varpi} u\sum_{m=1}^{\infty} e^{-2\pi mn}\\=&\log\frac {\varpi}{\pi}+\log\sin {\frac {\pi}{\varpi} u}-2\sum_{n=1}^{\infty}\frac {1}{n (e^{2n\pi}-1)}\cos 2n\frac {\pi}{\varpi} u.\end {align}$$

为了数值计算三角级数各项的系数，需要算出两个关键常数，也就是 $\frac {\varpi}{\pi},e^{-\pi}$ 。Gauss 是怎么进行这两个常数的数值计算的呢？

我们在这里已经提过，直接用幂级数展开来计算 $\varpi=2\int_{0}^{1}\frac {\mathrm {d} x}{\sqrt {1-x^4}}$ 是不现实的。为了数值计算这个常数，有必要找到该常数收敛更快的表达式。我们知道 $\begin {align}\varpi&=2\int_{0}^{\pi/2}\frac {\mathrm {d}\varphi}{\sqrt {1+\sin^2\varphi}}\\&=\sqrt {2}\int_{0}^{\pi/2}\frac {\mathrm {d}\varphi}{\sqrt {1-\frac {1}{2}\sin^2\varphi}}\\&=2\sqrt {\frac {2}{3}}\int_{0}^{\pi/2}\frac {\mathrm {d}\varphi}{\sqrt {1-\frac {1}{3}\cos 2\varphi}}\end {align}$

将被积函数展成幂级数，并且逐项积分，我们就可以得到 Gauss 在 Leiste 笔记中 (全集第十卷第一册，169 页，下图选用的是 Gauss 1799 年 11 月的记录 [同上，184 页]，这一纪录更加清晰完整) 的表达式：

最后一个表达式收敛最快，Gauss 利用这个表达式，算出 (全集第十卷第一册，169 页)

$$\frac {\varpi}{\pi}\approx0.834626841674030.$$

Gauss 手稿中 (全集第三卷，426-432 页) 有对 $e^{-\pi},e^{-\pi/2},e^{-\pi/4}$ 的数值计算 (精确到小数点后 50 位)。我们依照高木贞治《近世数学史谈》第八章和 Gauss 全集第三卷的内容来解说 Gauss 的数值计算。

Gauss 首先给出了 $e^{-\pi}$ 精确到小数点后 14 位的数值 $e^{-\pi}\approx0.04321391826377$ 。Gauss 没有写他本人是怎么得到这个初始值的。不过我们可以猜一猜 Gauss 的方法 [注意，这是作者本人的猜测，Gauss 关于这个初值没有写任何计算动机和方法]。从题图的对数表可以看到 $\log 23\approx 3.135$ 。这个值还不够接近 $\pi$ 。我们继续观察 $\log 529$ 附近的值。从 Wolfram 的表中可以得到 $\log 535\approx 6.2822667,\log 536\approx 6.2841342$ 。平均值与 $2\pi$ 的误差大约在 $2\times 10^{-5}$ 的量级。所以我们可以用 $\sqrt {\frac {2}{1071}}$ 来近似 $e^{-\pi}$ 。通过查 Wolfram 的表我们可以得到 $(\log1071-\log2)/2-\pi\approx0.00000779135411$ 。用幂级数展开三项就可以算到小数点后 15 位精度。 $\sqrt {\frac {2}{1071}}$ 可以通过连分数展开 / Newton 方法计算其任意精度近似值。这样我们就得到了 $e^{-\pi}$ 的初始近似值。

上面的方法可以说是 Gauss 第三卷中计算 $e^{-\pi}$ 到 50 位精度算法的雏形。Gauss 本人接下来的计算就是找到自然数的对数的线性组合$\log N$ ，使之与 $\pi$ 的差 $\delta$ 的绝对值足够小。自然数如果都落在 [1,10000] 的区间当中，那么 Gauss 就可以用 Wolfram 的结果计算 $\delta$ 的值。显然 $e^{-\pi}=\frac {1}{N} e^{\log N-\pi}=\frac {e^{\delta}}{N}$ 。 $e^{\delta}$ 的值可以由幂级数展开得到。问题就在于如何找到合适的线性组合。Gauss 首先计算了 $11\times10^{10}\times0.04321391826377$ ，它的整数部分应当是 $4753531009$ 。 $4753531008=128\times243\times152827$ 。如何判定 152827 是否可以分解质因数呢？Gauss 计算了 $\sqrt {152827\times2}\approx 552,\sqrt {152827\times3}\approx 677.$ 而 $\begin {align} 552^2&\equiv950\mod 152827\\677^2&\equiv152\mod 152827\end {align}$

所以 $x^2\equiv3800 \mod 152827$ 有两个不同的解 1104 和 3385。那么 152827 一定不是质数 (为什么？)。如果 p 是 152827 的质因数，那么 p 一定整除 $3385\pm1104$ 中的其中一个 (为什么？)。据此我们立刻可以得到 $152827=67\times 2281$ 。那么我们上面提到的 $\delta$ 的值为 $\delta=\log11+10\log 10-\log 2281-\log 2144-\log 972-\pi$ 。根据 Wolfram 的表格 [$\pi$ 的 100 位小数早在 1706 年由John Machin得到]，我们可以得到 $\delta=2.135\cdots\times10^{-10}$ 。这一 $\delta$ 的值仍然偏大。Gauss 作了一些修正，令 $\delta^{\prime}=\delta-\log (1+2135\times10^{-13}),\delta^{\prime\prime}=\delta^{\prime}-\log (1+1443\times10^{-17})$ 。修正项每项都可以用幂级数计算。此时 $|\delta^{\prime\prime}|< 10^{-17}$ 。计算 $e^{\delta^{\prime\prime}}$ 到小数点后 50 位只用展开幂级数前三项即可。遵循以上的过程，Gauss 得到了

$$e^{-\pi}=0.04321391826377224977441773717172801127572810981063$$

Gauss 又对 $59\times10^{9}\times0.04321391826377$ 作了同样的处理，结果完全一致。Gauss 借此也间接验证了 Wolfram 表中相关的项直到小数点最末一位都是准确的。