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