[PD] haversine formula in Pd

Claude Heiland-Allen claude at mathr.co.uk
Sun Jun 7 05:27:24 CEST 2015


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>


Claude
-- 
http://mathr.co.uk




More information about the Pd-list mailing list