浮動小数点数のパーサを書いてみた

先日、こういうツイートを見かけた。で、さいきん浮動小数点数づいているのもあって自分も浮動小数点数の文字列表記をパースして単精度および倍精度浮動小数点数に変換するパーサを書いてみた。

非常に苦労したので書いた
もっといいやり方がありそうだし、もっとうまい説明ができるかもしれないけどこれが限界でした・・・

文字列少数点数表記を IEEE754 倍精度浮動小数点数にエンコードする方法｜Sukesan1984 #note https://t.co/2v5f1eMzea

— Sukesan1984 (@sukesan1984) October 22, 2020

この記事では、自分が実装したパーサを説明する。説明はいいからコードを見せろ、という人はこちらをどうぞ。

https://github.com/hkurokawa/FloatingPointNumberParser

ゴール

今回は浮動小数点数の理解を深めるために自分でパーサを書くのが目的である。したがって、パフォーマンスは気にしない。また、シンタックスエラーも気にしないことにした¹。さいわい先ほど引用したツイートへの補足ツイートにもあるがGo言語の標準ライブラリのテストケースがちょうどよさそうなので、このテストが通ることを目標とする。

https://golang.org/src/strconv/atof_test.go

実装方針

基本的にはWikipediaの記事にあるのと同じことをすればよい。

入力を10進数表記の文字列としよう（2進数表記でも同じ議論が成り立つ）。たとえば 0.15625 といった文字列だ。これを2進数表記にしたい。最初に気付くのは、無限精度のままこの10進数表記の値を扱う必要があるということだ²。一見、floatなら仮数部が23桁なので2進数表記で小数点以下23桁、切り上げや切り下げを考えても24桁まで見ればよさそうに思えるが、そうではない。最終的に最近接丸めをするには小数点以下24位が0だったとしても、そのさらに先で0以外の値があればそれは結果に影響する。したがって10進数表記の文字列を1文字目から順に見ていったとして、どこかでこれ以上読まなくてよいという場所は存在しない（ずっと0が続くのか、それともどこかで0以外が登場するのかは知る必要がある）。

以上のことから、まず10進数表記の文字列はそのまま無限精度の10進数の値として格納することにする。

つづいてこれを2進数表記に変換することを考えよう。整数部分については既知なので以降では小数部分についてだけ考える。

もっともナイーブな方法は変換したい値が\({0.5}\)より大きいか調べ、大きければ小数点以下第1位のビットを立てて\({0.5}\)を引く。つぎにその引いた結果と\({0.25}\)を比較し、それより大きければ小数点以下第2位のビットを立てる、という方法だ。この操作を行うと、たとえば10進表記の\({0.15625}\)という値は2進表記で\({0.00101}\)となることがわかる。

これは2進数の定義そのものなので理解はしやすいが、無限精度の10進数で引き算をやるのはあまりぞっとしない。ちょっと考えると、この\({0.5}\)を引いて、つぎは\({0.25}\)を引いて、という操作は、対象の値を\({2}\)倍して\({1}\)を引いて、さらに\({2}\)倍して\({1}\)を引くという操作と等価なことがわかる。じっさい、さきほどの\({0.15625}\)にこの操作を行うと同じ結果が得られる³。

さて、以上で無限精度の10進数を2進表記する方法がわかった。また、浮動小数点数の定義上、この値を\(a \times 2^{k}\)という形式にする必要がある。ただし\({a}\)は範囲\({[1, 2)}\)に収まる実数。すなわち、対象の数に\({2}\)もしくは\({\frac{1}{2}}\)を掛けて\({[1, 2)}\)の範囲に収める必要がある。このことから、無限精度の10進数が\({2}\)倍および\({\frac{1}{2}}\)倍の演算をサポートしている必要がある。

parseFloatのおおまかな実装

以下では、さきほど述べた\({2}\)倍および\({\frac{1}{2}}\)倍の演算をサポートした無限精度の10進数を表すクラス BigDecimal があるとする。これを使って値を浮動小数点数に変換する。

まずは\({[1, 2)}\)の範囲に収めよう。\({2}\)倍あるいは\({\frac{1}{2}}\)倍した場合はそれに合わせて指数部を増減させる。

    BigDeciaml d = new BigDecimal(s); // 10進数文字列表記を無限精度の10進数に保存する
    int exponent = 0;
    while (d.isEqualToOrGreaterThanTwo()) { // 2以上なら2未満になるまで1/2を掛ける
      d.divideByTwo();
      exponent++;
    }
    while (d.isLessThanOne()) { // 1未満なら1以上になるまで2を掛ける
      d.multiplyByTwo();
      exponent--;
    }
    exponent += 127; // バイアスを足す

