Implement normalization and simplify power table

This commit is contained in:
Victor Zverovich
2018-05-27 08:04:30 -07:00
parent 6a5bb6e268
commit 2f257b7291
4 changed files with 90 additions and 52 deletions

View File

@ -294,35 +294,28 @@ const uint64_t basic_data<T>::POWERS_OF_10_64[] = {
// These are generated by support/compute-powers.py. // These are generated by support/compute-powers.py.
template <typename T> template <typename T>
const uint64_t basic_data<T>::POW10_SIGNIFICANDS[] = { const uint64_t basic_data<T>::POW10_SIGNIFICANDS[] = {
0xfa8fd5a0081c0288ull, 0xbaaee17fa23ebf76ull, 0x8b16fb203055ac76ull, 0xfa8fd5a0081c0288, 0xbaaee17fa23ebf76, 0x8b16fb203055ac76, 0xcf42894a5dce35ea,
0xcf42894a5dce35eaull, 0x9a6bb0aa55653b2dull, 0xe61acf033d1a45dfull, 0x9a6bb0aa55653b2d, 0xe61acf033d1a45df, 0xab70fe17c79ac6ca, 0xff77b1fcbebcdc4f,
0xab70fe17c79ac6caull, 0xff77b1fcbebcdc4full, 0xbe5691ef416bd60cull, 0xbe5691ef416bd60c, 0x8dd01fad907ffc3c, 0xd3515c2831559a83, 0x9d71ac8fada6c9b5,
0x8dd01fad907ffc3cull, 0xd3515c2831559a83ull, 0x9d71ac8fada6c9b5ull, 0xea9c227723ee8bcb, 0xaecc49914078536d, 0x823c12795db6ce57, 0xc21094364dfb5637,
0xea9c227723ee8bcbull, 0xaecc49914078536dull, 0x823c12795db6ce57ull, 0x9096ea6f3848984f, 0xd77485cb25823ac7, 0xa086cfcd97bf97f4, 0xef340a98172aace5,
0xc21094364dfb5637ull, 0x9096ea6f3848984full, 0xd77485cb25823ac7ull, 0xb23867fb2a35b28e, 0x84c8d4dfd2c63f3b, 0xc5dd44271ad3cdba, 0x936b9fcebb25c996,
0xa086cfcd97bf97f4ull, 0xef340a98172aace5ull, 0xb23867fb2a35b28eull, 0xdbac6c247d62a584, 0xa3ab66580d5fdaf6, 0xf3e2f893dec3f126, 0xb5b5ada8aaff80b8,
0x84c8d4dfd2c63f3bull, 0xc5dd44271ad3cdbaull, 0x936b9fcebb25c996ull, 0x87625f056c7c4a8b, 0xc9bcff6034c13053, 0x964e858c91ba2655, 0xdff9772470297ebd,
0xdbac6c247d62a584ull, 0xa3ab66580d5fdaf6ull, 0xf3e2f893dec3f126ull, 0xa6dfbd9fb8e5b88f, 0xf8a95fcf88747d94, 0xb94470938fa89bcf, 0x8a08f0f8bf0f156b,
0xb5b5ada8aaff80b8ull, 0x87625f056c7c4a8bull, 0xc9bcff6034c13053ull, 0xcdb02555653131b6, 0x993fe2c6d07b7fac, 0xe45c10c42a2b3b06, 0xaa242499697392d3,
0x964e858c91ba2655ull, 0xdff9772470297ebdull, 0xa6dfbd9fb8e5b88full, 0xfd87b5f28300ca0e, 0xbce5086492111aeb, 0x8cbccc096f5088cc, 0xd1b71758e219652c,
0xf8a95fcf88747d94ull, 0xb94470938fa89bcfull, 0x8a08f0f8bf0f156bull, 0x9c40000000000000, 0xe8d4a51000000000, 0xad78ebc5ac620000, 0x813f3978f8940984,
0xcdb02555653131b6ull, 0x993fe2c6d07b7facull, 0xe45c10c42a2b3b06ull, 0xc097ce7bc90715b3, 0x8f7e32ce7bea5c70, 0xd5d238a4abe98068, 0x9f4f2726179a2245,
0xaa242499697392d3ull, 0xfd87b5f28300ca0eull, 0xbce5086492111aebull, 0xed63a231d4c4fb27, 0xb0de65388cc8ada8, 0x83c7088e1aab65db, 0xc45d1df942711d9a,
0x8cbccc096f5088ccull, 0xd1b71758e219652cull, 0x9c40000000000000ull, 0x924d692ca61be758, 0xda01ee641a708dea, 0xa26da3999aef774a, 0xf209787bb47d6b85,
0xe8d4a51000000000ull, 0xad78ebc5ac620000ull, 0x813f3978f8940984ull, 0xb454e4a179dd1877, 0x865b86925b9bc5c2, 0xc83553c5c8965d3d, 0x952ab45cfa97a0b3,
0xc097ce7bc90715b3ull, 0x8f7e32ce7bea5c70ull, 0xd5d238a4abe98068ull, 0xde469fbd99a05fe3, 0xa59bc234db398c25, 0xf6c69a72a3989f5c, 0xb7dcbf5354e9bece,
0x9f4f2726179a2245ull, 0xed63a231d4c4fb27ull, 0xb0de65388cc8ada8ull, 0x88fcf317f22241e2, 0xcc20ce9bd35c78a5, 0x98165af37b2153df, 0xe2a0b5dc971f303a,
0x83c7088e1aab65dbull, 0xc45d1df942711d9aull, 0x924d692ca61be758ull, 0xa8d9d1535ce3b396, 0xfb9b7cd9a4a7443c, 0xbb764c4ca7a44410, 0x8bab8eefb6409c1a,
0xda01ee641a708deaull, 0xa26da3999aef774aull, 0xf209787bb47d6b85ull, 0xd01fef10a657842c, 0x9b10a4e5e9913129, 0xe7109bfba19c0c9d, 0xac2820d9623bf429,
0xb454e4a179dd1877ull, 0x865b86925b9bc5c2ull, 0xc83553c5c8965d3dull, 0x80444b5e7aa7cf85, 0xbf21e44003acdd2d, 0x8e679c2f5e44ff8f, 0xd433179d9c8cb841,
0x952ab45cfa97a0b3ull, 0xde469fbd99a05fe3ull, 0xa59bc234db398c25ull, 0x9e19db92b4e31ba9, 0xeb96bf6ebadf77d9, 0xaf87023b9bf0ee6b
0xf6c69a72a3989f5cull, 0xb7dcbf5354e9beceull, 0x88fcf317f22241e2ull,
0xcc20ce9bd35c78a5ull, 0x98165af37b2153dfull, 0xe2a0b5dc971f303aull,
0xa8d9d1535ce3b396ull, 0xfb9b7cd9a4a7443cull, 0xbb764c4ca7a44410ull,
0x8bab8eefb6409c1aull, 0xd01fef10a657842cull, 0x9b10a4e5e9913129ull,
0xe7109bfba19c0c9dull, 0xac2820d9623bf429ull, 0x80444b5e7aa7cf85ull,
0xbf21e44003acdd2dull, 0x8e679c2f5e44ff8full, 0xd433179d9c8cb841ull,
0x9e19db92b4e31ba9ull, 0xeb96bf6ebadf77d9ull, 0xaf87023b9bf0ee6bull
}; };
// Binary exponents of pow(10, k), for k = -348, -340, ..., 340, corresponding // Binary exponents of pow(10, k), for k = -348, -340, ..., 340, corresponding

View File

@ -258,7 +258,17 @@ inline dummy_int isnan(...) { return dummy_int(); }
inline dummy_int _isnan(...) { return dummy_int(); } inline dummy_int _isnan(...) { return dummy_int(); }
// A handmade floating-point number f * pow(2, e). // A handmade floating-point number f * pow(2, e).
struct fp { class fp {
private:
// All sizes are in bits.
static constexpr int char_size = std::numeric_limits<unsigned char>::digits;
// Subtract 1 to account for an implicit most significant bit in the
// normalized form.
static constexpr int double_significand_size =
std::numeric_limits<double>::digits - 1;
static constexpr uint64_t implicit_bit = 1ull << double_significand_size;
public:
uint64_t f; uint64_t f;
int e; int e;
@ -270,24 +280,36 @@ struct fp {
explicit fp(Double d) { explicit fp(Double d) {
// Assume double is in the format [sign][exponent][significand]. // Assume double is in the format [sign][exponent][significand].
typedef std::numeric_limits<Double> limits; typedef std::numeric_limits<Double> limits;
const int double_size = const int double_size = sizeof(Double) * char_size;
sizeof(Double) * std::numeric_limits<unsigned char>::digits; const int exponent_size =
// Subtract 1 to account for an implicit most significant bit in the double_size - double_significand_size - 1; // -1 for sign
// normalized form.
const int significand_size = limits::digits - 1;
const int exponent_size = double_size - significand_size - 1; // -1 for sign
const uint64_t implicit_bit = 1ull << significand_size;
const uint64_t significand_mask = implicit_bit - 1; const uint64_t significand_mask = implicit_bit - 1;
const uint64_t exponent_mask = (~0ull >> 1) & ~significand_mask; const uint64_t exponent_mask = (~0ull >> 1) & ~significand_mask;
const int exponent_bias = (1 << exponent_size) - limits::max_exponent - 1; const int exponent_bias = (1 << exponent_size) - limits::max_exponent - 1;
auto u = bit_cast<uint64_t>(d); auto u = bit_cast<uint64_t>(d);
auto biased_e = (u & exponent_mask) >> significand_size; auto biased_e = (u & exponent_mask) >> double_significand_size;
f = u & significand_mask; f = u & significand_mask;
if (biased_e != 0) if (biased_e != 0)
f += implicit_bit; f += implicit_bit;
else else
biased_e = 1; // Subnormals use biased exponent 1 (min exponent). biased_e = 1; // Subnormals use biased exponent 1 (min exponent).
e = biased_e - exponent_bias - significand_size; e = biased_e - exponent_bias - double_significand_size;
}
// Normalizes the value converted from double and multiplied by (1 << SHIFT).
template <int SHIFT = 0>
void normalize() {
// Handle subnormals.
auto shifted_implicit_bit = implicit_bit << SHIFT;
while ((f & shifted_implicit_bit) == 0) {
f <<= 1;
--e;
}
auto fp_significand_size = sizeof(f) * char_size;
// Subtract 1 to account for hidden bit.
auto offset = fp_significand_size - double_significand_size - SHIFT - 1;
f <<= offset;
e -= offset;
} }
}; };

View File

@ -39,9 +39,9 @@ for i, exp in enumerate(range(min_exponent, max_exponent + 1, step)):
print('Significands:', end='') print('Significands:', end='')
for i, fp in enumerate(powers): for i, fp in enumerate(powers):
if i % 3 == 0: if i % 4 == 0:
print(end='\n ') print(end='\n ')
print(' {:0<#16x}ull'.format(fp.f, ), end=',') print(' {:0<#16x}'.format(fp.f, ), end=',')
print('\n\nExponents:', end='') print('\n\nExponents:', end='')
for i, fp in enumerate(powers): for i, fp in enumerate(powers):

View File

@ -881,17 +881,40 @@ TEST(UtilTest, ParseNonnegativeInt) {
fmt::format_error, "number is too big"); fmt::format_error, "number is too big");
} }
TEST(UtilTest, FPSubtract) { template <bool is_iec559>
auto r = fp(123, 1) - fp(102, 1); void test_construct_from_double() {
EXPECT_EQ(r.f, 21u); fmt::print("warning: double is not IEC559, skipping FP tests\n");
EXPECT_EQ(r.e, 1);
} }
TEST(UtilTest, FPMultiply) { template <>
auto r = fp(123ULL << 32, 4) * fp(56ULL << 32, 7); void test_construct_from_double<true>() {
EXPECT_EQ(r.f, 123u * 56u); auto v = fp(1.23);
EXPECT_EQ(r.e, 4 + 7 + 64); EXPECT_EQ(v.f, 0x13ae147ae147aeu);
r = fp(123ULL << 32, 4) * fp(567ULL << 31, 8); EXPECT_EQ(v.e, -52);
EXPECT_EQ(r.f, (123 * 567 + 1u) / 2); }
EXPECT_EQ(r.e, 4 + 8 + 64);
TEST(FPTest, ConstructFromDouble) {
test_construct_from_double<std::numeric_limits<double>::is_iec559>();
}
TEST(FPTest, Normalize) {
auto v = fp(0xbeef, 42);
v.normalize();
EXPECT_EQ(0xbeef000000000000, v.f);
EXPECT_EQ(-6, v.e);
}
TEST(FPTest, Subtract) {
auto v = fp(123, 1) - fp(102, 1);
EXPECT_EQ(v.f, 21u);
EXPECT_EQ(v.e, 1);
}
TEST(FPTest, Multiply) {
auto v = fp(123ULL << 32, 4) * fp(56ULL << 32, 7);
EXPECT_EQ(v.f, 123u * 56u);
EXPECT_EQ(v.e, 4 + 7 + 64);
v = fp(123ULL << 32, 4) * fp(567ULL << 31, 8);
EXPECT_EQ(v.f, (123 * 567 + 1u) / 2);
EXPECT_EQ(v.e, 4 + 8 + 64);
} }