文字列少数点数表記を IEEE754 倍精度浮動小数点数にエンコードする方法

Sukesan1984

2020年10月22日 22:30

42tokyoというパリ発のエンジニア養成機関で学んでいるのですが、表題の関数が課題で必要になったため、絶対に将来忘れる気がしてならないので自分への備忘録として書き記しておきたいと思います

恥ずかしながら、今まであまりちゃんと浮動小数点数について考えたことがなく、かなり初心者なので、間違っていることがあるかもしれないのですが、そのときはご指摘いただけると嬉しいです

こちらの記事を多分に参考にさせていただいています

何度か読み直してみましたが、大変わかりにくくてすみません・・・
自分でも読んでてよくわからなくなりそうになります
都度、参考実装へのリンクを貼っているのでもしか、興味ある方いればそちらを見てもらったほうが早いかもしれません

ソースコードはこちら

動機

まず、表題の処理が何故必要なのか？というと42の課題で、Raytracingに関する課題があり、その中でファイルから物体の位置やら色やらを読み込んで表示させる必要があり、その座標が少数点数を扱うからです

C言語なので、stdlib.hなどをincludeできればatofが使えるのですが、課題の条件として、stdlibの使用は禁じられているため自作する必要があります

#include <stdlib.h>
int main()
{
       double d;
       d = atof("123.456");
       return (0);
}

やりたいのはこれのatofを作りたいということになります

前提

言語: C言語
インプット: 10進数文字列 (char *)
ex) 123.456, 1.4e+15, -1.4e-15など
アウトプット: 倍精度浮動小数点数 (double)
double atof(char *str);

を実装したい。

まず思いつく方法

同じような関数であるatoi（文字列からintへの変換をする関数）は、オーバーフローに気をつけて、数値を1文字ずつ読み取って、10倍して足していく方法で比較的簡単に実装できます


while (*str != '\0')
    i = i * 10 + (*str++ - '0');

イメージはこんな感じですね

それと同じように整数部分を処理し、"."が現れたらそこからは、値を1/10して足していく方法が考えられます

途中まで読み込んだときの結果をどう処理するか？
次第ではあるのですが、上記と同じようにdoubleに格納していく方法ではうまく行かないことがわかります

＊）実は、あとの処理と合わせて一緒にやれば一文字ずつ読み込む方法でもできるかもしれないです。試してないですが・・

例）
18014398509482012という数字を考えてみると、この数字自体はdoubleで正確に表現できます
そして、上記のロジックで、一つ手前まで読み込んで最後の数字を読み込むときのことを考えると

18014398509482012 = 1801439850948201 * 10 + 2

と書き表されますが、
10倍した18014398509482010という数字は、doubleで正確に表現できず
18014398509482008という値に丸められてしまい、その値に2を足した最終結果も
18014398509482008となり、 18014398509482012にはなりません・・なので、1文字ずつ読んでいく方法ではうまく行かないことがわかります

浮動小数点のビット列表現の読み方

まず最初に倍精度浮動小数点数のビット列表記の読み方を見ておきます

64bitのbit列で、レイアウトは以下のようになります

1bit目が符号部で0なら正の数, 1なら負の数 b(1bit)
2bitから12bit目までが指数部 e (11bit)
13bit目から64bit目が仮数部 m (52bit)

を表します

これは10進数表記でみると以下のように表せます

それぞれの値の範囲は
指数部に関しては11bitあるので、0 - 2047まで表現力があるわけですが
すべて0の場合を非正規化数、すべて1の場合を無限大やNaNに利用するため
実際には、1 - 2046の表現力があり、
小数点以下も表すために-1022から1023になります

一方仮数部に関しては、52bitあるので、0 - 4503599627370495(2^52 - 1)の表現力があり、かつ常に右から53bit目に1がたっているとするケチ表現をとっています

例えば、以下のようなbit列であれば

10進2進混じった書き方をすると

となり

と読めます

stringからdoubleへの変換のステップ概要

では、どのように読み込んで上記表現に変換するかを説明します

最終ゴールとしては

・符号部 b (0 or 1)
・2進数指数 e (-1022 - 1023)
・仮数部となる整数 m (0 - 4503599627370495)

を得ることです

この3つがあれば、IEEE754の倍精度浮動小数点数にエンコードできます

大きく分けると4つのステップがあります

1) 1文字ずつを読み込んでstructに格納する

1文字ずつ読んで十進数の状態のまま符号部、指数部、仮数部のstructとして中間状態を表現する

typedef struct  s_prep_number {
       int             negative;
       int32_t         exponent;
       uint64_t        mantissa;
}               t_prep_number;

2) m * 10^eをm' * 2^e'に変換する（後ほど詳しく）

mを10倍したり10で割ったりしながら、上記のように変換します

3) m'を一番ケチ表現に当たる部分まで左ビットシフトする

ケチ表現に当たる部分にビットを立てるために左シフトします
（その一つ右のビットから52ビットを仮数部としてエンコードするため）

