注意：この記事は懐古趣味、歴史的興味の側面が強い。現代の環境に当てはまるかどうかは、本文を参照して判断されたい。

C言語やJavaなどのプログラミング言語では、float と double など、精度の異なる複数の浮動小数点数型を提供している。
プログラムを実行するCPUでも何らかの方法で複数の浮動小数点型に対応する必要があるわけだが、その際のやり方として2通りが考えられる：

「大は小を兼ねる」方式：高精度な浮動小数点数型用の演算器を一個用意し、 double にも float にもその演算器（一種類の演算命令）を使う方式。
「ぴったりサイズ」方式：それぞれの型に応じた演算器（演算命令）を用意する方式。

Intelのx87 FPUやモトローラのFPUは前者を採用した。それ以外のアーキテクチャーは後者を採用するものが多い。

x87 FPUは「高精度な浮動小数点型」として精度64ビット、指数部15ビット（全体の幅が80ビット）のフォーマットを採用し、float も double もそれを使って計算することにした。このことがこの記事で解説する種々の問題、この記事のタイトルで言う「呪い」を生むことになる。

x86_64はこの呪いから解放されているので、面倒に首を突っ込みたくない人は32ビット環境を捨てて64ビット環境のみをターゲットとするべきである。

（残念なことに、MinGWのインストール方法でググると（-w64ではない）「32ビット用のMinGW」の導入方法を紹介するページが複数ヒットする。そういうページは現代においてはもはや有害なのでさっさと削除されるべきである。）

x87 FPU

x87 FPUのレジスターは80ビットの幅があり、符号1ビット、指数部15ビット、仮数部64ビットから成る二進浮動小数点数を表現できる。

x87 FPUの挙動はControl Wordと呼ばれる16ビットのレジスターをいじることによってある程度変更できる。16ビットの内訳は次の通りである：

名前	Reserved	X	RC	PC	Reserved	PM	UM	OM	ZM	DM	IM
ビット	15 / 14 / 13	12	11 / 10	9 / 8	7 / 6	5	4	3	2	1	0

それぞれの説明はこんな感じである：

X: Infinity Control
- 無限大の符号の扱いを変える設定だったらしい。IEEE 754が標準化される前の痕跡で、387以降は意味を持たない。
RC: Rounding Control
- 丸め方法を指定する。
- 00B: Round to nearest (even)
- 01B: Round down (toward $-\infty$)
- 10B: Round up (toward $+\infty$)
- 11B: Round toward zero (Truncate)
PC: Precision Control
- 演算結果の仮数部の精度を指定する。
- 00B: Signle Precision (24 bits)
- 01B: Reserved
- 10B: Double Precision (53 bits)
- 11B: Double Extended Precision (64 bits)
Exception Mask Bits
- 浮動小数点例外がトラップするかどうかを指定する。例外に対応するマスクビットが立っていたらトラップせずに実行を継続し、マスクビットがクリアされていたらトラップする。
- PM: Precision Mask
- UM: Underflow Mask
- OM: Overflow Mask
- ZM: Zero Divide Mask
- DM: Denormal Operand Mask
  - これはIEEE 754では規定されていない種類の例外で、入力に非正規化数が含まれていたことを表す。
- IM: Invalid Operation Mask

Rounding ControlとException Mask Bitsは他のFPU（x86のSSEや、ARMのFPU）にも見られる、ありがちなフィールドである。一方、Infinity ControlとPrecision Controlはx87に特徴的である。

この記事の話で重要なのはPrecision Control (PC)で、演算結果の仮数部の精度を3種類から選択できる。Precision Controlを変えても指数部の幅は15ビットで変わらない。

ちなみに、Cで書かれたプログラムからFPUの設定をいじる方法としては、C99で規定されたいくつかの関数（fesetround, feholdexceptなど）と、それ以外の環境依存な方法がある。後者の例として、MSVCでは _control87, _controlfp, __control87_2 などの関数でControl Wordをいじることができる。

_control87, _controlfp, __control87_2 | Microsoft Docs

SSEとの関係

