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

Question

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?

There might be a numeric or large-value library for your language, but asking for such is off-topic (as you should know). Instead, what have you search for? What have you found? What have you tried? And what is the language you're programming in? — Some programmer dude, Aug 25 '20 at 17:05
@Someprogrammerdude Why off-topic? I am programming in Julia, but the question should apply to any language which follows floating point. I am not asking for a big float library solution. I am asking for a way of rewriting this for better accuracy. — a06e, Aug 25 '20 at 17:09
If you want a language-agnostic solution then please add that tag. If you want a Julia-specific solution then please add that tag. — Some programmer dude, Aug 25 '20 at 17:15
@Someprogrammerdude Done. I think the solution can be discussed at pseudo-code level — a06e, Aug 25 '20 at 17:29

chtz · Accepted Answer · 2020-08-25T17:15:18.333

6

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

edited Aug 25 '20 at 17:15

answered Aug 25 '20 at 17:08

chtz

17,329
4
26
56

Maybe I can just use this expression for all `x`. – a06e Aug 25 '20 at 17:09
I should branch at `x = 1`? – a06e Aug 25 '20 at 17:10
I added a suggestion how to cover all floats (no need for division) – chtz Aug 25 '20 at 17:15
@BrettHale `(4/x)` could overflow for very small `x`, so at least one should factor out another `sqrt(4)`: `sqrt(x) * sqrt((1/x) + x/4) + (x/2)` -- but you'd still need to branch for 0 (and possibly for sub-normals) *and* you need to calculate two `sqrt` (instead of one) as well as one inverse. – chtz Aug 26 '20 at 02:17
@chtz - you're right. And that `(4/x)` would propagate `+Inf` into the calculation for say, a small sub-normal value - the question explicitly only deals with `(x > 0)`. – Brett Hale Aug 26 '20 at 02:40

njuffa · Answer 2 · 2020-08-25T20:30:36.083

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.

alias · Answer 3 · 2020-08-27T17:26:00.650

3

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:

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.

edited Aug 27 '20 at 17:26

answered Aug 25 '20 at 18:39

alias

28,120
2
23
40

1

The split at `7.23e-11` looks spurious -- and indeed trying an argument like `1e-9` using the last branch gives something near `-1e27`, while it should be a value close to 1. When providing a `x>=0` precondition, it gives a result which looks more plausible, though: [link (not sure if this is stable)](https://herbie.uwplse.org/demo/8953c6a594c9835cdc2e99e57ee53948dc86e0ef.bcfc386df7dce61b9961ec7b0a896641c066a0aa/graph.html). – chtz Aug 25 '20 at 22:57
@chtz Interesting, If you think this is a bug, I'm sure they'd appreciate a ticket: https://github.com/uwplse/herbie/issues – alias Aug 26 '20 at 00:41
I've yet to see any case of herbie's output that's not either (a) unhelpful, or (b) plain wrong. Even something simple like `sqrt(1+x^2) - x` (which has numerical issues that are resolved by using `1 / (sqrt(1+x^2) + x)` for positive `x`) gives crazy output. – Mark Dickinson Aug 26 '20 at 10:38
That is really unfortunate. Did you try reporting the issues you found to them? In my somewhat limited experience, they're quite responsive to bug reports, and recording the issues you saw would definitely help. – alias Aug 26 '20 at 13:16
@alias: I'm not sure there's much constructive I could put into a bug report, that hasn't already been reported by others. For example, one of the problems I'm observing is this: https://github.com/uwplse/herbie/issues/261, which is essentially closed as a "won't fix". I may leave a comment on that issue. – Mark Dickinson Aug 26 '20 at 13:50
I've opened https://github.com/uwplse/herbie/issues/333, too. – Mark Dickinson Aug 27 '20 at 07:57
@MarkDickinson Thanks. Your concerns are valid. I've edited the answer to point out the Herbie output needs to be separately validated. – alias Aug 27 '20 at 15:49
It seems that's what really going on here is not really a bug: Herbie is working as designed, but the metrics that Herbie is using aren't necessarily a good match for what a typical user might want. – Mark Dickinson Aug 27 '20 at 16:32

score 1 · Answer 4 · answered Aug 25 '20 at 20:12

1

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;

answered Aug 25 '20 at 20:12

chux - Reinstate Monica

143,097
13
135
256

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

4 Answers4