x86系CPU限定floatとdoubleに関する意外と知られていないかもしれない事実.
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版よりメモリアクセス回数がを減らしてる分若干優秀な感じ. まあ, この程度は完全にキャッシュに入るので実行速度的には変わりませんが.
挙動の異なる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の計算結果が異なってるのはおかしくありませんか.
; 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
演算精度を変更してみましょう. 変更するには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 // 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などの処理時間は変わりません.
ちなみにこの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をつけよということです.
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になりますが, 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だぜ, を実現するためなんでしょうね.
// 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の場合加筆