昔のx86系プロセッサーでは浮動小数点演算を行いたかったらx87 FPUを使うかソフトウェアエミュレーションを使うしかなかったが、Pentium IIIで変化が訪れた。

Pentium IIIでは、SIMDのために、x87 FPUとは別系統の浮動小数点命令 (SSE; Streaming SIMD Extensions) が追加された。SSEはベクトル演算 (SIMD) だけではなくスカラー演算（通常の浮動小数点演算）にも使うことができる。この記事では、SSEについてはもっぱらスカラー演算としての側面を取り扱う。

SSE系の命令はプロセッサーの世代を追うごとに拡張されている。スカラー演算に関して言うと、重要なポイントは次の3つである：

無印SSE (Pentium III以降): 単精度 (float) 演算ができるようになった。
SSE2 (Pentium 4, Xeon以降): 倍精度 (double) 演算ができるようになった。
- 今日のパソコン向けのx86 CPUには必ずと言っていいほどSSE2が搭載されているはずである。
  パソコン向けに限らなければ、数年前までIntel QuarkというSSE非搭載のプロセッサーが製造されていたらしい。
x86_64はSSE2を必須としている。つまり、
- 64ビットに対応するx86系プロセッサーはSSE2を実装しなければならないし、
- 64ビット向けのx86コードは何の断りもなくSSE2命令を使うことができる。
- x86_64向けにプログラムをコンパイルすると、float や double の演算にはSSE2が使われる。

というわけで、プログラムを64ビット環境（x86_64）向けにコンパイルする場合は何の問題もない。問題は、ターゲットが32ビット環境の場合である。

現代のx86 CPUではほぼ必ずSSE2が使えるとはいえ、32ビット環境をターゲットとする場合に無条件でSSEを使ってしまうと、（20年以上前の）古いプロセッサーでプログラムが動作しなくなる。そのため、保守的なコンパイラー（GCC等）ではターゲットが32ビット環境の場合にはSSEを使わず、従来のx87 FPUを使うコードを出力する（詳しくは後述）。

一方、いくつかのコンパイラーでは32ビットがターゲットであってもデフォルトでSSE2を使ったコードを出力する。具体的に筆者が知っているものをいくつか挙げる：

今のMSVCは32ビットコードを出力する際にデフォルトでSSE2を使う。コンパイルオプション /arch:IA32 を明示的に指定するとx87 FPUを使うようになる（詳しくは後述）。
HaskellコンパイラーのGHCも最近のバージョン (8.10.1以降) では常にSSE2を使うようになった（x87 FPUを使うモードを削除した）。
- Release notes for version 8.10.1より：(x86) Native code generator support for legacy x87 floating point coprocessor has been removed. From this point forth GHC will only support floating point via SSE2.
Goも1.16でx87 FPUのサポートを削除するらしい。
- Go 1.16 Release Notes - The Go Programming Language

x87 FPUの呪い

float や double の演算にSSE2ではなくx87 FPUを使うと具体的にどういう問題が発生するのか解説する。

精度をどうするか

x87 FPUでは演算の精度を変更できるとは言っても、演算のたびに精度を変更するのは煩雑である。現実的なのは、「プログラムの最初にx87 FPUの精度を決めて、以後ずっとそれを使う」である。3種類の精度のうち、どれを使うべきか考えてみよう。

PC=24：普通は使われない。…が、昔のDirectXでは float の演算を高速化するため、初期化時に精度を24ビットにするという挙動をしていたらしい。
PC=53：double の挙動をIEEE 754準拠に近づけたい場合はこれを使う。その場合 long double も double 相当の精度となる。
PC=64：幅80ビットの long double をフル活用したい場合はこれを使うしかない。代わりに、一部の入力に対して double 演算が少し怪しくなる。→「一つ目の呪い」

なお、x87 FPUの精度を変えても指数部の範囲は変わらないままである。なので、Precision Controlを24ビットや53ビットとしてもIEEEのbinary32やbinary64とは同一の挙動にはならない。したがって、PC=24やPC=53の状態を「幅32ビット」「幅64ビット」と称するのは間違いである。

一つ目の呪い：PC=64と倍精度演算

