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.

You might want to read [What Every Computer Scientist Should Know About Floating-Point Arithmetic](http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html). In short, it's not possible to get exact representations of some numbers, so you will always have rounding problems. You can alleviate some of them by using a larger precision type, like `long double`. — Some programmer dude, Dec 31 '13 at 12:02
This isn't a problem with `acos`; the problem is that `arg` is not exactly `1.0`. — Oliver Charlesworth, Dec 31 '13 at 12:07
Side note: `` contains definitions for various mathematical constants, so you can define `#define DEG_TO_RAD M_PI/180.0` instead. — Martin R, Dec 31 '13 at 12:09
A format that displays enough digits to distinguish all double-precision numbers is `%.16e`. When you use `%.12f`, you may get the impression that `arg` is zero when it isn't. — Pascal Cuoq, Dec 31 '13 at 12:10
When the points are close together the curvature of the earth can be ignored and you can apply Pythagoras. You just need to apply a scaling factor to longitude, based on latitude. — Hot Licks, Dec 31 '13 at 13:26

Joni · Accepted Answer · 2013-12-31T12:24:02.190

6

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.

edited Dec 31 '13 at 12:24

answered Dec 31 '13 at 12:22

Joni

108,737
14
143
193

This one seems to be obtaining better precision even using only double and not long double... – Jesús Ros Dec 31 '13 at 12:52
2

@gunbl4d3 Yes, knowledge of math usually beats knowledge of floating-point programming. – Pascal Cuoq Dec 31 '13 at 13:04
I'll take this solution. Only one has to make sure that the values within the square root are bounded between 0 and 1 and do not cross those boundaries due to floating point errors. – Jesús Ros Dec 31 '13 at 16:12

Pascal Cuoq · Answer 2 · 2013-12-31T13:03:43.913

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

In windows I do not seem to have this `long double`... I am using MinGW. — Jesús Ros, Dec 31 '13 at 12:39
@gunbl4d3 I was going to say something about the definition of `long double` as `double` on Windows. GCC has a software implementation of quad-precision. I will try to use it and update my answer with an example if I manage to: http://gcc.gnu.org/onlinedocs/libquadmath.pdf — Pascal Cuoq, Dec 31 '13 at 12:43
@gunbl4d3 Done. The quadmath library is only part of recent enough GCC versions. — Pascal Cuoq, Dec 31 '13 at 12:59
Thanks for the effort but I could not implement it in my machine. I will take it into account when I am making code for Unix systems, where it worked properly or if I get a newer gcc version on windows. — Jesús Ros, Dec 31 '13 at 16:15

score 0 · Answer 3 · answered Feb 11 '14 at 13:47

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.

acos() bad results with two close points in C

3 Answers3

EDIT: