Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Merge pull request #390 from magumagu/fp-reciprocal 
Add accurate implementations of frsqrte and fres
  • Loading branch information
@neobrain
neobrain committed May 24, 2014
2 parents ca30029 + 6955fef commit cffa848
Show file tree
Hide file tree
Showing 13 changed files with 153 additions and 161 deletions.
136 changes: 136 additions & 0 deletions Source/Core/Core/PowerPC/Interpreter/Interpreter_FPUtils.h
View file
Expand Up @@ -4,6 +4,8 @@

#pragma once

#include <limits>

#include "Common/CPUDetect.h"
#include "Common/MathUtil.h"
#include "Core/PowerPC/Interpreter/Interpreter.h"
Expand Down Expand Up @@ -271,3 +273,137 @@ inline u64 ConvertToDouble(u32 _x)
return ((x & 0xc0000000) << 32) | z | ((x & 0x3fffffff) << 29);
}
}

// Used by fres and ps_res.
inline double ApproximateReciprocal(double val)
{
static const int expected_base[] = {
0x7ff800, 0x783800, 0x70ea00, 0x6a0800,
0x638800, 0x5d6200, 0x579000, 0x520800,
0x4cc800, 0x47ca00, 0x430800, 0x3e8000,
0x3a2c00, 0x360800, 0x321400, 0x2e4a00,
0x2aa800, 0x272c00, 0x23d600, 0x209e00,
0x1d8800, 0x1a9000, 0x17ae00, 0x14f800,
0x124400, 0x0fbe00, 0x0d3800, 0x0ade00,
0x088400, 0x065000, 0x041c00, 0x020c00,
};
static const int expected_dec[] = {
0x3e1, 0x3a7, 0x371, 0x340,
0x313, 0x2ea, 0x2c4, 0x2a0,
0x27f, 0x261, 0x245, 0x22a,
0x212, 0x1fb, 0x1e5, 0x1d1,
0x1be, 0x1ac, 0x19b, 0x18b,
0x17c, 0x16e, 0x15b, 0x15b,
0x143, 0x143, 0x12d, 0x12d,
0x11a, 0x11a, 0x108, 0x106,
};

union
{
double valf;
s64 vali;
};
valf = val;
s64 mantissa = vali & ((1LL << 52) - 1);
s64 sign = vali & (1ULL << 63);
s64 exponent = vali & (0x7FFLL << 52);

// Special case 0
if (mantissa == 0 && exponent == 0)
return sign ? -std::numeric_limits<double>::infinity() :
std::numeric_limits<double>::infinity();
// Special case NaN-ish numbers
if (exponent == (0x7FFLL << 52))
{
if (mantissa == 0)
return sign ? -0.0 : 0.0;
return 0.0 + valf;
}
// Special case small inputs
if (exponent < (895LL << 52))
return sign ? -FLT_MAX : FLT_MAX;
// Special case large inputs
if (exponent >= (1149LL << 52))
return sign ? -0.0f : 0.0f;

exponent = (0x7FDLL << 52) - exponent;

int i = (int)(mantissa >> 37);
vali = sign | exponent;
vali |= (s64)(expected_base[i / 1024] - (expected_dec[i / 1024] * (i % 1024) + 1) / 2) << 29;
return valf;
}

