浮動小数点数の足し算と掛け算は可換か

読むのが面倒な人向けの結論：**可換です**。

「可換です」以外の答えを知りたい人はこの記事を最後まで読んでください。

# 結合法則と交換法則

実数の足し算や掛け算については結合法則 $x+(y+z)=(x+y)+z$ が成り立ちます。これに対し、浮動小数点数の足し算・掛け算が結合的でないことはとても有名な話です。

例えば、倍精度で `(0x1p-200 + 1) + (-1)` を計算すると、結合法則が成り立てば答えは `0x1p-200` となるはずですが、実際には `0` が返ってきます。

浮動小数点演算が結合的でないことは有名な話なので、ここではこれ以上詳しくは取りあげません。

一方で、交換法則（可換性）はどうでしょうか？「浮動小数点演算はこういう法則を満たさない！クソ！」みたいな話題で槍玉に上がるのはほとんどの場合結合法則で、交換法則に言及するものはあまり見かけない気がします。

## 交換法則が成り立つとどういう場面で嬉しいか

結合法則は並列化とかでガッツリ役に立ちますが、交換法則はそれに比べて地味に思えます。しかし、役に立たない場面がないわけではありません。

例えば、ベクトルとスカラーの積を取る関数を作っているとしましょう。

```typescript
// 配列 xs で表されるベクトルを x 倍する
function scale(x: number, xs: number[]): number[]
{
    let ys: number[] = [];
    for (let i = 0; i < xs.length; ++i) {
        ys[i] = x * xs[i]
    }
    return ys;
}
```

仮に交換法則が成り立たなかったら、「左からスカラー倍する」`scaleLeft` と「右からスカラー倍する」`scaleRight` という関数を2種類用意しなければなりません。

一方で、交換法則が成り立てば、用意する関数は一つだけで済むか、あるいは `scaleRight(xs, x) := scaleLeft(x, xs)` という風に一方の実装をもう一方で使いまわせて、経済的です。

メリットが地味ですか？地味ですね。

もっと地味な例を挙げると、コンパイラーの最適化（共通部分式除去, CSE）で利用できるかもしれない、というのがあります。例えば、次のようなCコードがあった時に、足し算が可換であれば `x` と `y` の足し算を1回で済ませることができます。

```c
void func(double x, double y)
{
    double u = x + y;
    double v = y + x;
    ...
}
```

# 浮動小数点数が表している物は何か

浮動小数点演算を考える前に、そもそも浮動小数点数とは何か確認しておきましょう。

有限な（無限大やNaNではない）浮動小数点数は、数直線上の点を表しています。つまり実数です。

注意して欲しいのは、メモリやレジスタ上のビット列として表される浮動小数点数それ自体には、「誤差」という情報は含まれません。各種数学関数は与えられた浮動小数点数が無限に正確な物として取り扱います。「誤差」の概念があるとしたら、コードを書くプログラマーの頭の中です。

有限な浮動小数点数のなす集合は、実数の離散的な部分集合と考えることができます。この集合はもちろん、浮動小数点数のフォーマット（基数、精度と、指数部の範囲）によって異なります。

最もポピュラーなフォーマットである倍精度浮動小数点数について、有限なものを集合として書き下せば、

```math
\mathrm{Float64}_{\mathit{fin}}=
\begin{gathered}\{-(2^{1024}-2^{971}),-(2^{1024}-2\times 2^{971}),\ldots,-(1+2^{52}),-1,-(1-2^{53}),\ldots,\\
\quad-2^{-1074},0,2^{-1074},\ldots,1-2^{53},1,1+2^{52},\ldots,2^{1024}-2\times 2^{971},2^{1024}-2^{971}\}
\end{gathered}
```

となるでしょう。

注意深い読者は「ゼロの符号を考慮していない」ことに気が付かれたかもしれませんが、ゼロの符号については後で考えることにします。[^many-to-one]

[^many-to-one]: 複数の浮動小数点数が同じ実数を表す状況は、符号付きのゼロ以外にもあり得ます。具体的には、10進浮動小数点数です。この節では、「浮動小数点数を表す集合」は「同じ実数を表す浮動小数点数」を同一視した、実数の部分集合として表せる物を表します。

さて、浮動小数点数演算は通常の実数の演算を基に定義されます。演算結果の実数が浮動小数点数のなす集合の元であれば良いのですが、残念ながらそうとは限りません。その場合、何らかの**丸め**を行って代わりとなる浮動小数点数を返します。

丸めのやり方はいろいろありますが、IEEE 754で規定されている物には

* 最近接偶数丸め (roundTiesToEven)
* 「四捨五入」[^ties-to-away] (roundTiesToAway)
* 正の無限大方向への丸め（切り上げ） (roundTowardPositive)
* 負の無限大方向への丸め（切り下げ） (roundTowardNegative)
* ゼロ方向への丸め（切り捨て） (roundTowardZero)

