float型とdouble型

x86系CPU限定floatとdoubleに関する意外と知られていないかもしれない事実.


注意
以下はSSEを使わない場合の話です. 32-bit環境ではSSEを使うようにコンパイルオプションを設定した場合(cf SSE用オプション), あるいは64-bit環境ではデフォルトでSSEを使うので以下の話は適用されません. その場合はSSEの場合を参照してください.
クイズ
唐突ですが, 下記のコードを実行したとき速いのはどちらでしょう. また演算結果はどうなると予想されますか.
double calcd(int n)
{
    double ret = 0;
    for (int i = 1; i < n; i++) {
        ret += 1 / double(i);
    }
    return ret;
}

float calcf(int n)
{
    float ret = 0;
    for (int i = 1; i < n; i++) {
        ret += 1 / float(i);
    }
    return ret;
}

int main()
{
    printf("calcf=%.15f\n", calcf(100000000));
    printf("calcd=%.15f\n", calcd(100000000));
}

float型は遅い? いえいえ, 結果はVC6だとcalcd, calcfともに同じ時間で同じ結果でした.

    18.997896403852554

mingw/gcc on Linuxでもcalcd, calcfともに同じ時間で同じ結果でした.

    18.997896403853908

VC2005でのみcalcd, calcfは異なる結果を出しました.

    ; calcf
    15.403682708740234
    ; calcd
    18.997896403852554

デバッグオプションでは他の要因が入ってくるので異なる結果になるかもしれませんが, -O2などの普通の最適化オプションをつけておけば同様の結果が得られると思います. そうなる理由はアセンブリ言語出力をしてみれば一目瞭然. VC2005以外ではfloatとdoubleで(ほぼ)同じコードが生成されているからです.

// calcdの内側のループ
.lp_d:
    fild   DWORD PTR _i$[esp+4]      // 変数iをFPUレジスタに
    inc    eax                       // i++
    cmp    eax, ecx                  // while (i < n) {
    mov    DWORD PTR _i$[esp+4], eax // i保存
    fdivr  QWORD PTR __real@3ff0000000000000 // 1/i
    faddp  ST(1), ST(0)              // ret += (1/i)
    jl     .lp_d                     // }

// calcfの内側のループ
.lp_f:
    fild   DWORD PTR _i$[esp+4]
    inc    eax
    cmp    eax, ecx
    mov    DWORD PTR _i$[esp+4], eax
    fdivr  DWORD PTR __real@3f800000
    faddp  ST(1), ST(0)
    jl     SHORT .lp_f

違うのはfdivrの部分ですが1はfloatであろうとdoubleであろうと同じ値なので(処理時間と結果は)変わりません. gccでは完全同一のコードを生成してました.

// stに1が代入されている
L14:
    pushl    %eax         // i保存
    incl     %eax         // i++
    fildl    (%esp)       // 変数iをFPUレジスタに
    addl     $4, %esp     // stack開放
    cmpl     %edx, %eax
    fdivr    %st(1), %st  // 1/i
    faddp    %st, %st(2)  // ret += 1/i
    jl    L14

gcc版はVC版よりメモリアクセス回数がを減らしてる分若干優秀な感じ. まあ, この程度は完全にキャッシュに入るので実行速度的には変わりませんが.

あれ, pushとpopはペアじゃないと遅いCPUって無かったっけ. P4はいいのかな
.

挙動の異なるVC2005のcalcfではretをレジスタではなくメモリ上に置いていました.

    fild    DWORD PTR _i$3677[esp+16]      // 変数iをFPUレジスタに
    add     ecx, 1                         // i++
    cmp     ecx, esi                       // while (i < n) {
    mov     DWORD PTR _i$3677[esp+16], ecx // i保存
    fdivr   ST(0), ST(1)                   // tmp = 1/i
    fadd    DWORD PTR _ret$[esp+16]        // st0 += tmp
    fstp    DWORD PTR _ret$[esp+16]        // ret = st0 ; *** ここ ***
    jl      SHORT $LC3@calcf

このようにVC2005を除くコンパイラは演算をレジスタ上ですませられる場合 doubleとfloatの区別をしないようです. なおVC2005には/fp:fastオプションがあり, これをつけると今までのVCと同じ振る舞いをします.

    fild    DWORD PTR _i$3677[esp+12]
    add     ecx, 1
    cmp     ecx, esi
    mov     DWORD PTR _i$3677[esp+12], ecx
    fdivr   ST(0), ST(1)
    faddp   ST(2), ST(0)
    jl      SHORT $LC3@calcf

「そうか, floatも常にdoubleとして計算されるのか」と思われた方, それは少々早合点です. 上記WindowsとLinuxの計算結果が異なってるのはおかしくありませんか.

