浮動小数点数の2進10進変換のこと

浮動小数点数のわかりにくさの一因

浮動小数点数は難しいし何だか怪しげな感じがする。たぶん主な理由は

だと思う。でもこれらとは別に、

という理由があるような気がする。
浮動小数点数は2進小数×2のベキの形で表現されているけど、浮動小数点数を出力したり入力したりする場合は10進小数×10のベキの形で表記される。つまり

  • 浮動小数点数演算のシステム自体(コンピュータ側の世界)は2進数だけで完結している。
  • 逆に、浮動小数点数の表示・入力(人間側の世界)では通常、10進数だけで完結している。
  • この二つの世界を行き来する場合に2進10進変換がおこなわれる。

ということになる。さらにこれは、整数での2進10進変換とは全く違っている。整数の場合、2進と10進の間で変換しても表現が違うだけで全く同じ数を表している。
浮動小数点数では、10進→2進変換は多くの数で正確にはおこなえない。例えば1/5は10進数では0.2と正確に書けるけど、2進数では1.100110011……と循環小数になってしまう。
逆向きの2進→10進変換は一応正確におこなうことは可能だけど、浮動小数点数の有効桁数を考慮すると「正確な」10進数にしてもしょうがない。例えば2-100は、2進小数として正確に表すことができる。この数を正確な10進数で表すと

7.888609052210118054117285652827862296732064351090230047702789306640625e-31

となる。でも倍精度浮動小数点数の精度が10進数表記で15桁ちょっとぐらいなので、こんなに書き出してもわりとしょうがない。結局、浮動小数点数では10進→2進変換の場合でも2進→10進変換の場合でも誤差が混入することになる。
浮動小数点数は10進数表記を通してしか値を見たり与えたりしないので結局、浮動小数点数を扱う場合は常に誤差を通して見ている、ということになる。

このように2進10進変換というのは浮動小数点数を扱う場合に考慮すべき事柄のひとつのように思うのだけど、実際には浮動小数点数の説明・解説で、2進10進変換についての説明はあまりおこなわれないような気がする。
例えばDavid Goldberg「What Every Computer Scientist Should Know About Floating-Point Arithmetic」の中で2進10進変換に触れているのは、

あたりしかない。しかも書かれているのは、単精度の数を10進数で出力して再び読み戻したときに同じ数に戻るようにするには9桁出力すれば良い、倍精度の場合は17桁出力すれば良い、ということぐらい。
具体的なアルゴリズムの説明までは必要ないにしても、(多くの処理系では)10進小数を入力したとき必ずしもそれに一番近い浮動小数点数に変換するとは限らないことなんかは説明があっても良いと思う。
2進10進変換のやり方としては、R5RS 6.2.6節、R6RS 11.7.4.4節にあるように、

  1. 外部表現文字列に変換してから浮動小数点数に戻したとき同じ数に戻るようにする。
  2. 1.を満たす最小桁数(最小文字数)で出力する。

が簡潔でわかりやすいと思う。必ずしもこれが最適なやり方とは限らないけど。あと問題はこのやり方(とこれを満たすアルゴリズム)がLisp Scheme界隈以外ではあまり普及していないこと。

浮動小数点数の2進指数表記

2進10進変換によるわかりにくさを避けるひとつのやり方として、そもそも2進10進変換しないというやり方が考えられる。つまり浮動小数点数をそのまま2進表現で出力する、というやり方。通常の10進小数+指数で表記する代わりに2進小数+指数で表記する。例えば(1.1001001)2×2-34を「1.1001001e-34]と書く。
倍精度浮動小数点のいくつかの値は次のように表記される(精度がわかるように、小数部末尾の0は省略しないことにした)。

2進指数表記 10進指数表記
最大の正規化数 1.1111111111111111111111111111111111111111111111111111e1023 1.7976931348623157e308
最小の正規化数 1.0000000000000000000000000000000000000000000000000000e-1022 2.2250738585072014e-308
最大の非正規化数 0.1111111111111111111111111111111111111111111111111111e-1022 2.225073858507201e-308
最小の非正規化数 0.0000000000000000000000000000000000000000000000000001e-1022 5.0e-324

2進指数表記が役に立つような場合も多少はあるかもしれない。

2進指数表記の書き出し・読み込みをおこなうプログラム

Ruby版とGauche版を作ってみた。ただし、すごく長い表記(1000字を越えるぐらい)を与えると、正しくない結果を返すバグあり(理由は浮動小数点演算の途中でオーバーフローが起こるため)。
使ってみるとこんな感じ。

