1

I read that the expm1 function is suited for computing 1 - exp(x) for small x, without losing precision from truncation of 1.0 at ~15 digits (for doubles). Is there such a function for 1 - sqrt(x)? At the moment I am simply using a very large Taylor expansion, because I require as many (preferably all) digits of precision as doubles/long doubles can offer.

Edit: I badly confused my intention here: I am looking to compute 1 - sqrt(1-x) for x anywhere between 10^-12 and 1.

zjw518
  • 13
  • 6
  • 4
    *"for computing 1 - exp(x) "*... then you read backwards :) it is for computing `exp(x) - 1` – Cory Kramer Jun 20 '17 at 19:46
  • A search of SO would have helped you: https://stackoverflow.com/questions/8721022/how-to-improve-fixed-point-square-root-for-small-values – duffymo Jun 20 '17 at 20:21
  • How is `double foo(double x) { return 1 - sqrt(x); }` not sufficient for your needs with small `x`? – chux - Reinstate Monica Jun 20 '17 at 21:47
  • 1
    @duffymo - that would be because I mis-stated the question. I meant 1 - sqrt(1-x). Should I repost the question or edit this one? – zjw518 Jun 20 '17 at 22:22
  • @zjw518 Some libraries (notably Boost) offer a function `sqrt1pm1()` that computes sqrt(x+1)-1, which should be suitable for your purposes, since 1 - sqrt(1-x) = - sqrt1pm1(-x). You could always create a reasonable implementation yourself: `double sqrt1pm1 (double a) { return a / (1.0 + sqrt (a + 1.0)); }` – njuffa Sep 06 '17 at 04:40

2 Answers2

5

The question here seems badly motivated. While exp(x) converges to 1 as x goes to 0, meaning that given the same floating point precision exp(x)-1 has more significant figures than exp(x) for small x, this is not true for sqrt(x), which converges to 0 as x goes to 0. In other words exp(x)-1 can be made fractionally more precise than exp(x) for small x, but the same is not true for 1-sqrt(x) -- which would in fact get worse, since you're taking it from something near 0 (1e-6) to something near 1 (0.999999).

If on the other hand you instead wanted to calculate sqrt(1+x) for very small x (as an accurate measurement of sqrt(x) very near x=1), sqrt(1+x)-1 would be a more accurate floating point computation. And its Taylor series would work very well; I find that for |x| < 1e-9, x/2 - x^2/8 + x^3/16 is a good approximation of sqrt(1+x)-1 to within an RMS fractional error of 3e-29 (with a maximum of 8e-29 on the edges) -- twice as many digits as are accurate in a double. Even the quadratic approximation is probably good enough (with roughly 20 digits of accuracy)

jwimberley
  • 1,696
  • 11
  • 24
  • 6
    Or, you could use `x / (sqrt(1+x) + 1)` which is algebraically equivalent to `sqrt(1+x) - 1` but doesn't have the cancellation artifacts for small `x`. – Daniel Schepler Jun 20 '17 at 20:47
  • 1
    You are correct, @jwimberley, I meant `1 - sqrt(1-x)`. My intent is to have a function which works for x as small as 10^-12 up to 1. I did a poor job phrasing the question (I confused myself after reading about `expm1`). I currently have ~20 terms in the expansion to be safe, so that I can switch to using `sqrt` at x~.01 - it would be great if a function existed that encapsulates this (and is faster, of course). Though Daniel's suggestion looks promising. – zjw518 Jun 20 '17 at 22:29
  • You can modify either of our suggestions to your purpose. @DanielSchepler's is probably better: just compute -x/(1+sqrt(1-x)) in that order. This is equal to x*(1-sqrt(1-x))/(1-(1-x)) = 1-sqrt(1-x), but doesn't rely on cancellation to achieve accuracy. – jwimberley Jun 20 '17 at 22:34
  • (My -x in the above comment a typo) – jwimberley Jun 20 '17 at 22:40
1

the expm1 function is suited for computing 1 - exp(x).

As can read in the ref:

For small magnitude values of x, expm1 may be more accurate than exp(x)-1.

Is there such a function for 1 - sqrt(x)?

No, at least not in the standard headers.

gsamaras
  • 71,951
  • 46
  • 188
  • 305
  • Note that 1 - exp(x) = - expm1(x), so expm1() is in fact well suited for this computation. – njuffa Sep 06 '17 at 04:30
  • You stated: "It's the other way around", which I took as a *contradiction* to "expm1 .. is suited for computing 1 - exp(x)." – njuffa Sep 06 '17 at 05:08
  • 1
    I am familiar with the problem: posting before the first cup of coffee is not advised :-) – njuffa Sep 06 '17 at 05:15