On 07/06/15 03:48, Max wrote:
In the test case the Pd implementation is 8917.74 km off the proper result (2887.26). However I need a precision of about 1m.
Circumference of the Earth in meters: 40,075,000 Accuracy of single precision (24bit): 16,777,216
So your input values will already be inaccurate before you even start messing with rounding errors etc.
So I assume the haversine formula is not implementable in Pd at all? (unless double precision will be there that is)
I guess so.
Or is there a workaround?
Not easy at all, but:
Maybe use two floats at different scales 'a' (close to the true value), 'b' (the residual difference from the true value) to express each coordinate 'c', where:
c = a + b
This would give around 48bits, which should be enough. Actually performing the addition would give just 'a', due to limited precision. You have to work with the unpleasant properties of floating point numbers, like "(a + b) - b" not always being equal to "a".
There are a few blog posts out there about using it on GPU with GLSL ("emulated double", "double-single"), and there is the QD package for "double-double" and "quad-double" in C++ and FORTRAN (maybe with C wrappers), there's a Haskell library called "compensated" which might give some ideas also.
You'll probably also need to do some algebraic manipulations with trigonometric identities etc to get accurate results.
See: https://www.thasler.com/blog/blog/glsl-part2-emu http://crd-legacy.lbl.gov/~dhbailey/mpdist/ http://hackage.haskell.org/package/compensated https://en.wikipedia.org/wiki/List_of_trigonometric_identities