4) IEEE754の倍精度浮動小数点数形式にエンコードする

ここまでくれば
符号部は1)のときに0か1かにできますし
e'は-1022から1023に収まっていればそのままビット列に変換できますし（範囲外の場合はそのまま、INFあるいは0にする）
m'もケチ表現にあたるビット列を除いてその右のビットから52ビット列分よみとってエンコードすればうまくいきます

これを一つずつ説明していきます

1) 1文字ずつを読み込んでstructに格納する

参考にさせていただいた記事の通りstate machineで表現できるのでそれに則って状態を遷移しながら最後の文字まで読み取っていきます

0000032
0000.73
3.474650000
+32.746
.43
3.6E00000004
7E+2

こういった表現を左から順番に読んでいって現れた文字に合わせて状態を変えていきます

ここは結構シンプルなのでそんなもんかと思って貰えれば次のステップ2)まで飛ばしてください

a) 空白をスキップする

空白以外の文字が来るまで、文字列を進めていきます

空白が来たらb)に進みます

b) 符号を読み取る

- が来た場合にはnegativeフラグを立てて次の文字に進み、c)へ

+が来た場合は、次の文字へ進み、c)へ

数字が来た場合は、そのままc)へ

それ以外の文字が来た場合は、STOPへ

c) .の前に0が続く場合の処理をする

0が来る限り次の文字に進めてc)に戻る　

.が来たら次の文字へ進めてd)へ

数字の場合は、e)へ

d) 0.のあとに0が続く場合の処理をする

0が来たら、exponentを-1して次の文字へ進めてd)へ

それ以外の文字が来たらf)へ

e) 仮数部の整数部分を読み取る(123.456の場合の123の部分)

数字が来たら、整数部分としてatoiと同じような処理で10倍してmantissaを更新していきます e)へ

.が来たら整数パートが終わりなので、文字を進めてf)へ

それ以外の文字（e,Eとか）の場合はf)へ

f) 仮数部の小数部分を読み取る(123.456の場合の456の部分)

数字が来たら、整数部分としてmantissaを更新して、exponentを-1します（123.4の場合整数部分を1234にして指数部分を10 ^ -1とするイメージです）

e,Eが来た場合は、一文字進めてg)へ

それ以外の文字が来た場合はg)へ

g) 指数部の符号を読み取ります (例えば、123.456e-10の-の部分）

+が来たら、一文字進めて 8)へ

-が来たら、指数部の符号のフラグを立てて、一文字すすめて8)へ

それ以外の文字が来たら、8)へ

h) 指数部の後ろに0が続く場合スキップする(3.6e00004は3.6e4にする）

0が来たら、一文字進めてh)へ

それ以外が来たらi)へ

i) 数字が来たら、指数部を整数として読み込んでいく

数字が来たら、指数部としてatoiのように整数として読み込んでi)へ

それ以外が来たらSTOPへ

どういう状態になるかというと

123.456という少数点数を読み込むと

sign: 0
exponent: -3
mantissa:123456

という形で表されます

ちなみに、途中でオーバーフローの対応をする必要がありますが、仮数部としては何桁分読み取れば十分なのでしょうか？

仮数部は最大でケチ表現を入れて53ビットすべてに1が立つパターンが最も桁数が多くなります
この場合、9007199254740991の16桁ですがWikipediaによると、binary64の場合、17桁以上20桁以下を保持するようにしておけば、十分であるので、今回は18桁分保持するようにしています

また指数部に関しても

2^1024が 10^308のオーダーになるので、それを超える場合
非正規化数の2^-1074を仮数部のビットシフト分(51bit左シフト)を考慮して2^(-1074-51)が10^-339のオーダーになるので、それを下回る場合をINFあるいは、0にしてしまっています

前準備）96ビットの演算を用意する

次のステップに移る前に、64ビットで仮数部を処理すると情報を落としてしまうので96ビットで演算できる準備をしておきます

具体的には、96ビットを32ビット* 3に分割しておき、既存の64ビットの数値を下位64ビットにマッピングして諸々の処理を行えるようにします

 typedef struct  s_bit96 {
       uint32_t        s2;
       uint32_t        s1;
       uint32_t        s0;
}                t_bit96;

このような構造体を準備します

t_bit96 add96(t_bit96 s, t_bit96 d);

t_bit96 sub96(t_bit96 s, t_bit96 d);

t_bit96 lsr96(t_bit96 s);

t_bit96 lsl96(t_bit96 s);

このように、96ビット同士の足し算、引き算、左シフト、右シフトを用意しておきます

この実装自体はそこまで複雑ではないので割愛します

**2) m * 10^eをm' * 2^e'に変換する**

ここまでで1.23456e+8という数字であれば

123456 * 10 ^ 3という表現として読み込むことができています

sign: 0
exponent: 3
mantissa:123456

これを2進数表現にする必要があります

ここでポイントとなるのが仮数部の123456はuint64として読み込んでいて、これは内部的にはすでに2進数表現と同等であるといえます