[^ties-to-away]: 2進法の場合「四捨五入」というのは変な言い方です（なのでカッコ書きしました）が、他に適当な言い方を思いつかないのでこうしておきます。

があります。場面によってはここに挙げた方法以外で丸めが行われることがあります。

これらの方法によって丸めを行う関数を $R\colon \mathbf{R}\rightarrow \mathrm{Float}$ で書くことにすると、浮動小数点演算は

```math
\begin{aligned}
x\oplus y&:=R(x+y), \\
x\ominus y&:=R(x-y), \\
x\otimes y&:=R(x\times y), \\
x\oslash y&:=R(x/y)
\end{aligned}
```

（ただし$x$, $y$は有限な浮動小数点数）

という風に定義できます。ただし、左辺の丸で囲まれた演算は浮動小数点演算を表し、右辺の丸で囲まれていない演算は実数としての正確な演算を表します。

この定義を見れば、**浮動小数点数の足し算と掛け算について交換法則が成り立つのは自然なことだと言えます**。式変形として書けば

```math
    \begin{aligned}
     x\oplus y&=R(x+y) & &(\text{演算の定義}) \\
     &=R(y+x) & &(\text{実数についての演算の可換性}) \\
     &=y\oplus x & &(\text{演算の定義})
    \end{aligned}
```

となります。

ですが、浮動小数点数については「有限な実数」以外の諸々が含まれるので、可換か否かの結論を出すにはもう少し議論が必要です。

# 符号付きのゼロ

IEEE 754で規定される浮動小数点数には**正負2種類のゼロ**があります。これらは通常の比較演算（C系の言語の `==` など）では等価として扱われますが、文字列化、割り算（の分母）、atan2などの一部の状況では区別されます。

符号付きのゼロが絡んだ場合に交換法則が成り立つか、ですが、その前に符号付きのゼロに関する浮動小数点数の足し算と掛け算の定義を確認しておきましょう。

まずは足し算です。話を簡単にするために、丸め方法は最近接偶数丸めとします。

```math
    \begin{aligned}
    \mathtt{+0}\oplus\mathtt{+0}&=\mathtt{+0}, \\
    \mathtt{+0}\oplus\mathtt{-0}&=\mathtt{+0}, \\
    \mathtt{-0}\oplus\mathtt{+0}&=\mathtt{+0}, \\
    \mathtt{-0}\oplus\mathtt{-0}&=\mathtt{-0}, \\
    x\oplus (-x)&=\mathtt{+0} & (x\ne 0)
    \end{aligned}
```

足し算は、同符号のゼロを足した場合は同じ符号のゼロが返ります。それ以外の場合は正のゼロが返ります[^zero-sum]。

[^zero-sum]: 丸め方法が負の無限大方向への丸め (roundTowardNegative) の場合は、負のゼロが返ります。

掛け算は単純で、積の符号はオペランドの符号のXORです。例としてゼロどうしの積がどうなるかを挙げてみましょう。

```math
    \begin{aligned}
    \mathtt{+0}\otimes\mathtt{+0}&=\mathtt{+0}, \\
    \mathtt{+0}\otimes\mathtt{-0}&=\mathtt{-0}, \\
    \mathtt{-0}\otimes\mathtt{+0}&=\mathtt{-0}, \\
    \mathtt{-0}\otimes\mathtt{-0}&=\mathtt{+0}
    \end{aligned}
```

浮動小数点数の掛け算は実数と違って、ゼロでない二数の積がゼロになることがあります（アンダーフロー）。その場合も符号はオペランドのそれのXORです。

交換法則、という点に着目して今書いた定義を確認すると、浮動小数点数の足し算、掛け算のいずれも、**ゼロの符号を考慮しても交換法則が成り立ちます**。

# 無限大

IEEE 754で規定される浮動小数点数には正負2種類の**無限大**も含まれます。これらはオーバーフローの結果やゼロ除算の結果として生成されます。

無限大に関する浮動小数点演算、特に足し算と掛け算は、

```math
    (\pm\infty) \oplus y =\begin{cases}
    \pm\infty & \text{$y$が有限な場合} \\
    \pm\infty & \text{$y$が同符号の無限大な場合} \\
    \mathrm{NaN} & \text{$y$が異符号の無限大またはNaNの場合}
    \end{cases}
```

```math
    (\pm\infty) \otimes y =\begin{cases}
    \pm\infty & \text{$y$が非0で有限または無限大な場合。符号は$y$に依存する} \\
    \mathrm{NaN} & \text{$y$が0またはNaNの場合}
    \end{cases}
```

という風に定義されます。$y$の方が無限大な場合は交換法則が成り立つように定義します。

……はい。というわけで、浮動小数点数の足し算と掛け算について、**無限大を考慮しても交換法則が成り立ちます**。

# NaN

浮動小数点数のとりうる値としては、**NaN** (not a number; 訳すなら**非数**)というのもあり得ます。典型的には、浮動小数点演算の中で「不定」が発生した時、あるいは数学関数の入力として定義域の外の入力が与えられた時にNaNが返ってきます。それから、入力のいずれかがNaNの場合も原則としてNaNが返ります（NaNは伝播する）。

