NaNてことだ　SIMD化したら結果が変わるなんて

処理をSIMD化してたらまれに計算結果が一致しないケースがあったので調査した。

**注：SIMD化はあくまでこの問題を発見したきっかけであり、本稿のテーマとはさほど関係ありません。**

## 結論

公式規格書は一応読んだけどクソザコイングリッシュ故読み違いあるかも。間違ってたら指摘してください。

- NaNにはQuiet NaN (QNaN)とSignaling NaN (SNaN)がある
  - 仮数部は0以外何でも選べるし符号ビットも自由なのでビットパターンはそれぞれ$2^{52}$個くらいある。
    『エラーコードっぽい使い方できるよ！』とのことらしい
- IEEE 754では浮動小数点数に対するMin/Max演算を定義しているが、引数の一方がNaNであった場合の振る舞いが2008年版と2019年版で定義が異なる
  - 2008年版ではQNaNとSNaNで異なる動作をする
  - 2019年版では2通りのmin/maxが定義された

- そうは言ってもこの仕様に対して律儀に従うかどうかは完全に実装依存であり、使うハードウェアやライブラリを切り替えると結果が変わったりする
  - 交換法則さえ成り立たないのもあってつらい

自分で試した範囲のまとめを表にまとめた[^1][^2][^3]。