inline double ApproximateReciprocalSquareRoot(double val)
{
static const int expected_base[] = {
0x3ffa000, 0x3c29000, 0x38aa000, 0x3572000,
0x3279000, 0x2fb7000, 0x2d26000, 0x2ac0000,
0x2881000, 0x2665000, 0x2468000, 0x2287000,
0x20c1000, 0x1f12000, 0x1d79000, 0x1bf4000,
0x1a7e800, 0x17cb800, 0x1552800, 0x130c000,
0x10f2000, 0x0eff000, 0x0d2e000, 0x0b7c000,
0x09e5000, 0x0867000, 0x06ff000, 0x05ab800,
0x046a000, 0x0339800, 0x0218800, 0x0105800,
};
static const int expected_dec[] = {
0x7a4, 0x700, 0x670, 0x5f2,
0x584, 0x524, 0x4cc, 0x47e,
0x43a, 0x3fa, 0x3c2, 0x38e,
0x35e, 0x332, 0x30a, 0x2e6,
0x568, 0x4f3, 0x48d, 0x435,
0x3e7, 0x3a2, 0x365, 0x32e,
0x2fc, 0x2d0, 0x2a8, 0x283,
0x261, 0x243, 0x226, 0x20b,
};

union
{
double valf;
s64 vali;
};
valf = val;
s64 mantissa = vali & ((1LL << 52) - 1);
s64 sign = vali & (1ULL << 63);
s64 exponent = vali & (0x7FFLL << 52);

// Special case 0
if (mantissa == 0 && exponent == 0)
return sign ? -std::numeric_limits<double>::infinity() :
std::numeric_limits<double>::infinity();
// Special case NaN-ish numbers
if (exponent == (0x7FFLL << 52))
{
if (mantissa == 0)
{
if (sign)
return std::numeric_limits<double>::quiet_NaN();
return 0.0;
}
return 0.0 + valf;
}
// Negative numbers return NaN
if (sign)
return std::numeric_limits<double>::quiet_NaN();

if (!exponent)
{
// "Normalize" denormal values
do
{
exponent -= 1LL << 52;
mantissa <<= 1;
} while (!(mantissa & (1LL << 52)));
mantissa &= (1LL << 52) - 1;
exponent += 1LL << 52;
}

bool odd_exponent = !(exponent & (1LL << 52));
exponent = ((0x3FFLL << 52) - ((exponent - (0x3FELL << 52)) / 2)) & (0x7FFLL << 52);

int i = (int)(mantissa >> 37);
vali = sign | exponent;
int index = i / 2048 + (odd_exponent ? 16 : 0);
vali |= (s64)(expected_base[index] - expected_dec[index] * (i % 2048)) << 26;
return valf;
}
48 changes: 4 additions & 44 deletions Source/Core/Core/PowerPC/Interpreter/Interpreter_FloatingPoint.cpp
View file
Expand Up @@ -18,7 +18,6 @@
#endif

#include "Common/MathUtil.h"
#include "Core/PowerPC/LUT_frsqrtex.h"
#include "Core/PowerPC/Interpreter/Interpreter.h"
#include "Core/PowerPC/Interpreter/Interpreter_FPUtils.h"

Expand Down Expand Up @@ -388,19 +387,7 @@ void Interpreter::fdivsx(UGeckoInstruction _inst)
void Interpreter::fresx(UGeckoInstruction _inst)
{
double b = rPS0(_inst.FB);
double one_over = ForceSingle(1.0 / b);
// this is based on the real hardware tests
if (b != 0.0 && IsINF(one_over))
{
if (one_over > 0)
riPS0(_inst.FD) = riPS1(_inst.FD) = MAX_SINGLE;
else
riPS0(_inst.FD) = riPS1(_inst.FD) = MIN_SINGLE;
}
else
{
rPS0(_inst.FD) = rPS1(_inst.FD) = one_over;
}
rPS0(_inst.FD) = rPS1(_inst.FD) = ApproximateReciprocal(b);
if (b == 0.0)
{
SetFPException(FPSCR_ZX);
Expand All @@ -415,39 +402,12 @@ void Interpreter::frsqrtex(UGeckoInstruction _inst)
if (b < 0.0)
{
SetFPException(FPSCR_VXSQRT);
rPS0(_inst.FD) = PPC_NAN;
}
else
else if (b == 0.0)
{
#ifdef VERY_ACCURATE_FP
if (b == 0.0)
{
SetFPException(FPSCR_ZX);
riPS0(_inst.FD) = 0x7ff0000000000000;
}
else
{
u32 fsa = Common::swap32(Common::swap64(riPS0(_inst.FB)));
u32 fsb = Common::swap32(Common::swap64(riPS0(_inst.FB)) >> 32);
u32 idx=(fsa >> 5) % (sizeof(frsqrtex_lut) / sizeof(frsqrtex_lut[0]));

s32 e = fsa >> (32-12);
e &= 2047;
e -= 1023;
s32 oe =- ((e + 1) / 2);
oe -= ((e + 1) & 1);

u32 outb = frsqrtex_lut[idx] << 20;
u32 outa = ((oe + 1023) & 2047) << 20;
outa |= frsqrtex_lut[idx] >> 12;
riPS0(_inst.FD) = ((u64)outa << 32) + (u64)outb;
}
#else
if (b == 0.0)
SetFPException(FPSCR_ZX);
rPS0(_inst.FD) = ForceDouble(1.0 / sqrt(b));
#endif
SetFPException(FPSCR_ZX);
}
rPS0(_inst.FD) = ApproximateReciprocalSquareRoot(b);
UpdateFPRF(rPS0(_inst.FD));
if (_inst.Rc) Helper_UpdateCR1();
}
Expand Down
45 changes: 7 additions & 38 deletions Source/Core/Core/PowerPC/Interpreter/Interpreter_Paired.cpp
View file
Expand Up @@ -178,58 +178,27 @@ void Interpreter::ps_res(UGeckoInstruction _inst)
{
SetFPException(FPSCR_ZX);
}
rPS0(_inst.FD) = ForceSingle(1.0 / a);
if (a != 0.0 && IsINF(rPS0(_inst.FD)))
{
if (rPS0(_inst.FD) > 0)
riPS0(_inst.FD) = MAX_SINGLE; // largest finite single
else
riPS0(_inst.FD) = MIN_SINGLE; // most negative finite single
}
rPS1(_inst.FD) = ForceSingle(1.0 / b);
if (b != 0.0 && IsINF(rPS1(_inst.FD)))
{
if (rPS1(_inst.FD) > 0)
riPS1(_inst.FD) = MAX_SINGLE;
else
riPS1(_inst.FD) = MIN_SINGLE;
}
rPS0(_inst.FD) = ApproximateReciprocal(a);
rPS1(_inst.FD) = ApproximateReciprocal(b);
UpdateFPRF(rPS0(_inst.FD));
if (_inst.Rc) Helper_UpdateCR1();
}