同じコードなのに
「演算順序が違ってたまたま異なる結果が出ただけじゃないの」と思われるかもしれません. 試しにcalcをアセンブリ言語で書いてWindows/Linux上で実行してみましょう.
どうでもいいけどFPU命令使ったの久しぶり. 忘れてる, 忘れてる
; for nasm
proc calc
        xor     eax, eax
        mov     ecx, [esp + 4]
        fldz
        inc     eax
        fld1                    ; 1, 0
        push    eax
.lp:
        fild    dword [esp]     ; i, 1, ret
        inc     eax
        cmp     eax, ecx
        mov     [esp], eax
        fdivr   st0, st1        ; 1/i, 1, ret
        faddp   st2, st0        ; 1, ret += 1/i
        jl      .lp
        pop     ecx
        fstp    st0             ; ret
        ret

同一ソースなんだから同一結果と思いきや, やはりWindowsとLinuxとでは結果が違ってます. 正確にはVC ≠ (mingw@Windows = gcc@Linux). しかも計算時間はgcc系の方が掛かってます. VC2005もVC6と同じ結果です. これは一体何故.

// VC + calc
    ret=18.997896403852554
    time=2403
// ming + calc
    ret=18.997896403853908
    time=2713
精度が違う
実はx86系FPUでは演算精度を変更することができます. 24bit/53bit/64bitの設定ができるのですが, デフォルトではWindowsは53bit, gcc@Linuxは64bitの精度を利用しているのでした.

演算精度を変更してみましょう. 変更するにはfnstcw命令を使いますが, Cの関数(非ANSI)でも可能です.

#if defined(_WIN32)
    #include <float.h>
#elif defined(__linux__)
    #include <fpu_control.h>
#endif

/*!
    演算精度を変更する
    @param len [in] 精度長(24, 53, 64)
*/
void ChangeFpuPrecision(int len)
{
#ifdef _WIN32
    unsigned int flag;
    switch (len) {
    case 24:
        flag = _PC_24;
        break;
    case 53:
        flag = _PC_53;
        break;
    case 64:
        flag = _PC_64;
        break;
    default:
        return;
    }
    _controlfp(flag, _MCW_PC);
#else
    unsigned int flag;
    switch (len) {
    case 24:
        flag = _FPU_SINGLE;
        break;
    case 53:
        flag = _FPU_DOUBLE;
        break;
    case 64:
        flag = _FPU_EXTENDED;
        break;
    default:
        return;
    }
    unsigned short cw;
    _FPU_GETCW(cw);
    cw = (cw & ~0x300) | flag;
    _FPU_SETCW(cw);
#endif
}

このChangeFpuPrecision()で設定してからcalcd/calcfを呼ぶようにしました. 下記の結果は設定無し+asm, 設定無し+float, 設定無し+double, 設定あり+doubleです.

ちなみに最初VC6では最適化によりcalcdを一度しか計算しなくて焦りました. しょうもない副作用のあるコードを入れて対処.
//VC6
// default
asm
ret=18.997896403852554 ; 53bit
time=2403
float
ret=18.997896403852554 ; 53bit
time=2433
double
ret=18.997896403852554 ; 53bit
time=2424

// ChangeFpuPrecision()
53bit
state=9001F
ret=18.997896403852554 ; 53bit
time=2423
24bit
state=A001F
ret=15.403682708740234 ; 24bit
time=1452
64bit
state=8001F
ret=18.997896403853908 ; 64bit
time=2764

// mingw
// default
asm
ret=18.997896403853908 ; 64bit
time=2743
float
ret=18.997896403853908 ; 64bit
time=2754
double
ret=18.997896403853908 ; 64bit
time=2734

// ChangeFpuPrecision()
53bit
state=9001F
ret=18.997896403852554 ; 53bit
time=2413
24bit
state=A001F
ret=15.403682708740234 ; 24bit
time=1462
64bit
state=8001F
ret=18.997896403853908 ; 64bit
time=2744

ChangeFpuPrecision()で設定すると環境によらず同一の結果が得られているのが分かります. 演算結果もぴったり同じ(ま, たまたまかもしれませんが). 処理時間が24bit < 53bit < 64bitであることも分かります. ただこれはdivが入ってるからで, addやsubなどの処理時間は変わりません.

結局答えは
「float/doubleともにFPUに設定された演算精度に依存する」でした.

ちなみにこのbit数は仮数のことなので24bitはfloat(32bit), 53bitはdouble(64bit), 64bitはlong double(80bit:gccのみ)相当です. いずれにせよ配列を扱うのでなければ, コンパイラにとってはfloatでもdoubleでも同じということです.

