Implement more FP formatting options

This commit is contained in:
Victor Zverovich 2018-08-29 09:34:57 -07:00
parent 46484da711
commit dd8c5ce442
4 changed files with 237 additions and 209 deletions

View File

@ -340,6 +340,96 @@ const int16_t basic_data<T>::POW10_EXPONENTS[] = {
template <typename T> const char basic_data<T>::RESET_COLOR[] = "\x1b[0m";
template <typename T> const wchar_t basic_data<T>::WRESET_COLOR[] = L"\x1b[0m";
// A handmade floating-point number f * pow(2, e).
class fp {
private:
typedef uint64_t significand_type;
// All sizes are in bits.
static FMT_CONSTEXPR_DECL const int char_size =
std::numeric_limits<unsigned char>::digits;
// Subtract 1 to account for an implicit most significant bit in the
// normalized form.
static FMT_CONSTEXPR_DECL const int double_significand_size =
std::numeric_limits<double>::digits - 1;
static FMT_CONSTEXPR_DECL const uint64_t implicit_bit =
1ull << double_significand_size;
public:
significand_type f;
int e;
static FMT_CONSTEXPR_DECL const int significand_size =
sizeof(significand_type) * char_size;
fp(): f(0), e(0) {}
fp(uint64_t f, int e): f(f), e(e) {}
// Constructs fp from an IEEE754 double. It is a template to prevent compile
// errors on platforms where double is not IEEE754.
template <typename Double>
explicit fp(Double d) {
// Assume double is in the format [sign][exponent][significand].
typedef std::numeric_limits<Double> limits;
const int double_size = sizeof(Double) * char_size;
const int exponent_size =
double_size - double_significand_size - 1; // -1 for sign
const uint64_t significand_mask = implicit_bit - 1;
const uint64_t exponent_mask = (~0ull >> 1) & ~significand_mask;
const int exponent_bias = (1 << exponent_size) - limits::max_exponent - 1;
auto u = bit_cast<uint64_t>(d);
auto biased_e = (u & exponent_mask) >> double_significand_size;
f = u & significand_mask;
if (biased_e != 0)
f += implicit_bit;
else
biased_e = 1; // Subnormals use biased exponent 1 (min exponent).
e = static_cast<int>(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;
}
// Subtract 1 to account for hidden bit.
auto offset = significand_size - double_significand_size - SHIFT - 1;
f <<= offset;
e -= offset;
}
// Compute lower and upper boundaries (m^- and m^+ in the Grisu paper), where
// a boundary is a value half way between the number and its predecessor
// (lower) or successor (upper). The upper boundary is normalized and lower
// has the same exponent but may be not normalized.
void compute_boundaries(fp &lower, fp &upper) const {
lower = f == implicit_bit ?
fp((f << 2) - 1, e - 2) : fp((f << 1) - 1, e - 1);
upper = fp((f << 1) + 1, e - 1);
upper.normalize<1>(); // 1 is to account for the exponent shift above.
lower.f <<= lower.e - upper.e;
lower.e = upper.e;
}
};
// Returns an fp number representing x - y. Result may not be normalized.
inline fp operator-(fp x, fp y) {
FMT_ASSERT(x.f >= y.f && x.e == y.e, "invalid operands");
return fp(x.f - y.f, x.e);
}
// Computes an fp number r with r.f = x.f * y.f / pow(2, 64) rounded to nearest
// with half-up tie breaking, r.e = x.e + y.e + 64. Result may not be normalized.
FMT_API fp operator*(fp x, fp y);
// Returns cached power (of 10) c_k = c_k.f * pow(2, c_k.e) such that its
// (binary) exponent satisfies min_exponent <= c_k.e <= min_exponent + 3.
FMT_API fp get_cached_power(int min_exponent, int &pow10_exponent);
FMT_FUNC fp operator*(fp x, fp y) {
// Multiply 32-bit parts of significands.
uint64_t mask = (1ULL << 32) - 1;
@ -446,31 +536,52 @@ FMT_FUNC void grisu2_gen_digits(
}
}
FMT_FUNC void grisu2_format_positive(double value, char *buffer, size_t &size,
int &dec_exp) {
FMT_ASSERT(value > 0, "value is nonpositive");
fp fp_value(value);
fp lower, upper; // w^- and w^+ in the Grisu paper.
fp_value.compute_boundaries(lower, upper);
// Find a cached power of 10 close to 1 / upper.
const int min_exp = -60; // alpha in Grisu.
auto dec_pow = get_cached_power( // \tilde{c}_{-k} in Grisu.
min_exp - (upper.e + fp::significand_size), dec_exp);
dec_exp = -dec_exp;
fp_value.normalize();
fp scaled_value = fp_value * dec_pow;
fp scaled_lower = lower * dec_pow; // \tilde{M}^- in Grisu.
fp scaled_upper = upper * dec_pow; // \tilde{M}^+ in Grisu.
++scaled_lower.f; // \tilde{M}^- + 1 ulp -> M^-_{\uparrow}.
--scaled_upper.f; // \tilde{M}^+ - 1 ulp -> M^+_{\downarrow}.
uint64_t delta = scaled_upper.f - scaled_lower.f;
grisu2_gen_digits(scaled_value, scaled_upper, delta, buffer, size, dec_exp);
}
// Prettifies the output of the Grisu2 algorithm.
// The number is given as v = buffer * 10^exp.
FMT_FUNC void grisu2_prettify(char *buffer, size_t &size, int exp, char type,
size_t precision, bool print_decimal_point) {
FMT_FUNC void grisu2_prettify(char *buffer, size_t &size, int exp,
int precision) {
// pow(10, full_exp - 1) <= v <= pow(10, full_exp).
int full_exp = static_cast<int>(size) + exp;
int int_size = static_cast<int>(size);
// 10^(full_exp - 1) <= v <= 10^full_exp.
int full_exp = int_size + exp;
if (int_size <= full_exp && full_exp <= 21) {
// 1234e7 -> 12340000000
const int exp_threshold = 21;
if (int_size <= full_exp && full_exp <= exp_threshold) {
// 1234e7 -> 12340000000[.0+]
std::uninitialized_fill_n(buffer + int_size, full_exp - int_size, '0');
char *p = buffer + full_exp;
if (print_decimal_point && size < precision) {
if (precision > 0) {
*p++ = '.';
auto fill_size = precision - size;
std::uninitialized_fill_n(p, fill_size, '0');
p += fill_size;
std::uninitialized_fill_n(p, precision, '0');
p += precision;
}
size = to_unsigned(p - buffer);
} else if (0 < full_exp && full_exp <= 21) {
// 1234e-2 -> 12.34
size_t fractional_size = to_unsigned(int_size - full_exp);
} else if (0 < full_exp && full_exp <= exp_threshold) {
// 1234e-2 -> 12.34[0+]
int fractional_size = -exp;
std::memmove(buffer + full_exp + 1, buffer + full_exp, fractional_size);
buffer[full_exp] = '.';
if (type == 'f' && fractional_size < precision) {
size_t num_zeros = precision - fractional_size;
if (fractional_size < precision) {
int num_zeros = precision - fractional_size;
std::uninitialized_fill_n(buffer + size + 1, num_zeros, '0');
size += num_zeros;
}
@ -493,29 +604,14 @@ FMT_FUNC void grisu2_prettify(char *buffer, size_t &size, int exp, char type,
}
}
FMT_FUNC void grisu2_format_positive(double value, char *buffer, size_t &size,
int &dec_exp) {
FMT_ASSERT(value > 0, "value is nonpositive");
fp fp_value(value);
fp lower, upper; // w^- and w^+ in the Grisu paper.
fp_value.compute_boundaries(lower, upper);
// Find a cached power of 10 close to 1 / upper.
const int min_exp = -60; // alpha in Grisu.
auto dec_pow = get_cached_power( // \tilde{c}_{-k} in Grisu.
min_exp - (upper.e + fp::significand_size), dec_exp);
dec_exp = -dec_exp;
fp_value.normalize();
fp scaled_value = fp_value * dec_pow;
fp scaled_lower = lower * dec_pow; // \tilde{M}^- in Grisu.
fp scaled_upper = upper * dec_pow; // \tilde{M}^+ in Grisu.
++scaled_lower.f; // \tilde{M}^- + 1 ulp -> M^-_{\uparrow}.
--scaled_upper.f; // \tilde{M}^+ - 1 ulp -> M^+_{\downarrow}.
uint64_t delta = scaled_upper.f - scaled_lower.f;
grisu2_gen_digits(scaled_value, scaled_upper, delta, buffer, size, dec_exp);
inline void round(char *, size_t &size, int &exp, int diff) {
// TODO: round instead of truncating
size -= to_unsigned(diff);
exp += diff;
}
// Formats value using Grisu2 algorithm. Grisu2 doesn't give any guarantees on
// the shortness of the result.
// Formats a nonnegative value using Grisu2 algorithm. Grisu2 doesn't give any
// guarantees on the shortness of the result.
FMT_FUNC void grisu2_format(double value, char *buffer, size_t &size, char type,
int precision, bool print_decimal_point) {
FMT_ASSERT(value >= 0, "value is negative");
@ -526,14 +622,27 @@ FMT_FUNC void grisu2_format(double value, char *buffer, size_t &size, char type,
*buffer = '0';
size = 1;
}
size_t unsigned_precision = precision >= 0 ? precision : 6;
if (size > unsigned_precision) {
// TODO: round instead of truncating
dec_exp += static_cast<int>(size - unsigned_precision);
size = unsigned_precision;
const int default_precision = 6;
if (precision < 0)
precision = default_precision;
if (!type || type == 'g' || type == 'G') {
int extra_digits = static_cast<int>(size) - precision;
if (extra_digits > 0)
round(buffer, size, dec_exp, extra_digits);
precision = 0;
} else if (type == 'f' || type == 'F') {
if (precision > 0)
print_decimal_point = true;
int extra_digits = -dec_exp - precision;
if (extra_digits > 0) {
if (extra_digits >= static_cast<int>(size))
extra_digits = static_cast<int>(size) - 1;
round(buffer, size, dec_exp, extra_digits);
}
}
grisu2_prettify(buffer, size, dec_exp, type, unsigned_precision,
print_decimal_point);
if (print_decimal_point && precision < 1)
precision = 1;
grisu2_prettify(buffer, size, dec_exp, precision);
}
} // namespace internal

View File

@ -150,7 +150,7 @@ FMT_END_NAMESPACE
#endif
#ifndef FMT_USE_GRISU
# define FMT_USE_GRISU 0
# define FMT_USE_GRISU 1
#endif
// __builtin_clz is broken in clang with Microsoft CodeGen:
@ -275,96 +275,10 @@ inline dummy_int _finite(...) { 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).
class fp {
private:
typedef uint64_t significand_type;
// All sizes are in bits.
static FMT_CONSTEXPR_DECL const int char_size =
std::numeric_limits<unsigned char>::digits;
// Subtract 1 to account for an implicit most significant bit in the
// normalized form.
static FMT_CONSTEXPR_DECL const int double_significand_size =
std::numeric_limits<double>::digits - 1;
static FMT_CONSTEXPR_DECL const uint64_t implicit_bit =
1ull << double_significand_size;
public:
significand_type f;
int e;
static FMT_CONSTEXPR_DECL const int significand_size =
sizeof(significand_type) * char_size;
fp(): f(0), e(0) {}
fp(uint64_t f, int e): f(f), e(e) {}
// Constructs fp from an IEEE754 double. It is a template to prevent compile
// errors on platforms where double is not IEEE754.
template <typename Double>
explicit fp(Double d) {
// Assume double is in the format [sign][exponent][significand].
typedef std::numeric_limits<Double> limits;
const int double_size = sizeof(Double) * char_size;
const int exponent_size =
double_size - double_significand_size - 1; // -1 for sign
const uint64_t significand_mask = implicit_bit - 1;
const uint64_t exponent_mask = (~0ull >> 1) & ~significand_mask;
const int exponent_bias = (1 << exponent_size) - limits::max_exponent - 1;
auto u = bit_cast<uint64_t>(d);
auto biased_e = (u & exponent_mask) >> double_significand_size;
f = u & significand_mask;
if (biased_e != 0)
f += implicit_bit;
else
biased_e = 1; // Subnormals use biased exponent 1 (min exponent).
e = static_cast<int>(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;
}
// Subtract 1 to account for hidden bit.
auto offset = significand_size - double_significand_size - SHIFT - 1;
f <<= offset;
e -= offset;
}
// Compute lower and upper boundaries (m^- and m^+ in the Grisu paper), where
// a boundary is a value half way between the number and its predecessor
// (lower) or successor (upper). The upper boundary is normalized and lower
// has the same exponent but may be not normalized.
void compute_boundaries(fp &lower, fp &upper) const {
lower = f == implicit_bit ?
fp((f << 2) - 1, e - 2) : fp((f << 1) - 1, e - 1);
upper = fp((f << 1) + 1, e - 1);
upper.normalize<1>(); // 1 is to account for the exponent shift above.
lower.f <<= lower.e - upper.e;
lower.e = upper.e;
}
};
// Returns an fp number representing x - y. Result may not be normalized.
inline fp operator-(fp x, fp y) {
FMT_ASSERT(x.f >= y.f && x.e == y.e, "invalid operands");
return fp(x.f - y.f, x.e);
inline bool use_grisu() {
return FMT_USE_GRISU && std::numeric_limits<double>::is_iec559;
}
// Computes an fp number r with r.f = x.f * y.f / pow(2, 64) rounded to nearest
// with half-up tie breaking, r.e = x.e + y.e + 64. Result may not be normalized.
FMT_API fp operator*(fp x, fp y);
// Returns cached power (of 10) c_k = c_k.f * pow(2, c_k.e) such that its
// (binary) exponent satisfies min_exponent <= c_k.e <= min_exponent + 3.
FMT_API fp get_cached_power(int min_exponent, int &pow10_exponent);
// Formats value using Grisu2 algorithm:
// https://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf
FMT_API void grisu2_format(double value, char *buffer, size_t &size, char type,
@ -2948,13 +2862,14 @@ void basic_writer<Range>::write_double(T value, const format_specs &spec) {
return write_inf_or_nan(handler.upper ? "INF" : "inf");
basic_memory_buffer<char_type> buffer;
if (internal::const_check(FMT_USE_GRISU && sizeof(T) <= sizeof(double) &&
std::numeric_limits<double>::is_iec559)) {
char type = static_cast<char>(spec.type());
if (internal::const_check(
internal::use_grisu() && sizeof(T) <= sizeof(double)) &&
type != 'a' && type != 'A') {
char buf[100]; // TODO: correct buffer size
size_t size = 0;
internal::grisu2_format(
static_cast<double>(value), buf, size, static_cast<char>(spec.type()),
spec.precision(), spec.flag(HASH_FLAG));
internal::grisu2_format(static_cast<double>(value), buf, size, type,
spec.precision(), spec.flag(HASH_FLAG));
FMT_ASSERT(size <= 100, "buffer overflow");
buffer.append(buf, buf + size); // TODO: avoid extra copy
} else {

View File

@ -24,6 +24,78 @@
#undef min
#undef max
using fmt::internal::fp;
template <bool is_iec559>
void test_construct_from_double() {
fmt::print("warning: double is not IEC559, skipping FP tests\n");
}
template <>
void test_construct_from_double<true>() {
auto v = fp(1.23);
EXPECT_EQ(v.f, 0x13ae147ae147aeu);
EXPECT_EQ(v.e, -52);
}
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, ComputeBoundariesSubnormal) {
auto v = fp(0xbeef, 42);
fp lower, upper;
v.compute_boundaries(lower, upper);
EXPECT_EQ(0xbeee800000000000, lower.f);
EXPECT_EQ(-6, lower.e);
EXPECT_EQ(0xbeef800000000000, upper.f);
EXPECT_EQ(-6, upper.e);
}
TEST(FPTest, ComputeBoundaries) {
auto v = fp(0x10000000000000, 42);
fp lower, upper;
v.compute_boundaries(lower, upper);
EXPECT_EQ(0x7ffffffffffffe00, lower.f);
EXPECT_EQ(31, lower.e);
EXPECT_EQ(0x8000000000000400, upper.f);
EXPECT_EQ(31, upper.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);
}
TEST(FPTest, GetCachedPower) {
typedef std::numeric_limits<double> limits;
for (auto exp = limits::min_exponent; exp <= limits::max_exponent; ++exp) {
int dec_exp = 0;
auto fp = fmt::internal::get_cached_power(exp, dec_exp);
EXPECT_LE(exp, fp.e);
int dec_exp_step = 8;
EXPECT_LE(fp.e, exp + dec_exp_step * log2(10));
EXPECT_DOUBLE_EQ(pow(10, dec_exp), ldexp(fp.f, fp.e));
}
}
template <typename T>
struct ValueExtractor: fmt::internal::function<T> {
T operator()(T value) {

View File

@ -31,77 +31,6 @@ using fmt::format_error;
using fmt::string_view;
using fmt::memory_buffer;
using fmt::wmemory_buffer;
using fmt::internal::fp;
template <bool is_iec559>
void test_construct_from_double() {
fmt::print("warning: double is not IEC559, skipping FP tests\n");
}
template <>
void test_construct_from_double<true>() {
auto v = fp(1.23);
EXPECT_EQ(v.f, 0x13ae147ae147aeu);
EXPECT_EQ(v.e, -52);
}
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, ComputeBoundariesSubnormal) {
auto v = fp(0xbeef, 42);
fp lower, upper;
v.compute_boundaries(lower, upper);
EXPECT_EQ(0xbeee800000000000, lower.f);
EXPECT_EQ(-6, lower.e);
EXPECT_EQ(0xbeef800000000000, upper.f);
EXPECT_EQ(-6, upper.e);
}
TEST(FPTest, ComputeBoundaries) {
auto v = fp(0x10000000000000, 42);
fp lower, upper;
v.compute_boundaries(lower, upper);
EXPECT_EQ(0x7ffffffffffffe00, lower.f);
EXPECT_EQ(31, lower.e);
EXPECT_EQ(0x8000000000000400, upper.f);
EXPECT_EQ(31, upper.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);
}
TEST(FPTest, GetCachedPower) {
typedef std::numeric_limits<double> limits;
for (auto exp = limits::min_exponent; exp <= limits::max_exponent; ++exp) {
int dec_exp = 0;
auto fp = fmt::internal::get_cached_power(exp, dec_exp);
EXPECT_LE(exp, fp.e);
int dec_exp_step = 8;
EXPECT_LE(fp.e, exp + dec_exp_step * log2(10));
EXPECT_DOUBLE_EQ(pow(10, dec_exp), ldexp(fp.f, fp.e));
}
}
namespace {
@ -626,7 +555,10 @@ TEST(FormatterTest, HashFlag) {
EXPECT_EQ("0x42", format("{0:#x}", 0x42ull));
EXPECT_EQ("042", format("{0:#o}", 042ull));
EXPECT_EQ("-42.0000", format("{0:#}", -42.0));
if (fmt::internal::use_grisu())
EXPECT_EQ("-42.0", format("{0:#}", -42.0));
else
EXPECT_EQ("-42.0000", format("{0:#}", -42.0));
EXPECT_EQ("-42.0000", format("{0:#}", -42.0l));
EXPECT_THROW_MSG(format("{0:#", 'c'),
format_error, "missing '}' in format string");