気まま研究所ブログ

ITとバイク、思ったことをてきとーに書きます。

浮動小数点数演算を模倣して誤差を見てみた

f:id:AonaSuzutsuki:20200423120142p:plain

浮動小数点数が値によっては誤差が発生することは広く有名な話だと思います。
例えば、「0.1 + 0.2 == 0.3」は多くの言語処理系において等値ではありません。
というのは私も知っていましたが、ではどうして等値ではないのかと聞かれるとフワーッとしたことしか言えないのでこの機会に勉強してみることにしました。
誤差を計算で明らかにするはずが書いてたら浮動小数点数の解説になってしまったのは良かったのか良くなかったのか。

この手の分野はあまり詳しくないので間違いがあればご指摘いただけると幸いです。

浮動小数点数とは

プログラミングにおいて非常に大きな数から小さな数まで扱いたいことが多々あります。
しかしながら、馬鹿正直に「0.0000000000000000000000000000009281」みたいに表現していてはメモリ空間がいくらあってもたりません。
そこで、実数を1.2345 x 10^-5といったように、仮数と指数で表現することで効率よくこのような極端な数を扱うことができます。
このように、小数点が指数により移動することから浮動小数点数と呼ばれます。

浮動小数点数の表現は現在ではIEEE 754により規程されているため多くの処理系ではこちらが持ちられています。
C#の公式リファレンスでも「IEEE バイナリ浮動小数点数」と表記されているほどで、IEEE規格で表現されていることが伺えます。
かつてはIBM方式があったりいろいろな方式があったためにカオスな状況だったらしいです。
現在でも専門分野では独自の方式が取られているらしい。

IEEE 754では32bit浮動小数点数を基本に、16bitを半精度32bitを単精度64bitを倍精度128bitを四倍精度と呼ばれます。
一般的な型を明示する言語ではfloatを単精度doubleを倍精度として扱っていると思います。

逆に小数点を固定する表現を用いた数を固定小数点数と言います。
浮動小数点数とは異なり、小数点を特定の位に固定されるため、固定小数点数と呼ばれます。
小数ではありませんが、整数型も固定小数点数と言えます。
固定小数点数は通常の四則演算命令で処理できるため比較的高速に処理できますが、値の表現できる範囲が狭くなります。

一般的なコンピュータにはCPUにFPUが内蔵されていますが、それでもやはり固定小数点数演算の方がサイクルが少ないようです。
演算処理速度について

単精度浮動小数点数

f:id:AonaSuzutsuki:20200423120142p:plain
IEEE 754 単精度浮動小数点数は上図のようなメモリ構造で表現されます。
例として今回は0.1の浮動小数点数を入れてみました。

符号は正数なら0、負数なら1となります。

指数部は8bit符号なし整数(倍精度は11bit)で表され、1-254の範囲で数値の表現ができます。
ん?ちょっと待て。
小数を表現する場合はE-5といったように負数が出てくるけど、この範囲って負数表現できなくない??

実は指数をビット列に変換する際にバイアス値を加算することで正の値にしています。
負の値のままだと大小の比較の際に手間がかかるためにこのような処置が取られていいます。
バイアス値は後にも出てきますが、単精度では127倍精度では1023を加算します。
符号付き8bit整数なら最小値が-126-127であるため、127を加算すれば1-254と必ず整数になるわけです。
また、0や255は特殊な意味を持つため正規化数では使いません。

仮数部は23bit符号なし整数(倍精度は52bit)で表されます。
後にも出てきますが、仮数をオーバーする値は表現できないので丸めることで近似値として扱います。
逆に言えばここに入りきる値であれば誤差なく計算ができることになります。
例えば、「0.75 + 0.75 = 1.5」とか。
この数は循環しないので=が成り立ちます。
ただし、仮数は正規化により1bit省いているので実際は24bitまで正確に現す事ができます。

以下は指数と仮数の表現値についてWikipediaより

