# Alternate expression for the following function

So if the following function is evaluated with the floating-point arithmetic, we get poor results for certain range of values of $x$. Therefore, I need to provide an alternate function that can be used for those values of $x$. The function is:

$f(x)= \sqrt{1+x}-\sqrt{1-x}$

I have found the range for this and it is: $-\sqrt{2} \le y \le \sqrt{2}$

So how would I make an alternate expression or function for this.

Do I just multiply by its conjugate.

-
Multiplying by its conjugate is a good idea. –  lhf Feb 12 '12 at 23:37
On a random sample of $10^7$ numbers in $[0,1)$, I got a maximum error of $3.33\times 10^{-16}$ between $f(x)$ as given and after multiplying by its conjugate. So I guess the two expressions have the same performance in floating-point arithmetic. I used double precision. –  lhf Feb 12 '12 at 23:48
@lhf: Did you use uniformly distributed samples? The stability problems with the naive expression ought to show up for very small $x$, so it's not likely for $10^7$ uniform samples to contain any small enough points. In fact a typical pseudorandom generator that produces uniformly distributed numbers in $[0,1)$ will return only numbers that can be subtracted from $1$ without loss of precision. –  Henning Makholm Feb 12 '12 at 23:57
Just wondering what do you mean by saying "we get poor results"? –  Emmad Kareem Feb 13 '12 at 0:03
@Emmad: if $x$ is tiny, the expression given is exceedingly prone to subtractive cancellation. –  Ｊ. Ｍ. Feb 13 '12 at 3:33

The problem with this function is that for small $|x|$ the value $f(x)$ is the difference of two almost equal numbers, as has been noted by other readers. Very often this poses serious problems (as, e.g., in numerical differentiation), but here is an easy way out (you have hinted at it yourself): Use the formula $$a-b ={a^2-b^2\over a+b}$$ in order to obtain an alternative expression for $f$: $$f(x)\ =\ {2x\over \sqrt{1+x}+\sqrt{1-x}}\ .$$

-
what if used the L'hopital's rule would that work....??? –  A_1_615 Feb 14 '12 at 0:49
@A_1_615: De l'Hopital's rule gives you a limit value in certain situations, but you cannot use it to derive an expression. –  Christian Blatter Feb 14 '12 at 9:17
This is more likely to be the expected solution than my attempt. I had considered something like this briefly but discarded it because I was afraid that rounding errors in the $1+x$ and $1-x$ subexpressions might do harm for small $x$. With a bit more thought I see that this isn't likely to be the case. (And an arcsine computation would likely be based on a hardware arc_tangent_ and have a $\sqrt{1-x^2}$ intermediate anyway). –  Henning Makholm Feb 15 '12 at 1:06

A good first step in getting a handle on a problem of this kind is to graph the function. We see (or would see, if Wolfram Alpha didn't insist on having the x and y axis scale independently) that $f(x)\approx x$ for small $x$. So the risk is that for very small $x$, adding to and subtracting from $1$ would throw away many bits of $x$, which we then won't get again by subtracting the square roots.

HINT: Write down an expression for $f(x)^2$. Does anything look trigonometrical? Suppose $x$ is the sine of something? Remember the half-angle formulas.

Edit: Here's the solution alluded to in the above hint: Let's see if we can do something about those square roots. The natural thing to try in order to make them go away is to square the entire thing (and worry about signs later): $$f(x)^2 = (1+x)+(1-x)-2\sqrt{(1+x)(1-x)} = 2 - 2\sqrt{1-x^2}$$ That doesn't immediately seem to help, but the square root that remains looks distinctly trigonometrical. Let's see if it helps to define $x=\sin \theta$. Then we can compute $$f(x)^2 = 2 - 2\sqrt{1-\sin^2\theta} = 2 - 2\cos\theta = 4\sin^2\left(\frac12\theta\right)$$ So we get (checking now that the signs are right): $$f(x) = 2\sin\left(\frac12\arcsin x\right)$$ A self-respecting floating-point implementation ought to be able to compute sines and arcsines of small arguments without significant loss of precision. Closer to $x=\pm 1$ you probably want to switch over to the naive formula, though.

-
For $x$ close to zero perhaps we can just use the series expansion: $x + \frac{x^3}{8} + \frac{7x^5}{128} + \mathcal{O}(x^7)$? –  Aryabhata Feb 13 '12 at 0:01
Im not sure I follow...so I shouldn't multiply by its conjugate....instead i should just solve it using the Taylor series...??? –  A_1_615 Feb 13 '12 at 2:35
@A_1: My suggestion is to square the function (multiply it by itself, not its conjugate). Aryabhata has a different (unrelated) suggestion. –  Henning Makholm Feb 13 '12 at 2:48
@Aryabhata: or a Padé approximant, for that matter.... –  Ｊ. Ｍ. Feb 15 '12 at 1:04
@J.M.: I am not too familiar with them. Thanks for mentioning those! –  Aryabhata Feb 15 '12 at 1:15