[libc-commits] [libc] [libc] Alternative algorithm for decimal FP printf (PR #123643)
Simon Tatham via libc-commits
libc-commits at lists.llvm.org
Wed Jan 22 09:01:44 PST 2025
================
@@ -0,0 +1,593 @@
+//===-- Decimal Float Converter for printf (320-bit float) ------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+//
+// This file implements an alternative to the Ryū printf algorithm in
+// float_dec_converter.h. Instead of generating output digits 9 at a time on
+// demand, in this implementation, a float is converted to decimal by computing
+// just one power of 10 and multiplying/dividing the entire input by it,
+// generating the whole string of decimal output digits in one go.
+//
+// This avoids the large constant lookup table of Ryū, making it more suitable
+// for low-memory embedded contexts; but it's also faster than the fallback
+// version of Ryū which computes table entries on demand using DyadicFloat,
+// because those must calculate a potentially large power of 10 per 9-digit
+// output block, whereas this computes just one, which does the whole job.
+//
+// (OK, not quite. It computes one power of 10 per _attempt_. It must start by
+// guessing the decimal exponent of the output, and if it guesses too high or
+// too low then it adjusts the guess and tries again from scratch.)
+//
+// The calculation is done in 320-bit DyadicFloat, which provides enough
+// precision to generate 39 correct digits of output from any floating-point
+// size up to and including 128-bit long double, because the rounding errors in
+// computing the largest necessary power of 10 are still smaller than the
+// distance (in the 320-bit float format) between adjacent 39-decimal-digit
+// outputs.
+//
+// No further digits beyond the 39th are generated: if the printf format string
+// asks for more precision than that, the answer is padded with 0s. This is a
+// permitted option in IEEE 754-2019 (section 5.12.2): you're allowed to define
+// a limit H on the number of decimal digits you can generate, and pad with 0s
+// if asked for more than that, subject to the constraint that H must be
+// consistent across all float formats you support (you can't use a smaller H
+// for single precision than double or long double), and must be large enough
+// that even in the largest supported precision the only numbers misrounded are
+// ones extremely close to a rounding boundary. 39 digits is the smallest
+// permitted value for an implementation supporting binary128.
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_STDIO_PRINTF_CORE_FLOAT_DEC_CONVERTER_LIMITED_H
+#define LLVM_LIBC_SRC_STDIO_PRINTF_CORE_FLOAT_DEC_CONVERTER_LIMITED_H
+
+#include "src/__support/CPP/algorithm.h"
+#include "src/__support/CPP/string.h"
+#include "src/__support/CPP/string_view.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/dyadic_float.h"
+#include "src/__support/FPUtil/rounding_mode.h"
+#include "src/__support/integer_to_string.h"
+#include "src/__support/macros/config.h"
+#include "src/stdio/printf_core/core_structs.h"
+#include "src/stdio/printf_core/float_inf_nan_converter.h"
+#include "src/stdio/printf_core/writer.h"
+
+#include <inttypes.h>
+#include <stddef.h>
+#include <string.h>
+
+namespace LIBC_NAMESPACE_DECL {
+namespace printf_core {
+
+enum class ConversionType { E, F, G };
+using StorageType = fputil::FPBits<long double>::StorageType;
+
+constexpr unsigned MAX_DIGITS = 39;
+constexpr size_t DF_BITS = 320;
+constexpr char DECIMAL_POINT = '.';
+
+struct DigitsInput {
+ // Input mantissa, stored with the explicit leading 1 bit (if any) at the
+ // top. So either it has a value in the range [2^127,2^128) representing a
+ // real number in [1,2), or it has the value 0, representing 0.
+ UInt128 mantissa;
+
+ // Input exponent, as a power of 2 to multiply into mantissa.
+ int exponent;
+
+ // Input sign.
+ Sign sign;
+
+ // Constructor which accepts a mantissa direct from a floating-point format,
+ // and shifts it up to the top of the UInt128 so that a function consuming
+ // this struct afterwards doesn't have to remember which format it came from.
+ DigitsInput(int32_t fraction_len, StorageType mantissa_, int exponent_,
+ Sign sign)
+ : mantissa(UInt128(mantissa_) << (127 - fraction_len)),
+ exponent(exponent_), sign(sign) {
+ if (!(mantissa & (UInt128(1) << 127)) && mantissa != 0) {
+ // Normalize a denormalized input.
+ int shift = cpp::countl_zero(mantissa);
+ mantissa <<= shift;
+ exponent -= shift;
+ }
+ }
+};
+
+struct DigitsOutput {
+ // Output from decimal_digits().
+ //
+ // `digits` is a buffer containing nothing but ASCII digits. Even if the
+ // decimal point needs to appear somewhere in the final output string, it
+ // isn't represented in _this_ string; the client of this object will insert
+ // it in an appropriate place. `ndigits` gives the buffer size.
+ //
+ // `exponent` represents the exponent you would display if the decimal point
+ // comes after the first digit of decimal_digits, e.g. if digits == "1234"
+ // and exponent = 3 then this represents 1.234e3, or just the integer 1234.
+ size_t ndigits;
+ int exponent;
+ char digits[MAX_DIGITS + 1];
+};
+
+// Calculate the actual digits of a decimal representation of an FP number.
+//
+// If `e_mode` is true, then `precision` indicates the desired number of output
+// decimal digits. On return, `decimal_digits` will be a string of length
+// exactly `precision` starting with a nonzero digit; `decimal_exponent` will
+// be filled in to indicate the exponent as shown above.
+//
+// If `e_mode` is false, then `precision` indicates the desired number of
+// digits after the decimal point. On return, the last digit in the string
+// `decimal_digits` has a place value of _at least_ 10^-precision. But also, at
+// most `MAX_DIGITS` digits are returned, so the caller may need to pad it at
+// the end with the appropriate number of extra 0s.
+LIBC_INLINE
+DigitsOutput decimal_digits(DigitsInput input, int precision, bool e_mode) {
+ if (input.mantissa == 0) {
+ // Special-case zero, by manually generating the right number of zero
+ // digits and setting an appropriate exponent.
+ DigitsOutput output;
+ if (!e_mode) {
+ // In F mode, it's enough to return an empty string of digits. That's the
+ // same thing we do when given a nonzero number that rounds down to 0.
+ output.ndigits = 0;
+ output.exponent = -precision - 1;
+ } else {
+ // In E mode, generate a string containing the expected number of 0s.
+ memset(output.digits, '0', precision);
+ output.ndigits = precision;
+ output.exponent = 0;
+ }
+ return output;
+ }
+
+ // Estimate log10 of the input value, by multiplying its binary exponent by
+ // 19728/65536, which is a good approximation to log10(2).
+ //
+ // If this turns out to be wrong, we'll find out when we get the wrong number
+ // of digits from the conversion in the loop below, and we can adjust the
+ // value and try again.
+ int log10_input = ((input.exponent - 1) * 19728) >> 16;
+
+ // Make a DyadicFloat containing the value 10, to use as the base for
+ // exponentation inside the following loop, potentially more than once if we
+ // need to iterate.
+ fputil::DyadicFloat<DF_BITS> ten(Sign::POS, 1, 5);
+
+ while (true) {
+ // Compute the exponent of the lowest-order digit we want as output. In F
+ // mode this depends only on the desired precision. In E mode it's based on
+ // log10_input, which is (an estimate of) the exponent corresponding to the
+ // _high_-order decimal digit of the number.
+ int log10_low_digit = e_mode ? log10_input + 1 - precision : -precision;
+
+ // Calculate an integer whose decimal representation is precisely the
+ // string of output digits, by doing a DyadicFloat computation of
+ // (input_mantissa / 10^(log10_low_digit)) and then rounding that to an
+ // integer.
+ //
+ // If log10_low_digit < 0 then we calculate 10^(-log10_low_digit) and
+ // multiply by it instead, so that the exponent is non-negative in all
+ // cases. This ensures that the power of 10 is always mathematically
+ // speaking an integer, so that it can be represented exactly in binary
+ // (without a recurring fraction), and when it's small enough to fit in
+ // DF_BITS, fputil::pow_n should return the exact answer, and then
+ // fputil::rounded_{div,mul} will introduce only the unavoidable rounding
+ // error of up to 1/2 ULP.
+ //
+ // Beyond that point, pow_n will be imprecise. But DF_BITS is set high
+ // enough that even for the most difficult cases in 128-bit long double,
+ // the extra precision in the calculation is enough to ensure we still get
+ // the right answer.
+ //
+ // If the output integer doesn't fit in DF_BITS, we set the `overflow`
+ // flag.
+ UInt<DF_BITS> integer;
+ bool overflow;
+ {
+ // Calculate the power of 10 to divide or multiply by.
+ fputil::DyadicFloat<DF_BITS> power_of_10 =
+ fputil::pow_n(ten, cpp::abs(log10_low_digit));
+
+ // Convert the mantissa into a DyadicFloat, making sure it has the right
+ // sign, so that directed rounding will go in the right direction, if
+ // enabled.
+ fputil::DyadicFloat<DF_BITS> flt_mantissa(
+ input.sign, input.exponent - 127, input.mantissa);
+
+ // Divide or multiply, depending on whether log10_low_digit was positive
+ // or negative.
+ fputil::DyadicFloat<DF_BITS> flt_quotient =
+ log10_low_digit > 0 ? fputil::rounded_div(flt_mantissa, power_of_10)
+ : fputil::rounded_mul(flt_mantissa, power_of_10);
+
+ // Convert to an integer.
+ integer = flt_quotient.as_mantissa_type_rounded(&overflow);
+
+ // And take the absolute value.
+ if (flt_quotient.sign.is_neg())
+ integer = -integer;
+ }
+
+ // If the conversion to an integer overflowed, then we certainly can't use
+ // it.
+ if (overflow) {
+ if (e_mode) {
+ // In E mode, just try again with the exponent one higher.
+ log10_input++;
+ } else {
+ // In F mode, if the number would generate too many digits to fit, then
+ // reset to E mode, so that we generate the most significant digits of
+ // the number, and rely on our caller to append zeroes at the end.
+ precision = MAX_DIGITS;
+ e_mode = true;
----------------
statham-arm wrote:
I don't think that formula can be quite right, because of the thing I mention in the comment above, where some binary exponents can go with two possible decimal exponents depending on the mantissa, because a power of 10 lies in their interval. So any test that doesn't even look at the mantissa _must_ be potentially inaccurate.
But I'm sure you're right that _most_ cases of "too many output digits, go to e mode" can be detected before the conversion to decimal digits, and we can narrow it down so that the remaining ones only generate one digit too many, which can be dealt with by re-rounding the same float320 output.
https://github.com/llvm/llvm-project/pull/123643
More information about the libc-commits
mailing list