准确的 sqrt(1 + (x/2)^2) + x/2

Question

我需要计算 sqrt(1 + (x/2)^2) + x/2 的数值，以获得正数 x。对于非常大的 x 值，直接使用此表达式会失败。我怎样才能重写它以获得更准确的评估？

Answer 1

对于非常大的 x 你可以分解出一个 x/2:

sqrt(1 + (x/2)^2) + x/2
 = (x/2) * sqrt( 1/(x/2)^2 + (x/2)^2/(x/2)^2) + x/2
 = (x/2) * sqrt( (2/x)^2 + 1 ) + x/2

对于 x > 2/sqrt(eps)，平方根实际上计算为 1，您的整个表达式将简化为 x。 假设你需要覆盖整个范围 [0, infinity]，我建议只在那个点分支，在这种情况下 return x 和你的原始公式，否则：

if x > 2/sqrt(eps)  // eps is the machine epsilon of your float type
    return x
else
    return sqrt(1 + (x/2)^2) + x/2

Answer 2

许多编程语言提供了一个函数 hypot(x,y) 来计算 sqrt (x*x + y*y)，同时避免中间计算中的上溢和下溢。 hypot 的许多实现也比朴素表达式更准确地计算结果。这些优势是以适度增加 运行 时间为代价的。

有了这个函数，给定的表达式可以写成hypot (1.0, 0.5*x) + 0.5*x。如果您选择的编程语言不支持 hypot 或等效函数，您可以调整我在 this answer.

中提供的实现

Answer 3

注意 有人指出，Herbie 生成的表达式可能不适合所有上下文。特别是，Herbie 用于“改进”表达式的指标可能会生成对您的特定场景表现更差的表达式。因此，将其输出与谷物盐一起使用。我认为您仍然可以咨询 Herbie 以获得想法，但是不要将其用作drop-in替代品。

Herbie (https://herbie.uwplse.org/) 推荐以下表达式替换：

或者，在 C 中：

double code(double x) {
    return ((double) (((double) sqrt(((double) (1.0 + ((double) pow((x / 2.0), 2.0)))))) + (x / 2.0)));
}

变为：

double code(double x) {
    double VAR;
    if (((x / 2.0) <= -8569.643649604539)) {
        VAR = (1.0 / ((double) ((1.0 / ((double) pow(x, 3.0))) - ((double) (x + (1.0 / x))))));
    } else {
        double VAR_1;
        if (((x / 2.0) <= 7.229769585372425e-11)) {
            VAR_1 = ((double) ((x / 2.0) + ((double) sqrt(((double) (1.0 + ((double) pow((x / 2.0), 2.0))))))));
        } else {
            VAR_1 = ((double) ((x / 2.0) + ((double) (((double) ((1.0 / x) + ((double) (x * 0.5)))) - (1.0 / ((double) pow(x, 3.0)))))));
        }
        VAR = VAR_1;
    }
    return VAR;
}

它会生成一份详细报告，说明将其分成三个区域的原因。 Herbie 的输出可能相当难以阅读，据报道它可能不会更好，但也许它可以提供另一种观点。

Answer 4

hypot()

The hypot functions compute the square root of the sum of the squares of x and y, without undue overflow or underflow. A range error may occur.

代码可能会通过 hypot(1,x/2) + x/2;

获得更好的结果