New frsqrte implementation; verified accurate.

This is similar to the old implementation, but it uses smaller tables, and
handles more edge cases correctly.  (hwtest coming soon.)
This commit is contained in:
magumagu 2014-05-22 06:12:29 -07:00
parent 129e76e60d
commit a9a2d3d98d
1 changed files with 55 additions and 19 deletions

View File

@ -6,7 +6,6 @@
#include "Common/CPUDetect.h"
#include "Common/MathUtil.h"
#include "Core/PowerPC/LUT_frsqrtex.h"
#include "Core/PowerPC/Interpreter/Interpreter.h"
using namespace MathUtil;
@ -333,31 +332,68 @@ inline double ApproximateReciprocal(double val)
inline double ApproximateReciprocalSquareRoot(double val)
{
if (val < 0)
return PPC_NAN;
if (val == 0.0)
return INFINITY;
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
{
union {
double valf;
u64 vali;
long long vali;
};
valf = val;
long long mantissa = vali & ((1LL << 52) - 1);
long long sign = vali & (1ULL << 63);
long long exponent = vali & (0x7FFLL << 52);
u32 fsa = vali >> 32;
u32 idx = (fsa >> 5) % (sizeof(frsqrtex_lut) / sizeof(frsqrtex_lut[0]));
// Special case 0
if (mantissa == 0 && exponent == 0)
return sign ? -INFINITY : INFINITY;
// Special case NaN-ish numbers
if (exponent == (0x7FFLL << 52))
{
if (mantissa == 0)
return sign ? NAN : 0.0;
return 0.0 + valf;
}
// Negative numbers return NaN
if (sign)
return NAN;
s32 e = fsa >> (32 - 12);
e &= 2047;
e -= 1023;
s32 oe = -((e + 1) / 2);
oe -= ((e + 1) & 1);
if (!exponent)
{
// "Normalize" denormal values
do
{
exponent -= 1LL << 52;
mantissa <<= 1;
} while (!(mantissa & (1LL << 52)));
mantissa &= (1LL << 52) - 1;
exponent += 1LL << 52;
}
u32 outb = frsqrtex_lut[idx] << 20;
u32 outa = ((oe + 1023) & 2047) << 20;
outa |= frsqrtex_lut[idx] >> 12;
bool odd_exponent = !(exponent & (1LL << 52));
exponent = ((0x3FFLL << 52) - ((exponent - (0x3FELL << 52)) / 2)) & (0x7FFLL << 52);
vali = ((u64)outa << 32) + (u64)outb;
int i = (int)(mantissa >> 37);
vali = sign | exponent;
int index = i / 2048 + (odd_exponent ? 16 : 0);
vali |= (long long)(expected_base[index] - expected_dec[index] * (i % 2048)) << 26;
return valf;
}