3

i have a function f(x)=$\sqrt{x+1/x}$ - $\sqrt{x-1/x}$ (x $\ge{1}$) and want to calculate $f(2.0004 * 10^{18})$ in floating point arithmetic with a mantissa t = 4.

First, I rounded the input to 2 * $10^{18}$ (since the 4 at the end would require a mantissa of at least 5).

Second, calculating 1/(2*$10^{18}$) we get 5.0 * $10^{-19}$.

Adding/Subtracting this value to/from 2*10$^{18}$ inside the square root however has to be rounded to 2*10$^{18}$ since $10^{-19}$ is extremely small and would require a mantissa of at least 19 to affect the result.

This way I get an result of $\sqrt{2*10^{18}}$ - $\sqrt{2*10^{18}}$ = 0, which is not the expected result according to my notes.

Also the problem seems to be numericaly well conditioned for large input values, so I would expect only a small error in the end.

Is this problem caused by the given algorithm being "bad", e.g. is there a way to write the same function avoiding such errors (e.g. by avoiding addition which is bad conditioned in case x1 ~= -x2 or is there simply an error in my calculations.

Thank you in advance for any help and suggestions

Philipp
  • 155

1 Answers1

2

Yes your formula suffers from 'catastrophic cancellation'. You can transform it to a formula with better properties: $$\sqrt{x+1/x} - \sqrt{x-1/x}= \sqrt{\frac{x^2+1}{x}}-\sqrt{\frac{x^2-1}{x}}$$ $$= \frac{1}{\sqrt{x}}\left(\sqrt{x^2+1}-\sqrt{x^2-1}\right)$$ $$= \frac{1}{\sqrt{x}}\left(\sqrt{x^2+1}-\sqrt{x^2-1}\right) \frac{\left(\sqrt{x^2+1}+\sqrt{x^2-1}\right)}{\left(\sqrt{x^2+1}+\sqrt{x^2-1}\right)}$$ $$=\frac{1}{\sqrt{x}}\frac{(x^2+1) - (x^2-1) }{\sqrt{x^2+1}+\sqrt{x^2-1}}$$ $$f(x)=\frac{1}{\sqrt{x}}\frac{2}{\sqrt{x^2+1}+\sqrt{x^2-1}}$$ With this form of $f(x)$ you get $f(2.0004\cdot10^{18}) = 0.35344\cdot10^{-27}$.

And if you are working with large numbers for which $x^2 \pm 1 \approx x^2$ (e.g. if $x>1000$ in your case), you can simplify the function even further: $$f(x) \approx \frac{1}{x \sqrt{x}}, \quad \text{if}\; x^2-1 \approx x^2$$

gammatester
  • 18,827
  • Thank you so much, its exactly the kind of answer I was looking for! – Philipp Sep 12 '13 at 09:44
  • @Philipp: I justed corrected a typo in the first expression. And if the answer is useful please consider an up-vote. – gammatester Sep 12 '13 at 10:06
  • Thank you! I would love to give you an upvote, however I am new to math.stackexchange and don´t have the 15 required rep-points. I will give you an upvote as soon as I get the final two points :) – Philipp Sep 12 '13 at 10:44
  • @Philipp: I see. I upvoted your question from Tuesday (I have not seen it until now) and which is related but not identical to this one. – gammatester Sep 12 '13 at 10:54
  • Thank you! Just upvoted your answer :) – Philipp Sep 12 '13 at 12:41