400多年前的常用对数表是如何制作的

hbghlyj

hbghlyj 发表于 2025-1-12 10:39
另一种用于极高精度计算的替代方法是公式
$ {\displaystyle \ln x\approx {\frac {\pi }{2M(1,4/s)}}-m\ln 2,} $
其中 M 表示 1 和 4/s 的算术-几何平均数，并且
$ {\displaystyle s=x2^{m}>2^{p/2},} $
选择 m 以达到 p 位精度。

为了证明该公式，我们首先总结椭圆函数理论中的必要公式。模数 $k$ 的第一类完全椭圆积分 $K$ 和 $K'$ 定义为
\begin{align*}
& K \equiv K(k)=\int_0^{\pi / 2} \frac{d \theta}{\sqrt{1-k^2 \sin ^2 \theta}} \\
& K' \equiv K\left(k'\right), \quad k'=\sqrt{1-k^2}
\end{align*}
我们使用椭圆 theta 函数 $\theta_0(u)$ 和 $\theta_3(u)$ 的著名 $q$ 展开公式。我们只需要 $u=0$ 时的值：
\[
\begin{aligned}
& \theta_0^0 \equiv \theta_0(0)=1+2 \sum_{n=1}^{\infty}(-1)^n q^{n^2}, \\
& \theta_3^0 \equiv \theta_3(0)=1+2 \sum_{n=1}^{\infty} q^{n^2} .
\end{aligned}
\]
如果 $|q| \ll 1$，这些级数收敛得非常快。完全椭圆积分 $K$ 和 $K'$ 与 $q$ 的关系公式为
\begin{aligned}
q=\exp \left(-\pi K' / K\right) \\ \pi K' / K=\log (1 / q)
\end{aligned}此外，$k$（或 $k'$）和 $K$ 可以通过 $\theta_0^0$ 和 $\theta_3^0$ 表示为
\[
\begin{aligned}
& k=\sqrt{1-\left(\theta_0^0 / \theta_3^0\right)^4} \quad\left(\text { 或 } k'=\left(\theta_0^0 / \theta_3^0\right)^2\right), \\
& K=\pi\left(\theta_3^0\right)^2 / 2
\end{aligned}
\]
我们需要的另一个公式是著名的高斯算术-几何平均。设 $a_0=1$ 和 $b_0=k, 0<k<1$，并通过公式计算第 $n$ 项 $a_n$ 和 $b_n$
\[
a_n=\left(a_{n-1}+b_{n-1}\right) / 2, \quad b_n=\sqrt{a_{n-1} b_{n-1}}, \quad n \geqq 1 .
\]
然后，序列 $\left\{a_n\right\}$ 和 $\left\{b_n\right\}$ 以二次速度收敛到一个共同的极限。该极限称为算术-几何平均
\[
\lim _{n \rightarrow \infty} a_n=\lim _{n \rightarrow \infty} b_n=\pi / 2 K'
\]

hbghlyj

hbghlyj 发表于 2025-1-12 10:39
另一种用于极高精度计算的替代方法是公式
$ {\displaystyle \ln x\approx {\frac {\pi }{2M(1,4/s)}}-m\ln 2,} $
其中 M 表示 1 和 4/s 的算术-几何平均数，并且
$ {\displaystyle s=x2^{m}>2^{p/2},} $
选择 m 以达到 p 位精度。

首先，把$m\ln2$移到等式左边，$ {\displaystyle \ln x+m\ln2\approx {\frac {\pi }{2M(1,4/s)}},} $
使用$s=x2^{m}$得$\ln x+m\ln2=\ln s$，化为$ {\displaystyle \ln s\approx {\frac {\pi }{2M(1,4/s)}},} $
使用$ {\frac {\pi }{2M(1,4/s)}}=K'(4/s)$，化为$K'(4/s)\approx\ln s$
见《Pi and the AGM》第7.2节theta function algorithms for Log

hbghlyj

欲计算$ a_{0}=1$和$ g_{0}=10000$的算术-几何平均数，首先算出它们的算术平均数和几何平均数：

$ a_{1}={\frac {1+10000}{2}}=5000.5, $
$ g_{1}={\sqrt {1\times 10000}}=\, $$ 100 $

然后进行迭代：

$ a_{2}={\frac {a_1+g_1}{2}}=2550.25, $
$ g_{2}={\sqrt {a_1g_1}}=\, $$ 707.142136 $

继续计算，可得出以下的值：

$n$	$a_n$	$g_n$
0	1	10000
1	5000.50000000000	100.000000000000
2	2550.25000000000	707.142135641768
3	1628.69606782088	1342.90328446259
4	1485.79967614174	1478.91220120329
5	1482.35593867252	1482.35193850513
6	1482.35393858882	1482.35393858747
7	1482.35393858815	1482.35393858815

1和10000的算术-几何平均数是两个数列的公共极限，约为1482.35393858815
因此，1和1/10000的算术-几何平均数约为0.148235393858815
将$s=40000$代入${\displaystyle \ln s\approx {\frac {\pi }{2M(1,4/s)}},}$得 $$\ln 40000\approx {\frac {\pi }{2M(1,1/10000)}}=\frac{\pi}{2\times 0.148235393858815}=10.5966347570877 $$ 而$\ln 40000$的精确值为$10.5966347330961\dots$

hbghlyj

hbghlyj 发表于 2025-1-12 15:22
化为$K'(4/s)\approx\ln s$

这个的证明见《Pi and the AGM》第7.1节“$\pi$ and Log”定理 7.2 式(7.1.1)\[\tag{7.1.1}
\left|K'(k)-\log \left(\frac{4}{k}\right)\right|=O\left(\left|k^2 \log k\right|\right) \quad \operatorname{re}(k)>0
\]书上留作习题了，没有证明
@青青子衿 请问这个(7.1.1)如何证明呢