CODEWITHSUNDEEP
floating-pointlanguage-agnosticnumericnumerical-computing
a06e
a06e

Reputation: 20774

Accurate sqrt(1 + (x/2)^2) + x/2

I need to compute sqrt(1 + (x/2)^2) + x/2 numerically, for positive x. Using this expression directly fails for very large values of x. How can I rewrite it to obtain a more accurate evaluation?

Upvotes: 3

Views: 874

Answers (4)

alias
alias

Reputation: 30460

Note It's been pointed out that the Herbie generated expression may not be suitable in all contexts. In particular, metrics used to "improve" the expression by Herbie may generate expressions that perform worse for your particular scenario. So, take its output with a grain salt. I think you can still consult with Herbie to get an idea, but do not use it as a drop-in replacement.

Herbie (https://herbie.uwplse.org/) recommends the following replacement for your expression:

enter image description here

Or, in C:

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

becomes:

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;
}

It generates a detailed report on why it splits it into three regions. The output of Herbie can be rather hard to read, and it's been reported that it may not be better, but perhaps it can provide an alternative view.

Upvotes: 5

njuffa
njuffa

Reputation: 26175

Many programming languages offer a function hypot(x,y) that computes sqrt (x*x + y*y) while avoiding overflow and underflow in intermediate computation. Many implementations of hypot also compute the result more accurately than the naive expression. These advantages come at the expense of a moderate increase in run time.

With this function, the given expression can be written as hypot (1.0, 0.5*x) + 0.5*x. If your programming language of choice does not support hypot or an equivalent function, you may be able to adapt the implementation I provided in this answer.

Upvotes: 6

chux
chux

Reputation: 154075

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.

Code may then get a better result with hypot(1,x/2) + x/2;

Upvotes: 1

chtz
chtz

Reputation: 18827

For very large x you can factor out an 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

For x > 2/sqrt(eps) the square root will actually evaluate to 1 and your whole expression will simplify to just x. Assuming you need to cover the entire range [0, infinity], I would suggest just branching at that point and return x in this case and your original formula, otherwise:

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

Upvotes: 7

Related Questions