3

我正在制作一个函数，该函数使用它们的纬度/经度（以度而不是弧度为单位）和余弦球面定律来计算两点的距离。我遇到的问题是，由于函数中的舍入误差，acos()当两点彼此非常接近时，我获得的结果远非良好。

#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;
}

输出是

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

正如你所看到的，当它计算时acos()它应该给出零，但它确实给出了一些误差，这些误差在乘以地球半径后会被极大地放大。当两点不相等但非常接近时也会发生这种情况。如果有任何用处，我关于纬度和经度的数据最多有 7 位小数。

4

3 回答 3

6

你得到的结果acos和它得到的一样好：问题是计算arg总是会有一个小错误，并返回一个稍微偏离的值。当两点相等或非常接近时，结果小于一，例如 1-10 -16。如果您查看图表，acos(x)您会发现它在 x=1 处几乎是垂直的，这意味着即使是最轻微的误差也会arg对相对误差产生巨大影响。换句话说，该算法在数值上是不稳定的。

您可以使用半正弦公式来获得更好的结果。

于 2013-12-31T12:22:05.490 回答
4

这正是扩展精度double的用途：以比用于参数和结果的精度更高的精度计算中间结果。

我因此更改了您的程序：

#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;
}

通过这些更改，在为 提供扩展精度的编译器上long double，结果符合预期：

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

编辑：

这是一个使用 GCC 的 quadmath 库获得中间结果的版本。它需要最近的 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;
}

我编译并运行：

$ gcc-206231/bin/gcc t.c -lquadmath && LD_LIBRARY_PATH=gcc-206231/lib64 ./a.out
于 2013-12-31T12:25:02.187 回答
0

在这种情况下，通常最好尝试重新表述问题，以避免像 acos( x )这样的不良条件公式x小。在球体上的距离计算情况下，这是一个屡获殊荣的问题，并且在维基百科的Great-circle 页面上提供了更好的公式。这些使用普通的双精度算术为短距离提供高精度。

于 2014-02-11T13:47:03.193 回答