ちら裏メモなのでよく考えてないです

先日、本屋で Ruby 本を立ち読みしていたのですが、実数どうしの合同の比較で、画面の出力と実行結果を一致させるために、数値を文字列に変換して文字列どうしで比較するという項目があって、面白く奇妙な発想だと思い、メモっておくことにしました。

FORTRAN では、初心者に x == y の代わりに ABS(x - y) < eps のように書けと教えるわけですが、この時 eps をどのくらいにすればいいかは教えないので、別の地獄に落ちたりします。

Fortran90 では EPSILON(x), TINY(x), SPACING(x) などという精度にかかわる組み込み関数が用意された訳ですが、これも色々問題があります。これらの関数は正規数しか返しませんが、非正規数が有効だと二つの実数間で差をとった時に TINY(x) よりも小さい数が出てくる可能性があります。（また NEAREST(x, sign) 関数では SPACING(x) よりも小さい間隔で再隣接の数が得られることがあります。）

ここで TINY(x) は引数の型の実数のとりうる正規数の中の最小の数を返す関数です。SPACING(x) は正規数を返す関数なので、取りうる値は TINY(x) 以上の数になります。EPSILON(x) は 1.0 近傍での仮数部の最下位ビットの表す数を返します。これは、有効桁の目安になります。

正規数であっても仮数部の下の方のビットは TINY より小さい量を表し得ます。（数直線上の浮動小数の図があればわかり易いのですが・・）よって EPSILON(x) や TINY(x) や SPACING(x) が、吟味せずに使える微小量を教えてくれるわけでもないことになります。

実行結果

倍精度変数に単精度の数を代入すると、仮数部52ビット中に24ビット分のみが代入されるので、十進表現で見たとき倍精度変数と単精度の数とは異なる値を表示することがあります。比較する時は仮数部のビット列が同じなのでイコールが成立します。

Fortranでは、定数をプログラムに書く時、精度をきっちり合わせて書けとうるさく教えられますが、確かに重要なことです。たった１個、倍精度にし忘れた単精度定数があったばかりに、精度が足りなくてとんでもない結果が返ってきたりします。

Fortran90で導入された再隣接数を返す関数 NEAREST(x, sign) で隣の数と比較すると、１０進表現で見た時同じ数に見えても（自由フォーマットで）、仮数部の最小ビットが変化しているのでイコールが成り立たちません。

ソースコード

PROGRAM test
  IMPLICIT NONE
  REAL(4) :: x, y
  REAL(8) :: d
  CHARACTER(20) :: cx, cy
  CHARACTER(32) :: tmp1
  CHARACTER(64) :: tmp2
  
  x = 0.1
  d = x
  
  PRINT *, 'x = 0.1, d = x ' 
  IF (x == d) THEN 
   PRINT *, 'x == d:', x, '-', d, '=', x - d
  ELSE
   PRINT *, 'x /= d:', x, '-', d, '=', x - d
  END IF 
  WRITE(tmp1, '(b32.32)') x
  PRINT *, tmp1(1:1), ':   ', tmp1(2:9), ':', tmp1(10:)
  WRITE(tmp2, '(b64.64)') d
  PRINT *, tmp2(1:1), ':', tmp2(2:12), ':', tmp2(13:)

  PRINT *
  PRINT *, 'x == y ; y = NEAREST(x, +1.0) '
  y = NEAREST(x, +1.0)
  IF (x == y) THEN 
   PRINT *, 'x == y:', x, '-', y, '=', x - y
  ELSE
   PRINT *, 'x /= y:', x, '-', y, '=', x - y
  END IF 

  WRITE(tmp1, '(b32.32)') x
  PRINT *, tmp1(1:1), ':', tmp1(2:9), ':', tmp1(10:)
  WRITE(tmp1, '(b32.32)') y
  PRINT *, tmp1(1:1), ':', tmp1(2:9), ':', tmp1(10:)

  PRINT *
  PRINT *, 'ABS(x - y) <= SPACING(x)'
  IF (ABS(x - y) <= SPACING(x)) THEN 
   PRINT *, 'x == y:', x, '-', y, '=', x - y
  ELSE
   PRINT *, 'x /= y:', x, '-', y, '=', x - y
  END IF 
  
  WRITE(cx, *) x
  WRITE(cy, *) y
  PRINT *
  PRINT *, 'ON SCREEN CHARACTERS FOR X AND Y:', TRIM(cx), '      ', TRIM(cy)
  IF (cx == cy) THEN 
   PRINT *, 'x == y:', x, '-', y, '=', x - y
  ELSE
   PRINT *, 'x /= y:', x, '-', y, '=', x - y
  END IF 

  STOP
END PROGRAM test