高精度自然对数函数计算
1 背景
有时需要进行高精度的对数函数计算，这时可以使用连分式的方式进行处理。
2 自然对数函数计算
由于所有对数函数均可以通过自然对数函数来计算，因此这里考虑自然对数函数的计算方式。
2.1 函数展开
对于自然对数，有如下展开式
(1)
$$
\begin{aligned}
\ln(x)&=\ln(\dfrac{1+z}{1-z}) =\dfrac{2z}{1-\dfrac{1^2z^2}{3-\dfrac{2^2z^2}{5-\dfrac{3^2z^2}{7-\dfrac{4^2z^2}{9-\dfrac{5^2z^2}{11-\dfrac{6^2z^2}{13-…}}}}}}}
\\
\\
z &= \dfrac{x-1}{x+1}
\end{aligned}
$$
从上面公式可以发现，只有\(x\)趋近1(即\(z\)趋近0)时，收敛速度才快一些。
2.2 性质
总存在整数\(a\),以及非负整数\(b\)可以将正数\(x\)的值\(\dfrac{x}{10^a2^b}\)变换到区间\([1,2)\)。这时求对数变成
(2)
$$
\begin{aligned}
\ln(x)&=\ln(10^a2^b\dfrac{x}{10^a2^b})
\\
&=a\ln(10)+b\ln(2)+\ln(\dfrac{x}{10^a2^b})
\\
&=a\ln(5)+(a+b)\ln(2)+\ln(\dfrac{x}{10^a2^b})
\end{aligned}
$$
因此这时\(\ln(x)\)参数区间变成了\([1,2)\),这时利用1，2边界，计算这个区间\(|z|\)的最大值。这里设\(|z|\)取最大值时参数为\(x'\)，则
(3)
$$
|z|=\dfrac{x'-1}{x'+1} = \dfrac{1-\dfrac{x'}{2}}{1+\dfrac{x'}{2}} \Rightarrow x'=\sqrt{2},|z|=3-2\sqrt{2}
$$
从(3)可以发现，在\(\sqrt{2}\)处有最大值\(3-2\sqrt{2}\)，这个值依然不小。但是，如果在这个区间内，将大于等于\(\sqrt{2}\)的值除以\(\sqrt{2}\)，这时可以将\(\ln(x)\)参数区间变成\([1,\sqrt{2})\)。同理，反复多次，最终，可以将\(\ln(x)\)参数区间变成\([1,2^\frac{1}{2^n})\)，这时有
(4)
$$
0 \leq z < \dfrac{2^\frac{1}{2^n}-1}{2^\frac{1}{2^n}+1}
$$
最终，按如下算式计算对数函数
(5)
$$
\ln(x)=\begin{cases}
a\ln(5)+(a+b)\ln(2)+\ln(\dfrac{x}{10^a2^b}) &, \dfrac{x}{10^a2^b} \leq \delta
\\
\\
a\ln(5)+(a+b+\sum_{i=1}^n\dfrac{1}{2^{k_i}})\ln(2)+\ln(\dfrac{x}{10^a2^b\prod_{i=1}^n2^\frac{1}{2^{k_i}}})&, \dfrac{x}{10^a2^b} > \delta
\end{cases}
$$
公式\(k_i\)表示满足条件的整数索引。另外公式里\(\delta\)为参数值上限，这里约定\(\delta = 2^\frac{1}{2^8}\)，则\(n\)按如下方式取值
(6)
$$
\dfrac{x}{10^a2^b\prod_{i=1}^n2^\frac{1}{2^{k_i}}} < \dfrac{\delta-1}{\delta+1}<0.0013538023
$$
2.3 截断
这里假设有效位数为\(p\)，则通过如下经验公式计算截断项数\(m\)
(7)$$
m=\lceil -0.30831792741497+0.14992079831328p \rceil
$$
2.4 算法
对于\(\ln(x)\)函数，使用如下流程进行计算
① 按公式(5)(6)确定出参数\(a,b,n\)。
② 根据\(a,b,n\)计算参数\(z\)。
③ 根据公式(7)计算出截断项数\(m\)。
④ 设\(y_m=2m-1\)。
⑤ 循环计算\(y_i = 2i-1-\dfrac{i^2z^2}{y_{i+1}}\)，这里\(i = m-1,m-2,…,1\)。
⑥ 然后 \(\ln(x')=\dfrac{2z}{y_1}\)。
⑦ 按公式(5)计算最终\(\ln(x)\).
PS:\(\ln(2)、\ln(5)\)预制数据，也可以使用本算法先提前计算出来。