さて、ここで d は\({[1, 2)}\)の範囲に入っているはずである。あとはこの小数部を2進数表記すればよい。なお d.discardNumberPart() は整数部分を0クリアするメソッドである。

    d.discardNumberPart();
    int mantissa = 0;
    for (int i = 22; i >= 0; i--) {
      d.multiplyByTwo();
      if (!d.isLessThanOne()) {
        mantissa |= 1 << i;
        d.discardNumberPart();
      }
    }

dを2倍して\({1.0}\)以上であればビットを立て、\({1}\)を引いている。

以上でほぼ終わりである。あとは符号部を組み合わせれば単精度浮動小数点数になる。

    int sign = d.isNegative() ? 1 : 0;
    int bits = mantissa | (exponent << 23) | (sign << 31);
    return Float.intBitsToFloat(bits);

いくつかのエッジケースの処理

ここまでで “1.0” や “3.1415” といった文字列はパースできるはずである。しかし、いくつかのエッジケースについて考慮が漏れている。

まずは “0.0” が与えられた場合。この場合はいくら2倍しても\({1}\)以上になることはないので無限ループになる。したがって “0” に等しい場合はearly returnする必要がある。

    if (d.isZero()) {
      return Float.intBitsToFloat(sign << 31);
    }

続いて、オーバーフローの場合。単精度浮動小数点数では\({2^{127}}\)を超える数を表すことはできない。したがって +Infinity もしくは -Infinity を返す必要がある。

    if (exponent > 127) return d.isNegative() ? Float.NEGATIVE_INFINITY : Float.POSITIVE_INFINITY;

最近接偶数丸め

さらに最近接偶数丸めも対応する必要がある。まず仮数部で表せる桁の最下位の1つ下の桁が2進数表記で\({1}\)かどうかを調べる。これはいままでと同様に\({2}\)倍して\({1}\)以上になるか調べればよい。ここまでで対象の数が単精度浮動小数点数で表せる数のちょうど中間か中間よりも少しだけ0より遠い側にあることが分かる。

    // 仮数部はすべて詰め終わっており、dには0.xxxxという形でそれより下の小数部が格納されている
    // 最近接偶数丸め
    if (!d.isZero()) {
      d.multiplyByTwo();
      if (!d.isLessThanOne()) { // dは1.xxxxという形
        d.discardNumberPart();
        if (d.isZero()) {
          // ちょうど中間なので偶数になるように丸める
          if ((mantissa & 1) == 1) {
            mantissa++;
          }
        } else {
          // 0より遠い側なので最近接にするために仮数部を1増やす
          mantissa++;
        }
      }
    }

以上で最近接偶数丸めができた。めでたしめでたし、となりたいところだが、そうは問屋が卸さない。実は仮数部を1増やしたことによって仮数部がオーバーフローする可能性があるのだ。すなわち仮数部23ビットすべてに1が立っている状態で\({1}\)を足せばとうぜん仮数部は足りなくなる。この場合は指数部を\({1}\)増やして仮数部は\({0}\)にクリアする必要がある。さらに指数部を\({1}\)増やしたことによって全体がオーバーフローしないかチェックする必要もある。

   if (mantissa == 0x800000) {
     // 仮数部が23ビットを超えている
     mantissa = 0;
     exponent++;
     if (exponent > 127) {
       return d.isNegative() ? Float.NEGATIVE_INFINITY : Float.POSITIVE_INFINITY;
     }
   }

さて、ここまでで正規化数はすべてパースできるはずだ。最後に非正規化数をサポートする。

非正規化数

非正規化数は対象の値が\({2^{-126}}\)より小さい場合に必要になる。まず2進数表記で\({0.xx\ldots x \times 2^{-126}}\)という形にして、そのうえで指数部は\({0}\)にクリアしておく。

    if (exponent >= -126) {
      // 正規化数
      // バイアスを足す
      exponent += 127;
    } else {
      // 非正規化数
      while (exponent < -126) {
        d.divideByTwo();
        exponent++;
      }
      exponent = 0;
    }

完成版

