c++ - 为 [-1, 1] 中的 c 计算 sqrt((b²*c²)/(1-c²)) 的数值稳定方法

Question

对于一些实际值(value)b和 c在 [-1, 1] , 我需要计算sqrt( (b²*c²) / (1-c²) ) = (|b|*|c|) / sqrt((1-c)*(1+c))当 c 时，分母中出现灾难性抵消。接近 1 或 -1。平方根可能也无济于事。
我想知道是否可以在这里应用一个聪明的技巧来避免 c=1 和 c=-1 周围的困难区域？

Answer 1

这种稳定性方面最有趣的部分是分母，sqrt(1 - c*c) .为此，您只需将其扩展为 sqrt(1 - c) * sqrt(1 + c) .我不认为这真的有资格作为一个“聪明的把戏”，但这就是所需要的。
对于典型的二进制浮点格式(例如 IEEE 754 binary64，但其他常见格式应该表现得同样好，除了不愉快的东西，如 double-double 格式)，如果 c接近 1然后 1 - c将由 Sterbenz' Lemma 精确计算, 而 1 + c没有任何稳定性问题。同样，如果 c接近 -1然后 1 + c将被精确计算，1 - c将被准确计算。平方根和乘法运算不会引入重大的新误差。

这是一个数值演示，在具有 IEEE 754 binary64 浮点和正确舍入 sqrt 的机器上使用 Python手术。
让我们来个c接近(但小于)1 :

>>> c = float.fromhex('0x1.ffffffff24190p-1')
>>> c
0.9999999999

我们在这里要小心一点:注意显示的十进制值 0.999999999 , 是 c 的精确值的近似值.确切的值如十六进制字符串或小数形式的构造所示，562949953365017/562949953421312 ，这正是我们关心获得良好结果的确切值。
表达式 sqrt(1 - c*c) 的确切值，在点后四舍五入到小数点后 100 位，是:

0.0000141421362084401590649378320134409069878639187055610216016949959890888003204161068184484972504813

我使用 Python 的 decimal 计算了这个值模块，并使用 Pari/GP 仔细检查了结果.这是 Python 的计算:

>>> from decimal import Decimal, getcontext
>>> getcontext().prec = 1000
>>> good = (1 - Decimal(c) * Decimal(c)).sqrt().quantize(Decimal("1e-100"))
>>> print(good)
0.0000141421362084401590649378320134409069878639187055610216016949959890888003204161068184484972504813

如果我们天真地计算，我们会得到以下结果:

>>> from math import sqrt
>>> naive = sqrt(1 - c*c)
>>> naive
1.4142136208793713e-05

我们可以很容易地计算出 ulps 错误的近似数量(对于正在进行的类型转换量表示歉意 - float 和 Decimal 实例不能直接在算术运算中混合):

>>> from math import ulp
>>> float((Decimal(naive) - good) / Decimal(ulp(float(good))))
208701.28298527992

所以天真的结果是几十万 ulps - 粗略地说，我们已经失去了大约 5 个小数位的准确性。
现在让我们尝试扩展版本:

>>> better = sqrt(1 - c) * sqrt(1 + c)
>>> better
1.4142136208440158e-05
>>> float((Decimal(better) - good) / Decimal(ulp(float(good))))
-0.7170147200803595

所以在这里我们的准确度优于 1 ulp 错误。不完全正确地四舍五入，但下一个最好的事情。
再做一些工作，应该可以在表达式 sqrt(1 - c) * sqrt(1 + c) 中陈述和证明 ulps 错误数的绝对上限。 , 在域 -1 < c < 1 ，假设 IEEE 754 二进制浮点、舍入到偶数舍入模式以及自始至终正确舍入的操作。我还没有这样做，但如果这个上限超过 10 ulps，我会感到非常惊讶。