Why is law of cosines more preferable than haversine when calculating distance between two latitude-longitude points?

The problem is indicated by the word "well-conditioned." It's an issue of computer arithmetic, not mathematics.

Here are the basic facts to consider:

1. One radian on the earth spans almost 10^7 meters.

2. The cosine function for arguments x near 0 is approximately equal to 1 - x^2/2.

3. Double-precision floating point has about 15 decimal digits of precision.

Points (2) and (3) imply that when x is around one meter, or 10^-7 radians (point 1), almost all precision is lost: 1 - (10^-7)^2 = 1 - 10^-14 is a calculation in which the first 14 of the 15 significant digits all cancel, leaving just one digit to represent the result. Flipping this around (which is what the inverse cosine, "acos", does) means that computing acos for angles that correspond to meter-length distances cannot be done with any meaningful accuracy. (In certain bad cases the loss of precision gives a value where acos is not even defined, so the code will break down and give no answer, a nonsense answer, or crash the machine.) Similar considerations suggest you should avoid using the inverse cosine if distances less than a few hundred meters are involved, depending on how much precision you're willing to lose.

The role played by acos in the naive law-of-cosines formula is to convert an angle to a distance. That role is played by atan2 in the haversine formula. The tangent of a small angle x is approximately equal to x itself. Consequently the inverse tangent of a number, being approximately that number, is computed essentially with no loss in precision. This is why the haversine formula, although mathematically equivalent to the law of cosines formula, is far superior for small distances (on the order of 1 meter or less).

Here is a comparison of the two formulas using 100 random point-pairs on the globe (using Mathematica's double-precision calculations).

You can see that for distances less than about 0.5 meters, the two formulas diverge. Above 0.5 meters they tend to agree. To show how closely they agree, the next plot shows the ratios of the law of cosines:haversine results for another 100 random point pairs, with their latitudes and longitudes randomly differing by up to 5 meters.

This shows that the law of cosines formula is good to 3-4 decimal places once the distance exceeds 5-10 meters. The number of decimal places of accuracy increases quadratically; thus at 50-100 meters (one order of magnitude) you get 5-6 dp accuracy (two orders of magnitude); at 500-1000 meters you get 7-8 dp, etc.

A historical footnote:

The haversine was a way of avoiding large round-off errors in computations such as

1 - cos(x)

when x is small. In terms of the haversine we have

1 - cos(x) = 2*sin(x/2)^2
           = 2*haversin(x)

and 2*sin(x/2)^2 can be computed accurately even when x is small.

In the old days, the haversine formula had an additional advantage of avoiding an addition (which entailed an antilog lookup, the addition, and a log lookup). A trigonometic formula which entailed only multiplications was said to be in "logarithmic form".

Nowadays, use of the haversine formulas is slightly anachronistic. It might be that the angle x is expressed in terms sin(x) and cos(x) (and x might not be explicitly known). In that case, computing 1 - cos(x) via the haversine formula entails an arctangent (to get the angle x), halving (to get x/2), a sine (to get sin(x/2)), a square (to get sin(x/2)^2) and a final doubling. You are far better off using evaluating

1 - cos(x) = sin(x)^2/(1 + cos(x))

which entails no evaluations of trigonometric functions. (Obviously use the right hand side only if cos(x) > 0; otherwise, it's OK to use 1 - cos(x) directly.)