acos() bad results with two close points in C

Question

I am making a function that calculates the distance of two points using their latitude/longitude (in degrees not radians) and the spherical law of cosines. The problem I have is that due to rounding errors in the function acos() the results I am obtaining are far from good when the two points are very close to each other.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct{
 double lat;
 double lon;
} point;

#define DEG_TO_RAD 0.017453292519943295769236907684886
#define EARTH_RADIUS_IN_METERS 6372797.560856

double distance(point a, point b) {
  double arg=sin(a.lat * DEG_TO_RAD) * sin(b.lat * DEG_TO_RAD) + cos(a.lat * DEG_TO_RAD) * cos(b.lat * DEG_TO_RAD) * cos((a.lon-b.lon) * DEG_TO_RAD);
  if(arg>1) arg=1;
  else if (arg<-1) arg=-1;
  printf("arg=%.12f acos(arg)=%.12f\n",arg, acos(arg));    //to see the problem
  return acos(arg) * EARTH_RADIUS_IN_METERS;
}

int main(){
  point p1,p2;

  p1.lat=63.0;
  p1.lon=27.0;
  p2.lat=p1.lat;
  p2.lon=p1.lon;

  printf("dist=%.8f\n",distance(p1,p2));

  return 0;
}

The output is

arg=1.000000000000 acos(arg)=0.000000014901
dist=0.09496208

as you can see, when it computes the acos() it should give zero but it does give some errors that get enormously magnified after multiplying by the earth radius. It also happens when the two points are not equal but very close. If it is of any use my data about latitude and longitude has up to 7 decimal digits.

Answer 1

The result you get from acos is as good as it gets: the problem is that the calculation of arg will always have a small error and return a value that is slightly off. When the two points are equal or very close the result is less than one, for example 1-10^-16. If you look at the graph of acos(x) you'll see that it's nearly vertical at x=1, which means that even the slightest error in arg has a huge impact on the relative error. In other words, the algorithm is numerically unstable.

You can use the haversine formula to get better results.

Answer 2

This is exactly what extended precision is for: to compute intermediate results at a higher precision than the double precision used for arguments and results.

I changed your program thus:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>

typedef struct{
 double lat;
 double lon;
} point;

#define DEG_TO_RAD 0.017453292519943295769236907684886L
#define EARTH_RADIUS_IN_METERS 6372797.560856L

double distance(point a, point b) {
  long double arg=sinl(a.lat * DEG_TO_RAD) * sinl(b.lat * DEG_TO_RAD) + cosl(a.lat * DEG_TO_RAD) * cosl(b.lat * DEG_TO_RAD) * cosl((a.lon-b.lon) * DEG_TO_RAD);
  if(arg>1) arg=1;
  else if (arg<-1) arg=-1;
  printf("arg=%.20Le acos(arg)=%.20Le\n",arg, acosl(arg));
  return acosl(arg) * EARTH_RADIUS_IN_METERS;
}

int main(){
  point p1,p2;

  p1.lat=63.0;
  p1.lon=27.0;
  p2.lat=p1.lat;
  p2.lon=p1.lon;

  printf("precision of long double:%Le\n", LDBL_EPSILON);

  printf("dist=%.8f\n",distance(p1,p2));

  return 0;
}

With these changes, on a compiler that offers extended precision for long double, the result is as expected:

precision of long double:1.084202e-19
arg=1.00000000000000000000e+00 acos(arg)=0.00000000000000000000e+00
dist=0.00000000

EDIT:

Here is a version that uses GCC's quadmath library for intermediate results. It requires a recent GCC.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <quadmath.h>

typedef struct{
 double lat;
 double lon;
} point;

#define DEG_TO_RAD (M_PIq / 180)
#define EARTH_RADIUS_IN_METERS ((__float128)6372797.560856L)

double distance(point a, point b) {
  __float128 arg=sinq(a.lat * DEG_TO_RAD) * sinq(b.lat * DEG_TO_RAD) + cosq(a.lat * DEG_TO_RAD) * cosq(b.lat * DEG_TO_RAD) * cosq((a.lon-b.lon) * DEG_TO_RAD);
  if(arg>1) arg=1;
  else if (arg<-1) arg=-1;
  printf("arg=%.20Le acos(arg)=%.20Le\n",(long double)arg, (long double)acosq(arg));
  return acosq(arg) * EARTH_RADIUS_IN_METERS;
}

int main(){
  point p1,p2;

  p1.lat=63.0;
  p1.lon=27.0;
  p2.lat=p1.lat;
  p2.lon=p1.lon;

  printf("dist=%.8f\n",distance(p1,p2));

  return 0;
}

I compiled and ran with:

$ gcc-206231/bin/gcc t.c -lquadmath && LD_LIBRARY_PATH=gcc-206231/lib64 ./a.out 

Answer 3

In cases like this it's usually best to try to reformulate the problem to avoid badly conditioned formulas like acos(x) with x small. This is a well-beaten problem in the case of distance computations on a sphere and better formulations are provided on the Great-circle page on Wikipedia. These give high accuracy for short distances with ordinary double precision arithmetic.