以上をすべて行った文字列をパースして単精度浮動小数点数を返すコードが以下である。今後バグが見つかった場合はGitHubレポジトリを更新していくので、最新版のコードが見たい方はそちらを参照いただきたい。

  public float parseFloat(String s) {
    BigNumber d = BigNumber.parse(s); // 10進数文字列表記を無限精度の10進数に保存する
    int sign = d.isNegative() ? 1 : 0;
    if (d.isZero()) {
      return Float.intBitsToFloat(sign << 31);
    }
    int mantissa = 0;
    int exponent = 0;
    while (d.isEqualToOrGreaterThanTwo()) { // 2以上なら2未満になるまで1/2を掛ける
      d.divideByTwo();
      exponent++;
    }
    while (d.isLessThanOne()) { // 1未満なら1以上になるまで2を掛ける
      d.multiplyByTwo();
      exponent--;
    }
    if (exponent > 127) return d.isNegative() ? Float.NEGATIVE_INFINITY : Float.POSITIVE_INFINITY;
    if (exponent >= -126) {
      // 正規化数
      // バイアスを足す
      exponent += 127;
    } else {
      // 非正規化数
      while (exponent < -126) {
        d.divideByTwo();
        exponent++;
      }
      exponent = 0;
    }
    d.discardNumberPart();
    for (int i = 22; i >= 0; i--) {
      d.multiplyByTwo();
      if (!d.isLessThanOne()) {
        mantissa |= 1 << i;
        d.discardNumberPart();
      }
    }
    // 仮数部はすべて詰め終わっており、dには0.xxxxという形でそれより下の小数部が格納されている
    // 最近接偶数丸め
    if (!d.isZero()) {
      d.multiplyByTwo();
      if (!d.isLessThanOne()) { // dは1.xxxxという形
        d.discardNumberPart();
        if (d.isZero()) {
          // ちょうど中間なので偶数になるように丸める
          if ((mantissa & 1) == 1) {
            mantissa++;
          }
        } else {
          // 0より遠い側なので最近接にするために仮数部を1増やす
          mantissa++;
        }
        if (mantissa == 0x800000) {
          // 仮数部が23ビットを超えている
          mantissa = 0;
          exponent++;
          if (exponent > 127) {
            return d.isNegative() ? Float.NEGATIVE_INFINITY : Float.POSITIVE_INFINITY;
          }
        }
      }
    }

    int bits = mantissa | (exponent << 23) | (sign << 31);
    return Float.intBitsToFloat(bits);
  }

今後の課題

冒頭にも書いたが、この実装はパフォーマンスを無視している。たとえば BigDecimal の演算はその桁数を \({n}\) とすれば毎回 \({O(n)}\) かかってしまう。したがって\({2}\)で割って\({2}\)未満にする操作は \({O(n \log n)}\) かかるので “1e+400000” みたいな文字列を与えられると処理が現実的な時間で終わらない。

また、これは手抜きなのだが E 表記や p 表記の指数部は Integer.parseInt(String) を呼んでいるだけなのでそこでオーバーフローする。

今回はGoのテストケースをそのままコピーしたが、これらの指数部が大きすぎるテストケースはコメントアウトせざるを得なかった。

ひとまず考えられる改善としては指数部のパース時にfloatもしくはdoubleで表せる範囲を超えていることが確実になった時点でパースを打ち切ることだろう。そのうえで無限精度の10進数型の演算についてはいくつかアルゴリズムがあるようなので、それを試してみようと思う⁴。

感想

というわけで、自分で浮動小数点数のパーサを書いてみた。ハマりポイントは最近接偶数丸めで仮数部が桁あふれをするところで、これはGoのテストケースを実行するまで気付かなかったので、自分で実装してよかったと思うポイントだ。

あと、自分でテストケースを書こうとすると、たとえば無限精度で \({2^{-53}}\) の値が欲しくなることがある。適当な電卓では表示が打ち切られてしまうので https://keisan.casio.com/calculator のような有効桁数を指定できる電卓を使う必要があった。

最後に、面識はありませんが、きっかけとなった記事を書いてくれたSukesan1984に感謝します。この記事が読者のなにかしら参考になれば幸いです。

1. たとえばJavaの浮動小数点数リテラル表記としては正しくない 1.p といった表記も今回は受容している。 ↩

2. 厳密には無限精度である必要はないのだが、最終的に10進数表記でどれだけの有効桁数を考慮すればよいかというのはそれほど自明ではない。今回の実装では簡単のために無限精度で扱う。 ↩

3. 0.15625 → 0.3125 → 0.625 → 1.25 → 0.25 → 0.5 → 1.0 → 0.0 ↩

4. 有名なのは Will Clinger, “How to Read Floating Point Numbers Accurately”, ACM SIGPLAN ‘90, pp.92–101, 1990. らしい。 ↩