Add cached powers of 10

This commit is contained in:
Victor Zverovich 2018-04-29 06:33:05 -07:00
parent dd296e1de0
commit 2768af2388
3 changed files with 92 additions and 7 deletions

View File

@ -280,6 +280,55 @@ const uint64_t basic_data<T>::POWERS_OF_10_64[] = {
10000000000000000000ull 10000000000000000000ull
}; };
// Normalized 64-bit significands of pow(10, k), for k = -348, -340, ..., 340.
// These are generated by support/compute-powers.py.
template <typename T>
const uint64_t basic_data<T>::POW10_SIGNIFICANDS[] = {
0xfa8fd5a0081c0288ull, 0xbaaee17fa23ebf76ull, 0x8b16fb203055ac76ull,
0xcf42894a5dce35eaull, 0x9a6bb0aa55653b2dull, 0xe61acf033d1a45dfull,
0xab70fe17c79ac6caull, 0xff77b1fcbebcdc4full, 0xbe5691ef416bd60cull,
0x8dd01fad907ffc3cull, 0xd3515c2831559a83ull, 0x9d71ac8fada6c9b5ull,
0xea9c227723ee8bcbull, 0xaecc49914078536dull, 0x823c12795db6ce57ull,
0xc21094364dfb5637ull, 0x9096ea6f3848984full, 0xd77485cb25823ac7ull,
0xa086cfcd97bf97f4ull, 0xef340a98172aace5ull, 0xb23867fb2a35b28eull,
0x84c8d4dfd2c63f3bull, 0xc5dd44271ad3cdbaull, 0x936b9fcebb25c996ull,
0xdbac6c247d62a584ull, 0xa3ab66580d5fdaf6ull, 0xf3e2f893dec3f126ull,
0xb5b5ada8aaff80b8ull, 0x87625f056c7c4a8bull, 0xc9bcff6034c13053ull,
0x964e858c91ba2655ull, 0xdff9772470297ebdull, 0xa6dfbd9fb8e5b88full,
0xf8a95fcf88747d94ull, 0xb94470938fa89bcfull, 0x8a08f0f8bf0f156bull,
0xcdb02555653131b6ull, 0x993fe2c6d07b7facull, 0xe45c10c42a2b3b06ull,
0xaa242499697392d3ull, 0xfd87b5f28300ca0eull, 0xbce5086492111aebull,
0x8cbccc096f5088ccull, 0xd1b71758e219652cull, 0x9c40000000000000ull,
0xe8d4a51000000000ull, 0xad78ebc5ac620000ull, 0x813f3978f8940984ull,
0xc097ce7bc90715b3ull, 0x8f7e32ce7bea5c70ull, 0xd5d238a4abe98068ull,
0x9f4f2726179a2245ull, 0xed63a231d4c4fb27ull, 0xb0de65388cc8ada8ull,
0x83c7088e1aab65dbull, 0xc45d1df942711d9aull, 0x924d692ca61be758ull,
0xda01ee641a708deaull, 0xa26da3999aef774aull, 0xf209787bb47d6b85ull,
0xb454e4a179dd1877ull, 0x865b86925b9bc5c2ull, 0xc83553c5c8965d3dull,
0x952ab45cfa97a0b3ull, 0xde469fbd99a05fe3ull, 0xa59bc234db398c25ull,
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
// to significands above.
template <typename T>
const int16_t basic_data<T>::POW10_EXPONENTS[] = {
-1220, -1193, -1166, -1140, -1113, -1087, -1060, -1034, -1007, -980, -954,
-927, -901, -874, -847, -821, -794, -768, -741, -715, -688, -661,
-635, -608, -582, -555, -529, -502, -475, -449, -422, -396, -369,
-343, -316, -289, -263, -236, -210, -183, -157, -130, -103, -77,
-50, -24, 3, 30, 56, 83, 109, 136, 162, 189, 216,
242, 269, 295, 322, 348, 375, 402, 428, 455, 481, 508,
534, 561, 588, 614, 641, 667, 694, 720, 747, 774, 800,
827, 853, 880, 907, 933, 960, 986, 1013, 1039, 1066
};
FMT_FUNC fp operator*(fp x, fp y) { FMT_FUNC fp operator*(fp x, fp y) {
// Multiply 32-bit parts of significands. // Multiply 32-bit parts of significands.
uint64_t mask = (1ULL << 32) - 1; uint64_t mask = (1ULL << 32) - 1;

View File

@ -263,10 +263,10 @@ inline fp operator-(fp x, fp y) {
fp operator*(fp x, fp y); fp operator*(fp x, fp y);
// Compute k such that its cached power c_k = c_k.f * pow(2, c_k.e) satisfies // Compute k such that its cached power c_k = c_k.f * pow(2, c_k.e) satisfies
// alpha <= c_k.e + e <= alpha + 3. // min_exponent <= c_k.e + e <= min_exponent + 3.
inline int compute_cached_power_index(int e, int alpha) { inline int compute_cached_power_index(int e, int min_exponent) {
constexpr double one_over_log2_10 = 0.30102999566398114; // 1 / log2(10) constexpr double one_over_log2_10 = 0.30102999566398114; // 1 / log2(10)
return static_cast<int>(std::ceil((alpha - e + 63) * one_over_log2_10)); return static_cast<int>(std::ceil((min_exponent - e + 63) * one_over_log2_10));
} }
template <typename Allocator> template <typename Allocator>
@ -795,6 +795,8 @@ template <typename T = void>
struct FMT_API basic_data { struct FMT_API basic_data {
static const uint32_t POWERS_OF_10_32[]; static const uint32_t POWERS_OF_10_32[];
static const uint64_t POWERS_OF_10_64[]; static const uint64_t POWERS_OF_10_64[];
static const uint64_t POW10_SIGNIFICANDS[];
static const int16_t POW10_EXPONENTS[];
static const char DIGITS[]; static const char DIGITS[];
}; };

View File

@ -8,9 +8,43 @@ min_exponent = -348
max_exponent = 340 max_exponent = 340
step = 8 step = 8
significand_size = 64 significand_size = 64
for exp in range(min_exponent, max_exponent + 1, step): exp_offset = 2000
n = 10 ** exp if exp >= 0 else 2 ** 2000 / 10 ** -exp
class fp:
pass
powers = []
for i, exp in enumerate(range(min_exponent, max_exponent + 1, step)):
result = fp()
n = 10 ** exp if exp >= 0 else 2 ** exp_offset / 10 ** -exp
k = significand_size + 1 k = significand_size + 1
# Convert to binary and round. # Convert to binary and round.
n = (int('{:0<{}b}'.format(n, k)[:k], 2) + 1) / 2 binary = '{:b}'.format(n)
print('{:0<#16x}ull'.format(n)) result.f = (int('{:0<{}}'.format(binary[:k], k), 2) + 1) / 2
result.e = len(binary) - (exp_offset if exp < 0 else 0) - significand_size
powers.append(result)
# Sanity check.
exp_offset10 = 400
actual = result.f * 10 ** exp_offset10
if result.e > 0:
actual *= 2 ** result.e
else:
for j in range(-result.e):
actual /= 2
expected = 10 ** (exp_offset10 + exp)
precision = len('{}'.format(expected)) - len('{}'.format(actual - expected))
if precision < 19:
print('low precision:', precision)
exit(1)
print('Significands:', end='')
for i, fp in enumerate(powers):
if i % 3 == 0:
print(end='\n ')
print(' {:0<#16x}ull'.format(fp.f, ), end=',')
print('\n\nExponents:', end='')
for i, fp in enumerate(powers):
if i % 11 == 0:
print(end='\n ')
print(' {:5}'.format(fp.e), end=',')