(Cephes) log1p的近似

Question

log1p

的 Cephes 实现在 documentation 中被描述为

Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)

其中 P(x) 和 Q(x) 是两个多项式（Q(x) 是一元多项式。

这部分是

log(1 + x)

的麦克劳林级数，但目前尚不清楚他们如何获得 P(x) 和 Q(x)。

是否有这种形式的近似的名称或有关如何复制该过程以获得系数的资源？

Answer 1

这是有理逼近的一种特殊形式，有时称为 Cody & Waite 形式。 Cody & Waite 形式是一种约束近似形式，可用于多项式近似和有理近似，并且在数值方面具有优势。

当前导系数以给定的浮点格式精确表示时，近似精度会提高。因此，与其使用诸如 sinh(x) ≈ a₀x + a₁x³ + ... + a_nx²ⁿ⁺¹ 之类的近似值，不如强制前导系数为 1 : sinh(x) ≈ x + a₁x³ + ... + a_nx²ⁿ⁺¹，因为函数参数 x 对于构建数学库实现而言被认为是精确的。

这种安排被称为 Cody & Waite 形式，因为它在以下书籍中被广泛使用，这本书是 1980 年代和 1990 年代大部分时间数学图书馆作者的标准参考：

William J. Cody 和 William Waite，基本功能软件手册，Englewood Cliffs，新泽西州：Prentice Hall 1980。

然而，这种排列的最初使用可以追溯到 Hirondo Kuki，例如，他将它用于 asin() 和 sinh() 的实现。参见：

Hirondo Kuki，数学函数。中心的 7094 Fortran II 数学函数库的描述， 芝加哥大学计算中心报告，1966 年 2 月。

Cephes 库的作者 Stephen L. Moshier 几乎专门使用了 minimax 近似值。这些是 minimize the maximum error over the approximation interval 的近似值：根据 L_∞ 范数优化的近似值。对于 n 次多项式极小极大逼近 P_n(x)，误差函数 d(x) = P(x) - f(x) 的图形具有特征：它在 x 轴上来回摆动，并且在近似区间恰好在 n+2 个点达到最大值，其中所有最大值具有完全相同的大小但相邻最大值的符号交替。

这导致一组 n+2 个方程，其中 n+2 个未知数 P(x) - f(x) ± d = 0。1934 年，一位苏联数学家设计了一种迭代 numerical 算法来最小化 d，通常称为 Remez 交换算法，其中每次迭代都需要求解线性方程组：

Evgeny Yakovlevich Remez，“Sur un procédé convergent d'approximations successives pour determiner les polynômes d'approximation”，Compt。撕裂。学院。 Sc. 198, 2063 (1934)

该算法至今仍在使用，并已扩展用于极小极大有理逼近，其中误差函数 d(x) = P_m(x)/Q_n(x) - f(x) 具有 m +n+2 最大值，同样大小但符号交替。 Moshier 写了一本描述 Cephes 库底层设计的书，其中包括 Remez 算法的实现：

Stephen L. B. Moshier，数学函数的方法和程序， Chichester：Ellis Horwood 1989.

正如书中指出的，Cody & Waite 形式的有理逼近 f(x) ≈ x + x³ P(x²) / Q(x^²)，导致方程组形式为 x³ P(x²) - (f(x) - x ± d) Q(x²) = 0。问题中 log1p(x) 核心近似的具体形式扩展这是通过约束前两个系数来实现的。

虽然可以构建自己的 Remez 算法实现，但根据我的经验，创建稳健的实现并非易事。如果不想从头开始，可以考虑使用 Netlib 存储库 中的 Moshier 的 Remez 实现作为他们努力的潜在起点。免费开源Sollya 工具 具有生成极小极大多项式 逼近的优秀工具，但最后我检查它无法生成极小极大有理逼近。人们可能想为此寻找像 Mathematica 和 Maple 这样的商业产品（我都没有用过）。