irb(main):002:0> (3.14).to_s2
=> "1.1001000111101011100001010001111010111000010100011111e1"
irb(main):003:0> Float::EPSILON
=> 2.22044604925031e-16
irb(main):004:0> Float::EPSILON.to_s2
=> "1.0000000000000000000000000000000000000000000000000000e-52"
irb(main):005:0> Float.from_s2("-101100111000.11110000e-234")
=> -1.04065600054107e-67

Ruby

class Float
  def to_s2
    sign     = [self].pack("G").unpack("B*")[0][0,1]
    exponent = [self].pack("G").unpack("B*")[0][1,11]
    mantissa = [self].pack("G").unpack("B*")[0][12,52]
    sign_string = if (sign == "0") then "" else "-" end
    exponent_int = Integer("0b" + exponent)
    mantissa_int = Integer("0b" + mantissa)
    if exponent_int == 2047    # 無限or非数?
      self.to_s
    elsif exponent_int == 0
      if mantissa_int == 0     # ゼロ?
        self.to_s
      else                     # 非正規化数
        format("%s0.%se-1022", sign_string, mantissa)
      end
    else                       # 正規化数
        format("%s1.%se%d", sign_string, mantissa, exponent_int - 1023)
    end
  end

  def self.from_s2(x)
    re_float =
      /^([-+]?)(?:(?:([01]+)(?:\.([01]*))?)|\.([01]+))(?:[eE]([-+]?[0-9]+))?$/
    if re_float =~ x
      sign = if $1 == "-" then -1.0 else 1.0 end
      exponent = if $5 then Integer($5) else 0 end
      if $4
        convert_to_float(sign, "", $4, exponent)
      else
        convert_to_float(sign, $2, $3, exponent)
      end
    else
      0.0
    end
  end

  def self.convert_to_float(sign, i_part, f_part, exponent)
    m        = i_part.sub(/^0*/, '') + f_part
    length_m = m.length - 1
    mantissa = sign * Integer("0b" + m) * (2.0 ** ( - length_m))

    mantissa * (2.0 ** (exponent + length_m - f_part.length))
  end
  private_class_method :convert_to_float
end

Gauche

(use srfi-13) ; string-drop string-trim

(define (float->string2 x)
  (let* ((m&e&s    (decode-float x))
         (sign     (cond ((= 1 (vector-ref m&e&s 2)) "")
                         (else "-")))
         (mantissa (vector-ref m&e&s 0))
         (exponent (+ (vector-ref m&e&s 1) 52))
         )
    (case (type-of-float mantissa)
      ((normalized-float)
         (let1 fractional-part (string-drop (number->string mantissa 2) 1)
           (format #f "~a1.~ae~a" sign fractional-part exponent)))
      ((denormalized-float)
         (format #f "~a0.~52,'0be~a" sign mantissa exponent))
      (else x))))

(define (type-of-float m)
  (cond ((not (number? m)) 'inf-or-nan)
        ((zero? m)         'zero)
        ((< m (expt 2 52)) 'denormalized-float)
        (else              'normalized-float)))

(define (string2->float x)
  (define rx-float
    #/^([-+]?)(?:(?:([01]+)(?:\.([01]*))?)|\.([01]+))(?:[eE]([-+]?[0-9]+))?$/)

  (rxmatch-if (rxmatch rx-float x)
    (#f sign-string int-part frac-part frac-only exponent)
    (let ((sign (cond ((string= sign-string "-") -1.0)
                      (else 1.0)))
          (exponent (cond (exponent (string->number exponent))
                          (else 0))))
      (cond (frac-only (convert-to-float sign ""       frac-only exponent))
            (else      (convert-to-float sign int-part frac-part exponent))))
    #;not-matched #f))

(define (convert-to-float sign int-part frac-part exponent)
  (let* ((length-frac (string-length frac-part))
         (m (string-append (string-trim int-part #\0) frac-part))
         (length-m (- (string-length m) 1))
         (mantissa (* sign (string->number m 2) (expt 2.0 (- length-m)))))
         ;; mantissaは1以上2未満の浮動小数点数になる(結果が0.0の場合は除く)。
    (* mantissa (expt 2.0 (+ exponent length-m (- length-frac))))))

(define f->s2 float->string2)
(define s2->f string2->float)

(define (invariant? x)
  (let1 xx (string2->float (float->string2 x))
    (= x xx)))