[linux-users:78154] optimization and accuracy with gcc

[Ken Ebisawa <ebisawa _at_ xxxxxxxxxxxxxxxxxxxx>]

海老沢@グリーンベルトです。

# 以前に関連した話題で投稿したことがあるのですが。。

大きなプログラムをSolarisからLinuxに移植していて、
計算精度に由来するエラーが起きていることがわかりました。
問題になっている部分を取り出して、簡略したプログラムを
最後につけます。やっていることは単純で、3次元単位ベクトル
を極座標(alpha, delta)で与え、その長さを計算するということです。
当然長さは1でなくてはいけません。最適化パラメーターを変えてこの
プログラムをコンパイル、実行すると、結果が変わってきます。

gcc -O0 のとき：
  alpha = 135.000000, delta = -40.000000
  vector=(-0.541675, 0.541675, -0.642788)
  Norm = 1.000000
  This is OK.
  (Norm - 1.0)*1e20 = 0.000000 

  alpha = 135.000000, delta = 50.000000
  vector=(-0.454519, 0.454519, 0.766044)
  Norm = 1.000000
  This is OK.
  (Norm - 1.0)*1e20 = 0.000000 

gcc -O1のとき：
  alpha = 135.000000, delta = -40.000000
  vector=(-0.541675, 0.541675, -0.642788)
  Norm = 1.000000
  This is OK.
  (Norm - 1.0)*1e20 = 0.000000 

  alpha = 135.000000, delta = 50.000000
  vector=(-0.454519, 0.454519, 0.766044)
  Norm = 1.000000
  This is OK.
  (Norm - 1.0)*1e20 = -11102.230246 

gcc -O2または-O3のとき：
  alpha = 135.000000, delta = -40.000000
  vector=(-0.541675, 0.541675, -0.642788)
  Norm = 1.000000
  Your should not see this !
  (Norm - 1.0)*1e20 = 2255.140519 

  alpha = 135.000000, delta = 50.000000
  vector=(-0.454519, 0.454519, 0.766044)
  Norm = 1.000000
  This is OK.
  (Norm - 1.0)*1e20 = -3984.442984 

このように最適化によって計算精度が変わってくるのは、
そういうもんなんでしょうかねえ。問題となっているのは、
このベクトルの長さをacosにいれていて、引数が1.0より
大きいときに落ちる、ということです。上の例では最適化を
はずすしたら正しく動作していますが(ベクトルの長さが1)、
そうならない場合もあります。最適化をはずすだけでは
問題は解決しない、と。

こういう場合に賢い対処の方法というのは、どうしたもの
だろうか、と悩んでいます。開発はSolarisの上でやっていて、
"acos(単位ベクトルどうしの内積)"、という記述はいたる
ところに出てくるので、ソースはできるだけ変更したくないのです。

小数約16桁より下のほうで1.0からずれていて、これはおそらく
doubleで保障している精度以下なので、1.0と同一視してくれれば
いいのですが。そういうgccのオプションを探したのですが、ちょっと
見つかりませんでした。

何かコメントを頂ければ幸いです。環境はRedhat6.2, カーネル
2.2.16-3, egcs-2.91.66, CPUはAMD 1GHzです。

よろしくお願いします。

#include <math.h>
main(){
double alpha[2] = {135.0,135.0};
double delta[2] = {-40.0,50.0};
double vector[3][2];
double DD2R=0.017453292519943295769236907684886127134428718885417;
int i;
double norm;

for(i=0;i<=1;i++){
  printf("alpha = %lf, delta = %lf\n",alpha[i],delta[i]);
  vector[0][i] = cos(DD2R*alpha[i])*cos(DD2R*delta[i]);
  vector[1][i] = sin(DD2R*alpha[i])*cos(DD2R*delta[i]);
  vector[2][i] = sin(DD2R*delta[i]);
  printf("vector=(%lf, %lf, %lf)\n",vector[0][i],vector[1][i],vector[2][i]);
  norm = vector[0][i]*vector[0][i]+vector[1][i]*vector[1][i]+vector[2][i]*vector[2][i];
  printf("Norm = %lf\n",norm);
  if(norm > 1.0L){
    printf("Your should not see this !\n");
    printf("(Norm - 1.0)*1e20 = %lf \n",(norm-1.0)*1e20);	   
  }
  else
    {
      printf("This is OK.\n");
      printf("(Norm - 1.0)*1e20 = %lf \n",(norm-1.0)*1e20);	   
    }
  printf("\n");
}
return;
}