Fast Power Bug ?

FastPower(Single, Single) does not work with Base = 0,474733531475067 eg and Exponent = 150, should return zero.
The same checks as in rtl before Result := FastExp2(AExponent * FastLog2(ABase)); make it work as it should.

function FastPower(const ABase, AExponent: Single): Single;
begin
Result := FastExp2(AExponent * FastLog2(ABase));
end;

// 32 bit implementation

function FastExp2(const A: Single): Single; assembler;
var
OldFlags, NewFlags: UInt32;
asm
// Set rounding mode to Round Positive (=Round Down)
stmxcsr [OldFlags]
movss xmm0, [A]
mov ecx, [OldFlags]
xorps xmm1, xmm1
and ecx, SSE_ROUND_MASK
movss xmm3, xmm0
or ecx, SSE_ROUND_DOWN
movss xmm5, xmm0
mov [NewFlags], ecx

movss xmm1, DWORD [SSE_EXP2_F1]
ldmxcsr [NewFlags]

// Z := A - RoundDown(A)
cvtps2dq xmm3, xmm3
addss xmm1, xmm5 // A + 121.2740575
cvtdq2ps xmm3, xmm3
movss xmm2, DWORD [SSE_EXP2_F2]
subss xmm0, xmm3

movss xmm3, DWORD [SSE_EXP2_F3]
movss xmm4, DWORD [SSE_EXP2_F4]
subss xmm3, xmm0 // (4.84252568 - Z)
mulss xmm0, xmm4 // 1.49012907 * Z
divss xmm2, xmm3
movss xmm5, DWORD [SSE_EXP2_F5]
addss xmm1, xmm2 // A + 121.2740575 + 27.7280233 / (4.84252568 - Z)
subss xmm1, xmm0 // A + 121.2740575 + 27.7280233 / (4.84252568 - Z) - 1.49012907 * Z
mulss xmm1, xmm5 // (1 shl 23) * (A + 121.2740575 + 27.7280233 / (4.84252568 - Z) - 1.49012907 * Z)
cvtps2dq xmm1, xmm1

// Restore rounding mode
ldmxcsr [OldFlags]

movss [Result], xmm1
end;

function FastLog2(const A: Single): Single; assembler;
asm
movss xmm0, [A]
movss xmm2, DWORD [SSE_MASK_FRACTION]
movss xmm1, xmm0

// MX.I := (VX.I and $007FFFFF) or $3F000000
movss xmm3, DWORD [SSE_LOG2_I1]
pand xmm0, xmm2
cvtdq2ps xmm1, xmm1
movss xmm4, DWORD [SSE_LOG2_F1]
por xmm0, xmm3

movss xmm2, DWORD [SSE_LOG2_F2]
mulss xmm1, xmm4 // VX.I * 1.1920928955078125e-7
movss xmm3, DWORD [SSE_LOG2_F3]
subss xmm1, xmm2 // Result - 124.22551499
mulss xmm3, xmm0
movss xmm4, DWORD [SSE_LOG2_F5]
subss xmm1, xmm3 // Result - 124.22551499 - 1.498030302 * MX.S
movss xmm2, DWORD [SSE_LOG2_F4]
addss xmm0, xmm4
divss xmm2, xmm0
subss xmm1, xmm2 // Result - 124.22551499 - 1.498030302 * MX.S - 1.72587999 / (0.3520887068 + MX.S)

movss [Result], xmm1
end;


// 64 bit implementation

function FastExp2(const A: Single): Single; assembler;
var
OldFlags, NewFlags: UInt32;
asm
// Set rounding mode to Round Positive (=Round Down)
stmxcsr [OldFlags]
mov ecx, [OldFlags]
xorps xmm1, xmm1
and ecx, SSE_ROUND_MASK
movss xmm3, xmm0
or ecx, SSE_ROUND_DOWN
movss xmm5, xmm0
mov [NewFlags], ecx

movss xmm1, DWORD [SSE_EXP2_F1]
ldmxcsr [NewFlags]

// Z := A - RoundDown(A)
cvtps2dq xmm3, xmm3
addss xmm1, xmm5 // A + 121.2740575
cvtdq2ps xmm3, xmm3
movss xmm2, DWORD [SSE_EXP2_F2]
subss xmm0, xmm3

movss xmm3, DWORD [SSE_EXP2_F3]
movss xmm4, DWORD [SSE_EXP2_F4]
subss xmm3, xmm0 // (4.84252568 - Z)
mulss xmm0, xmm4 // 1.49012907 * Z
divss xmm2, xmm3
movss xmm5, DWORD [SSE_EXP2_F5]
addss xmm1, xmm2 // A + 121.2740575 + 27.7280233 / (4.84252568 - Z)
subss xmm1, xmm0 // A + 121.2740575 + 27.7280233 / (4.84252568 - Z) - 1.49012907 * Z
mulss xmm1, xmm5 // (1 shl 23) * (A + 121.2740575 + 27.7280233 / (4.84252568 - Z) - 1.49012907 * Z)
cvtps2dq xmm1, xmm1

// Restore rounding mode
ldmxcsr [OldFlags]

movss xmm0, xmm1
end;

function FastLog2(const A: Single): Single; assembler;
asm
movss xmm2, DWORD [SSE_MASK_FRACTION]
movss xmm1, xmm0

// MX.I := (VX.I and $007FFFFF) or $3F000000
movss xmm3, DWORD [SSE_LOG2_I1]
pand xmm0, xmm2
cvtdq2ps xmm1, xmm1
movss xmm4, DWORD [SSE_LOG2_F1]
por xmm0, xmm3

movss xmm2, DWORD [SSE_LOG2_F2]
mulss xmm1, xmm4 // VX.I * 1.1920928955078125e-7
movss xmm3, DWORD [SSE_LOG2_F3]
subss xmm1, xmm2 // Result - 124.22551499
mulss xmm3, xmm0
movss xmm4, DWORD [SSE_LOG2_F5]
subss xmm1, xmm3 // Result - 124.22551499 - 1.498030302 * MX.S
movss xmm2, DWORD [SSE_LOG2_F4]
addss xmm0, xmm4
divss xmm2, xmm0
subss xmm1, xmm2 // Result - 124.22551499 - 1.498030302 * MX.S - 1.72587999 / (0.3520887068 + MX.S)

movss xmm0, xmm1
end;


Any idea what bug could be ?

Bye,
Skybuck.