Mathematik ist einfach, sie auf Rechner zu bringen nicht…

Mathe ist ja eigentlich ganz einfach. Die nachfolgende Routine berechnet die kürzeste Distanz zwischen zwei Punkten auf der Erde unter der Berücksichtigung, dass die Erde eine perfekte Kugel ist.
Benutzt wird diese Routine um die Luftlinie zwischen zwei geokodierten Adressen zu ermitteln.

double DistanceBetweenCoordinates(
          double dLatitude1, double dLongitude1, 
          double dLatitude2, double dLongitude2)
{
  // Convert to radient
  dLatitude1 *= M_PI / 180;
  dLongitude1 *= M_PI / 180;
  dLatitude2 *= M_PI / 180;
  dLongitude2 *= M_PI / 180;

  // Get distance
  double dDistance;
  dDistance = sin(dLatitude1) * sin(dLatitude2) + cos(dLatitude1) * 
        cos(dLatitude2) * cos(dLongitude2-dLongitude1);
  // Ohne diese Zeile haben wir ein Problem... dDistance ist evtl. >1
  // dDistance = min(dDistance,1.0);
  dDistance = acos(dDistance) * EARTH_RADIUS;
  return dDistance;
}

Eigentlich ganz einfach. Aber dieser Code hat ein massives Problem, denn dDistance kann vor dem acos durch Rundungsfehler größer als 1 werden.

Für die folgenden Werte passiert etwas Übles:

dLatitude1 0.86208095750908087 
dLongitude1 0.15333764167657613 
dLatitude2 0.86208095750908087 
dLongitude2 0.15333764167657613

Eigentlich müsste glatt 1 herauskommen, denn die Koordinaten sind gleich, aber dDistance wird bei diesen Werten 1.0000000000000002. Der Sinus bzw. Cosinus ergibt eben doch leicht abweichende Werte. Der folgende acos scheitert dann aber und das Resultat wird NaN (Not a Number). Erlaubt sind natürlich nur Werte <=1.0 und >=-1.0. Die Folge ist, dass die Distanz nicht 0.0 ist.

Also muss hier eine kleine Sicherung eingebaut werden, die Aufgrund der Rundungs-Genauigkeitsfehler, das überschreiten von 1 verhindert.
Eigentlich ist Mathe ganz einfach, aber auf einem Rechner ist das alles manchmal ganz anders.

PS: Ich weiß schon warum ich noch nie Fließkommaarithmetik auf PCs mochte.