Question

我需要在数值上计算sqrt(1 + (x/2)^2) + x/2,对于正x。对于非常大的x值,使用此表达式直接失败。如何重写它以获得更准确的评估?

6

4

3

1

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]，我建议只在该点分支并返回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

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

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

Answer 3

注意有人指出，Herbie 生成的表达式可能不适用于所有上下文。特别是，Herbie 用于“改进”表达式的指标可能会生成在特定情况下性能较差的表达式。因此，请使用谷物盐获取其输出。我认为您仍然可以咨询 Herbie 以获得想法，但是不要将其用作直接替换。

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 的输出可能很难阅读，据报道它可能不会更好，但也许它可以提供另一种观点。