指数部 仮数部
ゼロ 0 0
非正規化数 0 0以外の任意の数
正規化数 1-254 任意の数
無限大 255 0
NaN 255 0以外の任意の数

正規化

IEEE 754では仮数部に入れる前に正規化が行われます。
正規化とは、一定のルールに則って値を変形することで、IEEE 754では1.Mとなるように小数点を移動させます

例えば、あとでも出てきますが、0.1なら二進数に変換すると「0.0001100110011」となるため、「1.100110011 * 2^-4」となります。

このルールにより、必ず1がくるため仮数部にはMの部分だけを入れれば1ビット分稼ぐ事ができます。
巷ではケチ表現とか言われてるらしい。
そのため実際には24bit分の仮数を表現する事ができます。

非正規化

対して0に極端に近いような小数の場合で指数部がオーバーするような際は正規化する事ができません
そのような数の場合は正規化を行わずに指数部を0とし、仮数部には正規化を行わずに二進数を入れます。

この場合、指数を2^-126(倍精度では-1022)とすることで指数部がオーバーする場合でも表現できるようにします。
具体的な計算については一番最後でやってみます。

浮動小数点数と誤差

上の方で「0.1 + 0.2 == 0.3」は等値ではないと述べましたが、C#でも漏れなくFalseが返ってきます。

Console.WriteLine((0.1 + 0.2) == 0.3);

こう言う場合は正確な等値を確認することはできません。
どうしても等値が必要ならdecimalを使うか許容範囲を定義するほかなさそう。
10進数浮動小数点数型であるdecimal型は10進法で計算されるためこのような二進法に起因する誤差はおきません。
といっても循環小数が出る場合(1m / 3m * 3m とか)は誤差が出るのでこちらでも確実に誤差がでないわけではありません。

では、なんでこんな現象が起こるのかと言えば、浮動小数点数を内部では0/1で表現されるために発生します。
後ほど計算しますが、0.1、0.2、0.3ともに二進数へ変換するとどれも循環するんですよ。
循環する場合は仮数の23bitに当然入りきらないため、どこかで丸なきゃいけないのでそこで誤差が発生します。

今回はそれをIEEE 754に則って演算を行って確認してみます。
なお、今回は単精度浮動小数点数で行うため、C#ではfloat型になります。
なので正確には「0.1f + 0.2f = 0.3f」です。

0.1を浮動小数点数へ

まずは0.1を十進数から二進数に変換します。
小数の二進数変換は小数に2をかけて1よりも小さければ0を、大きければ1となります。
これをかけた答えが1になるまで繰り返します。

0.1 * 2 = 0.2 ... 0
0.2 * 2 = 0.4 ... 0
0.4 * 2 = 0.8 ... 0
0.8 * 2 = 1.6 ... 1
0.6 * 2 = 1.2 ... 1
0.2 * 2 = 0.4 ... 0
0.4 * 2 = 0.8 ... 0
0.8 * 2 = 1.6 ... 1
0.6 * 2 = 1.2 ... 1
0.2 * 2 = 0.4 ... 0
0.4 * 2 = 0.8 ... 0
0.8 * 2 = 1.6 ... 1
0.6 * 2 = 1.2 ... 1
...

1になるまで繰り返したいところですが、循環するのでキリのいいところで止めておきます。
あとは上から結合していき、小数点を加えます。

0.0001100110011

次に、これを正規化により指数表現にします。
IEEE 754では先頭を必ず1.Mとなるように小数点をずらします。

1.100110011 * 2^-4

正規化より、必ず1.Mとなるため、1.を省略してMの部分を仮数部とします。
なお、二進数の数が少ないためMの部分が24bitになるように一番上の二進変換を続けます。

10011001100110011001100 1

ここで、仮数部は23bitであるため、最下位ビットを丸めます。
通常は偶数丸めみたいなので24bit目が1なら切り上げます。

10011001100110011001100 | 1
↓
10011001100110011001101