void Interpreter::ps_rsqrte(UGeckoInstruction _inst)
{
// this code is based on the real hardware tests
if (rPS0(_inst.FB) == 0.0 || rPS1(_inst.FB) == 0.0)
{
SetFPException(FPSCR_ZX);
}
// PS0
if (rPS0(_inst.FB) < 0.0)
{
SetFPException(FPSCR_VXSQRT);
rPS0(_inst.FD) = PPC_NAN;
}
else
{
rPS0(_inst.FD) = 1.0 / sqrt(rPS0(_inst.FB));
u32 t = ConvertToSingle(riPS0(_inst.FD));
rPS0(_inst.FD) = *(float*)&t;
}
// PS1
if (rPS1(_inst.FB) < 0.0)

if (rPS0(_inst.FB) < 0.0 || rPS1(_inst.FB) < 0.0)
{
SetFPException(FPSCR_VXSQRT);
rPS1(_inst.FD) = PPC_NAN;
}
else
{
rPS1(_inst.FD) = 1.0 / sqrt(rPS1(_inst.FB));
u32 t = ConvertToSingle(riPS1(_inst.FD));
rPS1(_inst.FD) = *(float*)&t;
}

rPS0(_inst.FD) = ApproximateReciprocalSquareRoot(rPS0(_inst.FB));
rPS1(_inst.FD) = ApproximateReciprocalSquareRoot(rPS1(_inst.FB));

UpdateFPRF(rPS0(_inst.FD));
if (_inst.Rc) Helper_UpdateCR1();
}
Expand Down
2 changes: 0 additions & 2 deletions Source/Core/Core/PowerPC/Jit64/Jit.h
View file
Expand Up @@ -168,12 +168,10 @@ class Jit64 : public Jitx86Base
void ps_arith(UGeckoInstruction inst); //aggregate
void ps_mergeXX(UGeckoInstruction inst);
void ps_maddXX(UGeckoInstruction inst);
void ps_recip(UGeckoInstruction inst);
void ps_sum(UGeckoInstruction inst);
void ps_muls(UGeckoInstruction inst);

void fp_arith(UGeckoInstruction inst);
void frsqrtex(UGeckoInstruction inst);

