Re: 浮動小数点を比較するには

comp sci

浮動小数点を比較するには[ruby-list:43103] からはじまるスレはスレ立て人が出てこなくなったし Ruby の話じゃないので隠れキャラのようにこっそりこっちに書く。

たいていはつぎのことが分かってれば十分なんじゃないでしょうか。

  • 浮動小数点数は2進法での有効数字の桁数が一定
  • 10進法で有限小数でも2進法では循環小数になることがあるので丸められる
  • 精度の考え方は理系の大学で実験の前にならう有効数字の扱いと大体一緒
  • たいていの場合全部の桁を気にすることはないし一致性を直接比較するのは間違い
  • IEEE754 に由来する浮動小数点の規格には丸め方などにオプションがあるので同じ入力同じプログラムでもプラットフォームや設定によって最下位に近い桁が一致しないことがある

しかし現実的な数値計算は、出来る限り既存のライブラリを使うべきです。難しいことが知られている問題、あるいはしばしば必要になる解法は Ruby での実装がなくても C か FORTRAN での実装はたいていあるので、拡張ライブラリを自分で作ればいい。多くの場合は、Ruby でラップするほうがずっと簡単で、問題固有の困難さの法がはるかにやっかいです。

ちなみに Goldberg の解説はほとんどの人には荷が重いんじゃないかと思います。伊理 正夫, 藤野 和建, 『数値計算の常識』, 共立あたりが入手も容易だし手ごろなんではないかと。僕の知識は古いのでもっといい本が出てるかもしれません。

もっとも、Goldberg の解説は確かにすばらしいので、数値計算のライブラリを作る人、メンテナンスする人、ばしばし使う人、それから本や雑誌記事を書く人、あたりには思い切り薦めます。

それと細かい議論を避けるためにどのメッセージも引用しなかったらみごとに成功したので、Ruby界のスルー力にちょっと感心したりしたのだけど、最初に書いたときに忘れたのは、ulp とか machine epsilon を等価性の検査に使うのはほとんどの場合にミスリードだということ。それがうまくいくのはリテラルが数値化される際の丸めとそこからせいぜい数度の四則演算のあとくらいだろう。Float#inspect を次のに差し替えて irb でいくつか計算してみればわかること。

# ieee754inspect.rb - overwrites Float#inspect.
# clipped from ieee754hack.rb:
#   http://www.notwork.org/~gotoken/ruby/p/as-is/ieee754hack.rb

class Float
  def human_readable
    if zero?
      format("%c0.0", sign < 0 ? ?- : ?+ )
    elsif normal?
      format("%+d * 0b1.%s * 2**%d",
             sign, mantissa_binary_string, exponent-1023)
    elsif denormal?
      format("%+d * 0b0.%s * 2**-1022", sign, mantissa_binary_string)
    elsif nan?
      format("NaN <%s, %s, %s>", *binary_string)
    else # infinity?
      format("%cInfinity", sign < 0 ? ?- : ?+)
    end
  end

  alias :inspect :human_readable

  def sign()
    if nan?
      [self].pack("G").unpack("B*").first[0] == ?1 ? -1 : 1
    else
      1/self > 0 ? 1 : -1
    end
  end

  FORM__ = "%b%011b%052b"

  def self::compile(s, e, m)
    [format(FORM__, (s<0 ? 1 : 0), e, m)].pack("B*").unpack("G").first
  end

  MINIMUM = Float::compile(+1, 1, 0)

  def normal?()
    finite? and abs >= MINIMUM
  end

  def denormal?()
    finite? and abs < MINIMUM and not zero?
  end

  def exponent
    exponent_binary_string.to_i(2)
  end

  def mantissa_binary_string()
    [self].pack("G").unpack("B*").first[12,52]
  end

  def exponent_binary_string
    [self].pack("G").unpack("B*").first[1,11]
  end

  def binary_string()
    x, = [self].pack("G").unpack("B*")
    [x[0,1], x[1,11], x[12,52]]
  end
end

使用例

require "rational"
    ==>true
def f(x) 4*x*(1-x) end
    ==>nil
x = Rational(1, 10) ; 5.times{ x = f(x) } ; a = x.to_f
    ==>+1 * 0b1.0010101110111100001111011010100010001010010010000010 * 2**-1
x = Rational(1, 10).to_f ; 5.times{ x = f(x) } ; b = x
    ==>+1 * 0b1.0010101110111100001111011010100010001010010001111101 * 2**-1

内容については、Cマガの2003年5月号に書いた記事とか、日経ソフトウェアの2001年11月号に書いた記事とか。これらの中で、numerical recipes を薦めてるのはミスリードだったかも。失敗した。

文字列に変換して調べるというのは sprintf でフォーマット指定しているならまっとうと思うけど高価すぎるんじゃないでしょうか。Ruby だったら気にならないのかな。引き算して絶対値をとって基準と比較すればよいと思います。フォーマット指定する場合は、絶対誤差か相対誤差かどっちか決めるべきで、%g ではなく %f か %e を使って精度まで指定するのが本来やりたいことのはず。精度を指定しない場合というのは、なんでもいいから似てるかどうか調べるだけ。精度を %.15e とかにしちゃうのは ulp だけ認めるということなので普通は厳しすぎるしそういうことがしたいのではないでしょう。

『Rubyレシピブック 268の技』は入手し損ねたので持ってないんだけど、どういうふうに書かれてるんだろうなあ。とかいうのはこわいほうにリンクしとくと気づいてもらえるのかな。

しらなかったけど、Rubyレシピブック 第2版 268の技がでるっぽい。

Rubyレシピブック 268の技

Rubyレシピブック 268の技