[libc-commits] [libc] [libc] Alternative algorithm for decimal FP printf (PR #123643)
Simon Tatham via libc-commits
libc-commits at lists.llvm.org
Tue Jan 28 01:58:58 PST 2025
================
@@ -0,0 +1,669 @@
+//===-- 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.
+//
+// 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"
+
+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.
+ __builtin_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
+ // 1292913986/2^32. That is a rounded-down approximation to log10(2),
+ // accurate enough that for any binary exponent in the range of float128 it
+ // will give the correct value of floor(log10(2^n)).
+ //
+ // This estimate is correct for deciding how many decimal digits we end up
+ // producing, unless a power of 10 falls in the interval corresponding to
+ // this binary exponent, in which case there might be one more decimal digit
+ // for larger mantissas. To detect this, we do the same computation for the
+ // next exponent up.
+ int log10_input_min = ((input.exponent - 1) * 1292913986LL) >> 32;
+ int log10_input_max = (input.exponent * 1292913986LL) >> 32;
+
+ // 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);
+
+ // 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_min + 1 - precision : -precision;
+
+ // The general plan is to 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.
+ //
+ // The number of output decimal digits (if the mathematical result of this
+ // operation were computed without overflow) will be one of these:
+ // (log10_input_min - log10_low_digit + 1)
+ // (log10_input_max - log10_low_digit + 1)
+ //
+ // In E mode, this means we'll either get the correct number of output digits
+ // immediately, or else one too many (in which case we can correct for that
+ // at the rounding stage). But in F mode, if the number is very large
+ // compared to the number of decimal places the user asked for, we might be
+ // about to generate far too many digits and overflow our float format. In
+ // that case, reset to E mode immediately, to avoid having to detect the
+ // overflow _after_ the multiplication and retry. So if even the smaller
+ // number of possible output digits is too many, we might as well change our
+ // mind right now and switch into E mode.
+ if (log10_input_max - log10_low_digit + 1 > MAX_DIGITS) {
+ precision = MAX_DIGITS;
+ e_mode = true;
+ log10_low_digit = log10_input_min + 1 - precision;
+ }
+
+ // Now actually calculate (input_mantissa / 10^(log10_low_digit)).
+ //
+ // 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.
+
+ // 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.
+ int round_dir;
+ UInt<DF_BITS> integer = flt_quotient.as_mantissa_type_rounded(&round_dir);
+
+ // And take the absolute value.
+ if (flt_quotient.sign.is_neg())
+ integer = -integer;
+
+ // Convert the mantissa integer into a string of decimal digits, and check
+ // to see if it's the right size.
+ const IntegerToString<decltype(integer), radix::Dec> buf{integer};
+ cpp::string_view view = buf.view();
+
+ // Start making the output struct, by copying in the digits from the above
+ // object. At this stage we may also have one digit too many (but that's OK,
+ // there's space for it in the DigitsOutput buffer).
+ DigitsOutput output;
+ output.ndigits = view.size();
+ __builtin_memcpy(output.digits, view.data(), output.ndigits);
+
+ // Set up the output exponent, which is done differently depending on mode.
+ // Also, figure out whether we have one digit too many, and if so, set the
+ // `need_reround` flag and adjust the exponent appropriately.
+ bool need_reround = false;
+ if (e_mode) {
+ // In E mode, the output exponent is the exponent of the first decimal
+ // digit, which we already calculated.
+ output.exponent = log10_input_min;
+
+ // In E mode, we're returning a fixed number of digits, given by
+ // `precision`, so if we have more than that, then we must shorten the
+ // buffer by one digit.
+ //
+ // If this happens, it's because the actual log10 of the input is
+ // log10_input_min + 1. Equivalently, we guessed we'd see something like
+ // X.YZe+NN and instead got WX.YZe+NN. So when we shorten the digit string
+ // by one, we'll also need to increment the output exponent.
+ if (output.ndigits > size_t(precision)) {
+ assert(output.ndigits == size_t(precision) + 1);
+ need_reround = true;
+ output.exponent++;
+ }
+ } else {
+ // In F mode, the output exponent is based on the place value of the _last_
+ // digit, so we must recover the exponent of the first digit by adding
+ // the number of digits.
+ //
+ // Because this takes the length of the buffer into account, it sets the
+ // correct decimal exponent even if this digit string is one too long. So
+ // we don't need to adjust the exponent if we reround.
+ output.exponent = int(output.ndigits) - precision - 1;
+
+ // In F mode, the number of returned digits isn't based on `precision`:
+ // it's variable, and we don't mind how many digits we get as long as it
+ // isn't beyond the limit MAX_DIGITS. If it is, we expect that it's only
+ // one digit too long, or else we'd have spotted the problem in advance and
+ // flipped into E mode already.
+ if (output.ndigits > MAX_DIGITS) {
+ assert(output.ndigits == MAX_DIGITS + 1);
+ need_reround = true;
+ }
+ }
+
+ if (need_reround) {
+ // If either of the branches above decided that we had one digit too many,
+ // we must now shorten the digit buffer by one. But we can't just truncate:
+ // we need to make sure the remaining n-1 digits are correctly rounded, as
+ // if we'd rounded just once from the original `flt_quotient`.
+ //
+ // In directed rounding modes this can't go wrong. If you had a real number
+ // x, and the first rounding produced floor(x), then the second rounding
+ // wants floor(x/10), and it doesn't matter if you actually compute
+ // floor(floor(x)/10): the result is the same, because each rounding
+ // boundary in the second rounding aligns with one in the first rounding,
+ // which nothing could have crossed. Similarly for rounding away from zero,
+ // with 'floor' replaced with 'ceil' throughout.
+ //
+ // In rounding to nearest, the danger is in the boundary case where the
+ // final digit of the original output is 5. Then if we just rerounded the
+ // digit string to remove the last digit, it would look like an exact
+ // halfway case, and we'd break the tie by choosing the even one of the two
+ // outputs. But if the original value before the first rounding was on one
+ // side or the other of 5, then that supersedes the 'round to even' tie
+ // break. So we need to consult `round_dir` from above, which tells us
+ // which way (if either) the value was adjusted during the first rounding.
+ // Effectively, we treat the last digit as 5+ε or 5-ε.
+ //
+ // To make this work in both directed modes and round-to-nearest mode
+ // without having to look up the rounding direction, a simple rule is: take
+ // account of round_dir if and only if the round digit (the one we're
+ // removing when shortening the buffer) is 5. In directed rounding modes
+ // this makes no difference.
+
+ // Extract the two relevant digits. round_digit is the one we're removing;
+ // new_low_digit is the last one we're keeping, so we need to know if it's
+ // even or odd to handle exact tie cases (when round_dir == 0).
+ int round_digit = output.digits[--output.ndigits] - '0';
+ int new_low_digit =
+ output.ndigits == 0 ? 0 : output.digits[output.ndigits - 1] - '0';
+
+ // Make a binary number that we can pass to `fputil::rounding_direction`.
+ // We put new_low_digit at bit 8, and imagine that we're rounding away the
+ // bottom 8 bits. Therefore round_digit must be "just below" bit 8, in the
+ // sense that we set the bottom 8 bits to (256/10 * round_digit) so that
+ // round_digit==5 corresponds to the binary half-way case of 0x80.
+ //
+ // Then we adjust by +1 or -1 based on round_dir if the round digit is 5,
+ // as described above.
+ LIBC_NAMESPACE::UInt<64> round_word = (new_low_digit * 256) +
+ ((round_digit * 0x19a) >> 4) +
+ (round_digit == 5 ? -round_dir : 0);
+
+ // Now we can call the existing binary rounding helper function, which
+ // takes account of the rounding mode.
+ if (fputil::rounding_direction(round_word, 8, flt_quotient.sign) > 0) {
+ // If that returned a positive answer, we must round the number up.
+ //
+ // The number is already in decimal, so we need to increment it one digit
+ // at a time. (A bit painful, but better than going back to the integer
+ // we made it from and doing the decimal conversion all over again.)
+ for (size_t i = output.ndigits; i-- > 0;) {
+ if (output.digits[i]++ != '9')
----------------
statham-arm wrote:
I'm confused – surely here I'm depending on exactly the same assumption about character encoding that I did in the earlier `char_value - '0'` expression? Namely, that the decimal digit characters appear consecutively, so that I can increment one to the next digit by writing `output.digits[i]++`.
Neither piece of code depends on the _specific_ ASCII encodings of the decimal digits. But both depend on them being consecutive. So either you've misunderstood, or I have: why is one of those OK and the other not?
https://github.com/llvm/llvm-project/pull/123643
More information about the libc-commits
mailing list