这是一个自我回答的问题，正如 Stackoverflow 上“鼓励”分享知识一样。因此，它是我之前关于精确计算 Lambert W 函数主分支 W0 的“自我回答问题”的后续内容。有关 Lambert W 函数的一些一般背景信息，请参阅此处。当我最初研究兰伯特 W 函数时，负分支 W-1 对我来说似乎并不是特别有用，但此后我遇到了许多实际使用它的出版物，例如 [1]， [2]、[3]、[4]. 实值 Lambert-W 函数的负分支 W-1 的输入域非常窄：
[-1/e, 0]。由于分支点位于 -1/e 且渐近行为接近于零，精确计算有些困难。如何使用 ISO-C99 标准数学库准确计算实值 Lambert W 函数 W-1 的负分支？
准确计算应定义为不超过 4 ulps 的最大误差，这与适用于英特尔数学库的 LA 配置文件的误差界限相同，并且我发现对于非平凡超越函数来说这是一个合理的误差界限实际使用。可以假定支持 IEEE-754 (2008) 二进制浮点算术以及对融合乘加 (FMA) 运算的功能正确支持，可通过 fma() 和
fmaf()标准库函数进行访问。
[1] A. C. Gormaz-Matamala 等人。 “使用兰伯特 W 函数解决热大质量恒星线驱动风的新流体动力学解决方案”，920:64（30 页），2021 年 10 月 [2] Santiago Pindado 等人，“应用于光伏系统建模时的简化兰伯特 w 函数数学方程。”
IEEE 行业应用汇刊，卷。 57，第 2 期，2021 年 1 月，第 1779-1788 页 [3] 华金·F·佩德雷耶斯。等人，“恒功率应用中基于兰伯特 W 函数的超级电容器电气变量的闭合形式表达式”。
能源，2021 年 3 月，218：119364 [4] Armengol Gasull 等。 al，“类威布尔分布的最大值和兰伯特 W 函数。”
测试24（2015）：714-733 我将首先介绍映射到 IEEE-754
lambert_wm1f()float。
与主分支的计算一样，函数迭代在-1/e处的分支点附近在数值上不稳定。因此，使用由 Remez 算法生成的多项式极小极大近似
p(t)，其中 t=√(x+1/e)
。 对于函数迭代，使用文献中最大相对误差为 1.488e-4 的起始近似值。只需进行一次二次收敛迭代即可达到完全的单精度精度。我通过实验确定牛顿迭代在接近零时效果最好，而使用 Iacono-Boyd 迭代比输入域的平衡更有优势。 当针对英特尔标准 C 数学库的实现构建时，在详尽的测试中，下面的 ISO C-99 实现的最大误差小于 2.5 ulp；标准 C 数学库的任何其他高质量实现都应该可以实现类似的误差范围。有关常用实现的准确性的概述，请参阅： Brian Gladman、Vincenzo Innocente、John Mather、Paul Zimmermann，“单精度、双精度、双扩展精度和四精度数学函数的准确性”。 HAL 预印本 hal-03141101v6
，2024 年 2 月
binary32映射到 IEEE-754
/*
Compute the negative branch of the real-valued Lambert W function, W_(-1). Maximum
error when built against the Intel MKL: 2.41426181 ulps @ -7.32652619e-2
*/
float lambert_wm1f (float x)
{
const float c0 = 1.68090820e-1f; // 0x1.584000p-3
const float c1 = -2.96497345e-3f; // -0x1.84a000p-9
const float c2 = -2.87322998e-2f; // -0x1.d6c000p-6
const float c3 = 7.07275391e-1f; // 0x1.6a2000p-1
const float ln2 = 0.69314718f;
const float l2e = 1.44269504f;
const float em1_fact_0 = 0.625529587f; // exp(-1)_factor_0
const float em1_fact_1 = 0.588108778f; // exp(-1)_factor_1
float redx, r, s, t, w, e, num, den, rden;
if (isnan (x)) return x + x;
if (x == 0.0f) return -INFINITY;
if (x > 0.0f) return INFINITY / INFINITY; // NaN
redx = fmaf (em1_fact_0, em1_fact_1, x); // x + exp(-1)
if (redx <= 0.09765625f) { // expansion at -(exp(-1))
r = sqrtf (redx);
w = -3.30250000e+2f; // -0x1.4a4000p+8
w = fmaf (w, r, 3.53563141e+2f); // 0x1.61902ap+8
w = fmaf (w, r, -1.91617889e+2f); // -0x1.7f3c5cp+7
w = fmaf (w, r, 4.94172478e+1f); // 0x1.8b5686p+5
w = fmaf (w, r, -1.23464909e+1f); // -0x1.8b1674p+3
w = fmaf (w, r, -1.38704872e+0f); // -0x1.6315a0p+0
w = fmaf (w, r, -1.99431837e+0f); // -0x1.fe8ba6p+0
w = fmaf (w, r, -1.81044364e+0f); // -0x1.cf793cp+0
w = fmaf (w, r, -2.33166337e+0f); // -0x1.2a73f2p+1
w = fmaf (w, r, -1.00000000e+0f); // -0x1.000000p+0
} else {
/* Initial approximation based on: D. A. Barry, L. Li, and D.-S. Jeng,
"Comments on 'Numerical Evaluation of the Lambert Function and
Application to Generation of Generalized Gaussian Noise with
Exponent 1/2'", IEEE Transactions on Signal Processing, Vol. 52,
No. 5, May 2004, pp. 1456-1457
*/
s = fmaf (log2f (-x), -ln2, -1.0f);
t = sqrtf (s);
w = -1.0f - s - (1.0f / fmaf (exp2f (c2 * t), c1 * t,
fmaf (1.0f / t, c3, c0)));
if (x > -0x1.0p-116f) {
/* Newton iteration, for arguments near zero */
e = exp2f (fmaf (w, l2e, 32));
num = fmaf (w, e, -0x1.0p32 * x);
den = fmaf (w, e, e);
rden = 1.0f / den;
w = fmaf (-num, rden, w);
} else {
/* Roberto Iacono and John Philip Boyd, "New approximations to the
principal real-valued branch of the Lambert W function",
Advances in Computational Mathematics, Vol. 43, No. 6,
December 2017, pp. 1403-1436
*/
t = w / (1.0f + w);
w = fmaf (logf (x / w), t, t);
}
}
return w;
}
的双精度实现 lambert_wm1()
double