一个迭代过程如何导致对数平均

hbghlyj

一个关于迭代如何导致意想不到结果的例子
取 $a_0=x,b_0=y$ 作为初值，然后迭代
$$a_1=\frac{a_0+\sqrt{a_0b_0}}{2},~~~~b_1=\frac{b_0+\sqrt{a_0b_0}}{2}$$
$$a_{n+1}=\frac{a_n+\sqrt{a_nb_n}}{2},~~~~b_{n+1}=\frac{b_n+\sqrt{a_nb_n}}{2}$$每次迭代后 \(a_n\) 和 \(b_n\) 之差减半，即 \(a_n - b_n = \frac{x - y}{2^n}\to0\)，则 $a_n$ 和 $b_n$ 都趋于同一个极限 $L$

求证 $L=\frac{x-y}{\ln x-\ln y}$？
证明：验证$$\frac{a_{n+1}-b_{n+1}}{\ln a_{n+1}-\ln b_{n+1}}=\frac{\frac{a_n-b_n}2}{\frac{\ln a_n-\ln b_n}2}=\frac{a_n-b_n}{\ln a_n-\ln b_n}$$
对于所有 $n$，$\frac{a_n-b_n}{\ln a_n-\ln b_n}$ 均不变，所以 $L=\frac{x-y}{\ln x-\ln y}$.
证毕

isee 发表于 2025-1-12 10:11
问：如何徒手算出10000以内任意整数的常用对数？

例如 $L(1,2)=\frac1{\ln2}$ 用 $a_n,b_n$ 可算到任意精度

1. FixedPoint[(# + Sqrt[#[[1]] #[[2]]])/2 &, {1., 2.}]

Copy the Code

ZCos666

置$x_n=\sqrt{\dfrac{a_n}{b_n}}$,则有

$x_{n+1}=\sqrt{\dfrac{a_n+\sqrt{a_nb_n}}{b_n+\sqrt{a_nb_n}}}=\sqrt{\dfrac{x_n^2+x_n}{1+x_n}}=\sqrt{x_n}$

于是有$\dfrac{a_n}{b_n}=\left(\dfrac{a_1}{b_1}\right)^{\frac{1}{2^{n-1}}},a_n-b_n=\dfrac{a_1-b_1}{2^{n-1}}$

解得$b_n=\dfrac{\dfrac{a_1-b_1}{2^{n-1}}}{\left(\dfrac{a_1}{b_1}\right)^{\frac{1}{2^{n-1}}}-1}$

取极限，洛必达得

$\displaystyle\lim_{n\to\infty}b_n=\lim_{n\to\infty}\dfrac{\dfrac{a_1-b_1}{2^{n-1}}\ln2}{\dfrac{1}{2^{n-1}}\left(\dfrac{a_1}{b_1}\right)^{\frac{1}{2^{n-1}}}\ln2\ln\dfrac{a_1}{b_1}}=\dfrac{a_1-b_1}{\ln{a_1}-\ln{b_1}}$

对数的近似计算本质上是利用了$\ln{x}\sim n\left(x^{\frac{1}{n}}-1\right)$

hbghlyj

是否与算术几何平均（AGM）的计算相关？不过这里的情况好像不同，因为AGM中的迭代是取算术平均和几何平均，而这里的每个a和b都是它们自己原来的值与几何平均的算术平均，所以极限不同

著名的高斯算术-几何平均。设 $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

欲计算$ a_{0}=1$和$ g_{0}=10000$的AGM，首先算出它们的AM和GM：

$ 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的AGM是两个数列的公共极限，约为1482.35393858815
因此，1和1/10000的AGM约为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

首先总结椭圆函数理论中的必要公式。模数 $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}
\]

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

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)呢