PC=64の設定下では、double の演算が本来の精度（仮数部53ビット）よりも高い精度で行われる。演算途中の精度が高くなるのは良いことのように思えるが、実際にはそうとは限らない。

例として、 x = 0x1.00002fff0p0 と y = 0x1.000000008p0 の積を考えよう。x * y の正確な値は 0x1.00002fff800017ff8p0 で、それを正しい丸めで double に変換すると 0x1.00002fff80001p0 となる。

しかし、PC=64なx87 FPUでこれを計算すると内部的には一旦64ビット精度に丸められ、その結果は 0x1.00002fff80002000p0 となる。さらにそれを double に丸めると 0x1.00002fff80002p0 となる。

コードで書くと、

# include <stdio.h>
int main(void)
{
    volatile double x = 0x1.00002fff0p0, y = 0x1.000000008p0;
    // x * yの正確な値は0x1.00002fff800017ff8p0で
    // それをdoubleへ正しく丸めると0x1.00002fff80001p0だが…
    double z = x * y;
    printf("x * y = %a\n", z);
}

このコードを32ビットのx86環境をターゲットとするGCCで（デフォルトの設定で）コンパイルすると、実行結果は

x * y = 0x1.00002fff80002p+0

となる。

この辺の話は、同人誌「浮動小数点数小話」に「二段階丸め」(double rounding)の問題として書いた。

要するに、x87 FPUでPC=64を設定すると、他の環境と演算結果が一致しなくなる（再現性が損なわれる）。

じゃあPC=53なら良いのか、という問題については次の節で扱う。

二つ目の呪い：指数部の範囲

x87 FPUのPrecision Controlフィールドをいじっても、変わるのは仮数部の精度だけで、指数部の精度は広いままである。

演算途中の指数部の精度が広いことにはメリットもある。例えば、「指数部のオーバーフローが要注意な関数」を素直に実装できる。hypot なら

double hypot(double x, double y)
{
    return sqrt(x * x + y * y);
}

となるだろう。

（hypot の「素直な実装」の問題を知らない方向けの説明：hypot(0x3p1000, 0x4p1000) は 0x5p1000 を返すべきだが、普通の環境で hypot(x, y) を sqrt(x * x + y * y) で実装してしまうと x * x = 0x3p1000 * 0x3p1000 = 0x9p2000 を倍精度で表せずに無限大が返ってしまう。）

この手の、指数部の範囲が広いことによるメリットがある処理には、複素数の割り算もある。詳しくは、筆者のブログ記事を参照されたい：

指数部の範囲が広いことのデメリットとしては、オーバーフローあるいはアンダーフローが起こるような場合に演算結果が他の環境と整合しなくなる（再現性が損なわれる）ことがある。

演算の度に（単精度あるいは倍精度として）メモリに書き出せば指数部の範囲に応じて適切にオーバーフロー・アンダーフローするようになるが、当然実行速度は遅くなる。

ただ、「演算の途中結果を都度メモリに書き出す」という涙ぐましい努力をしても、まだ問題が残っている。それは、アンダーフローの際の「二段階丸め」だ。

筆者の同人誌にも書いたが、浮動小数点数のアンダーフローは「結果が0になる」場合だけでなく「結果が非正規化数となる」場合も含む。その際に丸めが起こる。

例として、 x = 0x1.fffe0effffffep-51 と y = 0x1.0000000000001p-1000 の積を考えよう。「呪われていない」環境でこれを倍精度で計算すると

正確な演算結果は 0x1.fffe 0eff ffff ffff e0ef ffff fep-1051 である。
これを倍精度浮動小数点数に変換すると、結果は非正規化数となり、仮数部の精度は24ビットとなる。仮数部を24ビットに丸めると 0x1.fffe 0ep-1051 となる。

により、結果は 0x1.fffe 0ep-1051 となる。

一方、これをPC=53なx87 FPUで計算すると、

正確な演算結果は 0x1.fffe 0eff ffff ffff e0ef ffff fep-1051 である。
浮動小数点数レジスターの指数部の範囲は十分広いので、アンダーフローは起こらない。仮数部を53ビットで丸めて 0x1.fffe 0f00 0000 0p-1051 がFPUレジスターに格納される。
メモリに格納する際にアンダーフローが起こり、仮数部が24ビットに丸められる。0x1.fffe 1p-1051 がメモリに格納される。

