FloatUtils: Remove union type punning from ApproximateReciprocal functions

This form of type punning invokes undefined behavior in C++
This commit is contained in:
Lioncash 2018-05-09 09:40:40 -04:00
parent 46a4243d9a
commit fe218ea3f6
1 changed files with 37 additions and 34 deletions

View File

@ -5,6 +5,7 @@
#include "Common/FloatUtils.h" #include "Common/FloatUtils.h"
#include <cmath> #include <cmath>
#include <cstring>
namespace Common namespace Common
{ {
@ -99,20 +100,20 @@ const std::array<BaseAndDec, 32> frsqrte_expected = {{
double ApproximateReciprocalSquareRoot(double val) double ApproximateReciprocalSquareRoot(double val)
{ {
union s64 integral;
{ std::memcpy(&integral, &val, sizeof(integral));
double valf;
s64 vali; s64 mantissa = integral & ((1LL << 52) - 1);
}; const s64 sign = integral & (1ULL << 63);
valf = val; s64 exponent = integral & (0x7FFLL << 52);
s64 mantissa = vali & ((1LL << 52) - 1);
s64 sign = vali & (1ULL << 63);
s64 exponent = vali & (0x7FFLL << 52);
// Special case 0 // Special case 0
if (mantissa == 0 && exponent == 0) if (mantissa == 0 && exponent == 0)
{
return sign ? -std::numeric_limits<double>::infinity() : return sign ? -std::numeric_limits<double>::infinity() :
std::numeric_limits<double>::infinity(); std::numeric_limits<double>::infinity();
}
// Special case NaN-ish numbers // Special case NaN-ish numbers
if (exponent == (0x7FFLL << 52)) if (exponent == (0x7FFLL << 52))
{ {
@ -124,7 +125,7 @@ double ApproximateReciprocalSquareRoot(double val)
return 0.0; return 0.0;
} }
return 0.0 + valf; return 0.0 + val;
} }
// Negative numbers return NaN // Negative numbers return NaN
@ -143,15 +144,18 @@ double ApproximateReciprocalSquareRoot(double val)
exponent += 1LL << 52; exponent += 1LL << 52;
} }
bool odd_exponent = !(exponent & (1LL << 52)); const bool odd_exponent = !(exponent & (1LL << 52));
exponent = ((0x3FFLL << 52) - ((exponent - (0x3FELL << 52)) / 2)) & (0x7FFLL << 52); exponent = ((0x3FFLL << 52) - ((exponent - (0x3FELL << 52)) / 2)) & (0x7FFLL << 52);
integral = sign | exponent;
int i = (int)(mantissa >> 37); const int i = static_cast<int>(mantissa >> 37);
vali = sign | exponent; const int index = i / 2048 + (odd_exponent ? 16 : 0);
int index = i / 2048 + (odd_exponent ? 16 : 0);
const auto& entry = frsqrte_expected[index]; const auto& entry = frsqrte_expected[index];
vali |= (s64)(entry.m_base - entry.m_dec * (i % 2048)) << 26; integral |= static_cast<s64>(entry.m_base - entry.m_dec * (i % 2048)) << 26;
return valf;
double result;
std::memcpy(&result, &integral, sizeof(result));
return result;
} }
const std::array<BaseAndDec, 32> fres_expected = {{ const std::array<BaseAndDec, 32> fres_expected = {{
@ -167,44 +171,43 @@ const std::array<BaseAndDec, 32> fres_expected = {{
// Used by fres and ps_res. // Used by fres and ps_res.
double ApproximateReciprocal(double val) double ApproximateReciprocal(double val)
{ {
union s64 integral;
{ std::memcpy(&integral, &val, sizeof(integral));
double valf;
s64 vali;
};
valf = val; const s64 mantissa = integral & ((1LL << 52) - 1);
s64 mantissa = vali & ((1LL << 52) - 1); const s64 sign = integral & (1ULL << 63);
s64 sign = vali & (1ULL << 63); s64 exponent = integral & (0x7FFLL << 52);
s64 exponent = vali & (0x7FFLL << 52);
// Special case 0 // Special case 0
if (mantissa == 0 && exponent == 0) if (mantissa == 0 && exponent == 0)
return std::copysign(std::numeric_limits<double>::infinity(), valf); return std::copysign(std::numeric_limits<double>::infinity(), val);
// Special case NaN-ish numbers // Special case NaN-ish numbers
if (exponent == (0x7FFLL << 52)) if (exponent == (0x7FFLL << 52))
{ {
if (mantissa == 0) if (mantissa == 0)
return std::copysign(0.0, valf); return std::copysign(0.0, val);
return 0.0 + valf; return 0.0 + val;
} }
// Special case small inputs // Special case small inputs
if (exponent < (895LL << 52)) if (exponent < (895LL << 52))
return std::copysign(std::numeric_limits<float>::max(), valf); return std::copysign(std::numeric_limits<float>::max(), val);
// Special case large inputs // Special case large inputs
if (exponent >= (1149LL << 52)) if (exponent >= (1149LL << 52))
return std::copysign(0.0, valf); return std::copysign(0.0, val);
exponent = (0x7FDLL << 52) - exponent; exponent = (0x7FDLL << 52) - exponent;
int i = (int)(mantissa >> 37); const int i = static_cast<int>(mantissa >> 37);
const auto& entry = fres_expected[i / 1024]; const auto& entry = fres_expected[i / 1024];
vali = sign | exponent; integral = sign | exponent;
vali |= (s64)(entry.m_base - (entry.m_dec * (i % 1024) + 1) / 2) << 29; integral |= static_cast<s64>(entry.m_base - (entry.m_dec * (i % 1024) + 1) / 2) << 29;
return valf;
double result;
std::memcpy(&result, &integral, sizeof(result));
return result;
} }
} // namespace Common } // namespace Common