これで仮数部がもとまりました。

では、続いて指数部を求めます。
IEEE 754では指数部にバイアス値である127を足します。

-4 + 127 = 123

足したら二進数に変換します。
小数とは逆に2で割ったあまりをビット数とし、商が1になるまで続けます。
答えが出たらしたから順に結合すれば二進数になります。
なお、8bitにならない場合は8bitになるように最上位ビットを0で埋めます。

123 / 2 = 61 ... 1
61 / 2 = 30 ... 1
30 / 2 = 15 ... 0
15 / 2 = 7 ... 1
7 / 2 = 3 ... 1
3 / 2 = 1 ... 1

01111011

あとは符号の0/1を最上位ビットに追加し、順に指数と仮数を結合すれば完成です。

0 01111011 10011001100110011001101

0.2を浮動小数点数へ

0.2も同様に浮動小数点数へ変換します。
全く同じなのでサクサクっと答えを出します。

0 01111100 10011001100110011001101

0.3を浮動小数点数へ

ついでに0.3も同様に浮動小数点数へ変換しておきます。

0 01111101 00110011001100110011001

0.1 + 0.2を浮動小数点数演算で

さて、ここからが本番です。
先ほど求めた0.1と0.2を浮動小数点数演算で加算してみましょう。

まずは仮数部に省略した1.Mの1を追加します。

0 01111011 110011001100110011001101
0 01111100 110011001100110011001101

次に指数が大きい方へ桁数を合わせます。
0.1と0.2では0.1のほうが指数が小さいので指数の差の分だけ右へシフトします。
「124 - 123」なので1bit分だけシフトします。

0 01111011 110011001100110011001101
↓
0 00111101 111001100110011001100110

次に仮数部を加算します。

  111001100110011001100110
+ 110011001100110011001101
=1101100110011001100110011

すると桁上がりするので指数の大きい方へ1加算します。
今回は0.2のほうなので01111100へ加算します。
なお、桁上がりした仮数部の最上位ビットは無視します。

01111100 + 1 = 01111101

また、桁上がりしたため仮数部を1bit右へシフトします。
その前に追加した1.Mの1を省略してからシフトします。

01100110011001100110011
↓
001100110011001100110011

すると仮数が24bitになってしまったのでここでも浮動小数点数変換と同様に、23bitへ丸めます。

00110011001100110011001 | 1
↓
00110011001100110011010

これで加算完了です。
あとは符号と指数、仮数部を結合します。

0 01111101 00110011001100110011010

ところで、加算処理の際に省略した1.Mの1を追加する方法で計算してる人がいたんですが、わざわざ追加して計算するもんなのでしょうか。
ビット数は33bitになるし、最後にどうせ省略するんならいらない気がするんだけど、実際のところどうなんでしょ。
なにぶん情報がなくてこの計算はちょっと自信がない。
一応C#のメモリの内部値と一致したので結果はあってるらしい。

1追加しないと計算結果合わないものを見つけたのでやっぱり追加しないとダメです。
よくよく考えれば省略しただけで演算対象なんだから当たり前な話でした。
末尾にその話追加してます。

0.1 + 0.2と0.3を比較

ここまで計算できたところで「0.1 + 0.2」と「0.3」の浮動小数点数の内部ビットを比較してみましょう。

0.1 + 0.2

0 01111101 00110011001100110011010

0.3

0 01111101 00110011001100110011001

最下位から見て2bit分が違うことがわかると思います。
なので「0.1 + 0.2 == 0.3」は等値じゃないんですね。

このことからも分かる通り、浮動小数点数は循環しない値であれば正確な表現ができますが、循環する値は近似値で表されます。
なのでバイナリ浮動小数点数を扱う場合はそれが近似値であることを意識しましょう。
っていうかほとんど循環するので全部近似値と思った方がいいです。

加算した浮動小数点数を戻してみよう