となって、結果は 0x1.fffe 1p-1051 となる。

この挙動を回避してポータブルな結果を得るにはx87 FPUのフラグをいじるだけではどうにもならず、「結果が非正規化数だったら浮動小数点演算のソフトウェアエミュレーションを利用する（あるいは指数部をゴニョゴニョした上でx87 FPUを使う）」というような対策が必要となる。

二重の呪い

というわけで、x87 FPUは「（PC=64の場合）中途半端な精度」と「指数部の範囲」の二重に呪われており、倍精度演算でIEEE準拠の結果を得るのはなかなか厄介な問題となる。

幸い、2020年に生きる我々は、「SSE2が実装されていないx86系CPUのサポートを切り捨てる」という選択によりこの問題から逃れることができる。

C言語の規定： `FLT_EVAL_METHOD`

C言語は懐が広いので、x87 FPUのような厄介な挙動も許容している。つまり、浮動小数点演算の途中結果を名目上の型とは異なる形式、例えば long double 相当のフォーマットで格納することを許している。

C言語では、浮動小数点数の演算がどのような方式で行われるかプログラム側に教えるために、定数 FLT_EVAL_METHOD と型 float_t, double_t を提供している。

定数 FLT_EVAL_METHOD は浮動小数点演算の途中結果がどういう形式で格納されるかを表す。

-1: 不定
0: 全ての演算や定数を、名目上の型の範囲と精度で計算する。
1: float と double の演算と定数を double の範囲と精度で計算する。long double は long double 自身の範囲と精度で計算する。
2: 全ての浮動小数点演算と定数を、 long double の範囲と精度で計算する。

x86_64やAArch64など、float と double にそれぞれの演算命令が用意されている環境では FLT_EVAL_METHOD は 0 となる。一方、x87 FPUを使用せざるを得ない環境では FLT_EVAL_METHOD は 2 となる。

float_t と double_t は、それぞれ float, double と同じかそれ以上の幅を持つと規定されており、

FLT_EVAL_METHOD == 0: float_t は float で、 double_t は double である。
FLT_EVAL_METHOD == 1: float_t と double_t のいずれも double である。
FLT_EVAL_METHOD == 2: float_t と double_t のいずれも long double である。
それ以外：実装依存

となる。ざっくり言うと「float 型の演算の途中結果は float_t で、 double 型の演算の途中結果は double_t で行われる」と思って良さそうだが、FLT_EVAL_METHOD が 0, 1, 2 のいずれでもない場合は規定がない。

なお、C言語は、「浮動小数点演算の中間結果を変数に代入する（あるいはキャストする）か否かで結果が変わる」ことを許している。例えば、

double three_mul_1(double x, double y, double z)
{
    return x * y * z;
}

と

double three_mul_2(double x, double y, double z)
{
    double xy = x * y;
    return xy * z;
}

double three_mul_3(double x, double y, double z)
{
    return (double)(x * y) * z;
}

は異なる結果を返す可能性がある。

Javaの `strictfp`

Javaには long double がないので、実行時のx87 FPUの演算精度はPC=53として良い。float の演算は都度メモリに格納すれば良いだろう。というわけで「呪い」の一つ目の影響は受けない。

Javaで問題となるのは「呪い」の二つ目、指数部の範囲が広いことによる弊害である。浮動小数点演算の結果が環境に依存しないことを保証するには、x87 FPUを使う環境では浮動小数点数乗算の際に毎回アンダーフローの検査を入れる必要がある。

この問題のため、Javaでは strictfp という修飾子が導入された。strictfp が有効なコードでは浮動小数点演算の結果が環境に依存しないことが保証される代わりに、SSE2以前のx86で効率が落ちる。strictfp が有効でないコードではSSE2以前のx86でも浮動小数点演算を高速に行うことができるが、浮動小数点演算の結果が環境に依存するようになる。

strictfp の有無により挙動が変わるコードの実例はこちらを参照されたい：

Javaの strictfp が実際に意味を持つ環境を用意する - Qiita