要は演算途中の精度が(通常)高いのは勘弁してよというスタンスなんでしょうね. そうじゃないと(Linuxならたとえdoubleであっても)変換コードを入れる必要があり, 遅くなって仕方がないし. ループ内変数にvolatileをつけたのに近い速度低下です.

実際, 区別するようになったVC2005のデフォルト動作の変更に気がつかないとこのコンパイラでのみ激しく遅いという現象に悩むはめになります. 他のコンパイラと同程度でよい場合は/fp:fastオプションをつけた方がよいでしょう.

gcc-4.5のマニュアルによるとgcc 4.5のstrict C99モードでは同様に厳密にハンドリングするようになったようでデフォルトでは遅くなるようです(未確認). -fexcess-precision=fastをつけよということです.

蛇足
ただ無視しているコンパイラでもdoubleをfloatに明示的にキャストすればちゃんと精度を落としてくれます.
    float f = 0.1f;
    double d = 0.1;
    printf("f=%.15f\n", f);
    printf("d=%.15f\n", d);
    printf("d=%.15f\n", float(d));

このコードの結果

// mingw 3.2.3
    f=0.100000001490116
    d=0.100000000000000
    d=0.100000001490116 // floatにcastされた

// VC6
    f=0.100000001490116
    d=0.100000000000000
    d=0.100000000000000 // してないやん(^^;

// bcc5.5.1
    f=0.100000001490116
    d=0.100000000000000
    d=0.100000000000000 // 同類(^^;

// icc
    f=0.100000001490116
    d=0.100000000000000
    d=0.100000001490116 // OK

// VC7(8), VC2005
    f=0.100000001490116
    d=0.100000000000000
    d=0.100000001490116 // OK
蛇足の蛇足
上記答えで「FPUに設定された演算精度に依存する」とはいいましたが, 精度を53bitにすればdoubleと完全等価になるか, というとそうではないのがまたややこしい.

FPUの演算精度は53bitになりますが, FPUのレジスタ幅が80bitあることには変わりません. これにより, 数値がレジスタ上にある限り, doubleのつもりなのにdoubleより大きな数を扱えてしまう, という奇妙な現象を引き起こします(80bitレジスタの指数部の14ビットは変化しないため).

下記コードは与えられたデータ *x を2n倍したあと2nで割った値を返すという(オーバーフローしない限り)意味のないコードなのですが,

proc mulTest;(const double *x, int n);
        mov        eax, [esp + 4]
        fld        qword [q0p5]    ; 0.5の読み込み
        fld        qword [eax]
        mov        eax, [esp + 8]
.lp:
        fadd       st0, st0        ; x <= x + x
        dec        eax
        jnz        .lp

        mov        eax, [esp + 8]
.lp2:
        fmul       st0, st1        ; x <= x / 2
        dec        eax
        jnz        .lp2

        fxch       st1
        fstp       st0
        ret

これを

    int n = 3;
    double f = DBL_MAX, g;
    printf("f=%e(%08X:%08X)\n", f, *((int*)(&f) + 1), *(int*)&f);
    g = test2(&f, n);
    printf("g=%e(%08X:%08X)\n", g, *((int*)(&g) + 1), *(int*)&g);

というコードで利用するとgはinfiniteにならずに元のDBL_MAXに戻るんです. nはどこまで大きくできるかというと, 80bitの限界値16384(=214)からDBL_MAXの指数値1024(=210)を引いた数15360までいけます.

本来なら例外なりオーバーフローなりを出さないといけないはずが, そうなってないところに他プロセッサとの微妙で(悲しい?)非互換性があります. SSE/SSE2で積極的に例外などを扱うようにしたのは本当のIEEE754 64bit doubleだぜ, を実現するためなんでしょうね.

SSEの場合
floatやdoubleをSSEレジスタに割り当てて演算する場合, FPUの様な内部精度が大きいことによる微妙な差異は現れません. floatは23bit精度, doubleは53bit精度で厳密に計算されます. つまり32bit VC(2005以降)のFPUを使った場合のデフォルトの挙動と同一の結果を得ます. 演算速度もfloatとdoubleで変わりません.
// cl /arch:SSE2 /Ox
calcf=15.403682708740234
calcd=18.997896403852554
// gcc -msse2 -mfpmath=sse -O3
calcf=15.403682708740234
calcd=18.997896403852554

cf. FPUの例外マスク
2007/03/20 VC2005での振る舞い加筆
2009/08/15 SSEの場合加筆

1st:2005/08/08, last update:2016/08/29

御意見は光成滋生<herumi@nifty.com>までお願いします