浮動小数点数で計算するとNaNが返ってくる例は$\infty-\infty$, $0\times\infty$, $0/0$, $\infty/\infty$, $\sqrt{-1}$などです。

NaNを「一つの値」だと思えば、足し算と掛け算について

```math
    \begin{gathered}
    x\oplus\mathrm{NaN}=\mathrm{NaN} \oplus x =\mathrm{NaN}, \\
    x\otimes\mathrm{NaN}=\mathrm{NaN} \otimes x =\mathrm{NaN}
    \end{gathered}
```

が成り立ちます。**つまり、NaN同士を区別しないのならNaNを考慮しても交換法則が成り立ちます**。

## NaN同士の区別

実は、「NaN」と呼ばれる値は一つだけではありません。IEEE 754で規定されるNaNには、

* 符号ビット
* signaling / quiet の区別
* ペイロード 

などの情報が含まれます。

このうち、「signaling / quietの区別」は、演算結果としては常にquiet化されたNaNが返ってくるので交換法則には関係ありません。関係するのは符号ビットとペイロードです。

浮動小数点数の四則演算についてNaNの符号とペイロードがどうなるかですが、IEEE 754では

* 浮動小数点演算の結果がNaNの場合、その符号ビットは規定しない（入力に含まれるNaNがただ一つだったとしても） (6.3)
    * 特に、 `NaN * (-1)` が符号を変えるような実装と、`NaN * (-1)` の符号が変わらないような実装の両方が許容されます。
* ペイロードに関して、
    * 入力の中にNaNがただ一つの場合、結果のNaNのペイロードは入力に含まれるNaNの物と同一であるべき (6.2.3)
    * 入力の中にNaNが複数含まれる場合、結果のNaNのペイロードは入力のNaNのペイロードのいずれかと同一であるべき (6.2.3)

となっています。つまり、**具体的なことは規定されていません**。

NaNの符号とペイロードに関して浮動小数点数の足し算と掛け算が具体的にどう振る舞うかは、浮動小数点数の実装に依存します。ここではポピュラーな実装（CPUの命令セットアーキテクチャ）をいくつか見てみます。

x86系の場合、入力の双方がNaNだった場合は

* x87 FPU: ペイロードが大きい方。
    * 符号ビットに関しては明記されていないが、足し算と掛け算について試した感じではオペランドの符号ビットのANDを取っている？
* SSE/SSE2/SSE3/SSE4.1/AVX: 最初のNaN

となっています（参照：Intel SDM Volume 1, 4.8.3.5）。x87 FPUの足し算と掛け算はNaNの符号ビットとペイロードを考慮しても可換だと言えそうです。

ARM (AArch64)の場合は、通常は引数の中で最初のNaNが返されますが、FPCRのDN (Default NaN)ビットが設定されている場合はDefault NaN（符号ビットとペイロードが0なquiet NaN）が返されます。

（参照：Arm Architecture Reference Manual  Armv8, for Armv8-A architucture profile (A1.5.5, C5.2.7)）

結論としては「NaNの符号やペイロードまで考慮すると、実装によっては浮動小数点数の足し算や掛け算が可換とは限らない」となります。

## NaN同士を区別する意味があるのか

前節での結論は「可換とは限らない」でしたが、IEEE 754では入力の両方がNaNだった場合にどちらを返しても良いことになっているので、**可換性を使った最適化は許されます**。

例えば、最初の方で挙げた例

```c
void func(double x, double y)
{
    double u = x + y;
    double v = y + x;
    ...
}
```

では、 `x` と `y` の両方がNaNだった場合に

* `x + y` の計算では `x` を返しても良いし `y` を返しても良い
* `y + x` の計算では `y` を返しても良いし `x` を返しても良い

ので、足し算の計算を1回で済ませて `x + y` と `y + x` の両方が同じNaNを返す、というような最適化が可能です。

そもそもプログラミング言語によってはNaNの符号やペイロードを参照できないので、そういう言語では浮動小数点数の足し算と掛け算は厳密に可換です（例：ECMAScript）。

# まとめ

IEEE 754準拠の浮動小数点数の足し算と掛け算の可換性は、

* NaNを考慮しない範囲（通常の浮動小数点数、符号付きの0、無限大）では可換である。
* NaNを1つの値として扱った場合、可換である。
* NaNの持っている符号ビットやペイロードを考慮した場合、可換かどうかは実装依存となる。ただし可換性を使った最適化は行える。

となります。なんとも歯切れが悪い結論ですが、「**普通に扱う分には可換である。みんなが気にするような誤差は発生しない。**」と覚えておいてください。NaNの符号やペイロードを気にする必要があるのなんて多分CPUを作る人くらいです。

筆者がテストに使ったコードをここに置いておきます：[float-commutativity.c](https://gist.github.com/minoki/3f281f915b2421da2c414276d5496804)