void fcmpx(UGeckoInstruction inst);
void fmrx(UGeckoInstruction inst);
Expand Down
6 changes: 3 additions & 3 deletions Source/Core/Core/PowerPC/Jit64/Jit64_Tables.cpp
View file
Expand Up @@ -138,9 +138,9 @@ static GekkoOPTemplate table4_2[] =
{20, &Jit64::ps_arith}, //"ps_sub", OPTYPE_PS, 0}},
{21, &Jit64::ps_arith}, //"ps_add", OPTYPE_PS, 0}},
{23, &Jit64::ps_sel}, //"ps_sel", OPTYPE_PS, 0}},
{24, &Jit64::ps_recip}, //"ps_res", OPTYPE_PS, 0}},
{24, &Jit64::FallBackToInterpreter}, //"ps_res", OPTYPE_PS, 0}},
{25, &Jit64::ps_arith}, //"ps_mul", OPTYPE_PS, 0}},
{26, &Jit64::ps_recip}, //"ps_rsqrte", OPTYPE_PS, 0, 1}},
{26, &Jit64::FallBackToInterpreter}, //"ps_rsqrte", OPTYPE_PS, 0, 1}},
{28, &Jit64::ps_maddXX}, //"ps_msub", OPTYPE_PS, 0}},
{29, &Jit64::ps_maddXX}, //"ps_madd", OPTYPE_PS, 0}},
{30, &Jit64::ps_maddXX}, //"ps_nmsub", OPTYPE_PS, 0}},
Expand Down Expand Up @@ -360,7 +360,7 @@ static GekkoOPTemplate table63_2[] =
{22, &Jit64::FallBackToInterpreter}, //"fsqrtx", OPTYPE_FPU, FL_RC_BIT_F}},
{23, &Jit64::FallBackToInterpreter}, //"fselx", OPTYPE_FPU, FL_RC_BIT_F}},
{25, &Jit64::fp_arith}, //"fmulx", OPTYPE_FPU, FL_RC_BIT_F}},
{26, &Jit64::frsqrtex}, //"frsqrtex", OPTYPE_FPU, FL_RC_BIT_F}},
{26, &Jit64::FallBackToInterpreter}, //"frsqrtex", OPTYPE_FPU, FL_RC_BIT_F}},
{28, &Jit64::fmaddXX}, //"fmsubx", OPTYPE_FPU, FL_RC_BIT_F}},
{29, &Jit64::fmaddXX}, //"fmaddx", OPTYPE_FPU, FL_RC_BIT_F}},
{30, &Jit64::fmaddXX}, //"fnmsubx", OPTYPE_FPU, FL_RC_BIT_F}},
Expand Down
15 changes: 0 additions & 15 deletions Source/Core/Core/PowerPC/Jit64/Jit_FloatingPoint.cpp
View file
Expand Up @@ -101,21 +101,6 @@ void Jit64::fp_arith(UGeckoInstruction inst)
}
}

void Jit64::frsqrtex(UGeckoInstruction inst)
{
INSTRUCTION_START
JITDISABLE(bJITFloatingPointOff)
int d = inst.FD;
int b = inst.FB;
fpr.Lock(b, d);
fpr.BindToRegister(d, true, true);
MOVSD(XMM0, M((void *)&one_const));
SQRTSD(XMM1, fpr.R(b));
DIVSD(XMM0, R(XMM1));
MOVSD(fpr.R(d), XMM0);
fpr.UnlockAll();
}

void Jit64::fmaddXX(UGeckoInstruction inst)
{
INSTRUCTION_START
Expand Down
35 changes: 0 additions & 35 deletions Source/Core/Core/PowerPC/Jit64/Jit_Paired.cpp
View file
Expand Up @@ -112,41 +112,6 @@ void Jit64::ps_sign(UGeckoInstruction inst)
fpr.UnlockAll();
}

// ps_res and ps_rsqrte
void Jit64::ps_recip(UGeckoInstruction inst)
{
INSTRUCTION_START
JITDISABLE(bJITPairedOff)

if (inst.Rc)
{
FallBackToInterpreter(inst);
return;
}

OpArg divisor;
int d = inst.FD;
int b = inst.FB;
fpr.Lock(d, b);
fpr.BindToRegister(d, (d == b));
switch (inst.SUBOP5)
{
case 24:
// ps_res
divisor = fpr.R(b);
break;
case 26:
// ps_rsqrte
SQRTPD(XMM0, fpr.R(b));
divisor = R(XMM0);
break;
}
MOVAPD(XMM1, M((void*)&psOneOne));
DIVPD(XMM1, divisor);
MOVAPD(fpr.R(d), XMM1);
fpr.UnlockAll();
}

//add a, b, c

//mov a, b
Expand Down

0 comments on commit cffa848

Please sign in to comment.