x86でdoubleがfloatより速いかどうかを検証してみた



# 昔話
それは昔々のこと。
x86には浮動小数点演算を行う手段がなく、外付けの浮動小数点演算ユニットを接続するという手法で、浮動小数点演算を実現していたのであった。
x87と呼ばれたそれはとてもエクセレントなシステムで…という話はwikipediaに譲ろう。
https://ja.wikipedia.org/wiki/Intel_8087
重要なのは、x87が内部表現として80bitの拡張倍精度を使っている、ということ。
これのおかげで、x87においては、確かに~~**doubleのほうが速かった**
(floatだとdoubleへのキャストコストが発生するため)~~
嘘だろそれ。ASM見たら別にキャストとかしてなかったわ。
どっちかというと丸めの影響で精度が異なることのほうが重要だわ。
改めて調べてみると、doubleが速いとされている資料についてはあんまりないことに気付く。
(同等としている資料はherumiさんの http://homepage1.nifty.com/herumi/prog/prog90.html がある)
速いってのはそれがデマだったのでは…？

(2016/06/13追記: コメントにてご指摘いただきましたが、確かにfloatがdoubleより遅い条件があります。詳しくは後述)

しかしx87の登場からすでに30年(!)
x86も時代が移り変わり、x87は内部に統合され、MMX(こいつももはや死語),SSE,AVXと色々な命令が増えてきた。
この時代にあって、doubleがfloatより本当に速いのかどうか、調べてみよう。

# ソースコード
gistに置いてあるので適当にとっていってほしい。
https://gist.github.com/telmin/846e1ee9ecddd4fcc063

# 仮説
速くないと一刀両断するのは簡単なので、doubleとfloatでどういう性能差が発生するかまず考えてみる。

## キャッシュのはなし
floatとdoubleはご存知のとおり32bitと64bitで、使用するデータ量は倍になる。
とすると、当然キャッシュヒットしないような状態で計算するとあっさり性能は倍だろう。これはあんまり面白くないのと、いろいろと他でややこしくなるので、全データがキャッシュに乗るような状態を考えよう。

## SSE命令のはなし
現在、浮動小数点SIMD命令として、x86(_64)には、SSE, AVXが用意されている。
話を簡単にするためにSSEのみにフォーカスする。
最近のAuto-Vectorizationは結構優秀で、簡単な演算は割と勝手にベクタライズしてくれる。ありがたいことだ。それでもまだ手で書きたい。一日一回感謝のIntrinsics!!!

でまぁ、floatとdoubleのmulでどのようにレイテンシとスループットが違うかというと
Intel Intrinsics guide
https://software.intel.com/sites/landingpage/IntrinsicsGuide/
より

- \_\_m128 \_mm\_mul\_ps (\_\_m128 a, \_\_m128 b)

Architecture|Latency|Throughput
----|----|----
Haswell|5|0.5
Ivy Bridge|5|1
Sandy Bridge|5|1
Westmere|4|1
Nehalem|4|1

- \_\_m128d \_mm\_mul\_pd (\_\_m128d a, \_\_m128d b)

Architecture|Latency|Throughput
----|----|----
Haswell|5|0.5
Ivy Bridge|5|1
Sandy Bridge|5|1
Westmere|4|1
Nehalem|4|1

となり、全く一緒になる。
なのであとは純粋にbit長に依存して、4Flop, 2Flopとなるわけだ。

あと、SSEにはxmmレジスタを入力としながらも先頭32bit, 64bitのみを演算として使用する命令もあり、普段x86_64で浮動小数点を使うとこちらが使用される。あとベクタライズしろって言ってんのにこいつを吐くこともある。
これが、それぞれ
- \_\_m128 \_mm\_mul\_ss (\_\_m128 a, \_\_m128 b)
- \_\_m128 \_mm\_mul\_sd (\_\_m128 a, \_\_m128 b)

レイテンシ、スループットはps, pdと同じ(当たり前)

# x87命令のはなし
SSEのときと同様に、x87命令のレイテンシ、スループットを考えてみる。
x87の浮動小数点乗算命令はfmulである。
レイテンシ、スループットはhttp://www.agner.org/optimize/instruction_tables.pdf より

Architecture | Latency | Throughput|
----|----|----
Haswell | 5 | 1  
Ivy Bridge | 5 | 1
Sandy Bridge | 5 | 1

となっている。
なんか全体的に同じ感じするわな。

## つまり
ベクタライズしていればfloatとdoubleが約2倍の性能になるはず。
してない場合はx87、SSEのいずれにおいても、doubleとfloatは同等の性能になるはず。
x87の場合は、~~float -> doubleへのキャストコストが発生するからdoubleのほうが若干速いはず~~ 嘘。同じ速度ぐらいになるはず。

# 実験環境 & 実験条件
## 環境
MacbookPro Retina 15(mid 2012)
CPU : Corei7-3720QM IvyBridge
RAM : 8GB
コンパイラ : Apple LLVM version 6.0 (clang-600.0.54) (based on LLVM 3.5svn)

## 条件
float, double 4096要素に対してaxpyを実施
中間変数含め、使うのは4096 * 4要素
axpyを100回実行してかかった時間を合計する。

x87命令動かすのに32bitバイナリじゃないとちゃんと動かなかったので以下はm32オプションを付与しているものとする。
~~つかなんで動かないのか誰か教えてプリーズ。~~
ABI絡みだった気がする…？

# 実験結果
とりあえずSSE

```
$ clang++ -O2 double.cpp -std=c++11 -m32
```

```
float time 185000 nsec
double time 390000 nsec
```

x87

```
$ clang++ -O2 double.cpp -std=c++11 -m32 -mno-sse
```

```
float time 561000 nsec
double time 549000 nsec
```

SSEでベクタライズ止めてみる

```
$ clang++ -O2 double.cpp -std=c++11 -m32 -fno-vectorize
```

```
float time 549000 nsec
double time 553000 nsec
```

表にまとめると以下のようになる。

 | float | double|
----|----|----
SSE| 185000 | 390000 |
x87| 561000 | 549000 |
SSE(no-vectorize)| 549000 | 553000 |

x87でfloatが若干遅いのはなんだろうか。何度か実行してみると結果が若干揺れるからそのせいかな…

# まとめ
あんまり優位な差は出なかったので結論が微妙になりそうな気がしないでもないが、少なくともSSEとかがあるこんにちでは、精度的に大丈夫な場合はfloatを使うと、SIMDの恩恵にあずかれるからそっちのほうがいいんじゃないかとか思う。
あと、x87は頑張らないと出ないので、もう拡張倍精度のことは忘れたほうがいい。マジで。

# おまけ
ARMでの結果は以下
http://d.hatena.ne.jp/nakamura001/20090301/1235924387
(2009年なので、最新だとどうなってるかはわからないけどあんまり変わらんでしょ)

# おまけその2
ちゃんとSSE命令とか出てるぅ？
SSE版

```
LBB0_17:                               ## %vector.body116
                                        ##   Parent Loop BB0_13 Depth=1
                                        ## =>  This Inner Loop Header: Depth=2
        movups  (%eax,%ebx,4), %xmm0
        movups  (%ecx,%ebx,4), %xmm1
        mulps   %xmm0, %xmm1
        movups  (%edx,%ebx,4), %xmm0
        addps   %xmm1, %xmm0
        movups  %xmm0, (%esi,%ebx,4)
        addl    $4, %ebx
        cmpl    $4096, %ebx             ## imm = 0x1000
        jne     LBB0_17

```

x87版

```
LBB0_14:                                ##   Parent Loop BB0_13 Depth=1
                                        ## =>  This Inner Loop Header: Depth=2
        flds    (%eax,%ebx,4)
        fmuls   (%ecx,%ebx,4)
        fadds   (%edx,%ebx,4)
        fstps   (%esi,%ebx,4)
        incl    %ebx
        cmpl    $4096, %ebx             ## imm = 0x1000
        jne     LBB0_14
```

出てるっぽい。

# おまけその3
floatがdoubleより遅くなるシチュエーションは存在した…！
K&Rの頃のC言語は、floatはdoubleに昇格して演算を行うルールであった。
(fujitanozomu さんにコメントいただきました。ありがとうございました)

https://docs.oracle.com/cd/E19957-01/806-4836/compat.html
これでも記載されている通り、

| | K&R | ANSI/ISO C |
|----|----|----|
| 単精度計算と倍精度計算 | 浮動小数点式のオペランドを double に拡張する。float を返すように宣言された関数の戻り値は、常に double に拡張される。 | float の演算を単精度計算で行うことができる。このような関数に float の戻り型を使用できる。 |

となっているので、ANSIの規格化を以って、ちゃんとfloatが(言語仕様的には)使えるようになった模様。

あと、現在でも下手するとこの型昇格は踏むらしい。それが以下のような条件。

```cpp:hoge.c
extern float fa, fb, fc, fd;
extern double da, db, dc, dd;

extern float fcomp();
extern double dcomp();

hoge()
{
    fd = fcomp(fa, fb, fc);
}

piyo()
{
    dd = dcomp(da, db, dc);
}
```


```
$ clang -O2 -m32 -S hoge.c -o -
_hoge:                                  # @hoge
        subl    $24, %esp
        movss   _fa, %xmm0              # xmm0 = mem[0],zero,zero,zero
        cvtss2sd        %xmm0, %xmm0
        movss   _fb, %xmm1              # xmm1 = mem[0],zero,zero,zero
        cvtss2sd        %xmm1, %xmm1
        movss   _fc, %xmm2              # xmm2 = mem[0],zero,zero,zero
        cvtss2sd        %xmm2, %xmm2
        movsd   %xmm2, 16(%esp)
        movsd   %xmm1, 8(%esp)
        movsd   %xmm0, (%esp)
        calll   _fcomp
        fstps   _fd
        addl    $24, %esp
        retl

_piyo:                                  # @piyo
        subl    $24, %esp
        movsd   _da, %xmm0              # xmm0 = mem[0],zero
        movsd   _db, %xmm1              # xmm1 = mem[0],zero
        movsd   _dc, %xmm2              # xmm2 = mem[0],zero
        movsd   %xmm2, 16(%esp)
        movsd   %xmm1, 8(%esp)
        movsd   %xmm0, (%esp)
        calll   _dcomp
        fstpl   _dd
        addl    $24, %esp
        retl
```

```cpp:piyo.c
float fcomp(a, b, c)
float a;
float b;
float c;
{
    return a + b * c;
}

double dcomp(a, b, c)
double a;
double b;
double c;
{
    return a + b * c;
}
```

```
$ clang -O2 -m32 -S piyo.c -o -
_fcomp:                                 # @fcomp
        pushl   %eax
        movsd   24(%esp), %xmm0         # xmm0 = mem[0],zero
        movsd   16(%esp), %xmm1         # xmm1 = mem[0],zero
        movsd   8(%esp), %xmm2          # xmm2 = mem[0],zero
        cvtsd2ss        %xmm2, %xmm2
        cvtsd2ss        %xmm1, %xmm1
        cvtsd2ss        %xmm0, %xmm0
        mulss   %xmm1, %xmm0
        addss   %xmm2, %xmm0
        movss   %xmm0, (%esp)
        flds    (%esp)
        popl    %eax
        retl

_dcomp:                                 # @dcomp
        pushl   %ebp
        movl    %esp, %ebp
        andl    $-8, %esp
        subl    $8, %esp
        movsd   16(%ebp), %xmm0         # xmm0 = mem[0],zero
        mulsd   24(%ebp), %xmm0
        addsd   8(%ebp), %xmm0
        movsd   %xmm0, (%esp)
        fldl    (%esp)
        movl    %ebp, %esp
        popl    %ebp
        retl
```

> 現代のコンパイラでも、外部関数についてプロトタイプ宣言を行わず、昔の K&R スタイルの仮引数のない外部宣言のみで float の値を引数として与えて関数呼び出しを行うと、float の値は一旦 double に昇格されてから引数として渡されるコードが生成される

とご指摘いただいたので、この辺は気を付けないといけなさそう。