BUG #5711: input out of error with haversine formual

On 15 October 2010 05:58, Vince <vincecarney@gmail.com> wrote:

The following bug has been logged online:

Bug reference:      5711
Logged by:          Vince
Email address:      vincecarney@gmail.com
PostgreSQL version: 8.4
Operating system:   Linux
Description:        input out of error with haversine formual
Details:

The following will return an input out of error as the acos() function
cannot be -1 <= x <= 1.

SELECT * FROM
               (SELECT *, (3959 * acos(cos(radians(37.7438640)) *
cos(radians(37.7438640)) * cos(radians(-97.4631299) -
               radians(-97.4631299)) + sin(radians(37.7438640)) *
sin(radians(37.7438640))))
               AS distance
               FROM foo) AS distances
               WHERE distance < 10
               ORDER BY distance

If I break this down the following returns 1:
SELECT (cos(radians(37.7438640)) * cos(radians(37.7438640)) *
cos(radians(-97.4631299) - radians(-97.4631299)) + sin(radians(37.7438640))
* sin(radians(37.743864000)));

acos(1) would give me 0.

Thoughts?

I don't think this is a bug. It's a well known issue with the
arccos(..) form of the the Haversine formula that it suffers from
large rounding errors when the distance between the points is small.
In this case the intermediate value is a little over 1, due to these
rounding errors, but you can't see that, due to limited precision of
the output format.

Using the arcsin(..) form of the Haversine formula cures that -
http://en.wikipedia.org/wiki/Great-circle_distance

Regards,
Dean

"Vince" <vincecarney@gmail.com> writes:

If I break this down the following returns 1:
SELECT (cos(radians(37.7438640)) * cos(radians(37.7438640)) *
cos(radians(-97.4631299) - radians(-97.4631299)) + sin(radians(37.7438640))
* sin(radians(37.743864000)));

No, it doesn't return 1, it returns 1-plus-a-little-bit, which is hidden
by float8out's desire to not print things like 1.0000000000000002.
Try it after setting extra_float_digits to 2 or so, or do

regression=# SELECT (cos(radians(37.7438640)) * cos(radians(37.7438640)) *
cos(radians(-97.4631299) - radians(-97.4631299)) + sin(radians(37.7438640))
* sin(radians(37.743864000))) - 1;
?column?
----------------------
2.22044604925031e-16
(1 row)

As somebody already remarked, you need to use a version of that formula
that's less prone to roundoff error.

regards, tom lane