Poor results for sqrt(1+x) - sqrt(1-x)

[mdickinson]

If I enter sqrt(1+x) - sqrt(1-x) into the web demo, it suggests unconditionally transforming it to x+(x^3⋅0.125+0.0546875⋅x^5). That's fine for small x (e.g., for abs(x) < 1e-3), since the truncation error from truncating the Taylor series is on the order of x^7, but is much much worse than the original expression for larger x. A better rewrite, that's numerically stable across the whole of the [-1.0, 1.0] domain, would be 2x / (sqrt(1+x) + sqrt(1-x)).

Here's the screenshot:

The graph in the screenshot is confusing: if I'm understanding correctly, the red line is supposed to indicate the number of bits error for the original expression when computed directly using IEEE 754 binary64 arithmetic, and the blue line indicates the number of bits error for the rewritten expression. Is that correct? I'm not sure what "bits error" means here: is that essentially measuring the error in ulps, or is it counting how many bits of the result are incorrect?

Either way, at the extreme right-hand-side of the graph, there's an input 0.8993 shown. For that input (or rather, for the input 0.8992999999999999882760448599583469331264495849609375, which is the closest exactly-representable binary64 value), the exact result is 1.06081830207766257859364328547738223076248106920543..., a direct computation using binary64 floating-point (with a correctly-rounded sqrt operation) gives exactly 1.0608183020776624783110264615970663726329803466796875, and the transformed expression (again computed using binary64 floating-point) gives exactly 1.022379575763839643087749209371395409107208251953125.

So from my calculations, the error for the direct expression is approximately 0.45 ulps, while the error for the transformed expression is a whopping 173 million million ulps. This seems to contradict what's shown in the graph, which seems to say that the "bits error" for the transformed result at 0.8993 is approximately 4.

Meta: I'm seeing answers on Stack Overflow that suggest using Herbie to rewrite and then just using the result directly. For example, the author of this answer says:

The output of Herbie can be rather hard to read, but you can usually just use it as a drop-in replacement.

(and the replacement suggested in that answer is also rather bad).

[pavpanchekha]

Thanks for the bug report! Let me what I can answer...

First, bits of error is absolutely measuring ULPs of error and taking the logarithm, not strictly counting bits. More precisely, it is computing the number of floating point values between the correctly-rounded and the computed floating-point outputs. This is the same as ULPs of error for small errors, but when there are many orders of magnitude between the correct and incorrect answer it will differ. I think of it as "incremental" ULPs, in that it satisfies the triangle inequality.

The output is definitely inaccurate, as you point out. The flipped version you recommend (2x / ...) is better. It is likely better even in Herbie's internal measures, or maybe would be if only Herbie had more sampled points.

Herbie equally weights all floating-point input values in its error measure. This means that most of [-1, 1] is the very small values in that range (half are in [-1e-150, 1e-150] or so). So from Herbie's point of view, while the flipped version is better it is probably only a very little bit better.[1] The error in Herbie's expression is highest at the extreme value of +/-1, where it is something like 30% relative error, but Herbie doesn't weight that part of the input space very high.

The graph only seems to go up to 4 bits of error because the thick read line is a smoothing (using a window of 32 sample points, I think) of the underlying data, which are the pale red / blue dots. That smoothing can be deceptive at the edges, where the window stops being possible to make; I might try to fix that.

As for using Herbie's results directly—yeah, Herbie's results often do have either some kind of bizarre weirdness or some kind of obvious inaccuracy, and using its results directly isn't advisable if you have access to a numerical analyst. I suspect most people on SO do not. Doesn't mean we shouldn't make Herbie better though.

Anyway, I'll take a look internally and see if Herbie actually generates the flipped version and, if it does, why it doesn't go with it.

[1]: You might confirm this by adding a :herbie-target (/ (* 2 x) (+ (sqrt (+ 1 x)) (sqrt (- 1 x)))) field to the FPCore (you'd need to do this with file input, not the web demo, sadly). That will cause Herbie to draw a green error line for that expression; I expect it to be very accurate.

[pavpanchekha]

Herbie now regularly achieves the correct suggested result, so marking this fixed.

On the subtler issue of people using Herbie’s results directly, uhh, can’t really fix that but I think the new Pareto mode, by giving people multiple options, is going to help make it a bit less automatic.