なので、10 ^ 3を2進数表現に変換できれば良いわけです

一般化して、

と表現される10進数表現を考えます（mは仮数部のuint64です）

指数部のeが正か負かで処理が分かれますが、基本方針としてはeを0にして

と表現することを目指します

mはuint64の表現力があるのでそれを超えるまでは、mを10倍あるいは1/10に変更していって、超える場合にはビット演算でmを右シフトしてe'を+1していきます

具体的に見ていきます

eが正の場合

123456 * 10 ^ 3

を考えるとこの場合はシンプルで

まず10倍していくことを考えると

123456を10倍しても、1234560はuint64の範囲に収まる整数なので

1234560 * 10 ^ 2とできます

同様に12345600 * 10 → 123456000としてもすべてuin64に収まるのでこの場合は

123456000 * 10^(0) * 2^(0)と変換できたわけです

mがuint64の範囲を超える場合はどうするか？

まず、先程準備した96ビットの構造体としてmを表現します

そして上位4ビットの範囲を0にしておいて、最大92ビットまで表現することとします
上からmを10倍した時にもし仮に上位1から4bitにビットが立つようであれば
1から4bitが再び0になるまでmを右シフト(1/2倍)しつつe'を+1して辻褄を合わせます

ここで、10倍というのが実は(2 + 8)倍なので、m << 1 + m << 3とすることで、ビット演算だけで表現できるミニtipsがあります

eが負の数の場合は基本的に10で割り込んでいくので、できるだけ表現力をキープするために、まずmを一番左までビットシフトします
その後、ちょっと面倒ですが、eが0になるまで96bit表現を10で割る処理をした後e'をincrementしていきます
一回処理が終わるごとにまた、左に寄せることでできるだけ精度を保つようにします

3) m'を一番ケチ表現に当たる部分まで左ビットシフトする

ここまでの処理で、96bitで表された仮数部分が整数で得られているのですが
例えば2進数の111というのが仮数部だったとしてbit列としては
16進数で表すと
0x00000000 00000000 00000007
のように右寄せになります
つまり右から3bit(うち一番左がケチ表現分)を
倍精度浮動小数点数にエンコードしたときには、左から見た時に13bit目と14bit目にビットをたてることになります
0x_(1bit) ___________(11bit) 11000...
このように右寄せになっていると、右から何ビット分読めばいいかわからないので
左寄せにしてしまいます

つまり、96bitのうち左から4bit目までにビットが立つまで左シフト(2倍)していきその都度辻褄を合わせるために、e'をデクリメント(1/2)しておきます

そうすることで、111を
0x00000000 00000000 00000007を 90 左シフトして
0x1C000000 00000000 00000000 * e^(-90)としておきます(注 e'は指数部と合わせた値になる）

4) IEEE754の倍精度浮動小数点数形式にエンコードする

ここまでの処理で、

符号部は、0か1
e'は整数として -1022 から 1023の値の範囲あるいはその範囲に収まらない値
m'は96bit整数で、かつ左から4bit目がケチ表現のbitを表す整数

になっているので、エンコードしていきます

実際には、e'は非正規化数を考えると -1074まで取れるので（後述）

e'が-1074より小さい値の場合は符号部に合わせて+0あるいは-0にエンコードします

e'が1023より大きい場合は符号部に合わせて+INFあるいは-INFにエンコードします

e'が-1022から1023の場合は、
e'に1023を足して1から2046の値に変換したあと52bit分左シフトしたものとm'のうち5bit目から56bit目までと符号部が-の場合は1を63bit分左シフトしたものを|でつなぎます

これで完成です

非正規化数の対応

参考にさせていただいた記事でも非正規化数は対応していなかったのですが非正規化数を対応する場合についても書いておきます

例えば、このようなビット列で考えると

仮数部の先頭にだけビットが立っているのですがこの場合正規化数のときは

こういう形で、必ずケチ表現で表すこととしていましたが、この場合は

これが正しい表現になります

他にも、非正規化数の中で最も小さい値は

このようなビット列は

という値になります

これに対応するには、仮数部のエンコードの部分だけ修正すればよく
正規化数では、4bit目をケチ表現として落としていたものを52bitの最初のbitとしてみなしてエンコードしてあげます

最後に

いろいろ複雑な処理をしていてなんだかんだこの一ヶ月ぐらいずっと浮動小数点数とにらめっこしていた気がします

そして、もっとうまいやり方があるんじゃないか？と思っていろいろ調べたのですがなかなか情報がなく苦しみました

参考にさせていただいた記事がわかりやすく大変ありがたかったです

もっとこうやったら簡単にできるよというのがあればぜひ教えていただけると幸いです

改めて記事を読み直してみましたが、わかりにくいことこの上ないのですが、もっとわかりやすく説明できる様になりたいものです・・

#IEEE754

#浮動小数点数

#42Tokyo

この記事が気に入ったらサポートをしてみませんか？