浮動小数点誤差の自分メモ
普段適当にfloat32やfloat16使ってるけど、適当すぎてたまに困るので勉強のために机上での考察とプログラムを動かして検証してみた。
オーバーフロー、アンダーフロー
- 指数部が表現可能な範囲外を超えるケース
- float32なら指数が+127を超えるとオーバーフローしてinfになる
- 下限の-126を下回ると非正規化数
- 演算結果が非正規化数の最小値より小さな値になるとアンダーフローして0になる
確認コード
// 16進数浮動小数リテラルはC++17以降で使える std::cout << 0x1p+127f * 0x1p+127f << std::endl; // オーバーフロー（生の数） std::cout << -0x1p+127f * 0x1p+127f << std::endl; // オーバーフロー（負の数） std::cout << 0x1p-126f * 0x1p-25f << std::endl; // アンダーフロー
実行結果
inf -inf 0
float16
- float16は指数部が5bitしかない
- 指数部は-14〜+15
- 15乗と聞くと大きな数値に聞こえそうだが２進数なので215==32768
- 仮数部がall 1でも65504が最大値
- 普通にuint16_tの最大値の方が大きい
- 255×255はセーフだが、256×256はオーバーフローする
- 正規化数の最小値は ±2^-14
- こちらも 1/16384 と正直小さな値と言い難い
- 非正規化数の最小値は ±2^-24
- 指数部2^-14 仮数部 2^-10
- 1/16777216。微妙な数値。うっかりアンダーフローしそう
- 以前float16でディープラーニングモデルを学習させようとしたら速攻でNaNになって学習できなかった
- おそらくアンダーフローが発生して0になり0除算でNaNという流れ
実験用ルーチン
NaN、INFなどの特殊値や非正規化数、なんなら±0にすら非対応のなんちゃって実験用ルーチンのコード。文字列のオレオレバイナリ表現とfloat値の相互変換ができる。（実用には耐えられない）
#include <cstdint> #include <string> #include <algorithm> #include <iostream> std::string remove_char(const std::string &str, char c) { std::string result = str; auto new_end = std::remove(std::begin(result), std::end(result), c); result.erase(new_end, std::end(result)); return result; } std::string cleansing(const std::string &str) { auto s = remove_char(str, ' '); // 半角空白を削除 return remove_char(s, ','); // カンマを削除 } float str_to_fp32(const std::string &str) // なんちゃってバイナリ表現文字列 → float変換 { auto s = cleansing(str); // 区切り文字と空白を削除 std::uint32_t raw = 0; // バイナリ表現 // 符号ビット if (s[0] == '-') { raw = 0x80000000; } // 仮数部 std::uint32_t mantissa = 0; for (auto i=0;i<23;i++) { auto bit = (s[3 + i] == '1')? 1: 0; mantissa <<= 1; mantissa |= bit; } // 指数部 auto exponent = std::stoi( std::string(std::begin(s) + 3 + 23 + 4, std::end(s)) ); exponent += 127; raw += exponent << 23; raw += mantissa; return *reinterpret_cast<float *>(&raw); } std::string fp32_to_str(float f) // float → バイナリ表現文字列変換 { std::uint32_t v = *reinterpret_cast<std::uint32_t *>(&f); std::string result; result = (f > 0)? "+": "-"; result += "1."; // 仮数部 std::uint32_t mantissa = v & 0x7fffff; for (auto i=0;i<23;i++) { auto bit = (v >> (22 -i)) & 0x01; result += (bit == 1)? "1": "0"; if((i % 4) == 2 && i!=22) { result += ","; } } // 指数部 auto exponent = (v >> 23) & 0xff; result += " x 2 ** " + std::to_string(static_cast<int>(exponent) - 127); return result; }
情報落ち
- 加減算する時に２数の差が大きいと計算結果が変化しない
- 2仮数部のbit数 だけ差があると情報落ちする
- 計算結果は指数部の大きい方になる
- 例えば 0.5f は 1.000...×2^-1 で指数部-1なので仮数部は 2^-2 〜 2^-24 の 0/1 を表現できる
- 指数部の数値によって分解能は変化する
- 2.0f は 1.000...×21 指数部1 で仮数部は 20〜2^-22
- 1.0f は 2^-1〜2^-23
- 0.25f は 2^-3〜2^-25
- 指数部が小さい(-126に近づく方向)ほど仮数部の分解能が小さくなる
// 情報落ちしないギリギリのライン // fp32の場合、仮数部23bitなので2^23までの差は情報落ちしない std::cout << "1.0 + 2**-23 -> " << fp32_to_str(1.f + std::pow(2, -23)) << std::endl; std::cout << "0.5 + 2**-24 -> " << fp32_to_str(.5f + std::pow(2, -24)) << std::endl; // 情報落ちする。2^23を超える差がある std::cout << "1.0 + 2**-24 -> " << fp32_to_str(1.f + std::pow(2, -24)) << std::endl; std::cout << "0.5 + 2**-25 -> " << fp32_to_str(.5f + std::pow(2, -25)) << std::endl;
実行結果
1.0 + 2**-23 -> +1.000,0000,0000,0000,0000,0001 x 2 ** 0 0.5 + 2**-24 -> +1.000,0000,0000,0000,0000,0001 x 2 ** -1 1.0 + 2**-24 -> +1.000,0000,0000,0000,0000,0000 x 2 ** 0 0.5 + 2**-25 -> +1.000,0000,0000,0000,0000,0000 x 2 ** -1
情報落ちと桁落ちのコンボ
- 桁落ちそのものは「有効数字減って困るんだっけ？」くらいの認識だった
- 情報落ちと桁落ちが組み合わさって誤差を観測するケースを体験したので共有しておく
- ①で小さな値(指数部-8)と大きな値(指数部0)を加算する
- 情報落ちによって小さな値の仮数部後半が消失する
- ➁で先程加算した値を減算する
① +1.000,1000,1000,1000,1000,1000 x 2 ** -8 + +1.000,0000,0000,0000,0000,0000 x 2 ** 0 結果１ +1.000,0000,1000,1000,1000,1001 x 2 ** 0 ➁ +1.000,0000,1000,1000,1000,1001 x 2 ** 0 - +1.000,0000,0000,0000,0000,0000 x 2 ** 0 結果２ +1.000,1000,1000,1001,0000,0000 x 2 ** -8
乗算の誤差
- 個人的には指数部は単純な指数同士の加減算イメージ
- 仮数部も乗算
除算の誤差
お試しコード
for (auto i=3;i<15;i++) { auto lhs = 15.f; std::cout << "15 / " << i << " -> " << fp32_to_str(lhs / (i * 16)) << std::endl; }
実行結果
15 / 3 -> +1.010,0000,0000,0000,0000,0000 x 2 ** -2 15 / 4 -> +1.111,0000,0000,0000,0000,0000 x 2 ** -3 15 / 5 -> +1.100,0000,0000,0000,0000,0000 x 2 ** -3 15 / 6 -> +1.010,0000,0000,0000,0000,0000 x 2 ** -3 15 / 7 -> +1.000,1001,0010,0100,1001,0010 x 2 ** -3 15 / 8 -> +1.111,0000,0000,0000,0000,0000 x 2 ** -4 15 / 9 -> +1.101,0101,0101,0101,0101,0101 x 2 ** -4 15 / 10 -> +1.100,0000,0000,0000,0000,0000 x 2 ** -4 15 / 11 -> +1.010,1110,1000,1011,1010,0011 x 2 ** -4 15 / 12 -> +1.010,0000,0000,0000,0000,0000 x 2 ** -4 15 / 13 -> +1.001,0011,1011,0001,0011,1011 x 2 ** -4 15 / 14 -> +1.000,1001,0010,0100,1001,0010 x 2 ** -4
