これは何？
ruby の数値の等値 という記事の冒頭に軽い気持ちで書いた C 言語の話が面白かったのでここに記す。
基礎知識
Usual arithmetic conversions
C89 だか C99 以降、四則や比較の二項演算においては “usual arithmetic conversions”「通常の二項変換」と呼ばれる変換を事前に行う。
この「通常の二項変換」は、たとえば float と 整数の場合、整数を float に変換する。
このルールは、C11、C17 などでも変わっていない。と思う。
現象
たとえば
#include <stdint.h>
#include <stdio.h>
static inline void
show_float( char const * msg, float x ){ printf( "%s = (float)%f\n", msg, x ); }
static inline void
show_double( char const * msg, double x ){ printf( "%s = (double)%f\n", msg, x ); }
#define show(msg, X) _Generic((X), \
float: show_float, \
double: show_double)(msg, X)
int main()
{
float f2 = 99999999.0f;
int32_t i2 = 99999999;
printf( "f2==i2 is %s\n", f2==i2 ? "true" : "false" ); // [A]
show("f2-i2", f2-i2); // [B]
float x = 1e-30;
show("x*x", x*x); // [C]
printf( "x*x is %e, %s\n", x*x, (0.0f==x*x ? "zero" : "NOT zero")); // [D]
return 0;
}
普通の場合
A: i2 を float に変換する。99999999 は（丸めモードによるかもしれないけど） float にすると 100000000.0f になる。 f2 も 100000000.0f なので、等しい。したがって、 f2==i2 is true と出力される。
B: [A] と同じ変換のあとで減算なので、 0.0f になる。_Generic で show_float に到達する。0.0f は ゼロなので、f2-i2 = (float)0.000000, which is zero が出てくる。
C: x*x は、1e-60 ぐらいなので、 float の最小値を下回り、 0.0f になる。_Generic で show_float に到達する。0.0f は ゼロなので、x*x = (float)0.000000, which is zero となる。
D: 前述の通り xx は 0.0f なので、0.0f==xx は当然 真になる。なので x*x is 0.000000e+00, zero と出力される。
呪われた世界
手元で、 gcc-11 -std=c11 -pedantic -mfpmath=387 でコンパイルすると上記とは異なる結果となった。
| 記号 | 普通の世界 | 呪われた世界 |
|---|---|---|
| A | f2==i2 is true | f2==i2 is false |
| B | f2-i2 = (float)0.000000, which is zero | f2-i2 = (float)1.000000, which is NOT zero |
| C | x*x = (float)0.000000, which is zero | x*x = (float)0.000000, which is zero |
| D | x*x is 0.000000e+00, zero | x*x is 1.000000e-60, NOT zero |
A と B は、float より高い精度で比較・計算しているように見える。
呪われていると C が「zero」なのに D が「NOT zero」となる。
xx は float では表現できない 1e-60 であるにもかかわらず、 _Generic の show_float に飛んでいるので、 xx が float なのは間違いないらしい。
何が起こっているのか
起こっていることについては x87 FPUの呪い に全部書いてある。
要するに。
x87 で計算すると単精度はかえって高コストなので、全部 long double で計算したい。C / C++ にはそれを許す例外規定みたいなものがあって、x87 を使う場合はその規定が適用される。ということみたい。
なぜ困るのか
誤差のある入力に対して適切な誤差を含む計算結果がほしい、という用途なら全然困らない。
計算間違いするわけじゃないので。
こまるのは、数値実験。
スレッドのような不確定要素のないプロセスで、同じソースコード、同じ乱数、同じ入力で計算をしたら同じ結果になってほしい。
実際、x86-64 でも ARM でも PowerPC でも RISC-V でもちゃんと書けば同じ結果になるとおもうんだけど、呪われた世界はそうならない。クリーンな数値実験なのに再現性が疑われてしまう。困る。
呪われる条件
呪われた世界になる条件を調べてみた。
OS は Linux(debian)で、 amd64(64bit) と x86(32bit)。
コンパイラは、 gcc-10, g++-10, clang-12, clang++12。
コンパイラオプションは以下を試した。
-
-pedanticをあり、なし。 -
-fexcess-precisionは、なし、standard、fastの三種。 - FP の指定は、なし、
-mfpmath=sse -msse2、-mfpmath=387の三種。 -
-stdは、C 言語ではc89,c99,c11,c17。C++では、c++98,c++03,c++11,c++17。 - 最適化は
-O0と-O2
組み合わせによってはコンパイラがサポートしていなかったので、できる範囲で。
以下のコードを動かし、 true が出るか false が出るかを調べた。
#include <stdint.h>
#include <stdio.h>
int main()
{
float f2 = 99999999.0f;
int32_t i2 = 99999999;
printf( "f2==i2 is %s\n", f2==i2 ? "true" : "false" ); /* [A] */
return 0;
}
下表で、true が ✅、false が 💀 である。
-pedantic は結果に影響を与えなかったので表から外した。
下表は、1つのマスに -O0 の結果（左）と -O2 の結果（右）を入れてある。
つまり、同じマスに💀と✅が両方入っているマスは最適化オプションで結果が変わることを意味する。
32bit / gcc-10
| fpmath | exess_prec | c89 | c99 | c11 | c17 |
|---|---|---|---|---|---|
| - | - | ✅ / ✅ | ✅ / ✅ | 💀 / 💀 | 💀 / 💀 |
| - | fast | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
| - | standard | ✅ / ✅ | ✅ / ✅ | 💀 / 💀 | 💀 / 💀 |
| sse | - | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| sse | fast | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| sse | standard | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| 387 | - | ✅ / ✅ | ✅ / ✅ | 💀 / 💀 | 💀 / 💀 |
| 387 | fast | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
| 387 | standard | ✅ / ✅ | ✅ / ✅ | 💀 / 💀 | 💀 / 💀 |
32bit / clang-12
| fpmath | exess_prec | c89 | c99 | c11 | c17 |
|---|---|---|---|---|---|
| - | - | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
| sse | - | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| 387 | - | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
32bit / g++-10
| fpmath | exess_prec | c++98 | c++03 | c++11 | c++17 |
|---|---|---|---|---|---|
| - | - | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
| - | fast | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
| sse | - | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| sse | fast | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| 387 | - | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
| 387 | fast | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
32bit / clang++-12
| fpmath | exess_prec | c++98 | c++03 | c++11 | c++17 |
|---|---|---|---|---|---|
| - | - | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
| sse | - | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| 387 | - | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
64bit / gcc-10
| fpmath | exess_prec | c89 | c99 | c11 | c17 |
|---|---|---|---|---|---|
| - | - | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| - | fast | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| - | standard | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| sse | - | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| sse | fast | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| sse | standard | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| 387 | - | ✅ / ✅ | ✅ / ✅ | 💀 / 💀 | 💀 / 💀 |
| 387 | fast | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
| 387 | standard | ✅ / ✅ | ✅ / ✅ | 💀 / 💀 | 💀 / 💀 |
64bit / clang-12
| fpmath | exess_prec | c89 | c99 | c11 | c17 |
|---|---|---|---|---|---|
| - | - | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| sse | - | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
64bit / g++-10
| fpmath | exess_prec | c++98 | c++03 | c++11 | c++17 |
|---|---|---|---|---|---|
| - | - | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| - | fast | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| sse | - | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| sse | fast | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| 387 | - | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
| 387 | fast | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ | 💀 / ✅ |
64bit / clang++-12
| fpmath | exess_prec | c++98 | c++03 | c++11 | c++17 |
|---|---|---|---|---|---|
| - | - | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
| sse | - | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ | ✅ / ✅ |
まとめ
- いずれの環境も
-mfpmath=sse -msse2であれば平和な世界が約束される。 - 32bit 環境だと
-mfpmathの省略は呪われた世界を意味する。 - c99 と c11 の間に断絶があるのは
gccだけで、clangは影響がない。 -
-fexcess-precision=standardはあまり役に立たない。 - 呪われた世界では、最適化を変えると
==の結果が変わることがある。おそろしい。
結局のところ
結局のところ、c99 と c11 の間に断絶がある理由はわからないんだけど、32bit の clang の動きを鑑みるに
-
-std=c11をリリースする前に、呪われた世界も C 言語として正しいということに gcc 開発者チームが気づいた。 - そこで、デフォルトの動作を、より速い呪われた世界にしたいと思った。
- しかし、急に変えるとひどい非互換を入れることになるので「
std=c11以降では呪われる」という選択をした。
辺りではないかと想像する。想像しているだけで証拠はない。