On 17 July 2010 20:19, Tom Lane <tgl@sss.pgh.pa.us> wrote:
> ... For instance, if x is 1e8 * y, then y*y
> fails to affect the sum at all (given typical float8 arithmetic), and
> you'll get back sqrt(x*x) even though y should have been able to affect
> the result at the 8th place or so. In the patch's calculation, y/x is
> computed accurately but then we'll lose the same precision when we form
> 1 + yx*yx --- the result will be just 1 if y is lots smaller than x.
>
No. If x is 1e8 * y, then y will only affect the result in the 16th
place. You can see this if you do a simple series expansion:
sqrt(1+yx^2) = 1 + 1/2 yx^2 + O(yx^4)
> If we were feeling tense about this, we could look for an alternate way
> of calculating sqrt(1 + yx*yx) that doesn't lose so much accuracy.
> In principle I think that's doable since this expression is related to
> ln(1+x) which can be calculated accurately even for very small x.
This algorithm is about as accurate as it could possibly be. The point
with ln(1+x) is that for small x:
ln(1+x) = x + O(x^2)
so you would loose precision if x were much smaller than 1. This is
not the case with sqrt(1+x).
For most cases, the new algorithm is no more accurate than the old
one. The exception is when *both* x and y are very small. In this
case, it protects against incorrect underflows to 0.
Regards,
Dean