I really don't think you need double precision to do this.
The following is a very common and widely used formula for calculating distance between two sets of GPS coordinates.
dist = acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon2-lon1));
Unfortunately, this formula can not be relied upon with only 32-bit floats. The challenge is that when you combine small and large scale values in calculations, you loose precision. Triangulation as above relates narrow distances to global coordinates several magnitudes larger and this is the challenge. Whether your required output is a single bit (true/false, left/right) or a high precision number really doesn't matter. You still need the intermediate calculations to be precise in order to determine correct output.
The following is a rewrite that work well with 32-bit floats.
// alternative formulae (haversine)
const float two=2.0;
dist = two*asin(sqrt(square(sin((lat1-lat2)/two)) +
Debugging precision issues is difficult because an algorithm may work well for one set of coordinates, but fail miserably for another set. Realizing this and figuring out how to rewrite the formulas (if at all possible) to avoid rounding/precision issues is anything but trivial.
I'm no expert whatsoever in astronomical calculations, but I do respect that our planet earth is pretty small in comparison to our solar system at large. Even with double precision it is a formidable challenge to come up with formulas that avoid computational pitfalls due to rounding and loss of precision. GPS calculations in this context is like first grade math in comparison so I have no reason to doubt that double precision may indeed be required.