x86_64を含む現代のアーキテクチャーでは strictfp があってもなくても同じ（IEEE準拠の）挙動となる。なので、 strictfp があってもなくてもIEEE準拠の挙動としよう、という提案が出ている。

[JDK-8175916] JEP 306: Restore Always-Strict Floating-Point Semantics - Java Bug System

コンパイラーの対応

C言語の FLT_EVAL_METHOD 周りを検証できるプログラムを書いた：

floating-point-test/eval-method.c at master · minoki/floating-point-test

MSVC

MSVCでは long double は double と同じ幅・精度（つまり倍精度）である。

32ビットのx86をターゲットとする場合、デフォルトでSSE2が仮定される（/arch:SSE2 相当）。

/arch:IA32 を指定するとSSE2を仮定しないコード、つまりx87 FPUを使って演算を行うコードが出力される。

/arch (x86) | Microsoft Docs

MSVCは long double が80ビットではないので、実行時のFPUの設定はPC=53に設定される。

FLT_EVAL_METHOD は /arch:IA32 または /arch:SSE の時に 2 で、/arch:SSE2 以上の時に 0 となる。

GCC

x86 Options (Using the GNU Compiler Collection (GCC))

32ビットなx86をターゲットとするGCCがデフォルトでx87 FPUを使って演算することはすでに書いたが、この設定は -mfpmath オプションにより上書きできる。

-mfpmath=387|sse|both

-mfpmath= に　387 を指定するとx87のみを使い、 sse を使うとSSEを使うようになる。ただ、実際に float や double の演算をSSEを使って行うには、 -march= や -msse, -msse2 なども合わせて指定する必要がある。

というわけで、 float や double の演算にSSE2を使いたい場合、典型的なオプションの指定は -mfpmath=sse -msse2 となるだろう。

ちなみに、 -mfpmath=sse -msse を指定すると「float にはSSEを使うが double にはx87 FPUを使う」という中途半端な状態となる。この場合の FLT_EVAL_METHOD は -1（不定）となる。

x87 FPU使用時の設定を変えるオプションもいくつかあるので軽く紹介する。

-mpc32, -mpc64, -mpc80

これらはPrecision Controlの初期設定を変えるオプションのようだ。

Precision Controlをいじっても指数部の範囲の関係で「幅32ビット」「幅64ビット」のフォーマットとは異なるものが得られるので、このオプション名はおかしい。

最適化オプションについても、x87 FPU特有のものがある。

Optimize Options (Using the GNU Compiler Collection (GCC))
-ffloat-store

浮動小数点数の変数をFPUレジスターに格納しないようになる。つまり、変数への代入の際に必ずその型への丸めが起こるようになる。

-fexcess-precision=standard との違いは、キャストに関しては -ffloat-store は関知しないことのようだ。

-fexcess-precision=fast|standard

-fexcess-precision=standard が有効な場合、浮動小数点数の過剰な精度についてはC標準に従う。つまり、キャストや代入の際に過剰な精度は削ぎ落とされる。

-fexcess-precision=fast が有効な場合、x87 FPU使用時に

double three_mul_1(double x, double y, double z)
{
    return x * y * z;
}
double three_mul_2(double x, double y, double z)
{
    double xy = x * y;
    return xy * z;
}
double three_mul_3(double x, double y, double z)
{
    return (double)(x * y) * z;
}

/*
three_mul_1(0x1p+1000, 0x1p+1000, 0x1p-1000) => 0x1p+1000; C標準で許可された動作
three_mul_2(0x1p+1000, 0x1p+1000, 0x1p-1000) => 0x1p+1000; 過剰な精度で計算されている！
three_mul_3(0x1p+1000, 0x1p+1000, 0x1p-1000) => 0x1p+1000; 過剰な精度で計算されている！
*/

となる可能性がある。

このオプションの初期値は -std=c99 や -ffast-math の影響を受ける。

まとめ

浮動小数点演算をちゃんとIEEE準拠にしたいなら、32ビットのx86をターゲットとするのを今すぐやめろ。　

もしどうしても32ビットをターゲットとする必要がある場合は、SSE2を有効にしろ。