このままでは何が何だかわからないと思うので、念のため十進数へ戻してみましょう。
変換ができたら復元もできて初めて理解が深まります。

まずは指数を元に戻します。
二進数を十進数に変換する場合は最下位ビットから順に2のべき乗を行い、最後にそれらを加算します。
f:id:AonaSuzutsuki:20200423164006p:plain
関係ないけどdraw.ioのHelveticaは上付き文字が安定しない。

01111101
↓
2^6 + 2^5 + 2^4 + 2^3 + 2^2 + 1 = 125
↓
125 - 127 = -2

次に指数を用いて仮数部を二進数の小数点表示に変換します。
ここで注意したいのが、正規化により1.Mが省略されていたことを忘れないでください。

1.00110011001100110011001 * 2^-2
=0.0100110011001100110011010

あとは二進数から十進数へ変換すれば終わりです。

0.0100110011001100110011010
2,5,6,9,10,13,14,17,18,21,22,24
=2^-2+2^-5+2^-6+2^-9+2^-10+2^-13+2^-14+2^-17+2^-18+2^-21+2^-22+2^-24
=0.30000001192

なお、符号部が1なら-にします。

浮動小数点数の最大値と最小値

本当は浮動小数点数のところで述べたかったのですが、計算ありきなので一番最後にしました。
C#の公式リファレンスではfloatの範囲が「±1.5 * E−45 から ±3.4 * E38」となっています。
最後にこれを求めて終了しましょう。

正規化数の最大値と最小値

まずは「±3.4 x E38」から。
これは上述した方法で正規化した「0 11111110 11111111111111111111111」と「1 11111110 11111111111111111111111」を十進数に戻せば「±340282346638528859811704183484516925440」となり、「≒ ±3.4 * E38」と求められます。

非正規化数最大値と最小値

ではもう一つある「±1.5 * E−45」はなんなんだってところですが、これは非正規化数の最大値と最小値です。
非正規化数も同様に「0 00000000 00000000000000000000001」と「1 00000000 00000000000000000000001」を単純に十進数に戻せばいいのですが、正規化数とは異なり、指数部は2^-126となります。

これに則って計算するのですが、正規化していないので仮数部は0.Mとします。

0.00000000000000000000001 * 2^-126
↓
1 * 2^-149
↓
1.4012985e-45

というように±1.5 * E−45ではありませんが、近似値が得られます。

省略した1.Mを追加して演算する話

0.1 + 0.2であれば1.Mを追加しなくても結果があうんですが、0.131 + 0.131だと1追加するしないでズレが発生します。

まずは1.M追加するパターンから

0.131: 0 01111100 00001100110011001100110
0.131: 0 01111100 00001100110011001100110

# 1.M追加
0 01111100 100001100110011001100110
0 01111100 100001100110011001100110

  100001100110011001100110
+ 100001100110011001100110
=1000011001100110011001100

0 01111101 00001100110011001100110
= 0.262

1.Mを追加しないパターン

0.131: 0 01111100 00001100110011001100110
0.131: 0 01111100 00001100110011001100110

 00001100110011001100110
+00001100110011001100110
=00011001100110011001100

0 01111100 00011001100110011001100
= 0.137

追加した1が残らない為省略できないものの、どうやら追加しないとダメっぽい。
1ならまた省略、0なら省略しないとかそういう話なのか?
結局自分なりの万能な答えは見出せませんでした。

参考サイト

  1. IEEE 754 - Wikipedia
  2. 浮動小数点加算器 - 桁合わせと加算
  3. 浮動小数点数について本気出して考えてみた
  4. 【5分で覚えるIT基礎の基礎】ゼロから学ぶ2進数 第4回
  5. 浮動小数点数内部表現シミュレーター
  6. バイナリ計算機
  7. Binary to Decimal converter
  8. Decimal to Binary converter
  9. 浮動小数点数の内部表現を取得してみよう
  10. 固定小数点数について
  11. 整数と浮動小数点数