<table>
  <tr>
    <th colspan="2">&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;</th>
    <th>$\scriptsize{\mathrm{max(0   , QNaN)}}$</th>
    <th>$\scriptsize{\mathrm{max(0   , SNaN)}}$</th>
    <th>$\scriptsize{\mathrm{max(QNaN, 0   )}}$</th>
    <th>$\scriptsize{\mathrm{max(SNaN, 0   )}}$</th>
    <th>$\scriptsize{\mathrm{max(QNaN, QNaN)}}$</th>
    <th>$\scriptsize{\mathrm{max(QNaN, SNaN)}}$</th>
    <th>$\scriptsize{\mathrm{max(SNaN, QNaN)}}$</th>
    <th>$\scriptsize{\mathrm{max(SNaN, SNaN)}}$</th>
  </tr>
  <tr>
    <th rowspan="3">定<br/>義</th>
    <th>IEEE-754:2088<br/><code>maxNum</code></th>
    <th>$\mathrm{0}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{0}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
  </tr>
  <tr>
    <!---->
    <th>IEEE-754:2019<br/><code>maximum</code></th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
  </tr>
  <tr>
    <!---->
    <th>IEEE-754:2019<br/><code>maximumNumber</code></th>
    <th>$\mathrm{0}$</th>
    <th>$\mathrm{SNaN}$</th>
    <th>$\mathrm{0}$</th>
    <th>$\mathrm{SNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
  </tr>
  <tr>
    <th rowspan="2">H/W<br/>実<br/>装</th>
    <th>SSE2<br/><code>_mm_max_pd</code></th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
    <th>$\mathrm{0}$</th>
    <th>$\mathrm{0}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
  </tr>
  <tr>
    <th>AVX<br/><code>_mm256_max_pd</code></th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
    <th>$\mathrm{0}$</th>
    <th>$\mathrm{0}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
  </tr>
  <tr>
    <th rowspan="2">S/W<br/>実<br/>装</th>
    <th>C++<br/><code>std::max</code></th>
    <th>$\mathrm{0}$</th>
    <th>$\mathrm{0}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
  </tr>
  <tr>
    <th>C#<br/><code>Math.Max</code></th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{QNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
    <th>$\mathrm{SNaN}$</th>
  </tr>
</table>

他の環境の結果をご存知でしたら教えていただけると助かります。


[^1]: Intel Core i7-6700K CPU 4.00GHz (Skylake)
[^2]: C++: g++ 7.4.0。-O0から-O3まですべて同じ結果だった。
[^3]: C#: .Net Core 3.0。corert見に行ったら"IEEE-754:2019準拠するわ"っぽいコメントが書いてあったが、QNaNとSNaNを区別する気はないらしい。

## 補足：実験コード

<details><summary>C++</summary><div>

```C++:main.cpp
#include <algorithm>
#include <iostream>
#include <limits>
#include <math.h>
#include <emmintrin.h>
#include <immintrin.h>


std::string display(double x)
{
    if(issignaling(x))
        return "SNaN";
    if(isnan(x))
        return "QNaN";
    return std::to_string(x);
}


int main()
{
    constexpr auto qnan = std::numeric_limits<double>::quiet_NaN();
    constexpr auto snan = std::numeric_limits<double>::signaling_NaN();

    std::cout << "std : max(0.0 , QNaN) = " << display(std::max(0.0 , qnan)) << std::endl;
    std::cout << "std : max(0.0 , SNaN) = " << display(std::max(0.0 , snan)) << std::endl;
    std::cout << "std : max(QNaN, 0.0 ) = " << display(std::max(qnan, 0.0 )) << std::endl;
    std::cout << "std : max(SNaN, 0.0 ) = " << display(std::max(snan, 0.0 )) << std::endl;
    std::cout << "std : max(QNaN, QNaN) = " << display(std::max(qnan, qnan)) << std::endl;
    std::cout << "std : max(QNaN, SNaN) = " << display(std::max(qnan, snan)) << std::endl;
    std::cout << "std : max(SNaN, QNaN) = " << display(std::max(snan, qnan)) << std::endl;
    std::cout << "std : max(SNaN, SNaN) = " << display(std::max(snan, snan)) << std::endl;
    std::cout << std::endl;
    std::cout << "SSE2: max(0.0 , QNaN) = " << display(_mm_max_pd(_mm_set1_pd(0.0 ), _mm_set1_pd(qnan))[0]) << std::endl;
    std::cout << "SSE2: max(0.0 , SNaN) = " << display(_mm_max_pd(_mm_set1_pd(0.0 ), _mm_set1_pd(snan))[0]) << std::endl;
    std::cout << "SSE2: max(QNaN, 0.0 ) = " << display(_mm_max_pd(_mm_set1_pd(qnan), _mm_set1_pd(0.0 ))[0]) << std::endl;
    std::cout << "SSE2: max(SNaN, 0.0 ) = " << display(_mm_max_pd(_mm_set1_pd(snan), _mm_set1_pd(0.0 ))[0]) << std::endl;
    std::cout << "SSE2: max(QNaN, QNaN) = " << display(_mm_max_pd(_mm_set1_pd(qnan), _mm_set1_pd(qnan))[0]) << std::endl;
    std::cout << "SSE2: max(QNaN, SNaN) = " << display(_mm_max_pd(_mm_set1_pd(qnan), _mm_set1_pd(snan))[0]) << std::endl;
    std::cout << "SSE2: max(SNaN, QNaN) = " << display(_mm_max_pd(_mm_set1_pd(snan), _mm_set1_pd(qnan))[0]) << std::endl;
    std::cout << "SSE2: max(SNaN, SNaN) = " << display(_mm_max_pd(_mm_set1_pd(snan), _mm_set1_pd(snan))[0]) << std::endl;
    std::cout << std::endl;
    std::cout << "AVX : max(0.0 , QNaN) = " << display(_mm256_max_pd(_mm256_set1_pd(0.0 ), _mm256_set1_pd(qnan))[0]) << std::endl;
    std::cout << "AVX : max(0.0 , SNaN) = " << display(_mm256_max_pd(_mm256_set1_pd(0.0 ), _mm256_set1_pd(snan))[0]) << std::endl;
    std::cout << "AVX : max(QNaN, 0.0 ) = " << display(_mm256_max_pd(_mm256_set1_pd(qnan), _mm256_set1_pd(0.0 ))[0]) << std::endl;
    std::cout << "AVX : max(SNaN, 0.0 ) = " << display(_mm256_max_pd(_mm256_set1_pd(snan), _mm256_set1_pd(0.0 ))[0]) << std::endl;
    std::cout << "AVX : max(QNaN, QNaN) = " << display(_mm256_max_pd(_mm256_set1_pd(qnan), _mm256_set1_pd(qnan))[0]) << std::endl;
    std::cout << "AVX : max(QNaN, SNaN) = " << display(_mm256_max_pd(_mm256_set1_pd(qnan), _mm256_set1_pd(snan))[0]) << std::endl;
    std::cout << "AVX : max(SNaN, QNaN) = " << display(_mm256_max_pd(_mm256_set1_pd(snan), _mm256_set1_pd(qnan))[0]) << std::endl;
    std::cout << "AVX : max(SNaN, SNaN) = " << display(_mm256_max_pd(_mm256_set1_pd(snan), _mm256_set1_pd(snan))[0]) << std::endl;
}
```

</div></summary>

<details><summary>C#</summary><div>

```C#:Main.cs
using System;
using System.Runtime.CompilerServices;

namespace NaNTest
{
    class Program
    {
        private const ulong SignalingNaNBit = 0x7FF0_0000_0000_0001uL;
        private const ulong QuietNaNBit = 0x7FF8_0000_0000_0000uL;
        public static double SNaN { get; } = Unsafe.As<ulong, double>(ref Unsafe.AsRef(SignalingNaNBit));
        public static double QNaN { get; } = Unsafe.As<ulong, double>(ref Unsafe.AsRef(QuietNaNBit));

        public static bool IsSignalingNaN(double value) => double.IsNaN(value) & (Unsafe.As<double, ulong>(ref value) & 0x0008_0000_0000_0000uL) == 0;
        public static bool IsQuietNaN(double value) => double.IsNaN(value) & (Unsafe.As<double, ulong>(ref value) & 0x0008_0000_0000_0000uL) != 0;

        static string Display(double value)
        {
            if(IsSignalingNaN(value))
                return "SNaN";
            if(IsQuietNaN(value))
                return "QNaN";
            return value.ToString();
        }


        static void Main(string[] args)
        {
            Console.WriteLine($"double.NaN is {Display(double.NaN)}");
            Console.WriteLine($"System.Math: max(0.0 , qnan) = {Display(Math.Max(0.0 , QNaN))}");
            Console.WriteLine($"System.Math: max(0.0 , snan) = {Display(Math.Max(0.0 , SNaN))}");
            Console.WriteLine($"System.Math: max(qnan, 0.0 ) = {Display(Math.Max(QNaN, 0.0 ))}");
            Console.WriteLine($"System.Math: max(snan, 0.0 ) = {Display(Math.Max(SNaN, 0.0 ))}");
            Console.WriteLine($"System.Math: max(qnan, qnan) = {Display(Math.Max(QNaN, QNaN))}");
            Console.WriteLine($"System.Math: max(qnan, snan) = {Display(Math.Max(QNaN, SNaN))}");
            Console.WriteLine($"System.Math: max(snan, qnan) = {Display(Math.Max(SNaN, QNaN))}");
            Console.WriteLine($"System.Math: max(snan, snan) = {Display(Math.Max(SNaN, SNaN))}");
        }
    }
}

```

</div></summary>