3

緯度/経度 (ラジアンではなく度数) とコサインの球面法則を使用して 2 点の距離を計算する関数を作成しています。私が抱えている問題は、関数の丸め誤差のためacos()に、2 つの点が互いに非常に近い場合に得られる結果があまり良くないことです。

#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()とゼロになるはずですが、地球の半径を掛けた後に非常に拡大されるいくつかのエラーが発生します。また、2 つの点が等しくないが非常に近い場合にも発生します。緯度と経度に関する私のデータは、10 進数で 7 桁までです。

4

3 に答える 3

6

得られる結果はacos最高です。問題は、 の計算でarg常に小さなエラーが発生し、わずかに異なる値が返されることです。2 つの点が等しいか非常に近い場合、結果は 1 未満になります (例: 1-10 -16 ) 。のグラフをacos(x)見ると、x=1 でほぼ垂直になっていることがわかります。これは、わずかな誤差でもarg相対誤差に大きな影響を与えることを意味します。つまり、アルゴリズムは数値的に不安定です。

hasersine 式を使用して、より良い結果を得ることができます。

于 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が小さいものにしないようにするのが最善です。これは、球体上の距離計算の場合にはよく知られている問題であり、ウィキペディアの大圏ページでより適切な定式化が提供されています。これらは、通常の倍精度演算で短い距離に対して高い精度を提供します。

于 2014-02-11T13:47:03.193 に答える