べき乗 pow() xy

ふとべき乗計算を考えてみた。

double pow( double x, int y ){
return y > 0 ? pow( x, y - 1 ) * x : 1;
}

整数の時は簡単。

実数になると

__declspec(noinline) double pow( double x, double y ){
__asm{
fld y
fld x
fyl2x
fst ST(1)
frndint
fxch ST(1)
fsub ST(0), ST(1)
f2xm1
fld1
faddp ST(1), ST(0)
fscale
}
}

なんかinline assemblerを使ってしまった。うん、何やってるのかさっぱりです。命令とそのときのスタック状態を表にすると

命令	ST(0)	ST(1)	ST(2)
fld y	y
fld x	x	y
fyl2x	y×log2x
fst ST(1)	y×log2x	y×log2x
frndint	⌊y×log2x⌋	y×log2x
fxch ST(1)	y×log2x	⌊y×log2x⌋
fsub ST(0), ST(1)	y×log2x-⌊y×log2x⌋	⌊y×log2x⌋
f2xm1	2y×log2x-⌊y×log2x⌋-1	⌊y×log2x⌋
fld1	1	2y×log2x-⌊y×log2x⌋-1	⌊y×log2x⌋
faddp ST(1), ST(0)	2y×log2x-⌊y×log2x⌋	⌊y×log2x⌋
fscale	2y×log2x-⌊y×log2x⌋×2⌊y×log2x⌋

ここで2^{y×log₂x-⌊y×log₂x⌋}×2^{⌊y×log₂x⌋} = 2^{y×log₂x-⌊y×log₂x⌋+⌊y×log₂x⌋} = 2^y×log₂x = x^yです。
はい、さっぱりわかりませんでした。