-1

I have made a sample project using dot42 and have taken this C# code and placed it inside a library:

var gcd = GetGreatCircleDistanceKm(52.0, -1.90,  21.0, 39.0); // returns NaN should be ~4915

public double GetGreatCircleDistanceKm(double startLat, double startLong, double endLat, double endLong)
{
    var earthRadius = Constants.EarthRadiusKm;
    var φ1 = DegreesToRadians(startLat);
    var φ2 = DegreesToRadians(endLat);
    var Δφ = DegreesToRadians(endLat - startLat);
    var Δλ = DegreesToRadians(endLong - startLong);
    var a = Math.Sin(Δφ / 2) * Math.Sin(Δφ / 2) + Math.Cos(φ1) * Math.Cos(φ2) * Math.Sin(Δλ / 2) * Math.Sin(Δλ / 2);
    var c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
    return earthRadius * c;
}

public static double DegreesToRadians(double angle)
{
    return Math.PI * angle / 180.0;
}

problem is however that c always returns NaN and thus the function returns NaN. This same method works fine in a normal .net project.

I have tried many versions of the same haversine formula thinking I may have mistyped something, but no, all return NaN.

Source of formula : http://www.movable-type.co.uk/scripts/latlong.html

The solution was (strangely) to do :

var a1 = Δφ/2;
var a = Math.Sin(a1) * Math.Sin(a1) + Math.Cos(φ1) * Math.Cos(φ2) * Math.Sin(Δλ / 2) * Math.Sin(Δλ / 2);

No idea how it works, but it works. Please if someone could shed light.

Full fixed version:

/// <summary>
/// Gets the shortest possible distance between two points on earth
/// </summary>
/// <param name="startLat"></param>
/// <param name="startLong"></param>
/// <param name="endLat"></param>
/// <param name="endLong"></param>
/// <returns></returns>
public double GetGreatCircleDistanceKm(double startLat, double startLong, double endLat, double endLong)
{
    var startLatRad = DegreesToRadians(startLat);
    var endLatRad = DegreesToRadians(endLat);
    var latDiffRad = DegreesToRadians(endLat - startLat);
    var longDiffRad = DegreesToRadians(endLong - startLong);
    var halfLatDiff = latDiffRad/2;
    var a = Math.Sin(halfLatDiff) * Math.Sin(halfLatDiff) + Math.Cos(startLatRad) * Math.Cos(endLatRad) * Math.Sin(longDiffRad / 2) * Math.Sin(longDiffRad / 2);
    var c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
    return Constants.EarthRadiusKm * c;
}
sprocket12
  • 5,368
  • 18
  • 64
  • 133
  • Have you printed out what the values are during the computation? Or used a debugger to inspect the values? This should give you a clue to where the computation went off the rails. – Max Dec 24 '14 at 23:12
  • This reads a lot more like C# than it does Java, so I'm going to remove that tag. If it's incorrect, go ahead and add it in again. – Makoto Dec 24 '14 at 23:49
  • What are the full-precision values of `a` and `1-a` before the assignment statement, and what are the values of `Math.Sqrt(a)` and `Math.Sqrt(1 - a)`? (Possible values include NaN & INFs.) – philipxy Jan 06 '15 at 01:26

1 Answers1

0

The IEEE floating point standard is poorly supported by most languages and compilers. This is because there are many modes and flags for each operator call and variations on INF, NAN and even 0 and the languages do not make it easy for the programmer to set, check or respond to these.

Eg: If here a Math.Sin division by 2 underflows (is too small to be represented) or causes precision to be lost (because results nearest 0 are not accurate to the full number of bits) then this can set flags or return values (depending also on modes); if this is in the middle of an expression like your original assignment to variable a then the language/compiler can only handle one case, which might not be what you need. Whereas the language/compiler rule for assigning to an intermediate variable like halfLatDiff, depending on the result value, flags and modes of the right hand expression and including possibly some changes to those from the assignment, might just happen to give a different answer.

From Lecture Notes on the Status of IEEE Standard 754 for Binary Floating-Point Arithmetic by its main architect/instigator (circa 1997 between the original 1984 and current 2008 versions:

Those features lack support in programming languages and compilers, so those features are mishandled and/or practically unusable, so those features are little known and less in demand, and so those features lack support in programming languages and compilers.

Maybe of relevance here is:

Fused MACs generate anomalies when used to evaluate ab ± cd

(If you read the assembly code for your programs (which may be sufficiently intelligible even if you only know the basics of assembly code) you might be able to see how the expressions have been interpreted. But then many compilers do not properly take into account caching and pipelining of their hardware. And neither may your debugger. So you still might not be able to reconstruct the reasons for your observed outputs.)

philipxy
  • 14,867
  • 6
  • 39
  • 83