QWen3翻译的论文COMPUTING LOGARITHM INTERVALS WITH THE ARITHMETIC‐GEOMETRIC‐MEAN ITERATION - l1t1/note GitHub Wiki
使用算术-几何平均迭代计算对数区间 DANIEL J. BERNSTEIN
中文概述
这篇论文提出了一种快速算法，用于在给定一个正实数 x 的高精度区间的情况下，计算自然对数 log x 的高精度区间。该算法利用了完全椭圆积分和算术-几何平均（AGM）迭代的性质。
主要特点
- 高效性：对于
p位精度的需求，算法仅需约2 log p次平方根运算和5 log p次乘法。 - 额外收益：在计算对数的同时，还可以高效地计算圆周率
π。 - 区间计算：在输入为包含
x的紧致区间时，可以输出包含log x和π的紧致区间。
技术背景
- 椭圆积分：文中介绍了两类完全椭圆积分
I(a, b)和I1(a, b)的基本性质，以及它们与对数的关系。 - AGM 迭代：通过不断计算两个数的算术平均和几何平均，最终收敛到一个共同值，称为算术-几何平均（AGM）。
- 对数边界估计：利用椭圆积分的性质，推导出对数函数的上下界，从而实现对
log x的精确计算。
应用场景
- 适用于需要高精度计算对数的应用场景，如数值分析、科学计算和密码学。
- 支持多个对数的批量计算，效率更高。
结论
该论文详细探讨了如何利用 AGM 迭代和椭圆积分理论来高效计算自然对数，并提供了严格的数学证明和误差分析。它还比较了多种已有的对数计算方法，并指出当前方法在性能上的优势。
1. 引言
本文提出了一种快速算法，给定一个正实数 $ x $ 的高精度区间，可以计算出自然对数 $ \log x $ 的高精度区间。该算法具有以下几个有用的特性：
- 它非常快。对于 $ p $ 位的精度需求，在 $ x $ 处于（例如）2 到 $ 2^p $ 的范围内时，算法的时间主要由大约 7 $\log p$ 次 $ p $ 位数的操作决定，具体包括 2 $\log p$ 次平方根运算和 5 $\log p$ 次乘法。这里，$\log$ 表示自然对数，而 $\log_2$ 是以 2 为底的对数。
- 在计算 $ \log x $ 的过程中，它还能几乎免费地计算出高精度的圆周率 $ \pi $。
- 对后续对数的计算速度更快：从渐近角度看，每个对数的计算仅需约 2 $\log p$ 次平方根运算和约 2 $\log p$ 次乘法。
- 给定一个包含 $ x $ 的紧致区间，它可以计算包含 $ \log x $ 和 $ \pi $ 的紧致区间。参见 [4] 中的应用。对于大的 $ p $，处理区间的额外成本是微不足道的。
有关算法的详细信息，请参见第 5 节。
该算法依赖于完全椭圆积分和算术-几何平均 (AGM) 迭代的各种性质。本文提供了这些性质的简化证明。第 2 节介绍了文中使用的椭圆积分 $ I $ 和 $ I_1 $；第 3 节将 $ I $ 和 $ I_1 $ 与对数联系起来；第 4 节介绍了 AGM 迭代。
以往的工作
三十多年前，Salamin 就指出椭圆积分 $ I $ 可用于快速计算 $ \log x $ 和 $ \pi $ 的高精度近似值；参见 [2, Item 143]。Salamin 的算法主要由 $ \Theta(\log p) $ 次高精度操作组成，但其常数因子比本文中的算法更大。
随后对 Salamin 算法的改进几乎包含了所有最小化常数因子所需的思想，但这些思想从未被正确结合。有关文献的综述请参见第 6 节。
Brent 在 [9, 第 6 节] 提出了另一种计算高精度指数和对数的方法，不依赖于椭圆积分。有关进一步讨论，请参见 [3, 第 12–16 节]。对于较小的 $ p $，Brent 的方法可能比本文中的算法更快，但对于较大的 $ p $，它的速度会慢一个大致与 $ \log p $ 成正比的因子。
一些对数算法有时会提供适合区间计算的显式边界。参见 [14] 和 [13]。
从语言的角度以及显式不等式与数量级估计之间的差异来看，第 2、3 和 4 节的结果已有两个世纪的历史，并且只是 [8] 中收集的丰富材料的一小部分。然而，我尚未见过如此简短的关于 $ I_1 $ 性质的证明，特别是定理 2.6 和定理 4.3 的证明。
2. 椭圆积分：基本性质
对于正实数 $ a, b $，定义： $$ I(a, b) = \int_0^\infty (x^2 + a^2)^{-1/2}(x^2 + b^2)^{-1/2} dx $$
也定义： $$ I_1(a, b) = \frac{\partial I(a, b)}{\partial a} = \int_0^\infty -a(x^2 + a^2)^{-3/2}(x^2 + b^2)^{-1/2} dx $$
可积性、在积分号下求导等操作的合法性来源于本节中所有被积函数在 $ x $ 上是连续的，在 $ a, b $ 上是可微的，并且符号保持不变。
定理 2.1
如果 $ b \leq a $，则有： $$ \frac{\pi}{2a} \leq I(a, b) \leq \frac{\pi}{2b} $$
证明
由于 $ x^2 + b^2 \leq x^2 + a^2 $，所以有：
$$
\int_0^\infty (x^2 + a^2)^{-1} dx \leq I(a, b) \leq \int_0^\infty (x^2 + b^2)^{-1} dx
$$
而这两个积分的结果分别是 $ \frac{\pi}{2a} $ 和 $ \frac{\pi}{2b} $。证毕。
定理 2.2
$$ 0 \leq -a I_1(a, b) \leq I(a, b) $$
证明
因为 $ a^2 \geq 0 $，所以 $ x^2 + a^2 \geq x^2 $，从而有：
$$
-a(x^2 + a^2)^{-3/2} \leq 0
$$
又因为积分结果为负数，乘以 $ -a $ 后变为非负数，因此 $ -a I_1(a, b) \geq 0 $。
同时，原积分形式小于等于 $ I(a, b) $ 的形式，因此上界成立。证毕。
定理 2.3
$$ I(a, b) = \frac{1}{a} I\left(1, \frac{b}{a}\right) $$
证明
令 $ x = a u $，代入得：
$$
a I(a, b) = \int_0^\infty a^2 (a^2 u^2 + a^2)^{-1/2} (a^2 u^2 + b^2)^{-1/2} du
= \int_0^\infty (u^2 + 1)^{-1/2} \left(u^2 + \frac{b^2}{a^2}\right)^{-1/2} du = I\left(1, \frac{b}{a}\right)
$$
两边同除以 $ a $，即得结论。证毕。
定理 2.4
$$ I(a, b) + a \frac{\partial I}{\partial a}(a, b) + b \frac{\partial I}{\partial b}(a, b) = 0 $$
证明
对定理 2.3 中的恒等式两边分别对 $ a $ 和 $ b $ 求导，然后应用算子 $ a \frac{\partial}{\partial a} + b \frac{\partial}{\partial b} $，即可得出结论。证毕。
定理 2.5
$$ I(a, b) = I\left( \frac{a + b}{2}, \sqrt{ab} \right) $$
证明
令变换 $ u = \frac{x - ab/x}{2} $，利用恒等式：
$$
(2u)^2 + (a + b)^2 = \frac{(x^2 + a^2)(x^2 + b^2)}{x^2}
$$
以及：
$$
x du = (u^2 + ab)^{1/2} dx
$$
代入积分表达式得： $$ I(a, b) = \int_{-\infty}^\infty ((2u)^2 + (a + b)^2)^{-1/2}(u^2 + ab)^{-1/2} du = 2 \int_0^\infty (2^2(u^2 + (\frac{a + b}{2})^2))^{-1/2}(u^2 + ab)^{-1/2} du = I\left( \frac{a + b}{2}, \sqrt{ab} \right) $$
证毕。
定理 2.6
$$ I(a, b) + 2a I_1(a, b) = I_1\left( \frac{a + b}{2}, \sqrt{ab} \right) \cdot \frac{a - b}{2} $$
证明
对定理 2.5 应用算子 $ a \frac{\partial}{\partial a} - b \frac{\partial}{\partial b} $，得到：
$$
a I_1(a, b) - b \frac{\partial I(a, b)}{\partial b} = I\left( \frac{a + b}{2}, \sqrt{ab} \right) \cdot \frac{a - b}{2}
$$
再结合定理 2.4： $$ I(a, b) + a I_1(a, b) + b \frac{\partial I(a, b)}{\partial b} = 0 $$
联立两式即可消去 $ \frac{\partial I}{\partial b} $，最终得到结论。证毕。
定理 2.7
$$ I(1, b) = 2 \int_0^{\sqrt{b}} (x^2 + 1)^{-1/2}(x^2 + b^2)^{-1/2} dx $$
证明
将积分区间从 $ \sqrt{b} $ 到无穷大进行变量替换 $ x = \frac{b}{u} $，并重新整理后合并到 $ [0, \sqrt{b}] $ 区间，即可得证。证毕。
定理 2.8
$$ (1 + b)^{-1} + I(1, b) + I_1(1, b) = 2 \int_0^{\sqrt{b}} b^2 (x^2 + 1)^{-1/2}(x^2 + b^2)^{-3/2} dx $$
证明
由定理 2.4 可得：
$$
I(1, b) + I_1(1, b) + b \frac{\partial I(1, b)}{\partial b} = 0
$$
由定理 2.7 推导出 $ \frac{\partial I(1, b)}{\partial b} $ 的表达式，代入后即可化简得到结论。证毕。
定理 2.9
$$ a^2 I(a, b) + (a^2 - b^2)a I_1(a, b) = \int_0^\infty a^2 (x^2 + a^2)^{-3/2}(x^2 + b^2)^{1/2} dx $$
证明
观察到：
$$
x^2 + a^2 - (a^2 - b^2) = x^2 + b^2
$$
将其代入积分表达式，即可得到左边的表达式。这个定理解释了为什么这些积分被称为“椭圆积分”，以及为什么相关的代数结构被称为“椭圆曲线”——通过令 $ x = a \tan \theta $，可以发现该积分表示的是椭圆 $ \theta \mapsto (a \cos \theta, b \sin \theta) $ 的一象限弧长。此定理在本文其余部分不再使用。
3. 椭圆积分：对数边界
本节介绍如下形式的精确边界：“如果 $ b \approx 0 $，则 $ I(1, b) \approx \log(4/b) $，且 $ I(1, b) + I_1(1, b) \approx 1 $。”
定义：
$$
L(b) = \log\left(\sqrt{b^{-1}} + \sqrt{b^{-1} + 1}\right)
$$
这里使用的简单证明方法并非新创，这种技术在对数边界方面的明确估计也并不新颖，但我并未见过以前关于显式边界的简明证明。对于使用相同证明技术但未给出显式边界的文献，可参考 [18, 第522页] 和 [15]；对于带有更复杂证明的显式边界，可参考 [7, 第355–356页] 和 [8, 定理7.2]。
定理 3.1
如果 $ 0 < b \leq 1 $，则有以下不等式成立：
$$ \left(2 + \frac{1}{2} b^2\right)L(b) - \left(\frac{1}{2} b\right)(1 + b)^{1/2} \leq I(1, b) \leq \left(2 + \frac{1}{2} b^2 + \frac{9}{32} b^4\right)L(b) - \left(\frac{1}{2} b - \frac{3}{16} b^2 + \frac{9}{32} b^3\right)(1 + b)^{1/2} $$
上下界之间的差异大致在 $ b^2 $ 的数量级。通过相同的技术可以进一步收紧这些边界。
证明
根据定理 2.7，有：
$$
I(1, b) = 2 \int_0^{\sqrt{b}} (1 + x^2)^{-1/2}(x^2 + b^2)^{-1/2} dx
$$
注意到当 $ 0 \leq x \leq 1 $ 时，有： $$ (1 + x^2)^{-1/2} \geq 1 - \frac{x^2}{2} $$ 因为： $$ 1 \geq 1 - \frac{(3 - x^2)x^4}{4} = (1 - \frac{x^2}{2})^2(1 + x^2) $$
因此： $$ I(1, b) \geq \int_0^{\sqrt{b}} (2 - x^2)(x^2 + b^2)^{-1/2} dx = \left[(2 + \frac{1}{2} b^2)\log(x + (x^2 + b^2)^{1/2}) - \frac{1}{2} x (x^2 + b^2)^{1/2} \right]_0^{\sqrt{b}} = \left(2 + \frac{1}{2} b^2\right)L(b) - \left(\frac{1}{2} b\right)(1 + b)^{1/2} $$
同样地，由于： $$ (1 + x^2)^{-1/2} \leq 1 - \frac{1}{2} x^2 + \frac{3}{8} x^4 \quad \text{（对所有 } x \geq 0 \text{ 成立）} $$
因此： $$ I(1, b) \leq \int_0^{\sqrt{b}} (2 - x^2 + \frac{3}{4} x^4)(x^2 + b^2)^{-1/2} dx = \left[(2 + \frac{1}{2} b^2 + \frac{9}{32} b^4)\log(x + (x^2 + b^2)^{1/2}) - \left(\frac{1}{2} - \frac{3}{16} x^2 + \frac{9}{32} b^2\right)x(x^2 + b^2)^{1/2} \right]_0^{\sqrt{b}} = \left(2 + \frac{1}{2} b^2 + \frac{9}{32} b^4\right)L(b) - \left(\frac{1}{2} b - \frac{3}{16} b^2 + \frac{9}{32} b^3\right)(1 + b)^{1/2} $$
证毕。
定理 3.2
如果 $ 0 < b \leq 1 $，则有以下不等式成立：
$$ (2 + b^2)(1 + b)^{-1/2} - b^2 L(b) \leq (1 + b)^{-1} + I(1, b) + I_1(1, b) \leq (2 + b^2 + \frac{3}{8} b^3 + \frac{9}{8} b^4)(1 + b)^{-1/2} - (b^2 + \frac{9}{8} b^4)L(b) $$
证明
由定理 2.8 得：
$$
(1 + b)^{-1} + I(1, b) + I_1(1, b) = 2 \int_0^{\sqrt{b}} b^2 (x^2 + 1)^{-1/2}(x^2 + b^2)^{-3/2} dx
$$
再次利用不等式： $$ (1 + x^2)^{-1/2} \geq 1 - \frac{x^2}{2} \quad \text{（如定理 3.1 所述）} $$
因此： $$ (1 + b)^{-1} + I(1, b) + I_1(1, b) \geq \int_0^{\sqrt{b}} b^2 (2 - x^2)(x^2 + b^2)^{-3/2} dx = \left[(2 + b^2)x(x^2 + b^2)^{-1/2} - b^2 \log(x + (x^2 + b^2)^{1/2}) \right]_0^{\sqrt{b}} = (2 + b^2)(1 + b)^{-1/2} - b^2 L(b) $$
同样地，使用： $$ (1 + x^2)^{-1/2} \leq 1 - \frac{1}{2} x^2 + \frac{3}{8} x^4 \quad \text{（对所有 } x \geq 0 \text{ 成立）} $$
得到： $$ (1 + b)^{-1} + I(1, b) + I_1(1, b) \leq \int_0^{\sqrt{b}} b^2 (2 - x^2 + \frac{3}{4} x^4)(x^2 + b^2)^{-3/2} dx = \left[(2 + b^2 + \frac{3}{8} b^2 x^2 + \frac{9}{8} b^4)x(x^2 + b^2)^{-1/2} - (b^2 + \frac{9}{8} b^4)\log(x + (x^2 + b^2)^{1/2}) \right]_0^{\sqrt{b}} = (2 + b^2 + \frac{3}{8} b^3 + \frac{9}{8} b^4)(1 + b)^{-1/2} - (b^2 + \frac{9}{8} b^4)L(b) $$
证毕。
总结
本节建立了椭圆积分 $ I(1, b) $ 和 $ I_1(1, b) $ 与自然对数函数 $ \log x $ 的关系，并给出了在 $ b $ 较小时的精确上下界估计。这些结果为后续基于 AGM 迭代计算高精度对数区间提供了理论基础。
4. 椭圆积分：AGM 迭代
在本节中，$ a_0 $ 和 $ b_0 $ 是正实数；通过递归定义 $ a_{n+1} = \frac{a_n + b_n}{2} $ 和 $ b_{n+1} = \sqrt{a_n b_n} $ 来得到 $ a_n $ 和 $ b_n $。随着 $ n $ 的增加，$ a_n $ 和 $ b_n $ 都快速收敛到 $ \frac{\pi}{2I(a_0, b_0)} $，即 $ a_0 $ 和 $ b_0 $ 的算术-几何平均（AGM）。见定理 4.2 和 4.5。
定理 4.1
如果 $ n \geq 1 $，则 $ a_n \geq b_n $。
证明
几何平均不会超过算术平均：如果 $ n \geq 0 $，则
$$
a_{n+1}^2 - b_{n+1}^2 = \frac{(a_n - b_n)^2}{4} \geq 0，
$$
因此 $ a_{n+1} \geq b_{n+1} $。
定理 4.2
对于 $ n \geq 1 $，有
$$
\frac{\pi}{2a_n} \leq I(a_0, b_0) \leq \frac{\pi}{2b_n}.
$$
证明
根据定理 2.5，
$$
I(a_0, b_0) = I(a_1, b_1) = I(a_2, b_2) = \cdots = I(a_n, b_n).
$$
由定理 4.1，$ a_n \geq b_n $。应用定理 2.1 即可得证。
定理 4.3
定义
$$
\varepsilon_n = 2^n (a_n^2 - b_n^2) \left( -a_n \frac{I_1}{I}(a_n, b_n) \right).
$$
则对于 $ n \geq 1 $，有 $ 0 \leq \varepsilon_n \leq 2^n (a_n^2 - b_n^2) $，并且
$$
(a_0^2 - b_0^2) \left( -a_0 \frac{I_1}{I}(a_0, b_0) \right) = \varepsilon_n + \sum_{0 \leq i < n} 2^{i-1} (a_i^2 - b_i^2).
$$
证明
对于 $ n \geq 1 $，有 $ 0 < b_n \leq a_n $，因此根据定理 2.2，
$$
0 \leq -a_n \frac{I_1}{I}(a_n, b_n) \leq 1,
$$
即
$$
0 \leq \varepsilon_n \leq 2^n (a_n^2 - b_n^2).
$$
代入公式
$$
(a_{n+1}^2 - b_{n+1}^2)a_{n+1} = \frac{(a_n^2 - b_n^2)(a_n - b_n)}{8},
$$
可以得到
$$
\varepsilon_{n+1} = 2^n (a_n^2 - b_n^2) \left( -\frac{(a_n - b_n)}{4} \frac{I_1}{I}(a_{n+1}, b_{n+1}) \right).
$$
然后应用定理 2.6 和 2.5：
$$
\varepsilon_{n+1} = 2^n (a_n^2 - b_n^2) \left( -\frac{1}{2} - a_n \frac{I_1}{I}(a_n, b_n) \right) = -2^{n-1} (a_n^2 - b_n^2) + \varepsilon_n.
$$
通过归纳法，可以得出
$$
\varepsilon_0 = \varepsilon_n + \sum_{0 \leq i < n} 2^{i-1} (a_i^2 - b_i^2).
$$
定理 4.4
如果 $ 1 \leq \frac{a_0}{b_0} \leq 1 + 2^{2m} $，那么对于 $ n \geq 0 $，有
$$
1 \leq \frac{a_n}{b_n} \leq 1 + 2^{2m - n}.
$$
证明
如果 $ 1 \leq \frac{a_n}{b_n} \leq 1 + 2^{2m - n} $，那么
$$
\frac{a_{n+1}}{b_{n+1}} = \frac{\sqrt{\frac{a_n}{b_n}} + \sqrt{\frac{b_n}{a_n}}}{2} \leq \sqrt{\frac{a_n}{b_n}} \leq \sqrt{1 + 2^{2m - n}} \leq 1 + 2^{2m - n - 1}.
$$
定理 4.5
如果 $ m \geq 0 $ 且 $ 1 \leq \frac{a_0}{b_0} \leq 1 + 2^{2m} $，那么对于 $ n \geq m $，有
$$
1 \leq \frac{a_n}{b_n} \leq 1 + 2^{3 - 2n + 1 - m}.
$$
证明
当 $ n = m $ 时，由定理 4.4 可知
$$
2^{2m - m} = 2^{3 - 2m + 1 - m}.
$$
假设 $ \frac{a_n}{b_n} = 1 + \varepsilon $，其中 $ 0 \leq \varepsilon \leq 2^{3 - 2n + 1 - m} $，则
$$
4 + 4\varepsilon + \varepsilon^2 \leq (4 + \varepsilon^2 + \frac{\varepsilon^4}{16})(1 + \varepsilon),
$$
从而
$$
(1 + \varepsilon) + 2 + \frac{1}{1 + \varepsilon} \leq (2 + \frac{\varepsilon^2}{4})^2,
$$
因此
$$
\sqrt{\frac{a_n}{b_n}} + \sqrt{\frac{b_n}{a_n}} \leq 2 + \frac{\varepsilon^2}{4},
$$
进一步得到
$$
\frac{a_{n+1}}{b_{n+1}} \leq 1 + \frac{\varepsilon^2}{8} \leq 1 + 2^{3 - 2n + 2 - m}.
$$
5. 计算对数区间
本节介绍几种算法，给定一个包含正实数 $ x $ 的区间，计算出一个包含 $ \log x $ 的区间。关于快速计算区间运算（如加法、减法、乘法、除法和平方根）的子程序，请参见 [6]。
这些算法使用一个参数 $ p $ 来决定何时停止迭代。对于“任意 $ x $”的算法，如果输入区间的精度为 $ p $ 位，则输出区间的精度也大约为 $ p $ 位。而对于“超大 $ x $”的算法，即使 $ p $ 非常大，输出精度最多也只能达到约 $ 4 \lg x $ 位。此外，当 $ \lg \lg x $ 远大于 $ \lg p $ 时，“超大 $ x $”算法也会变慢；如果 $ \lg \lg x $ 明显大于 $ \lg p $，则应改用“任意 $ x $”算法。
需要注意的是，正如文献 [5] 所讨论的那样，常规的算术运算（如平方根）和一系列算术操作中通常包含许多冗余步骤，这些冗余可以被消除。例如，在 [11, Theorem 9.1] 中报告的 AGM 步骤序列可以被优化得比原报告快近三倍。同样地，Borwein 和 Borwein 在 [8, page 222] 中观察到，通过一种“四次 AGM”方法，即直接从 $ \sqrt{a_n} $ 和 $ \sqrt{b_n} $ 计算 $ \sqrt{a_{n+2}} $ 和 $ \sqrt{b_{n+2}} $，可以获得速度提升；但笔者认为，当平方根算法得到适当优化后，这些所谓的加速可能会反而变慢，因此实验留给读者进行。
计算超大 $ x $ 的对数
设 $ x > 4 $ 是一个实数。定义：
- $ a_0 = 1 $
- $ b_0 = b = \left( \frac{2x}{x^2 - 1} \right)^2 $
注意到： $$ \sqrt{b^{-1}} + \sqrt{b^{-1} + 1} = x $$ 所以根据第 3 节定义的函数 $ L(b) = \log x $。同时注意 $ 0 < b < \frac{1}{2} $。
按照第 4 节的方法，递归定义： $$ a_{n+1} = \frac{a_n + b_n}{2}, \quad b_{n+1} = \sqrt{a_n b_n} $$
定义： $$ c = 1 + \frac{I_1}{I}(1, b) $$ 根据定理 2.2，有 $ 0 \leq c \leq 1 $。
步骤如下：
- 计算 $ b $ 的区间。
- 依次计算 $ a_1, b_1, a_2, b_2, \dots, a_n, b_n $，直到 $ a_n - b_n $ 不再明显大于 $ \frac{1}{2^p} $。根据定理 4.5，步数 $ n $ 最多约为 $ \lg \lg x + \lg p $。
- 计算区间 $ [0, 2^n(a_n^2 - b_n^2)] + \sum_{0 \leq i < n} 2^{i-1}(a_i^2 - b_i^2) $。根据定理 4.3，该区间包含： $$ (1 - b^2)(1 - c) $$
- 将上述区间除以 $ 1 - b^2 $ 的区间，并从 1 中减去结果，得到 $ c $ 的区间。
- 最后，计算以下区间： $$ \left[ \frac{( \frac{1}{2} b - \frac{3}{16} b^2 + \frac{9}{32} b^3 ) c (1 + b)^{1/2} - (1 + b)^{-1} + (2 + b^2)(1 + b)^{-1/2} }{(2 + \frac{1}{2} b^2 + \frac{9}{32} b^4)c + b^2 }, \frac{( \frac{1}{2} b ) c (1 + b)^{1/2} - (1 + b)^{-1} + (2 + b^2 + \frac{3}{8} b^3 + \frac{9}{8} b^4)(1 + b)^{-1/2} }{(2 + \frac{1}{2} b^2)c + (b^2 + \frac{9}{8} b^4)} \right] $$
根据定理 3.1 和 3.2，该区间包含 $ L(b) = \log x $。注意像 $ \frac{9}{8} b^4 $ 这样的项对输出影响很小，可以用粗略估计代替。
此时还可以通过额外几步计算 $ \pi $ 的区间：
- 区间： $$ \left[ (2 + \frac{1}{2} b^2)L(b) - (\frac{1}{2} b)(1 + b)^{1/2},; (2 + \frac{1}{2} b^2 + \frac{9}{32} b^4)L(b) - (\frac{1}{2} b - \frac{3}{16} b^2 + \frac{9}{32} b^3)(1 + b)^{1/2} \right] $$ 包含 $ I(1, b) $（根据定理 3.1），而区间： $$ [2b_n I(1, b),; 2a_n I(1, b)] $$ 包含 $ \pi $（根据定理 4.2）。这种非常快速的 $ \pi $ 计算将在后续被利用。
数值稳定性分析
- 计算 $ c $：由于 $ c \approx 1 / \log(4/b) $，因此在计算 $ c = 1 - (1 - c) $ 时会损失大约 $ \lg \log(4/b) \approx \lg \lg x $ 位精度。
- 其他算术运算：每个操作仅损失有限数量的精度位，整体稳定。
- 定理 3.1 和 3.2 中的区间精度限制：受限于 $ \lg(1/b^2) \approx 4 \lg x $ 位精度。
计算多个超大 $ x $ 的对数
在计算一个超大 $ \log x $ 后，可以通过利用第一次计算中获得的 $ \pi $ 区间来更快地计算另一个超大 $ \log $ 值。
具体步骤如下：
- 从 $ x $ 出发，定义并计算 $ b, a_1, b_1, \dots, a_n, b_n $，跳过 $ \sum_i 2^{i-1}(a_i^2 - b_i^2) $ 的计算。
- 计算区间 $ [\pi/(2a_n), \pi/(2b_n)] $，根据定理 4.2，该区间包含 $ I(1, b) $。
- 计算区间： $$ \left[ \frac{I(1, b) + (\frac{1}{2} b - \frac{3}{16} b^2 + \frac{9}{32} b^3)(1 + b)^{1/2}}{2 + \frac{1}{2} b^2 + \frac{9}{32} b^4},; \frac{I(1, b) + (\frac{1}{2} b)(1 + b)^{1/2}}{2 + \frac{1}{2} b^2} \right] $$ 根据定理 3.1，该区间包含 $ L(b) = \log x $。
进一步优化：在某些应用中，只需计算对数之间的比值即可。例如，先用 $ I_1/I $ 计算 $ \log x $ 和 $ \pi $，然后用 $ I $ 计算 $ \pi/\log y $，最后相除得到 $ (\log y)/\log x $；但更高效的方式是两次使用 $ I $ 分别计算 $ \pi/\log x $ 和 $ \pi/\log y $，然后再相除得到 $ (\log y)/\log x $。
计算任意 $ x $ 的对数
给定一个包含正实数 $ x $ 的区间，找到一个整数 $ k $，使得 $ 2^k x $ 略大于 $ 2^{p/4} $，然后使用上述算法计算超大数 $ 2^k x $ 和 $ 2^{\lceil p/4 \rceil} $ 的对数区间。从中提取： $$ \log 2 = \frac{\log 2^{\lceil p/4 \rceil}}{\lceil p/4 \rceil}, \quad \log x = \log 2^k x - k \log 2 $$
该算法的时间几乎与 $ x $ 无关。但需注意，若 $ k $ 的位数较多，可能会影响性能。
当 $ x \in [2, 2^p] $ 时，可采用以下更快的算法：
- 选择一个整数 $ m > \lg p - 2 - \lg \lg x $
- 通过重复平方计算 $ x^{2^m} $ 的区间
- 计算 $ \log x^{2^m} $ 的区间
- 将结果除以 $ 2^m $ 得到 $ \log x $
注意：当 $ x $ 接近 1 时，这种方法效率较低。
计算多个任意 $ x $ 的对数
要计算 $ \ell $ 个数的对数问题，可以将其转化为计算 $ \ell + 1 $ 个超大数的对数问题（其中一个为 2 的幂），这与 $ \ell = 1 $ 的情况类似。
当所有数都处于合理范围内且 $ \ell $ 较小时，可以通过对每个数反复平方得到 $ \ell $ 个超大数，然后计算它们的对数。但对于较大的 $ \ell $，这种方法比使用 2 的幂的方法更慢。
6. 使用椭圆积分的早期对数算法
在以下综述中，“最优”一词的意思是“在我所知道的所有技术中最快”。这并不意味着未来不能有改进。
Salamin 的方法
如文献 [2, Item 143] 所报道，Salamin 提出使用 AGM（算术-几何平均）迭代来计算超大 $ x $ 对应的 $ \pi / \log x $，从而通过已知的 $ \pi $ 来计算 $ \log x $。当 $ x $ 是一个非常大的数且 $ \pi $ 已知时，这是最优策略。
Salamin 提出的计算 $ \pi $ 的方法如下：
- 先用其他方法计算 $ e = \exp 1 $，
- 然后通过对 $ e $ 进行反复平方得到 $ e^{2^m} = \exp(2^m) $，
- 最后使用 AGM 迭代计算 $ \pi / \log(\exp(2^m)) = \pi / 2^m $。
Schroeppel 在 [2, Item 144] 中提出改用 $ 2^{2^m} $ 和 $ 2^{2^m} \cdot e $ 来代替 $ e^{2^m} $。这两种方法都比最优方案慢很多。
对于处于合理范围内的 $ x $（例如在区间 $[2, 4]$ 内），Salamin 建议通过计算 $ \log(x^{2^m}) $ 来获得 $ \log x $。当只需要计算一个对数时，这是最优策略；但如果需要计算多个对数，则不是最优的。
Brent 的方法
Brent 在文献 [10] 中提出了另一种基于不完全椭圆积分的对数算法。该算法略慢于最优方案：
- 它通过专注于中等大小的 $ b $ 值节省了一半的 AGM 步骤，
- 但它处理了超过两个不同的 $ b $ 值，导致效率下降。
Legendre-Gauss 公式与 π 的计算
Brent 和 Salamin（文献 [16]）分别提出使用 Legendre-Gauss 公式： $$ I(1, 2^{-1/2})(I(1, 2^{-1/2}) + I_1(1, 2^{-1/2})) = \frac{\pi}{2} $$ 来计算 $ \pi $。这种方法比 Salamin 之前用于计算 $ \pi $ 的方法更快，但在对数计算的上下文中并不是最优的。后续一些计算 $ \pi $ 的方法也存在类似的问题（未在此引用）。
Brent 的任意 $ x $ 对数算法
Brent 在文献 [11, Section 9] 中提出了一种适用于任意 $ x $ 的对数算法：
- 选择合适的整数 $ k $，使得 $ 2^k x $ 是一个非常大的数；
- 然后计算 $ \log(2^k x) $。
当需要计算多个对数时，这是最优策略。
Borwein 和 Borwein 的方法
Borwein 和 Borwein 在文献 [7, Section 4] 中提出了一种新的对数算法，灵感来自 Brent 的方法，但仅依赖于完全椭圆积分（而不是不完全积分）。虽然这个方法有效，但其速度仍略逊于最优方案。
Newman 的方法
Newman 在文献 [15] 中提出了一种使用 AGM 迭代计算 $ \pi $ 和超大 $ x $ 的 $ \log x $ 的方法：
- 首先计算 $ (\log x)/\pi $，
- 然后再次计算 $ (\log(x+1))/\pi \approx (\log x + 1/x)/\pi $，
- 最后相减得到结果。
需要注意的是，这种减法会导致大约 $ \lg x $ 位精度的损失。一种更好的（但仍非最优）策略是用 $ x \cdot \exp 1 $ 替换 $ x + 1 $，其中 $ \exp 1 $ 用其他方法计算。
Newman 的目标是追求最大简化性而非最大速度。他避免使用 Legendre-Gauss 公式以及任何涉及 $ I_1 $ 的方法。然而，我认为将 Landen 变换法则（即本文中的定理 2.6）视为“大量使用椭圆函数理论”是错误的；我希望我已经说服读者：$ I_1 $ 的基本性质可以像 $ I $ 的基本性质一样轻松建立。（Legendre-Gauss 公式也不难多少。）
Borwein 和 Borwein 的另一种方法
在文献 [8, Algorithm 7.2] 中，Borwein 和 Borwein 提出了一个用于计算 $ \log x $ 的方法，例如当 $ x \in [2, 4] $ 时：
- 使用 $ I_1/I $ 计算两个超大对数，分别是 $ \log(2^k) $ 和 $ \log(2^k x) $，其中 $ k $ 是一个适当选取的整数。
虽然使用 $ I_1/I $ 来计算单个超大对数是最优的，但用来计算两个或更多对数则不是最优的。
总结
本节回顾了几种使用椭圆积分和 AGM 方法计算对数的经典算法，并比较了它们的效率。尽管这些方法在历史上具有重要意义，但作者指出目前提出的算法在速度上更接近最优，尤其是在处理多个对数或高精度区间计算时表现更好。