stdlib/public/runtime/SwiftDtoa.cpp

//===--- SwiftDtoa.c ---------------------------------------------*- c -*-===//
//
// This source file is part of the Swift.org open source project
//
// Copyright (c) 2018 Apple Inc. and the Swift project authors
// Licensed under Apache License v2.0 with Runtime Library Exception
//
// See https://swift.org/LICENSE.txt for license information
// See https://swift.org/CONTRIBUTORS.txt for the list of Swift project authors
//
//===---------------------------------------------------------------------===//
//
// Note: This is really a C file, but Swift's build system for Linux is
// partially allergic to C, so it's being compiled as ".cpp" for now.  Please
// don't infect it with C++-isms.
//
///
/// The core algorithm here (see `swift_decompose_double` below) is a
/// modified form of the Grisu2 algorithm from Florian Loitsch;
/// "Printing Floating-Point Numbers Quickly and Accurately with
/// Integers", 2010.  https://dl.acm.org/citation.cfm?id=1806623
///
/// This includes some improvements suggested by the "Errol paper":
/// Marc Andrysco, Ranjit Jhala, Sorin Lerner; "Printing
/// Floating-Point Numbers: A Faster, Always Correct Method", 2016.
/// https://dl.acm.org/citation.cfm?id=2837654
///
/// The following summary assumes you're familiar with Grisu-style
/// algorithms in general:
///
/// Loitsch' original Grisu2 implementation guarantees round-trip
/// accuracy but only generates the shortest decimal expansion about 99%
/// of the time.  Grisu3 is similar, but fails rather than producing
/// a result that is not the shortest possible.
///
/// The Errol paper provides a deeper analysis of the cases where
/// Grisu2 fails to find the shortest decimal expansion. There
/// are two root causes of such failures:
///
/// * Insufficient precision leads to scattered failures across the
///   entire range.  The enumeration technique described in the Errol
///   paper shows a way to construct a superset of the numbers subject
///   to such failures.  With this list, we can simply test whether we
///   have sufficient precision.
///
///   For Double, the Errol3 algorithm uses double-double arithmetic
///   with about 106 bits precision. This turns out to be not quite
///   sufficient, requiring Errol3 to pre-screen the input against a
///   list of exceptions culled from the larger list of possible
///   failures.  Using high-precision integers, we've discovered that
///   110 bit precision is sufficient to satisfy the Errol test cases
///   without requiring any pre-screening.
///
///   For Float and Float80, the same approach shows that we need 53
///   and 135 bits, respectively.  It is an interesting coincidence
///   that for all three cases, an n-bit significand can be formatted
///   optimally with no more than 2n+7 bits of intermediate precision.
///
/// * Sometimes, the shortest value might occur exactly at the
///   midpoint between two adjacent binary floating-point values.
///   When converted back to binary, this will round to the adjacent
///   even significand.  We handle this by widening the interval
///   whenever the significand is even in order to allow these
///   exact midpoints to be considered.
///
/// In addition to addressing the shortness failures characterized in the Errol
/// paper, the implementation here also incorporates final-digit corrections
/// that allow it to produce the optimal decimal decomposition in all cases.
///
/// In summary, this implementation is:
///
/// * Fast.  It uses only fixed-width integer arithmetic and has
///   constant memory requirements.
///
/// * Simple. It is only a little more complex than Loitsch' original
///   implementation of Grisu2.  The full digit decomposition for double
///   is less than 300 lines of standard C, including routine arithmetic
///   helper functions.
///
/// * Always Accurate. Converting the decimal form back to binary
///   will always yield exactly the same value.  For the IEEE 754
///   formats, the round-trip will produce exactly the same bit
///   pattern in memory.
///
/// * Always Short.  This always selects an accurate result with the
///   minimum number of decimal digits.
///
/// * Always Close.  Among all accurate, short results, this always
///   chooses the result that is closest to the exact floating-point
///   value. (In case of an exact tie, it rounds the last digit even.)
///
/// For single-precision Float, we can simply test all 2^32 values.
/// This requires only a few minutes on a mid-range modern laptop. The
/// Double and Float80 formatters rely on results from the Errol paper
/// to ensure correctness.  In addition, we have verified more than
/// 10^14 randomly-chosen Double values by comparing the results to the
/// output of Grisu3 (where Grisu3 fails, we've compared to Errol4).
///
// ----------------------------------------------------------------------------

#include <float.h>
#include <inttypes.h>
#include <math.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "swift/Runtime/SwiftDtoa.h"

#if defined(__SIZEOF_INT128__)
  // We get a significant speed boost if we can use the __uint128_t
  // type that's present in GCC and Clang on 64-bit architectures.  In
  // particular, we can do 128-bit arithmetic directly and can
  // represent 192-bit integers as a collection of 64-bit elements.
  #define HAVE_UINT128_T 1
#else
  // On 32-bit, we have to use slower code that manipulates 128-bit
  // and 192-bit integers as collections of 32-bit elements.
  #define HAVE_UINT128_T 0
#endif

// Try to verify that the system floating-point types really are what we
// expect.  Note that the code below is specific to these exact
// floating-point representations.

#if (FLT_RADIX != 2)
// Either you're using hexadecimal float format on S390, or you have a
// really broken C environment.
#error "This platform claims to not use binary floating point."
#endif

#if SWIFT_DTOA_FLOAT_SUPPORT
#if (FLT_MANT_DIG != 24) || (FLT_MIN_EXP != -125) || (FLT_MAX_EXP != 128)
#error "Are you certain `float` on this platform is really IEEE 754 single-precision binary32 format?"
#endif
#endif

#if SWIFT_DTOA_DOUBLE_SUPPORT
#if (DBL_MANT_DIG != 53) || (DBL_MIN_EXP != -1021) || (DBL_MAX_EXP != 1024)
#error "Are you certain `double` on this platform is really IEEE 754 double-precision binary64 format?"
#endif
#endif

#if SWIFT_DTOA_FLOAT80_SUPPORT
#if (LDBL_MANT_DIG != 64) || (LDBL_MIN_EXP != -16381) || (LDBL_MAX_EXP != 16384)
#error "Are you certain `long double` on this platform is really Intel 80-bit extended precision?"
#endif
#endif


// See the implementations at the bottom of this file for detailed explanations
// of the purpose of each function.

//
// Helper functions used by float, double, and float80 machinery.
//

#if SWIFT_DTOA_FLOAT_SUPPORT || SWIFT_DTOA_DOUBLE_SUPPORT || SWIFT_DTOA_FLOAT80_SUPPORT
static int binaryExponentFor10ToThe(int p);
static int decimalExponentFor2ToThe(int e);
#endif

//
// Helper functions used only by the single-precision float formatter
//

#if SWIFT_DTOA_FLOAT_SUPPORT
static uint64_t multiply64x32RoundingDown(uint64_t lhs, uint32_t rhs);
static uint64_t multiply64x32RoundingUp(uint64_t lhs, uint32_t rhs);
static uint64_t multiply64x64RoundingDown(uint64_t lhs, uint64_t rhs);
static uint64_t multiply64x64RoundingUp(uint64_t lhs, uint64_t rhs);
static void intervalContainingPowerOf10_Float(int p, uint64_t *lower, uint64_t *upper, int *exponent);
#endif

//
// Helpers used by both the double-precision and float80 formatters
//

#if SWIFT_DTOA_DOUBLE_SUPPORT || SWIFT_DTOA_FLOAT80_SUPPORT
#if HAVE_UINT128_T
typedef __uint128_t swift_uint128_t;
#define initialize128WithHighLow64(dest, high64, low64) ((dest) = ((__uint128_t)(high64) << 64) | (low64))
#else
typedef struct {
    uint32_t low, b, c, high;
} swift_uint128_t;
#define initialize128WithHighLow64(dest, high64, low64) \
    ((dest).low = (uint32_t)(low64),                    \
     (dest).b = (uint32_t)((low64) >> 32),              \
     (dest).c = (uint32_t)(high64),                     \
     (dest).high = (uint32_t)((high64) >> 32))
#endif
#endif


//
// Helper functions used only by the double-precision formatter
//

#if SWIFT_DTOA_DOUBLE_SUPPORT
#if HAVE_UINT128_T
#define increment128(dest) ((dest) += 1)
#define isLessThan128x128(lhs, rhs) ((lhs) < (rhs))
#define subtract128x128(lhs, rhs) (*(lhs) -= (rhs))
#define multiply128xi32(lhs, rhs) (*(lhs) *= (rhs))
#define initialize128WithHigh64(dest, value) ((dest) = (__uint128_t)(value) << 64)
#define extractHigh64From128(arg) ((uint64_t)((arg) >> 64))
static int extractIntegerPart128(__uint128_t *fixed128, int fractionBits) {
    return (int)(*fixed128 >> fractionBits);
}
static void clearIntegerPart128(__uint128_t *fixed128, int fractionBits) {
    const swift_uint128_t fixedPointMask = (((__uint128_t)1 << fractionBits) - 1);
    *fixed128 &= fixedPointMask;
}
#else
#define increment128(dest)             \
    do {                               \
        uint64_t t = (dest).low + 1;   \
        (dest).low = (uint32_t)t;      \
        t >>= 32;                      \
        t += (dest).b;                 \
        (dest).b = (uint32_t)t;        \
        t >>= 32;                      \
        t += (dest).c;                 \
        (dest).c = (uint32_t)t;        \
        t >>= 32;                      \
        (dest).high += (uint32_t)t;    \
    } while (0)
static int isLessThan128x128(swift_uint128_t lhs, swift_uint128_t rhs);
static void subtract128x128(swift_uint128_t *lhs, swift_uint128_t rhs);
static void multiply128xi32(swift_uint128_t *lhs, uint32_t rhs);
#define initialize128WithHigh64(dest, value)            \
    ((dest).low = (dest).b = 0,                         \
     (dest).c = (uint32_t)(value),                      \
     (dest).high = (uint32_t)((value) >> 32))
#define extractHigh64From128(arg) (((uint64_t)(arg).high << 32)|((arg).c))
static int extractIntegerPart128(swift_uint128_t *fixed128, int fractionBits) {
    const int highFractionBits = fractionBits % 32;
    return (int)(fixed128->high >> highFractionBits);
}
static void clearIntegerPart128(swift_uint128_t *fixed128, int fractionBits) {
    const int highFractionBits = fractionBits % 32;
    fixed128->high &= ((uint32_t)1 << highFractionBits) - 1;
}
#endif
static swift_uint128_t multiply128x64RoundingDown(swift_uint128_t lhs, uint64_t rhs);
static swift_uint128_t multiply128x64RoundingUp(swift_uint128_t lhs, uint64_t rhs);
static swift_uint128_t shiftRightRoundingDown128(swift_uint128_t lhs, int shift);
static swift_uint128_t shiftRightRoundingUp128(swift_uint128_t lhs, int shift);
static void intervalContainingPowerOf10_Double(int p, swift_uint128_t *lower, swift_uint128_t *upper, int *exponent);
#endif

//
// Helper functions used only by the extended-precision long double formatter
//

#if SWIFT_DTOA_FLOAT80_SUPPORT
#if HAVE_UINT128_T
// A 192-bit unsigned integer type stored as 3 64-bit words
typedef struct {uint64_t low, mid, high;} swift_uint192_t;
#define initialize192WithHighMidLow64(dest, high64, mid64, low64) \
    ((dest).low = (low64),                                        \
     (dest).mid = (mid64),                                        \
     (dest).high = (high64))
#else
// A 192-bit unsigned integer type stored as 6 32-bit words
typedef struct {uint32_t low, b, c, d, e, high;} swift_uint192_t;
#define initialize192WithHighMidLow64(dest, high64, mid64, low64) \
    ((dest).low = (uint64_t)(low64),                              \
     (dest).b = (uint64_t)(low64) >> 32,                          \
     (dest).c = (uint64_t)(mid64),                                \
     (dest).d = (uint64_t)(mid64) >> 32,                          \
     (dest).e = (uint64_t)(high64),                               \
     (dest).high = (uint64_t)(high64) >> 32)
#endif
static void multiply192x64RoundingDown(swift_uint192_t *lhs, uint64_t rhs);
static void multiply192x64RoundingUp(swift_uint192_t *lhs, uint64_t rhs);
static void multiply192xi32(swift_uint192_t *lhs, uint32_t rhs);
static void multiply192x128RoundingDown(swift_uint192_t *lhs, swift_uint128_t rhs);
static void multiply192x128RoundingUp(swift_uint192_t *lhs, swift_uint128_t rhs);
static void subtract192x192(swift_uint192_t *lhs, swift_uint192_t rhs);
static int isLessThan192x192(swift_uint192_t lhs, swift_uint192_t rhs);
static void shiftRightRoundingDown192(swift_uint192_t *lhs, int shift);
static void shiftRightRoundingUp192(swift_uint192_t *lhs, int shift);
static void intervalContainingPowerOf10_Float80(int p, swift_uint192_t *lower, swift_uint192_t *upper, int *exponent);
#endif

//
// --------------- Digit generation ---------------------
//

// This is the interesting part.

// These routines take a floating-point value and efficiently compute
// everything necessary to write an optimal base-10 representation of
// that value.  In particular, they compute the base-10 exponent and
// corresponding digits.

// swift_decompose_double is thoroughly commented; swift_decompose_float
// and swift_decompose_float80 are fundamentally the same algorithm, but
// adjusted to perform optimally for those types.

#if SWIFT_DTOA_DOUBLE_SUPPORT
// Return raw bits encoding the double
static uint64_t bitPatternForDouble(double d) {
    union { double d; uint64_t u; } converter;
    converter.d = d;
    return converter.u;
}

int swift_decompose_double(double d,
    int8_t *digits, size_t digits_length, int *decimalExponent)
{
    // Bits in raw significand (not including hidden bit, if present)
    static const int significandBitCount = DBL_MANT_DIG - 1;
    static const uint64_t significandMask
        = ((uint64_t)1 << significandBitCount) - 1;
    // Bits in raw exponent
    static const int exponentBitCount = 11;
    static const int exponentMask = (1 << exponentBitCount) - 1;
    // Note: IEEE 754 conventionally uses 1023 as the exponent
    // bias.  That's because they treat the significand as a
    // fixed-point number with one bit (the hidden bit) integer
    // portion.  The logic here reconstructs the significand as a
    // pure fraction, so we need to accomodate that when
    // reconstructing the binary exponent.
    static const int64_t exponentBias = (1 << (exponentBitCount - 1)) - 2; // 1022

    // Step 0: Deconstruct the target number
    // Note: this strongly assumes IEEE 754 binary64 format
    uint64_t raw = bitPatternForDouble(d);
    int exponentBitPattern = (raw >> significandBitCount) & exponentMask;
    uint64_t significandBitPattern = raw & significandMask;

    // Step 1: Handle the various input cases:
    int binaryExponent;
    uint64_t significand;
    if (digits_length < 17) {
        return 0;
    } else if (exponentBitPattern == exponentMask) { // NaN or Infinity
        // Return no digits
        return 0;
    } else if (exponentBitPattern == 0) {
        if (significandBitPattern == 0) { // Zero
            digits[0] = 0;
            *decimalExponent = 0;
            return 1;
        } else { // subnormal
            binaryExponent = 1 - exponentBias;
            significand = significandBitPattern
                          << (64 - significandBitCount - 1);
        }
    } else { // normal
        binaryExponent = exponentBitPattern - exponentBias;
        uint64_t hiddenBit = (uint64_t)1 << significandBitCount;
        uint64_t fullSignificand = significandBitPattern + hiddenBit;
        significand = fullSignificand << (64 - significandBitCount - 1);
    }

    // Step 2: Determine the exact unscaled target interval

    // Grisu-style algorithms construct the shortest decimal digit
    // sequence within a specific interval.  To build the appropriate
    // interval, we start by computing the exact midpoints between
    // this floating-point value and the adjacent ones.

    uint64_t halfUlp = (uint64_t)1 << (64 - significandBitCount - 2);
    uint64_t quarterUlp = halfUlp >> 1;
    uint64_t upperMidpointExact = significand + halfUlp;

    int isBoundary = significandBitPattern == 0;
    uint64_t lowerMidpointExact
        = significand - (isBoundary ? quarterUlp : halfUlp);

    // Step 3: Estimate the base 10 exponent

    // Grisu algorithms are based in part on a simple technique for
    // generating a base-10 form for a binary floating-point number.
    // Start with a binary floating-point number `f * 2^e` and then
    // estimate the decimal exponent `p`. You can then rewrite your
    // original number as:
    //
    // ```
    //     f * 2^e * 10^-p * 10^p
    // ```
    //
    // The last term is part of our output, and a good estimate for
    // `p` will ensure that `2^e * 10^-p` is going to be close to 1.
    // So multiplying the first three terms yields a fraction suitable
    // for producing the decimal digits.  So we need to estimate `p`:

    int base10Exponent = decimalExponentFor2ToThe(binaryExponent);

    // Step 4: Compute a power-of-10 scale factor
    // Compute `10^-p` to 128-bit precision.  We generate
    // both over- and under-estimates to ensure we can exactly
    // bound the later use of these values.
    swift_uint128_t powerOfTenRoundedDown;
    swift_uint128_t powerOfTenRoundedUp;
    int powerOfTenExponent = 0;
    intervalContainingPowerOf10_Double(-base10Exponent,
                                       &powerOfTenRoundedDown,
                                       &powerOfTenRoundedUp,
                                       &powerOfTenExponent);
    const int extraBits = binaryExponent + powerOfTenExponent;

    // Step 5: Scale the interval (with rounding)

    // As mentioned above, the final digit generation works
    // with an interval, so we actually apply the scaling
    // to the upper and lower midpoint values separately.

    // As part of the scaling here, we'll switch from a pure
    // fraction to a fixed-point form.  Using 14 bit integer portion
    // will allow us to compute four decimal digits at a time.
    static const int integerBits = 14;
    static const int fractionBits = 128 - integerBits;

    // We scale the interval in one of two different ways,
    // depending on whether the significand is even or odd...

    swift_uint128_t u, l;
    if (significandBitPattern & 1) {
        // Case A: Narrow the interval (odd significand)

        // Loitsch' original Grisu2 always narrows the interval.
        // Since our digit generation will select a value within
        // the scaled interval, narrowing the interval guarantees
        // that we will find a digit sequence that converts back
        // to the original value.

        // This ensures accuracy but, as explained in Loitsch' paper,
        // this carries a risk that there will be a shorter digit
        // sequence outside of our narrowed interval that we will
        // miss.  This risk obviously gets lower with increased
        // precision, but it wasn't until the Errol paper that anyone
        // had a good way to test whether a particular implementation
        // had sufficient precision.  That paper shows a way to enumerate
        // the worst-case numbers; those numbers that are extremely close
        // to the mid-points between adjacent floating-point values.
        // These are the values that might sit just outside of the
        // narrowed interval.  By testing these values, we can verify
        // the correctness of our implementation.

        // Multiply out the upper midpoint, rounding down...
        swift_uint128_t u1 = multiply128x64RoundingDown(powerOfTenRoundedDown,
                                                    upperMidpointExact);
        // Account for residual binary exponent and adjust
        // to the fixed-point format
        u = shiftRightRoundingDown128(u1, integerBits - extraBits);

        // Conversely for the lower midpoint...
        swift_uint128_t l1 = multiply128x64RoundingUp(powerOfTenRoundedUp,
                                                  lowerMidpointExact);
        l = shiftRightRoundingUp128(l1, integerBits - extraBits);

    } else {
        // Case B: Widen the interval (even significand)

        // As explained in Errol Theorem 6, in certain cases there is
        // a short decimal representation at the exact boundary of the
        // scaled interval.  When such a number is converted back to
        // binary, it will get rounded to the adjacent even
        // significand.

        // So when the significand is even, we widen the interval in
        // order to ensure that the exact midpoints are considered.
        // Of couse, this ensures that we find a short result but
        // carries a risk of selecting a result outside of the exact
        // scaled interval (which would be inaccurate).

        // The same testing approach described above also applies
        // to this case.

        swift_uint128_t u1 = multiply128x64RoundingUp(powerOfTenRoundedUp,
                                                  upperMidpointExact);
        u = shiftRightRoundingUp128(u1, integerBits - extraBits);

        swift_uint128_t l1 = multiply128x64RoundingDown(powerOfTenRoundedDown,
                                                    lowerMidpointExact);
        l = shiftRightRoundingDown128(l1, integerBits - extraBits);
    }

    // Step 6: Align first digit, adjust exponent

    // This preps for digit generation.  It just multiplies repeatedly
    // by 10 until we have exactly one decimal digit in the integer
    // part, adjusting the exponent as we go.

    // In particular, this prunes leading zeros from subnormals.
    // Generate digits for `t` with interval width `delta`
    swift_uint128_t t = u;
    swift_uint128_t delta = u;
    subtract128x128(&delta, l); // Explained below.
    int exponent = base10Exponent + 1;

    // Except for subnormals, this loop should never run more than once.
#if HAVE_UINT128_T
    static const swift_uint128_t fixedPointOne = (__uint128_t)1 << fractionBits;
    while (t < fixedPointOne)
#else
    // Because 1.0 in fixed point has a lot of zeros, it suffices
    // to only compare the high-order word here.  This is a minor
    // performance win.
    while (t.high < ((uint32_t)1 << (fractionBits % 32)))
#endif
    {
        exponent -= 1;
        multiply128xi32(&delta, 10);
        multiply128xi32(&t, 10);
    }

    // Step 7: Generate digits

    // This is a common part of Grisu-style algorithms.  The
    // underlying idea is to generate digits for the scaled upper and
    // lower boundaries, and stop when we hit the first different
    // digit (at which point, the digit for the upper midpoint is the
    // candidate final digit).  To understand this, just note that
    // 0.1234 is the shortest decimal between u = 0.123456 and l =
    // 0.123345.

    // Grisu uses a slightly optimized technique: it generates digits
    // for the upper bound (multiplying by 10 to isolate each digit)
    // and multiplies the interval width `delta` at the same time.
    // The `different digit` criteria above translates to a test for
    // `delta` being larger than the remainder.

    int8_t *digit_p = digits;

    // Adjustment above already set up the first digit:
    int nextDigit = extractIntegerPart128(&t, fractionBits);
    clearIntegerPart128(&t, fractionBits);

    // Further optimization: Generating four digits at a time reduces
    // the total arithmetic required per digit.  Note: The following
    // block can be entirely removed with no effect on the result.
    // If you're trying to understand this algorithm, just skip this
    // block on first reading.
    swift_uint128_t d0 = delta;
    multiply128xi32(&d0, 10000);
    swift_uint128_t t0 = t;
    multiply128xi32(&t0, 10000);
    int fourDigits = extractIntegerPart128(&t0, fractionBits); // 4 digits
    clearIntegerPart128(&t0, fractionBits);
    while (isLessThan128x128(d0, t0)) {
        *digit_p++ = nextDigit;
        int d = fourDigits / 100; // top 2 digits
        *digit_p++ = d / 10;
        *digit_p++ = d % 10;
        d = fourDigits % 100; // bottom 2 digits
        *digit_p++ = d / 10;
        nextDigit = d % 10;
        t = t0;
        delta = d0;
        multiply128xi32(&d0, 10000);
        multiply128xi32(&t0, 10000);
        fourDigits = extractIntegerPart128(&t0, fractionBits);
        clearIntegerPart128(&t0, fractionBits);
    }

    // Finish by generating one digit at a time.
    while (isLessThan128x128(delta, t)) {
        *digit_p++ = nextDigit;
        multiply128xi32(&delta, 10);
        multiply128xi32(&t, 10);
        nextDigit = extractIntegerPart128(&t, fractionBits);
        clearIntegerPart128(&t, fractionBits);
    }

    // Adjust the final digit to be closer to the original value.  It accounts
    // for the fact that sometimes there is more than one shortest digit
    // sequence.

    // For example, consider how the above would work if you had the
    // value 0.1234 and computed u = 0.1257, l = 0.1211.  The above
    // digit generation works with `u`, so produces 0.125.  But the
    // values 0.122, 0.123, and 0.124 are just as short and 0.123 is
    // the best choice, since it's closest to the original value.

    // If `delta <= t + 1.0`, then the interval is narrower than
    // one decimal digit, so there is no other option.

    // Note: We've already consumed most of our available precision,
    // so it's okay to just work in 64 bits here...
    uint64_t deltaHigh64 = extractHigh64From128(delta);
    uint64_t tHigh64 = extractHigh64From128(t);
    if (deltaHigh64 > tHigh64 + ((uint64_t)1 << (fractionBits % 64))) {
        // Note: 64-bit arithmetic is okay here
        uint64_t skew;
        if (isBoundary) {
            // If we're at the boundary where the exponent shifts,
            // then the original value is 1/3 of the way from
            // the bottom of the interval ...
            skew = deltaHigh64 - deltaHigh64 / 3 - tHigh64;
        } else {
            // ... otherwise it's exactly in the middle.
            skew = deltaHigh64 / 2 - tHigh64;
        }

        // The `skew` above is the difference between our
        // computed digits and the original exact value.
        // Use that to offset the final digit:
        uint64_t one = (uint64_t)(1) << (64 - integerBits);
        uint64_t fractionMask = one - 1;
        uint64_t oneHalf = one >> 1;
        if ((skew & fractionMask) == oneHalf) {
            int adjust = (int)(skew >> (64 - integerBits));
            // If the skew is exactly integer + 1/2, round the
            // last digit even after adjustment
            nextDigit = (nextDigit - adjust) & ~1;
        } else {
            // Else round to nearest...
            int adjust = (int)((skew + oneHalf) >> (64 - integerBits));
            nextDigit = (nextDigit - adjust);
        }
    }
    *digit_p++ = nextDigit; // Store the final digit.

    *decimalExponent = exponent;
    return digit_p - digits;
}
#endif

#if SWIFT_DTOA_FLOAT_SUPPORT
// Return raw bits encoding the float
static uint64_t bitPatternForFloat(float f) {
    union { float f; uint32_t u; } converter;
    converter.f = f;
    return converter.u;
}

// Decompose an IEEE 754 binary32 single-precision float
// into decimal digits and a corresponding decimal exponent.

// See swift_decompose_double for detailed comments on the algorithm here
int swift_decompose_float(float f,
    int8_t *digits, size_t digits_length, int *decimalExponent)
{
    static const int significandBitCount = FLT_MANT_DIG - 1;
    static const uint32_t significandMask
        = ((uint32_t)1 << significandBitCount) - 1;
    static const int exponentBitCount = 8;
    static const int exponentMask = (1 << exponentBitCount) - 1;
    // See comments in swift_decompose_double
    static const int64_t exponentBias = (1 << (exponentBitCount - 1)) - 2; // 125

    // Step 0: Deconstruct the target number
    // Note: this strongly assumes IEEE 754 binary32 format
    uint32_t raw = bitPatternForFloat(f);
    int exponentBitPattern = (raw >> significandBitCount) & exponentMask;
    uint32_t significandBitPattern = raw & significandMask;

    // Step 1: Handle the various input cases:
    int binaryExponent;
    uint32_t significand;
    if (digits_length < 9) {
        // Ensure we have space for 9 digits
        return 0;
    } else if (exponentBitPattern == exponentMask) { // NaN or Infinity
        // Return no digits
        return 0;
    } else if (exponentBitPattern == 0) {
        if (significandBitPattern == 0) { // Zero
            // Return one zero digit and decimalExponent = 0.
            digits[0] = 0;
            *decimalExponent = 0;
            return 1;
        } else { // Subnormal
            binaryExponent = 1 - exponentBias;
            significand = significandBitPattern << (32 - significandBitCount - 1);
        }
    } else { // normal
        binaryExponent = exponentBitPattern - exponentBias;
        uint32_t hiddenBit = (uint32_t)1 << (uint32_t)significandBitCount;
        uint32_t fullSignificand = significandBitPattern + hiddenBit;
        significand = fullSignificand << (32 - significandBitCount - 1);
    }

    // Step 2: Determine the exact unscaled target interval
    static const uint32_t halfUlp = (uint32_t)1 << (32 - significandBitCount - 2);
    uint32_t upperMidpointExact = significand + halfUlp;

    int isBoundary = significandBitPattern == 0;
    static const uint32_t quarterUlp = halfUlp >> 1;
    uint32_t lowerMidpointExact
        = significand - (isBoundary ? quarterUlp : halfUlp);

    // Step 3: Estimate the base 10 exponent
    int base10Exponent = decimalExponentFor2ToThe(binaryExponent);

    // Step 4: Compute a power-of-10 scale factor
    uint64_t powerOfTenRoundedDown = 0;
    uint64_t powerOfTenRoundedUp = 0;
    int powerOfTenExponent = 0;
    intervalContainingPowerOf10_Float(-base10Exponent,
                                      &powerOfTenRoundedDown,
                                      &powerOfTenRoundedUp,
                                      &powerOfTenExponent);
    const int extraBits = binaryExponent + powerOfTenExponent;

    // Step 5: Scale the interval (with rounding)
    static const int integerBits = 5;
    const int shift = integerBits - extraBits;
    const int roundUpBias = (1 << shift) - 1;
    static const int fractionBits = 64 - integerBits;
    uint64_t u, l;
    if (significandBitPattern & 1) {
        // Narrow the interval (odd significand)
        uint64_t u1 = multiply64x32RoundingDown(powerOfTenRoundedDown,
                                                upperMidpointExact);
        u = u1 >> shift; // Rounding down

        uint64_t l1 = multiply64x32RoundingUp(powerOfTenRoundedUp,
                                              lowerMidpointExact);
        l = (l1 + roundUpBias) >> shift; // Rounding Up
    } else {
        // Widen the interval (even significand)
        uint64_t u1 = multiply64x32RoundingUp(powerOfTenRoundedUp,
                                              upperMidpointExact);
        u = (u1 + roundUpBias) >> shift; // Rounding Up

        uint64_t l1 = multiply64x32RoundingDown(powerOfTenRoundedDown,
                                                lowerMidpointExact);
        l = l1 >> shift; // Rounding down
    }

    // Step 6: Align first digit, adjust exponent
    // In particular, this prunes leading zeros from subnormals
    static const uint64_t fixedPointOne = (uint64_t)1 << fractionBits;
    static const uint64_t fixedPointMask = fixedPointOne - 1;
    uint64_t t = u;
    uint64_t delta = u - l;
    int exponent = base10Exponent + 1;

    while (t < fixedPointOne) {
        exponent -= 1;
        delta *= 10;
        t *= 10;
    }

    // Step 7: Generate digits
    int8_t *digit_p = digits;
    int nextDigit = (int)(t >> fractionBits);
    t &= fixedPointMask;

    // Generate one digit at a time...
    while (t > delta) {
        *digit_p++ = nextDigit;
        delta *= 10;
        t *= 10;
        nextDigit = (int)(t >> fractionBits);
        t &= fixedPointMask;
    }

    // Adjust the final digit to be closer to the original value
    if (delta > t + fixedPointOne) {
        uint64_t skew;
        if (isBoundary) {
            skew = delta - delta / 3 - t;
        } else {
            skew = delta / 2 - t;
        }
        uint64_t one = (uint64_t)(1) << (64 - integerBits);
        uint64_t lastAccurateBit = 1ULL << 24;
        uint64_t fractionMask = (one - 1) & ~(lastAccurateBit - 1);
        uint64_t oneHalf = one >> 1;
        if (((skew + (lastAccurateBit >> 1)) & fractionMask) == oneHalf) {
            // If the skew is exactly integer + 1/2, round the last
            // digit even after adjustment
            int adjust = (int)(skew >> (64 - integerBits));
            nextDigit = (nextDigit - adjust) & ~1;
        } else {
            // Else round to nearest...
            int adjust = (int)((skew + oneHalf) >> (64 - integerBits));
            nextDigit = (nextDigit - adjust);
        }
    }
    *digit_p++ = nextDigit;

    *decimalExponent = exponent;
    return digit_p - digits;
}
#endif

#if SWIFT_DTOA_FLOAT80_SUPPORT
// See `swift_decompose_double` for detailed comments on this implementatoin.
int swift_decompose_float80(long double d,
    int8_t *digits, size_t digits_length, int *decimalExponent)
{
    static const int exponentBitCount = 15;
    static const int exponentMask = (1 << exponentBitCount) - 1;
    // See comments in swift_decompose_double to understand
    // why we use 16,382 instead of 16,383 here.
    static const int64_t exponentBias = (1 << (exponentBitCount - 1)) - 2; // 16,382

    // Step 0: Deconstruct the target number
    // Note: this strongly assumes Intel 80-bit extended format in LSB
    // byte order
    const uint64_t *raw_p = (const uint64_t *)&d;
    int exponentBitPattern = raw_p[1] & exponentMask;
    uint64_t significandBitPattern = raw_p[0];

    // Step 1: Handle the various input cases:
    int64_t binaryExponent;
    uint64_t significand;
    if (digits_length < 21) {
        return 0;
    } else if (exponentBitPattern == exponentMask) { // NaN or Infinity
        // Return no digits
        return 0;
    } else if (exponentBitPattern == 0) {
        if (significandBitPattern == 0) { // Zero
            digits[0] = 0;
            *decimalExponent = 0;
            return 1;
        } else { // subnormal
            binaryExponent = 1 - exponentBias;
            significand = significandBitPattern;
        }
    } else if (significandBitPattern >> 63) { // Normal
        binaryExponent = exponentBitPattern - exponentBias;
        significand = significandBitPattern;
    } else {
        // "Unnormal" values are rejected as invalid by 80387 and later.
        // Treat them the same as NaNs here.
        return 0;
    }

    // Step 2: Determine the exact unscaled target interval
    uint64_t halfUlp = (uint64_t)1 << 63;
    uint64_t quarterUlp = halfUlp >> 1;
    uint64_t threeQuarterUlp = halfUlp + quarterUlp;
    swift_uint128_t upperMidpointExact, lowerMidpointExact;
    initialize128WithHighLow64(upperMidpointExact, significand, halfUlp);
    int isBoundary = (significandBitPattern & 0x7fffffffffffffff) == 0;
    // Subtract 1/4 or 1/2 ULP by first subtracting 1 full ULP, then adding some back
    initialize128WithHighLow64(lowerMidpointExact, significand - 1, isBoundary ? threeQuarterUlp : halfUlp);

    // Step 3: Estimate the base 10 exponent
    int base10Exponent = decimalExponentFor2ToThe(binaryExponent);

    // Step 4: Compute a power-of-10 scale factor
    swift_uint192_t powerOfTenRoundedDown;
    swift_uint192_t powerOfTenRoundedUp;
    int powerOfTenExponent = 0;
    intervalContainingPowerOf10_Float80(-base10Exponent,
                                        &powerOfTenRoundedDown,
                                        &powerOfTenRoundedUp,
                                        &powerOfTenExponent);
    const int extraBits = binaryExponent + powerOfTenExponent;

    // Step 5: Scale the interval (with rounding)
    static const int integerBits = 14;
    static const int fractionBits = 192 - integerBits;
#if HAVE_UINT128_T
    static const int highFractionBits = fractionBits % 64;
#else
    static const int highFractionBits = fractionBits % 32;
#endif
    swift_uint192_t u, l;
    if (significandBitPattern & 1) {
        // Narrow the interval (odd significand)
        u = powerOfTenRoundedDown;
        multiply192x128RoundingDown(&u, upperMidpointExact);
        shiftRightRoundingDown192(&u, integerBits - extraBits);

        l = powerOfTenRoundedUp;
        multiply192x128RoundingUp(&l, lowerMidpointExact);
        shiftRightRoundingUp192(&l, integerBits - extraBits);
    } else {
        // Widen the interval (even significand)
        u = powerOfTenRoundedUp;
        multiply192x128RoundingUp(&u, upperMidpointExact);
        shiftRightRoundingUp192(&u, integerBits - extraBits);

        l = powerOfTenRoundedDown;
        multiply192x128RoundingDown(&l, lowerMidpointExact);
        shiftRightRoundingDown192(&l, integerBits - extraBits);
    }

    // Step 6: Align first digit, adjust exponent
    // In particular, this prunes leading zeros from subnormals
    static const uint64_t fixedPointOneHigh = (uint64_t)1 << highFractionBits;
    static const uint64_t fixedPointMaskHigh = fixedPointOneHigh - 1;
    swift_uint192_t t = u;
    swift_uint192_t delta = u;
    subtract192x192(&delta, l);
    int exponent = base10Exponent + 1;

    while (t.high < fixedPointOneHigh) {
        exponent -= 1;
        multiply192xi32(&delta, 10);
        multiply192xi32(&t, 10);
    }

    // Step 7: Generate digits
    int8_t *digit_p = digits;

    // Adjustment above already set up the first digit:
    int nextDigit = (int)(t.high >> highFractionBits);
    t.high &= fixedPointMaskHigh;

    // Generate four digits at a time ...
    swift_uint192_t d0 = delta;
    swift_uint192_t t0 = t;
    multiply192xi32(&d0, 10000);
    multiply192xi32(&t0, 10000);
    int fourDigits = (int)(t0.high >> highFractionBits);
    t0.high &= fixedPointMaskHigh;
    while (isLessThan192x192(d0, t0)) {
        *digit_p++ = nextDigit;
        int d = fourDigits / 100;
        *digit_p++ = d / 10;
        *digit_p++ = d % 10;
        d = fourDigits % 100;
        *digit_p++ = d / 10;
        nextDigit = d % 10;
        t = t0;
        delta = d0;
        multiply192xi32(&d0, 10000);
        multiply192xi32(&t0, 10000);
        fourDigits = (int)(t0.high >> highFractionBits);
        t0.high &= fixedPointMaskHigh;
    }

    // Generate one digit at a time...
    while (isLessThan192x192(delta, t)) {
        *digit_p++ = nextDigit;
        multiply192xi32(&delta, 10);
        multiply192xi32(&t, 10);
        nextDigit = (int)(t.high >> highFractionBits);
        t.high &= fixedPointMaskHigh;
    }

    // Adjust the final digit to be closer to the original value
    // We've already consumed most of our available precision, so it's
    // okay to just work in 64 bits here...
#if HAVE_UINT128_T
    uint64_t deltaHigh64 = delta.high;
    uint64_t tHigh64 = t.high;
#else
    uint64_t deltaHigh64 = ((uint64_t)delta.high << 32) + delta.e;
    uint64_t tHigh64 = ((uint64_t)t.high << 32) + t.e;
#endif
    if (deltaHigh64 > tHigh64 + ((uint64_t)1 << (fractionBits % 64))) {
        uint64_t skew;
        if (isBoundary) {
            skew = deltaHigh64 - deltaHigh64 / 3 - tHigh64;
        } else {
            skew = deltaHigh64 / 2 - tHigh64;
        }
        uint64_t one = (uint64_t)(1) << (64 - integerBits);
        uint64_t fractionMask = one - 1;
        uint64_t oneHalf = one >> 1;
        if ((skew & fractionMask) == oneHalf) {
            int adjust = (int)(skew >> (64 - integerBits));
            // If the skew is integer + 1/2, round the last digit even
            // after adjustment
            nextDigit = (nextDigit - adjust) & ~1;
        } else {
            // Else round to nearest...
            int adjust = (int)((skew + oneHalf) >> (64 - integerBits));
            nextDigit = (nextDigit - adjust);
        }
    }
    *digit_p++ = nextDigit;

    *decimalExponent = exponent;
    return digit_p - digits;
}
#endif

//
// ---------------- High-level API -----------------
//
// Functions that format a Float/Double/Float80
// directly into a buffer.
//
// These handle various exception cases (infinity, Nan, zero)
// before invoking the general base-10 conversion.

#if SWIFT_DTOA_FLOAT_SUPPORT || SWIFT_DTOA_DOUBLE_SUPPORT || SWIFT_DTOA_FLOAT80_SUPPORT
static size_t swift_format_constant(char *dest, size_t length, const char *s) {
    const size_t l = strlen(s);
    if (length <= l) {
        return 0;
    }
    strcpy(dest, s);
    return l;
}
#endif

#if SWIFT_DTOA_FLOAT_SUPPORT
size_t swift_format_float(float d, char *dest, size_t length)
{
    if (!isfinite(d)) {
        if (isinf(d)) {
            // Infinity
            if (signbit(d)) {
                return swift_format_constant(dest, length, "-inf");
            } else {
                return swift_format_constant(dest, length, "inf");
            }
        } else {
            // NaN
            static const int significandBitCount = 23;
            uint32_t raw = bitPatternForFloat(d);
            const char *sign = signbit(d) ? "-" : "";
            const char *signaling = ((raw >> (significandBitCount - 1)) & 1) ? "" : "s";
            uint32_t payload = raw & ((1L << (significandBitCount - 2)) - 1);
            char buff[32];
            if (payload != 0) {
                snprintf(buff, sizeof(buff), "%s%snan(0x%x)",
                         sign, signaling, payload);
            } else {
                snprintf(buff, sizeof(buff), "%s%snan",
                         sign, signaling);
            }
            return swift_format_constant(dest, length, buff);
        }
    }

    // zero
    if (d == 0.0) {
        if (signbit(d)) {
            return swift_format_constant(dest, length, "-0.0");
        } else {
            return swift_format_constant(dest, length, "0.0");
        }
    }

    // Decimal numeric formatting
    int decimalExponent;
    int8_t digits[9];
    int digitCount =
        swift_decompose_float(d, digits, sizeof(digits), &decimalExponent);
    // People use float to model integers <= 2^24, so we use that
    // as a cutoff for decimal vs. exponential format.
    if (decimalExponent < -3 || fabsf(d) > 0x1.0p24F) {
        return swift_format_exponential(dest, length, signbit(d),
                 digits, digitCount, decimalExponent);
    } else {
        return swift_format_decimal(dest, length, signbit(d),
                 digits, digitCount, decimalExponent);
    }
}
#endif

#if SWIFT_DTOA_DOUBLE_SUPPORT
size_t swift_format_double(double d, char *dest, size_t length)
{
    if (!isfinite(d)) {
        if (isinf(d)) {
            // Infinity
            if (signbit(d)) {
                return swift_format_constant(dest, length, "-inf");
            } else {
                return swift_format_constant(dest, length, "inf");
            }
        } else {
            // NaN
            static const int significandBitCount = 52;
            uint64_t raw = bitPatternForDouble(d);
            const char *sign = signbit(d) ? "-" : "";
            const char *signaling = ((raw >> (significandBitCount - 1)) & 1) ? "" : "s";
            uint64_t payload = raw & ((1ull << (significandBitCount - 2)) - 1);
            char buff[32];
            if (payload != 0) {
                snprintf(buff, sizeof(buff), "%s%snan(0x%" PRIx64 ")",
                         sign, signaling, payload);
            } else {
                snprintf(buff, sizeof(buff), "%s%snan",
                         sign, signaling);
            }
            return swift_format_constant(dest, length, buff);
        }
    }

    // zero
    if (d == 0.0) {
        if (signbit(d)) {
            return swift_format_constant(dest, length, "-0.0");
        } else {
            return swift_format_constant(dest, length, "0.0");
        }
    }

    // Decimal numeric formatting
    int decimalExponent;
    int8_t digits[17];
    int digitCount =
        swift_decompose_double(d, digits, sizeof(digits), &decimalExponent);
    // People use double to model integers <= 2^53, so we use that
    // as a cutoff for decimal vs. exponential format.
    if (decimalExponent < -3 || fabs(d) > 0x1.0p53) {
        return swift_format_exponential(dest, length, signbit(d),
                 digits, digitCount, decimalExponent);
    } else {
        return swift_format_decimal(dest, length, signbit(d),
                 digits, digitCount, decimalExponent);
    }
}
#endif

#if SWIFT_DTOA_FLOAT80_SUPPORT
size_t swift_format_float80(long double d, char *dest, size_t length)
{
    if (!isfinite(d)) {
        if (isinf(d)) {
            // Infinity
            if (signbit(d)) {
                return swift_format_constant(dest, length, "-inf");
            } else {
                return swift_format_constant(dest, length, "inf");
            }
        } else {
            // NaN
            // Assumes Intel 80-bit extended format in LSB byte order:
            uint64_t significandBitPattern = *(const uint64_t *)&d;
            const char *sign = signbit(d) ? "-" : "";
            const char *signaling = ((significandBitPattern >> 62) & 1) ? "" : "s";
            uint64_t payload = significandBitPattern & (((uint64_t)1 << 61) - 1);
            char buff[32];
            if (payload != 0) {
                snprintf(buff, sizeof(buff), "%s%snan(0x%" PRIx64 ")",
                         sign, signaling, payload);
            } else {
                snprintf(buff, sizeof(buff), "%s%snan",
                         sign, signaling);
            }
            return swift_format_constant(dest, length, buff);
        }
    }

    // zero
    if (d == 0.0) {
        if (signbit(d)) {
            return swift_format_constant(dest, length, "-0.0");
        } else {
            return swift_format_constant(dest, length, "0.0");
        }
    }

    // Decimal numeric formatting
    int decimalExponent;
    int8_t digits[21];
    int digitCount =
        swift_decompose_float80(d, digits, sizeof(digits), &decimalExponent);
    // People use long double to model integers <= 2^64, so we use that
    // as a cutoff for decimal vs. exponential format.
    // The constant is written out as a float80 (aka "long double") literal
    // here since it can't be expressed as a 64-bit integer.
    if (decimalExponent < -3 || fabsl(d) > 0x1.0p64L) {
        return swift_format_exponential(dest, length, signbit(d),
                 digits, digitCount, decimalExponent);
    } else {
        return swift_format_decimal(dest, length, signbit(d),
                 digits, digitCount, decimalExponent);
    }
}
#endif

/**
 * Routines to format a decomposed value into a standard string form.
 */

// Format into exponential format: "1.234e+56"
// Returns number of characters actually written to `dest`.
// Returns zero if buffer is too small.
size_t swift_format_exponential(char *dest, size_t length,
    bool negative, const int8_t *digits, int digit_count, int exponent)
{
    // Largest buffer we could possibly need:
    size_t maximum_size = digit_count + 9;
    if (length < maximum_size) {
        // We only do the detailed check if the size is borderline.
        size_t actual_size =
            + (negative ? 1 : 0) // Leading minus
            + digit_count // digits
            + (digit_count > 1 ? 1 : 0) // decimal
            + 1 // 'e'
            + 1 // sign
            + (exponent > 99 ? (exponent > 999 ? 4 : 3) : 2) // exponent
            + 1; // trailing zero byte
        if (length < actual_size) {
            if (length > 0) {
                dest[0] = 0;
            }
            return 0;
        }
    }
    char *p = dest;
    if (negative) {
        *p++ = '-';
    }

    *p++ = digits[0] + '0';
    exponent -= 1;
    if (digit_count > 1) {
        *p++ = '.';
        for (int i = 1; i < digit_count; i++) {
            *p++ = digits[i] + '0';
        }
    }
    *p++ = 'e';
    if (exponent < 0) {
        *p++ = '-';
        exponent = -exponent;
    } else {
        *p++ = '+';
    }
    if (exponent > 99) {
        if (exponent > 999) {
            *p++ = (exponent / 1000 % 10) + '0';
        }
        *p++ = (exponent / 100 % 10) + '0';
        exponent %= 100;
    }
    *p++ = (exponent / 10) + '0';
    *p++ = (exponent % 10) + '0';
    *p = '\0';
    return p - dest;
}

// Format into decimal form: "123456789000.0", "1234.5678", "0.0000001234"
// Returns number of bytes of `dest` actually used, or zero if
// provided buffer is too small.
size_t swift_format_decimal(char *dest, size_t length,
    bool negative, const int8_t *digits, int digit_count, int exponent)
{
    // Largest buffer we could possibly need:
    size_t maximum_size =
        digit_count // All the digits
        + (exponent > 0 ? exponent : -exponent) // Max # of extra zeros
        + 4; // Max # of other items
    if (length < maximum_size) {
        // We only do the detailed check if the size is borderline.
        if (exponent <= 0) { // "0.0000001234"
            size_t actual_size =
                (negative ? 1 : 0) // Leading minus
                + 2 // Leading "0."
                + (-exponent) // Leading zeros after decimal point
                + digit_count
                + 1; // Trailing zero byte
            if (length < actual_size) {
                if (length > 0) {
                    dest[0] = 0;
                }
                return 0;
            }
        } else if (exponent < digit_count) { // "123.45"
            size_t actual_size =
                (negative ? 1 : 0) // Leading minus
                + digit_count
                + 1 // embedded decimal point
                + 1; // Trailing zero byte
            if (length < actual_size) {
                if (length > 0) {
                    dest[0] = 0;
                }
                return 0;
            }
        } else { // "12345000.0"
            size_t actual_size =
                (negative ? 1 : 0) // Leading minus
                + digit_count
                + (exponent - digit_count) // trailing zeros
                + 2 // ".0" to mark this as floating point
                + 1; // Trailing zero byte
            if (length < actual_size) {
                if (length > 0) {
                    dest[0] = 0;
                }
                return 0;
            }
        }
    }

    char *p = dest;
    if (negative) {
        *p++ = '-';
    }

    if (exponent <= 0) {
        *p++ = '0';
        *p++ = '.';
        while (exponent < 0) {
            *p++ = '0';
            exponent += 1;
        }
        for (int i = 0; i < digit_count; ++i) {
            *p++ = digits[i] + '0';
        }
    } else if (exponent < digit_count) {
        for (int i = 0; i < digit_count; i++) {
            if (exponent == 0) {
                *p++ = '.';
            }
            *p++ = digits[i] + '0';
            exponent -= 1;
        }
    } else {
        for (int i = 0; i < digit_count; i++) {
            *p++ = digits[i] + '0';
            exponent -= 1;
        }
        while (exponent > 0) {
            *p++ = '0';
            exponent -= 1;
        }
        *p++ = '.';
        *p++ = '0';
    }
    *p = '\0';
    return p - dest;
}

//
// ------------  Arithmetic helpers ----------------
//

// The core algorithm relies heavily on fraction and fixed-point
// arithmetic with 64-bit, 128-bit, and 192-bit integer values. (For
// float, double, and float80, respectively.) They also need precise
// control over all rounding.
//
// Note that most arithmetic operations are the same for integers and
// fractions, so we can just use the normal integer operations in most
// places.  Multiplication however, is different for fixed-size
// fractions.  Integer multiplication preserves the low-order part and
// discards the high-order part (ignoring overflow).  Fraction
// multiplication preserves the high-order part and discards the
// low-order part (rounding).  So most of the arithmetic helpers here
// are for multiplication.

// Note: With 64-bit GCC and Clang, we get a noticable performance
// gain by using `__uint128_t`.  Otherwise, we have to break things
// down into 32-bit chunks so we don't overflow 64-bit temporaries.

#if SWIFT_DTOA_FLOAT_SUPPORT
// Multiply a 64-bit fraction by a 32-bit fraction, rounding down.
static uint64_t multiply64x32RoundingDown(uint64_t lhs, uint32_t rhs) {
    static const uint64_t mask32 = UINT32_MAX;
    uint64_t t = ((lhs & mask32) * rhs) >> 32;
    return t + (lhs >> 32) * rhs;
}

// Multiply a 64-bit fraction by a 32-bit fraction, rounding up.
static uint64_t multiply64x32RoundingUp(uint64_t lhs, uint32_t rhs) {
    static const uint64_t mask32 = UINT32_MAX;
    uint64_t t = (((lhs & mask32) * rhs) + mask32) >> 32;
    return t + (lhs >> 32) * rhs;
}

// Multiply a 64-bit fraction by a 64-bit fraction, rounding down.
static uint64_t multiply64x64RoundingDown(uint64_t lhs, uint64_t rhs) {
#if HAVE_UINT128_T
    __uint128_t full = (__uint128_t)lhs * rhs;
    return (uint64_t)(full >> 64);
#else
    static const uint64_t mask32 = UINT32_MAX;
    uint64_t t = (lhs & mask32) * (rhs & mask32);
    t >>= 32;
    uint64_t a = (lhs >> 32) * (rhs & mask32);
    uint64_t b = (lhs & mask32) * (rhs >> 32);
    // Useful: If w,x,y,z are all 32-bit values, then:
    // w * x + y + z
    //   <= (2^64 - 2^33 + 1) + (2^32 - 1) + (2^32 - 1)
    //   <= 2^64 - 1
    //
    // That is, a product of two 32-bit values plus two more 32-bit
    // values can't overflow 64 bits.  (But "three more" can, so be
    // careful!)
    //
    // Here: t + a + (b & mask32) <= 2^64 - 1
    t += a + (b & mask32);
    t >>= 32;
    t += (b >> 32);
    return t + (lhs >> 32) * (rhs >> 32);
#endif
}

// Multiply a 64-bit fraction by a 64-bit fraction, rounding up.
static uint64_t multiply64x64RoundingUp(uint64_t lhs, uint64_t rhs) {
#if HAVE_UINT128_T
    static const __uint128_t roundingUpBias = ((__uint128_t)1 << 64) - 1;
    __uint128_t full = (__uint128_t)lhs * rhs;
    return (uint64_t)((full + roundingUpBias) >> 64);
#else
    static const uint64_t mask32 = UINT32_MAX;
    uint64_t t = (lhs & mask32) * (rhs & mask32);
    t = (t + mask32) >> 32;
    uint64_t a = (lhs >> 32) * (rhs & mask32);
    uint64_t b = (lhs & mask32) * (rhs >> 32);
    t += (a & mask32) + (b & mask32) + mask32;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    return t + (lhs >> 32) * (rhs >> 32);
#endif
}

#endif

#if SWIFT_DTOA_DOUBLE_SUPPORT
// Multiply a 128-bit fraction by a 64-bit fraction, rounding down.
static swift_uint128_t multiply128x64RoundingDown(swift_uint128_t lhs, uint64_t rhs) {
#if HAVE_UINT128_T
    uint64_t lhsl = (uint64_t)lhs;
    uint64_t lhsh = (uint64_t)(lhs >> 64);
    swift_uint128_t h = (swift_uint128_t)lhsh * rhs;
    swift_uint128_t l = (swift_uint128_t)lhsl * rhs;
    return h + (l >> 64);
#else
    swift_uint128_t result;
    static const uint64_t mask32 = UINT32_MAX;
    uint64_t rhs0 = rhs & mask32;
    uint64_t rhs1 = rhs >> 32;
    uint64_t t = (lhs.low) * rhs0;
    t >>= 32;
    uint64_t a = (lhs.b) * rhs0;
    uint64_t b = (lhs.low) * rhs1;
    t += a + (b & mask32);
    t >>= 32;
    t += (b >> 32);
    a = lhs.c * rhs0;
    b = lhs.b * rhs1;
    t += (a & mask32) + (b & mask32);
    result.low = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    a = lhs.high * rhs0;
    b = lhs.c * rhs1;
    t += (a & mask32) + (b & mask32);
    result.b = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    t += lhs.high * rhs1;
    result.c = t;
    result.high = t >> 32;
    return result;
#endif
}

// Multiply a 128-bit fraction by a 64-bit fraction, rounding up.
static swift_uint128_t multiply128x64RoundingUp(swift_uint128_t lhs, uint64_t rhs) {
#if HAVE_UINT128_T
    uint64_t lhsl = (uint64_t)lhs;
    uint64_t lhsh = (uint64_t)(lhs >> 64);
    swift_uint128_t h = (swift_uint128_t)lhsh * rhs;
    swift_uint128_t l = (swift_uint128_t)lhsl * rhs;
    const static __uint128_t bias = ((__uint128_t)1 << 64) - 1;
    return h + ((l + bias) >> 64);
#else
    swift_uint128_t result;
    static const uint64_t mask32 = UINT32_MAX;
    uint64_t rhs0 = rhs & mask32;
    uint64_t rhs1 = rhs >> 32;
    uint64_t t = (lhs.low) * rhs0 + mask32;
    t >>= 32;
    uint64_t a = (lhs.b) * rhs0;
    uint64_t b = (lhs.low) * rhs1;
    t += (a & mask32) + (b & mask32) + mask32;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    a = lhs.c * rhs0;
    b = lhs.b * rhs1;
    t += (a & mask32) + (b & mask32);
    result.low = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    a = lhs.high * rhs0;
    b = lhs.c * rhs1;
    t += (a & mask32) + (b & mask32);
    result.b = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    t += lhs.high * rhs1;
    result.c = t;
    result.high = t >> 32;
    return result;
#endif
}

#if !HAVE_UINT128_T
// Multiply a 128-bit fraction by a 32-bit integer in a 32-bit environment.
// (On 64-bit, we use a fast inline macro.)
static void multiply128xi32(swift_uint128_t *lhs, uint32_t rhs) {
    uint64_t t = (uint64_t)(lhs->low) * rhs;
    lhs->low = (uint32_t)t;
    t = (t >> 32) + (uint64_t)(lhs->b) * rhs;
    lhs->b = (uint32_t)t;
    t = (t >> 32) + (uint64_t)(lhs->c) * rhs;
    lhs->c = (uint32_t)t;
    t = (t >> 32) + (uint64_t)(lhs->high) * rhs;
    lhs->high = (uint32_t)t;
}

// Compare two 128-bit integers in a 32-bit environment
// (On 64-bit, we use a fast inline macro.)
static int isLessThan128x128(swift_uint128_t lhs, swift_uint128_t rhs) {
    return ((lhs.high < rhs.high)
            || ((lhs.high == rhs.high)
                && ((lhs.c < rhs.c)
                    || ((lhs.c == rhs.c)
                        && ((lhs.b < rhs.b)
                            || ((lhs.b == rhs.b)
                                && (lhs.low < rhs.low)))))));
}

// Subtract 128-bit values in a 32-bit environment
static void subtract128x128(swift_uint128_t *lhs, swift_uint128_t rhs) {
    uint64_t t = (uint64_t)lhs->low + (~rhs.low) + 1;
    lhs->low = (uint32_t)t;
    t = (t >> 32) + lhs->b + (~rhs.b);
    lhs->b = (uint32_t)t;
    t = (t >> 32) + lhs->c + (~rhs.c);
    lhs->c = (uint32_t)t;
    t = (t >> 32) + lhs->high + (~rhs.high);
    lhs->high = (uint32_t)t;
}
#endif

// Shift a 128-bit integer right, rounding down.
static swift_uint128_t shiftRightRoundingDown128(swift_uint128_t lhs, int shift) {
#if HAVE_UINT128_T
    return lhs >> shift;
#else
    // Note: Shift is always less than 32
    swift_uint128_t result;
    uint64_t t = (uint64_t)lhs.low >> shift;
    t += ((uint64_t)lhs.b << (32 - shift));
    result.low = t;
    t >>= 32;
    t += ((uint64_t)lhs.c << (32 - shift));
    result.b = t;
    t >>= 32;
    t += ((uint64_t)lhs.high << (32 - shift));
    result.c = t;
    t >>= 32;
    result.high = t;
    return result;
#endif
}

// Shift a 128-bit integer right, rounding up.
static swift_uint128_t shiftRightRoundingUp128(swift_uint128_t lhs, int shift) {
#if HAVE_UINT128_T
    uint64_t bias = ((uint64_t)1 << shift) - 1;
    return ((lhs + bias) >> shift);
#else
    swift_uint128_t result;
    const uint64_t bias = (1 << shift) - 1;
    uint64_t t = ((uint64_t)lhs.low + bias) >> shift;
    t += ((uint64_t)lhs.b << (32 - shift));
    result.low = t;
    t >>= 32;
    t += ((uint64_t)lhs.c << (32 - shift));
    result.b = t;
    t >>= 32;
    t += ((uint64_t)lhs.high << (32 - shift));
    result.c = t;
    t >>= 32;
    result.high = t;
    return result;
#endif
}
#endif

#if SWIFT_DTOA_FLOAT80_SUPPORT
// Multiply a 192-bit fraction by a 64-bit fraction, rounding down.
static void multiply192x64RoundingDown(swift_uint192_t *lhs, uint64_t rhs) {
#if HAVE_UINT128_T
    // Compute the three 128-bit cross-products
    __uint128_t cd = (__uint128_t)lhs->low * rhs;
    __uint128_t bc = (__uint128_t)lhs->mid * rhs;
    __uint128_t ab = (__uint128_t)lhs->high * rhs;
    // Add up the three 64-bit outputs (including carries)
    __uint128_t c = (cd >> 64) + (uint64_t)bc;
    __uint128_t b = (bc >> 64) + (uint64_t)ab + (c >> 64);
    __uint128_t a = (ab >> 64) + (b >> 64);
    lhs->high = a;
    lhs->mid = b;
    lhs->low = c;
#else
    static const uint64_t mask32 = UINT32_MAX;
    uint64_t rhs0 = rhs & mask32;
    uint64_t rhs1 = rhs >> 32;
    uint64_t t = lhs->low * rhs0;
    t >>= 32;
    uint64_t a = lhs->low * rhs1;
    uint64_t b = lhs->b * rhs0;
    t += a + (b & mask32);
    t >>= 32;
    t += (b >> 32);
    a = lhs->b * rhs1;
    b = lhs->c * rhs0;
    t += (a & mask32) + (b & mask32);
    lhs->low = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    a = lhs->c * rhs1;
    b = lhs->d * rhs0;
    t += (a & mask32) + (b & mask32);
    lhs->b = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    a = lhs->d * rhs1;
    b = lhs->e * rhs0;
    t += (a & mask32) + (b & mask32);
    lhs->c = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    a = lhs->e * rhs1;
    b = lhs->high * rhs0;
    t += (a & mask32) + (b & mask32);
    lhs->d = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    t += lhs->high * rhs1;
    lhs->e = t;
    lhs->high = t >> 32;
#endif
}

// Multiply a 192-bit fraction by a 64-bit fraction, rounding up.
static void multiply192x64RoundingUp(swift_uint192_t *lhs, uint64_t rhs) {
#if HAVE_UINT128_T
    // Compute the three 128-bit cross-products
    __uint128_t cd = (__uint128_t)lhs->low * rhs + UINT64_MAX;
    __uint128_t bc = (__uint128_t)lhs->mid * rhs;
    __uint128_t ab = (__uint128_t)lhs->high * rhs;
    // Add up the three 64-bit outputs (including carries)
    __uint128_t c = (cd >> 64) + (uint64_t)bc;
    __uint128_t b = (bc >> 64) + (uint64_t)ab + (c >> 64);
    __uint128_t a = (ab >> 64) + (b >> 64);
    lhs->high = a;
    lhs->mid = b;
    lhs->low = c;
#else
    static const uint64_t mask32 = UINT32_MAX;
    static const uint64_t bias = mask32;
    uint64_t rhs0 = rhs & mask32;
    uint64_t rhs1 = rhs >> 32;
    uint64_t t = lhs->low * rhs0 + bias;
    t >>= 32;
    uint64_t a = lhs->low * rhs1;
    uint64_t b = lhs->b * rhs0;
    t += (a & mask32) + (b & mask32) + bias;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    a = lhs->b * rhs1;
    b = lhs->c * rhs0;
    t += (a & mask32) + (b & mask32);
    lhs->low = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    a = lhs->c * rhs1;
    b = lhs->d * rhs0;
    t += (a & mask32) + (b & mask32);
    lhs->b = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    a = lhs->d * rhs1;
    b = lhs->e * rhs0;
    t += (a & mask32) + (b & mask32);
    lhs->c = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    a = lhs->e * rhs1;
    b = lhs->high * rhs0;
    t += (a & mask32) + (b & mask32);
    lhs->d = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    t += lhs->high * rhs1;
    lhs->e = t;
    lhs->high = t >> 32;
#endif
}

// Multiply a 192-bit fraction by a 64-bit integer.
// This is used in the digit generation to multiply by ten or
// 10,000. Note that rounding never appliles here.
// As used below, this will never overflow.
static void multiply192xi32(swift_uint192_t *lhs, uint32_t rhs) {
#if HAVE_UINT128_T
    __uint128_t t = (__uint128_t)lhs->low * rhs;
    lhs->low = (uint64_t)t;
    t = (t >> 64) + (__uint128_t)lhs->mid * rhs;
    lhs->mid = (uint64_t)t;
    t = (t >> 64) + (__uint128_t)lhs->high * rhs;
    lhs->high = (uint64_t)t;
#else
    uint64_t t = (uint64_t)lhs->low * rhs;
    lhs->low = t;
    t = (t >> 32) + (uint64_t)lhs->b * rhs;
    lhs->b = t;
    t = (t >> 32) + (uint64_t)lhs->c * rhs;
    lhs->c = t;
    t = (t >> 32) + (uint64_t)lhs->d * rhs;
    lhs->d = t;
    t = (t >> 32) + (uint64_t)lhs->e * rhs;
    lhs->e = t;
    t = (t >> 32) + (uint64_t)lhs->high * rhs;
    lhs->high = t;
#endif
}

// Multiply a 192-bit fraction by a 128-bit fraction, rounding down.
static void multiply192x128RoundingDown(swift_uint192_t *lhs, swift_uint128_t rhs) {
#if HAVE_UINT128_T
    // A full multiply of three 64-bit values by two 64-bit values
    // yields five such components.  We discard the bottom two (except
    // for carries) to get a rounded-down three-element result.
    __uint128_t current = (__uint128_t)lhs->low * (uint64_t)rhs;

    current = (current >> 64);
    __uint128_t t = (__uint128_t)lhs->low * (rhs >> 64);
    current += (uint64_t)t;
    __uint128_t next = t >> 64;
    t = (__uint128_t)lhs->mid * (uint64_t)rhs;
    current += (uint64_t)t;
    next += t >> 64;

    current = next + (current >> 64);
    t = (__uint128_t)lhs->mid * (rhs >> 64);
    current += (uint64_t)t;
    next = t >> 64;
    t = (__uint128_t)lhs->high * (uint64_t)rhs;
    current += (uint64_t)t;
    next += t >> 64;
    lhs->low = (uint64_t)current;

    current = next + (current >> 64);
    t = (__uint128_t)lhs->high * (rhs >> 64);
    current += t;
    lhs->mid = (uint64_t)current;
    lhs->high = (uint64_t)(current >> 64);
#else
    uint64_t a, b, c, d; // temporaries
    // Six 32-bit values multiplied by 4 32-bit values.  Oh my.
    static const uint64_t mask32 = UINT32_MAX;
    uint64_t t = lhs->low * rhs.low;
    t >>= 32;
    a = (uint64_t)lhs->low * rhs.b;
    b = (uint64_t)lhs->b * rhs.low;
    t += a + (b & mask32);
    t >>= 32;
    t += (b >> 32);
    a = (uint64_t)lhs->low * rhs.c;
    b = (uint64_t)lhs->b * rhs.b;
    c = (uint64_t)lhs->c * rhs.low;
    t += (a & mask32) + (b & mask32) + (c & mask32);
    t >>= 32;
    t += (a >> 32) + (b >> 32) + (c >> 32);
    a = (uint64_t)lhs->low * rhs.high;
    b = (uint64_t)lhs->b * rhs.c;
    c = (uint64_t)lhs->c * rhs.b;
    d = (uint64_t)lhs->d * rhs.low;
    t += (a & mask32) + (b & mask32) + (c & mask32) + (d & mask32);
    t >>= 32;
    t += (a >> 32) + (b >> 32) + (c >> 32) + (d >> 32);
    a = (uint64_t)lhs->b * rhs.high;
    b = (uint64_t)lhs->c * rhs.c;
    c = (uint64_t)lhs->d * rhs.b;
    d = (uint64_t)lhs->e * rhs.low;
    t += (a & mask32) + (b & mask32) + (c & mask32) + (d & mask32);
    lhs->low = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32) + (c >> 32) + (d >> 32);
    a = (uint64_t)lhs->c * rhs.high;
    b = (uint64_t)lhs->d * rhs.c;
    c = (uint64_t)lhs->e * rhs.b;
    d = (uint64_t)lhs->high * rhs.low;
    t += (a & mask32) + (b & mask32) + (c & mask32) + (d & mask32);
    lhs->b = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32) + (c >> 32) + (d >> 32);
    a = (uint64_t)lhs->d * rhs.high;
    b = (uint64_t)lhs->e * rhs.c;
    c = (uint64_t)lhs->high * rhs.b;
    t += (a & mask32) + (b & mask32) + (c & mask32);
    lhs->c = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32) + (c >> 32);
    a = (uint64_t)lhs->e * rhs.high;
    b = (uint64_t)lhs->high * rhs.c;
    t += (a & mask32) + (b & mask32);
    lhs->d = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    t += (uint64_t)lhs->high * rhs.high;
    lhs->e = t;
    lhs->high = t >> 32;
#endif
}

// Multiply a 192-bit fraction by a 128-bit fraction, rounding up.
static void multiply192x128RoundingUp(swift_uint192_t *lhs, swift_uint128_t rhs) {
#if HAVE_UINT128_T
    // Same as the rounding-down version, but we add
    // UINT128_MAX to the bottom two to force an extra
    // carry if they are non-zero.
    swift_uint128_t current = (swift_uint128_t)lhs->low * (uint64_t)rhs;
    current += UINT64_MAX;

    current = (current >> 64);
    swift_uint128_t t = (swift_uint128_t)lhs->low * (rhs >> 64);
    current += (uint64_t)t;
    swift_uint128_t next = t >> 64;
    t = (swift_uint128_t)lhs->mid * (uint64_t)rhs;
    current += (uint64_t)t;
    next += t >> 64;
    // Round up by adding UINT128_MAX (upper half)
    current += UINT64_MAX;

    current = next + (current >> 64);
    t = (swift_uint128_t)lhs->mid * (rhs >> 64);
    current += (uint64_t)t;
    next = t >> 64;
    t = (swift_uint128_t)lhs->high * (uint64_t)rhs;
    current += (uint64_t)t;
    next += t >> 64;
    lhs->low = (uint64_t)current;

    current = next + (current >> 64);
    t = (swift_uint128_t)lhs->high * (rhs >> 64);
    current += t;
    lhs->mid = (uint64_t)current;
    lhs->high = (uint64_t)(current >> 64);
#else
    uint64_t a, b, c, d; // temporaries
    // Six 32-bit values multiplied by 4 32-bit values.  Oh my.
    static const uint64_t mask32 = UINT32_MAX;
    uint64_t t = (uint64_t)lhs->low * rhs.low + mask32;
    t >>= 32;
    a = (uint64_t)lhs->low * rhs.b;
    b = (uint64_t)lhs->b * rhs.low;
    t += (a & mask32) + (b & mask32) + mask32;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    a = (uint64_t)lhs->low * rhs.c;
    b = (uint64_t)lhs->b * rhs.b;
    c = (uint64_t)lhs->c * rhs.low;
    t += (a & mask32) + (b & mask32) + (c & mask32) + mask32;
    t >>= 32;
    t += (a >> 32) + (b >> 32) + (c >> 32);
    a = (uint64_t)lhs->low * rhs.high;
    b = (uint64_t)lhs->b * rhs.c;
    c = (uint64_t)lhs->c * rhs.b;
    d = (uint64_t)lhs->d * rhs.low;
    t += (a & mask32) + (b & mask32) + (c & mask32) + (d & mask32) + mask32;
    t >>= 32;
    t += (a >> 32) + (b >> 32) + (c >> 32) + (d >> 32);
    a = (uint64_t)lhs->b * rhs.high;
    b = (uint64_t)lhs->c * rhs.c;
    c = (uint64_t)lhs->d * rhs.b;
    d = (uint64_t)lhs->e * rhs.low;
    t += (a & mask32) + (b & mask32) + (c & mask32) + (d & mask32);
    lhs->low = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32) + (c >> 32) + (d >> 32);
    a = (uint64_t)lhs->c * rhs.high;
    b = (uint64_t)lhs->d * rhs.c;
    c = (uint64_t)lhs->e * rhs.b;
    d = (uint64_t)lhs->high * rhs.low;
    t += (a & mask32) + (b & mask32) + (c & mask32) + (d & mask32);
    lhs->b = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32) + (c >> 32) + (d >> 32);
    a = (uint64_t)lhs->d * rhs.high;
    b = (uint64_t)lhs->e * rhs.c;
    c = (uint64_t)lhs->high * rhs.b;
    t += (a & mask32) + (b & mask32) + (c & mask32);
    lhs->c = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32) + (c >> 32);
    a = (uint64_t)lhs->e * rhs.high;
    b = (uint64_t)lhs->high * rhs.c;
    t += (a & mask32) + (b & mask32);
    lhs->d = t;
    t >>= 32;
    t += (a >> 32) + (b >> 32);
    t += (uint64_t)lhs->high * rhs.high;
    lhs->e = t;
    lhs->high = t >> 32;
#endif
}

// Subtract two 192-bit integers or fractions.
static void subtract192x192(swift_uint192_t *lhs, swift_uint192_t rhs) {
#if HAVE_UINT128_T
    swift_uint128_t t = (swift_uint128_t)lhs->low + (~rhs.low) + 1;
    lhs->low = t;
    t = (t >> 64) + lhs->mid + (~rhs.mid);
    lhs->mid = t;
    lhs->high += (t >> 64) + (~rhs.high);
#else
    uint64_t t = (uint64_t)lhs->low + (~rhs.low) + 1;
    lhs->low = t;
    t = (t >> 32) + lhs->b + (~rhs.b);
    lhs->b = t;
    t = (t >> 32) + lhs->c + (~rhs.c);
    lhs->c = t;
    t = (t >> 32) + lhs->d + (~rhs.d);
    lhs->d = t;
    t = (t >> 32) + lhs->e + (~rhs.e);
    lhs->e = t;
    lhs->high += (t >> 32) + (~rhs.high);
#endif
}

// Compare two 192-bit integers or fractions.
static int isLessThan192x192(swift_uint192_t lhs, swift_uint192_t rhs) {
#if HAVE_UINT128_T
    return (lhs.high < rhs.high)
        || (lhs.high == rhs.high
            && (lhs.mid < rhs.mid
                || (lhs.mid == rhs.mid
                    && lhs.low < rhs.low)));
#else
    return (lhs.high < rhs.high
            || (lhs.high == rhs.high
                && (lhs.e < rhs.e
                    || (lhs.e == rhs.e
                        && (lhs.d < rhs.d
                            || (lhs.d == rhs.d
                                && (lhs.c < rhs.c
                                    || (lhs.c == rhs.c
                                        && (lhs.b < rhs.b
                                            || (lhs.b == rhs.b
                                                && (lhs.low < rhs.low)))))))))));
#endif
}

// Shift a 192-bit integer right, rounding down.
static void shiftRightRoundingDown192(swift_uint192_t *lhs, int shift) {
#if HAVE_UINT128_T
    __uint128_t t = (__uint128_t)lhs->low >> shift;
    t += ((__uint128_t)lhs->mid << (64 - shift));
    lhs->low = t;
    t >>= 64;
    t += ((__uint128_t)lhs->high << (64 - shift));
    lhs->mid = t;
    t >>= 64;
    lhs->high = t;
#else
    uint64_t t = (uint64_t)lhs->low >> shift;
    t += ((uint64_t)lhs->b << (32 - shift));
    lhs->low = t;
    t >>= 32;
    t += ((uint64_t)lhs->c << (32 - shift));
    lhs->b = t;
    t >>= 32;
    t += ((uint64_t)lhs->d << (32 - shift));
    lhs->c = t;
    t >>= 32;
    t += ((uint64_t)lhs->e << (32 - shift));
    lhs->d = t;
    t >>= 32;
    t += ((uint64_t)lhs->high << (32 - shift));
    lhs->e = t;
    t >>= 32;
    lhs->high = t;
#endif
}

// Shift a 192-bit integer right, rounding up.
// Note: The shift will always be less than 20.  Someday, that
// might suggest a way to further optimize this.
static void shiftRightRoundingUp192(swift_uint192_t *lhs, int shift) {
#if HAVE_UINT128_T
    const uint64_t bias = (1 << shift) - 1;
    __uint128_t t = ((__uint128_t)lhs->low + bias) >> shift;
    t += ((__uint128_t)lhs->mid << (64 - shift));
    lhs->low = t;
    t >>= 64;
    t += ((__uint128_t)lhs->high << (64 - shift));
    lhs->mid = t;
    t >>= 64;
    lhs->high = t;
#else
    const uint64_t bias = (1 << shift) - 1;
    uint64_t t = ((uint64_t)lhs->low + bias) >> shift;
    t += ((uint64_t)lhs->b << (32 - shift));
    lhs->low = t;
    t >>= 32;
    t += ((uint64_t)lhs->c << (32 - shift));
    lhs->b = t;
    t >>= 32;
    t += ((uint64_t)lhs->d << (32 - shift));
    lhs->c = t;
    t >>= 32;
    t += ((uint64_t)lhs->e << (32 - shift));
    lhs->d = t;
    t >>= 32;
    t += ((uint64_t)lhs->high << (32 - shift));
    lhs->e = t;
    t >>= 32;
    lhs->high = t;
#endif
}
#endif

//
// ------------ Power of 10 calculation ----------------
//

//
// ------------  Power-of-10 tables. --------------------------
//
// Grisu-style algorithms rely on being able to rapidly
// find a high-precision approximation of any power of 10.
// These values were computed by a simple script that
// relied on Python's excellent variable-length
// integer support.

#if SWIFT_DTOA_FLOAT_SUPPORT || SWIFT_DTOA_DOUBLE_SUPPORT || SWIFT_DTOA_FLOAT80_SUPPORT
// Float table
//
// The constant powers of 10 here represent pure fractions
// with a binary point at the far left. (Each number in
// this first table is implicitly divided by 2^64.)
//
// Table size: 320 bytes
//
// A 64-bit significand allows us to exactly represent
// powers of 10 up to 10^27.  For larger powers, the
// value here is rounded DOWN from the exact value.
// For those powers, the value here is less than the
// exact power of 10; adding one gives a value greater
// than the exact power of 10.
//
// For single-precision Float, we use these directly
// for positive powers of 10.  For negative powers of
// ten, we multiply a value here by 10^-40.
//
// For Double and Float80, we use the 28 exact values
// here to help reduce the size of those tables.
static const uint64_t powersOf10_Float[40] = {
    0x8000000000000000, // x 2^1 == 10^0 exactly
    0xa000000000000000, // x 2^4 == 10^1 exactly
    0xc800000000000000, // x 2^7 == 10^2 exactly
    0xfa00000000000000, // x 2^10 == 10^3 exactly
    0x9c40000000000000, // x 2^14 == 10^4 exactly
    0xc350000000000000, // x 2^17 == 10^5 exactly
    0xf424000000000000, // x 2^20 == 10^6 exactly
    0x9896800000000000, // x 2^24 == 10^7 exactly
    0xbebc200000000000, // x 2^27 == 10^8 exactly
    0xee6b280000000000, // x 2^30 == 10^9 exactly
    0x9502f90000000000, // x 2^34 == 10^10 exactly
    0xba43b74000000000, // x 2^37 == 10^11 exactly
    0xe8d4a51000000000, // x 2^40 == 10^12 exactly
    0x9184e72a00000000, // x 2^44 == 10^13 exactly
    0xb5e620f480000000, // x 2^47 == 10^14 exactly
    0xe35fa931a0000000, // x 2^50 == 10^15 exactly
    0x8e1bc9bf04000000, // x 2^54 == 10^16 exactly
    0xb1a2bc2ec5000000, // x 2^57 == 10^17 exactly
    0xde0b6b3a76400000, // x 2^60 == 10^18 exactly
    0x8ac7230489e80000, // x 2^64 == 10^19 exactly
    0xad78ebc5ac620000, // x 2^67 == 10^20 exactly
    0xd8d726b7177a8000, // x 2^70 == 10^21 exactly
    0x878678326eac9000, // x 2^74 == 10^22 exactly
    0xa968163f0a57b400, // x 2^77 == 10^23 exactly
    0xd3c21bcecceda100, // x 2^80 == 10^24 exactly
    0x84595161401484a0, // x 2^84 == 10^25 exactly
    0xa56fa5b99019a5c8, // x 2^87 == 10^26 exactly
    0xcecb8f27f4200f3a, // x 2^90 == 10^27 exactly
    0x813f3978f8940984, // x 2^94 ~= 10^28
    0xa18f07d736b90be5, // x 2^97 ~= 10^29
    0xc9f2c9cd04674ede, // x 2^100 ~= 10^30
    0xfc6f7c4045812296, // x 2^103 ~= 10^31
    0x9dc5ada82b70b59d, // x 2^107 ~= 10^32
    0xc5371912364ce305, // x 2^110 ~= 10^33
    0xf684df56c3e01bc6, // x 2^113 ~= 10^34
    0x9a130b963a6c115c, // x 2^117 ~= 10^35
    0xc097ce7bc90715b3, // x 2^120 ~= 10^36
    0xf0bdc21abb48db20, // x 2^123 ~= 10^37
    0x96769950b50d88f4, // x 2^127 ~= 10^38
    0xbc143fa4e250eb31, // x 2^130 ~= 10^39
};
#endif

#if SWIFT_DTOA_DOUBLE_SUPPORT
// As above, but with 128-bit fractions.
//
// Table size: 464 bytes
//
// We only store every 28th power of ten here.
// We can multiply by an exact 64-bit power of
// ten from the table above to reconstruct the
// significand for any power of 10.
static const uint64_t powersOf10_Double[] = {
    // low-order half, high-order half
    0x3931b850df08e738, 0x95fe7e07c91efafa, // x 2^-1328 ~= 10^-400
    0xba954f8e758fecb3, 0x9774919ef68662a3, // x 2^-1235 ~= 10^-372
    0x9028bed2939a635c, 0x98ee4a22ecf3188b, // x 2^-1142 ~= 10^-344
    0x47b233c92125366e, 0x9a6bb0aa55653b2d, // x 2^-1049 ~= 10^-316
    0x4ee367f9430aec32, 0x9becce62836ac577, // x 2^-956 ~= 10^-288
    0x6f773fc3603db4a9, 0x9d71ac8fada6c9b5, // x 2^-863 ~= 10^-260
    0xc47bc5014a1a6daf, 0x9efa548d26e5a6e1, // x 2^-770 ~= 10^-232
    0x80e8a40eccd228a4, 0xa086cfcd97bf97f3, // x 2^-677 ~= 10^-204
    0xb8ada00e5a506a7c, 0xa21727db38cb002f, // x 2^-584 ~= 10^-176
    0xc13e60d0d2e0ebba, 0xa3ab66580d5fdaf5, // x 2^-491 ~= 10^-148
    0xc2974eb4ee658828, 0xa54394fe1eedb8fe, // x 2^-398 ~= 10^-120
    0xcb4ccd500f6bb952, 0xa6dfbd9fb8e5b88e, // x 2^-305 ~= 10^-92
    0x3f2398d747b36224, 0xa87fea27a539e9a5, // x 2^-212 ~= 10^-64
    0xdde50bd1d5d0b9e9, 0xaa242499697392d2, // x 2^-119 ~= 10^-36
    0xfdc20d2b36ba7c3d, 0xabcc77118461cefc, // x 2^-26 ~= 10^-8
    0x0000000000000000, 0xad78ebc5ac620000, // x 2^67 == 10^20 exactly
    0x9670b12b7f410000, 0xaf298d050e4395d6, // x 2^160 == 10^48 exactly
    0x3b25a55f43294bcb, 0xb0de65388cc8ada8, // x 2^253 ~= 10^76
    0x58edec91ec2cb657, 0xb2977ee300c50fe7, // x 2^346 ~= 10^104
    0x29babe4598c311fb, 0xb454e4a179dd1877, // x 2^439 ~= 10^132
    0x577b986b314d6009, 0xb616a12b7fe617aa, // x 2^532 ~= 10^160
    0x0c11ed6d538aeb2f, 0xb7dcbf5354e9bece, // x 2^625 ~= 10^188
    0x6d953e2bd7173692, 0xb9a74a0637ce2ee1, // x 2^718 ~= 10^216
    0x9d6d1ad41abe37f1, 0xbb764c4ca7a4440f, // x 2^811 ~= 10^244
    0x4b2d8644d8a74e18, 0xbd49d14aa79dbc82, // x 2^904 ~= 10^272
    0xe0470a63e6bd56c3, 0xbf21e44003acdd2c, // x 2^997 ~= 10^300
    0x505f522e53053ff2, 0xc0fe908895cf3b44, // x 2^1090 ~= 10^328
    0xcca845ab2beafa9a, 0xc2dfe19c8c055535, // x 2^1183 ~= 10^356
    0x1027fff56784f444, 0xc4c5e310aef8aa17, // x 2^1276 ~= 10^384
};
#endif

#if SWIFT_DTOA_FLOAT80_SUPPORT
// Every 83rd power of 10 across the range of Float80
//
// Table size: 2,928 bytes
//
// Note: We could cut this in half at the cost of one additional
// 192-bit multiply by only storing the positive values and
// multiplying by 10^-5063 to obtain the negative ones, similar
// to how we handle Float above.
static const uint64_t powersOf10_Float80[] = {
    // low 64 bits, middle 64 bits, high 64 bits
    0x56825ada98468526, 0x0abc23fcfda07e29, 0x871db786569ca2dd, // x 2^-16818 ~= 10^-5063
    0xe3885beff5ee930d, 0xd1e638f68c97f0e5, 0xde90d1b113564507, // x 2^-16543 ~= 10^-4980
    0x5a7d22b5236bcad4, 0xbab3a28a6f0a489c, 0xb74ea21ab2946479, // x 2^-16267 ~= 10^-4897
    0x7d1f78bf0f4e2878, 0xcf4aea39615ffe6e, 0x96f9351fcfdd2686, // x 2^-15991 ~= 10^-4814
    0xad725c0d6c214314, 0x5bd5c19f18bd2857, 0xf8afafd6ef185238, // x 2^-15716 ~= 10^-4731
    0xe418e9217ce83755, 0x801e38463183fc88, 0xccd1ffc6bba63e21, // x 2^-15440 ~= 10^-4648
    0x4dcd52747e029d0c, 0x867b3096b1619df8, 0xa8b11ff4721d92fb, // x 2^-15164 ~= 10^-4565
    0xed2903f7e1df2b78, 0xc846664fe1364ee8, 0x8aefaaae9060380f, // x 2^-14888 ~= 10^-4482
    0xed7a7f4e1e171498, 0x7da1a627b88527f1, 0xe4dbb751faa311b0, // x 2^-14613 ~= 10^-4399
    0x320796dc9b1a158c, 0x2a11a871597b8012, 0xbc7d620092481a7e, // x 2^-14337 ~= 10^-4316
    0x796014ec6f4c0dcb, 0xcfa99f62903708d7, 0x9b3dee433c1311e9, // x 2^-14061 ~= 10^-4233
    0x08920ae76bdb8282, 0x952b06c385a08ff6, 0xffb7a402531fd4c9, // x 2^-13786 ~= 10^-4150
    0x18faa162f2c4d6b9, 0x050be8c5d21c6db6, 0xd29c7528965ae5bd, // x 2^-13510 ~= 10^-4067
    0x576a6c1abab7f7e7, 0xc05fb0b5c550f28d, 0xad7617610634129e, // x 2^-13234 ~= 10^-3984
    0x6cc3b2fb9cae4875, 0xe1bd5b09b7202157, 0x8edd4417ae3dc210, // x 2^-12958 ~= 10^-3901
    0xfafc1dc8fd12a6f8, 0xc3f29d230036529b, 0xeb542860da9bc7d8, // x 2^-12683 ~= 10^-3818
    0x9d198c56bc799f35, 0x42960adf9591f02a, 0xc1d1a4b6bc2eafc8, // x 2^-12407 ~= 10^-3735
    0x1803d096e43ff4bc, 0xdef9759d6432dab7, 0x9fa18c5de65a16fb, // x 2^-12131 ~= 10^-3652
    0x74872779576a577c, 0x9150140eb5101a96, 0x83793dfc83fd2f9b, // x 2^-11855 ~= 10^-3569
    0x415d0667f0f88262, 0x4898d98d1314d99f, 0xd890d45a257f4644, // x 2^-11580 ~= 10^-3486
    0x30a5256a610c1c72, 0x6ebca4d0365d504d, 0xb25d93f98145bdab, // x 2^-11304 ~= 10^-3403
    0xfb149b1f86a46376, 0xb5143323a7f8e16c, 0x92e74be10524c389, // x 2^-11028 ~= 10^-3320
    0x7b7532e25fead4c8, 0x0df6ab8ac6a0ec2f, 0xf1fb6e82e4df4a77, // x 2^-10753 ~= 10^-3237
    0x3b738d6a3caae67a, 0x346b2dd31826cfed, 0xc74c79bce7fe949f, // x 2^-10477 ~= 10^-3154
    0x22d3a12777cee527, 0xe185ac46f6ef1993, 0xa424ef0bb5ad3129, // x 2^-10201 ~= 10^-3071
    0xb7b6bccb9f60adec, 0x30ad7df78fe30cc8, 0x8730d40821cd89f3, // x 2^-9925 ~= 10^-2988
    0x21049149d72d44d5, 0x1e86debc54dd290d, 0xdeb04cb82aec22cb, // x 2^-9650 ~= 10^-2905
    0xeb69d287f5e7f920, 0x8d7e84a13801d034, 0xb7688f9800d26b3a, // x 2^-9374 ~= 10^-2822
    0x47b3da9aeff7df71, 0x82678774d6c0ac59, 0x970e8fd25ead01e3, // x 2^-9098 ~= 10^-2739
    0x50d3e92e2e4f8210, 0x9492db3d978aaca8, 0xf8d2dcaf37504b51, // x 2^-8823 ~= 10^-2656
    0xf4ce72058f777a4c, 0xb9eb2c585e924efe, 0xcceef83fdedcc506, // x 2^-8547 ~= 10^-2573
    0xb0e624a35b791884, 0x7104a98dbbb38b94, 0xa8c8fc3b03c3d1ed, // x 2^-8271 ~= 10^-2490
    0x6c94ecd8291e2ac9, 0x978b03821014b68c, 0x8b03518396007c08, // x 2^-7995 ~= 10^-2407
    0x96475208daa03ee3, 0x10eaa1481b149e5a, 0xe4fc163319551441, // x 2^-7720 ~= 10^-2324
    0xd709a820ff4ac847, 0x75aab7cb5cb15414, 0xbc980b270680156a, // x 2^-7444 ~= 10^-2241
    0xf273bca6b4de9e24, 0xb5025d88d4b252e1, 0x9b53e384eccabcf7, // x 2^-7168 ~= 10^-2158
    0x6c55a0603a928f40, 0x9e20db16dfb6b461, 0xffdbcf7251089796, // x 2^-6893 ~= 10^-2075
    0x9006db4d43cffe42, 0x9b4bca4cd6cec2db, 0xd2ba3f510a3aa638, // x 2^-6617 ~= 10^-1992
    0xa6b3c457fd0cd4d6, 0x28a4de91ba868fbf, 0xad8ea05a5f27642a, // x 2^-6341 ~= 10^-1909
    0xe7b14ed140f8d98e, 0x7b2f7d61ce5d426c, 0x8ef179291b6f5424, // x 2^-6065 ~= 10^-1826
    0x4a964d052fd03e10, 0x06897060bf491e6e, 0xeb75718d285cd8bf, // x 2^-5790 ~= 10^-1743
    0x22b2270f0e8dd87c, 0xa8510fa2f5a9e4de, 0xc1ed0ed498f7c54c, // x 2^-5514 ~= 10^-1660
    0x09102915726a9905, 0x5a0eb896edc89b54, 0x9fb8208d65ea5eda, // x 2^-5238 ~= 10^-1577
    0xb80d5f481d01deb9, 0x673f2aa50486f5ba, 0x838bd699b7c539e6, // x 2^-4962 ~= 10^-1494
    0x668b62b20ec2633b, 0x8682604c7123f859, 0xd8af761f94d2db2c, // x 2^-4687 ~= 10^-1411
    0xb2adaed8559cc199, 0x712339ba54f12372, 0xb276ce87987995d5, // x 2^-4411 ~= 10^-1328
    0x6beb873308685711, 0xac1ce34246ed56ad, 0x92fc133455668c02, // x 2^-4135 ~= 10^-1245
    0x593293d68a2261bc, 0x3c368f9497ca075d, 0xf21da89a29fa1c61, // x 2^-3860 ~= 10^-1162
    0x854051f9f0e4ca66, 0x8c5d5a234eda57f7, 0xc768aa46d6d1b675, // x 2^-3584 ~= 10^-1079
    0x333d09b2299c5e6b, 0xcf1f49c33399c5ac, 0xa43c26a751d4f7e7, // x 2^-3308 ~= 10^-996
    0x25a440d8b1620532, 0x274ebc67c3e21943, 0x8743f33df0feed29, // x 2^-3032 ~= 10^-913
    0x3ca95e3deb5be648, 0x52d18ccca1c558c2, 0xdecfcc3329238dd8, // x 2^-2757 ~= 10^-830
    0x59d1a7704af3acd7, 0xfae7722c6af19467, 0xb78280c024488353, // x 2^-2481 ~= 10^-747
    0x78813f3e80148049, 0x73b2baf13aa1c233, 0x9723ed8a28baf5ac, // x 2^-2205 ~= 10^-664
    0xf296a8198aa40fb8, 0x235532b08487fe6a, 0xf8f60e812de0cd7d, // x 2^-1930 ~= 10^-581
    0xa7fbdcb40b4f648f, 0x4f20ba9a64a7f6e7, 0xcd0bf4d206072167, // x 2^-1654 ~= 10^-498
    0x8dbf63ea468c724f, 0xa0e25c08b5c189d6, 0xa8e0dbe18ffb82cf, // x 2^-1378 ~= 10^-415
    0x765995c6cfd406ce, 0x4c3bcb5021afcc31, 0x8b16fb203055ac76, // x 2^-1102 ~= 10^-332
    0xe1d09ab6fb409872, 0x82b7e12780e7401a, 0xe51c79a85916f484, // x 2^-827 ~= 10^-249
    0x89cbe2422f9e1df9, 0x7415d448f6b6f0e7, 0xbcb2b812db11a5de, // x 2^-551 ~= 10^-166
    0x605fe83842e4d290, 0xc986afbe3ee11aba, 0x9b69dbe1b548ce7c, // x 2^-275 ~= 10^-83
    0x0000000000000000, 0x0000000000000000, 0x8000000000000000, // x 2^1 == 10^0 exactly
    0x0d6953169e1c7a1e, 0xf50a3fa490c30190, 0xd2d80db02aabd62b, // x 2^276 ~= 10^83
    0x3720b80c7d8ee39d, 0xaf561aa79a10ae6a, 0xada72ccc20054ae9, // x 2^552 ~= 10^166
    0x7cb3f026a212df74, 0x29cb4d87f2a7400e, 0x8f05b1163ba6832d, // x 2^828 ~= 10^249
    0x7dda22f9451d28a4, 0xe41c5bd18c57e88f, 0xeb96bf6ebadf77d8, // x 2^1103 ~= 10^332
    0xd5da00e6e2d05e5d, 0x5e510c5a752f0f8e, 0xc2087cd3215a16ad, // x 2^1379 ~= 10^415
    0x5603ba353e0b2fac, 0x48bbddc4d7359e49, 0x9fceb7ee780436f0, // x 2^1655 ~= 10^498
    0x15ceea8df15e47c7, 0x6a83c85cf158c652, 0x839e71d847c1779e, // x 2^1931 ~= 10^581
    0x514478d1fcd48eea, 0x3a4181cdda0d6e24, 0xd8ce1c3a2fffaea7, // x 2^2206 ~= 10^664
    0xe8b634620f1062be, 0x7304c7fb8a2f8a8a, 0xb2900ca735bdf121, // x 2^2482 ~= 10^747
    0xc3ec2fd9302c9bda, 0x729a6a7e830e1cf2, 0x9310dd78089bd66f, // x 2^2758 ~= 10^830
    0x1750ef5f751be079, 0x52ccabc96fc88a23, 0xf23fe788c763dffa, // x 2^3033 ~= 10^913
    0xaa80925ec1c80b65, 0x97681c548ff6c12f, 0xc784decd820a6180, // x 2^3309 ~= 10^996
    0xb4212d4b435a2317, 0x8df0a55abbb2c99a, 0xa453618b9dfd92db, // x 2^3585 ~= 10^1079
    0x939c2fedd434642a, 0x7d7de34bf5aa96b4, 0x875715282612729b, // x 2^3861 ~= 10^1162
    0xb7d9a0e46bbebb36, 0x3b2057dea52d686b, 0xdeef5022af37f1f6, // x 2^4136 ~= 10^1245
    0x931a74148ea64e59, 0xfe7fe67bd1074d0c, 0xb79c7593a1c17df0, // x 2^4412 ~= 10^1328
    0xabd20df0f1f1ad54, 0x0fff83cc7fa6b77b, 0x97394e479b6573b1, // x 2^4688 ~= 10^1411
    0x25f2467421674b7a, 0x828ff55a248bc026, 0xf919454d86f16685, // x 2^4963 ~= 10^1494
    0xe32dbd7131e6ab7d, 0xf674dd4821982084, 0xcd28f57dc585d094, // x 2^5239 ~= 10^1577
    0x866ab816a532b07d, 0xdc567471f9639b4e, 0xa8f8bee890f905c7, // x 2^5515 ~= 10^1660
    0xaca3a975993a2626, 0xc41cf207a71d87e4, 0x8b2aa784c405e2f1, // x 2^5791 ~= 10^1743
    0x731f0b7d820918bd, 0x355fde18e8448607, 0xe53ce1b25fb31788, // x 2^6066 ~= 10^1826
    0x0be7c29568db3f20, 0xc1328a3f1bf4d2b8, 0xbccd68c49888be61, // x 2^6342 ~= 10^1909
    0x9c6e0b1b927b7d3f, 0xffc9b96619da642a, 0x9b7fd75a060350cd, // x 2^6618 ~= 10^1992
    0x2a84c8fb4bd2edc9, 0xa679df45d339389b, 0x80121ad60ca2c518, // x 2^6894 ~= 10^2075
    0x0d4648d0876cf1c3, 0x48c67661c087fb5a, 0xd2f5e0469040e0eb, // x 2^7169 ~= 10^2158
    0xfe2d99a281a011ac, 0x65c13361e6b2c078, 0xadbfbcb6c676a69b, // x 2^7445 ~= 10^2241
    0xf4ec157aa4147562, 0xac89bfa5e79484a6, 0x8f19ebdf7661e3e9, // x 2^7721 ~= 10^2324
    0xe945f8c80090be1f, 0x9b8672e64aadbed2, 0xebb812063c9e01db, // x 2^7996 ~= 10^2407
    0xca12d8ad3b36d2e4, 0xd252322ea50ad274, 0xc223eeb2e1bde452, // x 2^8272 ~= 10^2490
    0xf554fb41e5b3e384, 0x977acb4d4af624fc, 0x9fe55281904ba38b, // x 2^8548 ~= 10^2573
    0xf3f69093398e2573, 0x111ae5735ec0e878, 0x83b10fb893300cde, // x 2^8824 ~= 10^2656
    0x30f65a8da0d10429, 0x1eecf4cf8a0b25f5, 0xd8ecc6aa93e876fc, // x 2^9099 ~= 10^2739
    0xa577f5f0c9f1e5a7, 0x91a5430ed623abf0, 0xb2a94e58da4930c3, // x 2^9375 ~= 10^2822
    0x202d2f87585ec0d7, 0x7a63589863efd480, 0x9325aaac89304b57, // x 2^9651 ~= 10^2905
    0x049544fbba3c01a1, 0x53accb0f60ac6095, 0xf2622b4f6c68d6ce, // x 2^9926 ~= 10^2988
    0x268a0d5f9d5ce861, 0xfe40703e1a91de57, 0xc7a117517a09153b, // x 2^10202 ~= 10^3071
    0xb6cd470a2a3b1d63, 0x5f963916b20ea587, 0xa46a9fb9111003bc, // x 2^10478 ~= 10^3154
    0xa3101c09fd8e6e96, 0xa2c6328011db5211, 0x876a39c722f798a7, // x 2^10754 ~= 10^3237
    0x8507e5fdb0ec5d83, 0x7ce93cc7f8feeed4, 0xdf0ed8875e7b8914, // x 2^11029 ~= 10^3320
    0xe19b7ebe4c7bfbca, 0x40930d1129943838, 0xb7b66e12fe1af499, // x 2^11305 ~= 10^3403
    0xc90c4ec15c21a357, 0xd91c86512d147305, 0x974eb20b241a65f6, // x 2^11581 ~= 10^3486
    0xb90341199c02a4eb, 0x69e684f53db6e8ce, 0xf93c8114f6c31f8a, // x 2^11856 ~= 10^3569
    0xa873f1318cef91cb, 0xbf8718466b31a7ca, 0xcd45fa43b1ce4c8e, // x 2^12132 ~= 10^3652
    0xacfe0dcc5262e273, 0x6e1bbb68662fd27a, 0xa910a550810203f8, // x 2^12408 ~= 10^3735
    0x59e7c8921bbe3758, 0x834743c5eab7dcea, 0x8b3e56b1b5c57589, // x 2^12684 ~= 10^3818
    0x145eed64fda2e6af, 0x1c605bdcc764238f, 0xe55d4e51d30b5592, // x 2^12959 ~= 10^3901
    0xb747164b17268ea2, 0xd8aa19f1d85da07d, 0xbce81d3cc784a1ca, // x 2^13235 ~= 10^3984
    0x3666af1cb2f0356b, 0xc8dd55687a68bb70, 0x9b95d5ee4f80366d, // x 2^13511 ~= 10^4067
    0x70d67261b5bde1e9, 0x1d76f2d15166ec20, 0x8024383bab19730d, // x 2^13787 ~= 10^4150
    0x084a3ba0b748546a, 0xc67f9026f83dca47, 0xd313b714d3a1c65e, // x 2^14062 ~= 10^4233
    0x4411a8127eea085e, 0x441eb397ffcdab0d, 0xadd8501ad0361d15, // x 2^14338 ~= 10^4316
    0x7b62a54ed6233032, 0x75458a1c8300e014, 0x8f2e2985332eae98, // x 2^14614 ~= 10^4399
    0x162d5b51a1dd9594, 0x655bb1b7aa4e8196, 0xebd96954582af06f, // x 2^14889 ~= 10^4482
    0x55e6c62f920d3682, 0x79fd57cf7c37941c, 0xc23f6474669f4abe, // x 2^15165 ~= 10^4565
    0x19482fa0ac45669c, 0x803c1cd864033781, 0x9ffbf04722750449, // x 2^15441 ~= 10^4648
    0xa412d1f95f4624cd, 0xc95abe9ce589e048, 0x83c3b03af95c9674, // x 2^15717 ~= 10^4731
    0xc1207e487c57b4e1, 0xf93dd2c7669a8ed1, 0xd90b75715d861b38, // x 2^15992 ~= 10^4814
    0xeb20d9a25e0372bd, 0xb5073df6adc221b4, 0xb2c2939d0763fcac, // x 2^16268 ~= 10^4897
    0x1a648c339e28cc45, 0xbd14f0fa3e24b6ae, 0x933a7ad2419ea0b5, // x 2^16544 ~= 10^4980
};
#endif

#if SWIFT_DTOA_FLOAT_SUPPORT || SWIFT_DTOA_DOUBLE_SUPPORT || SWIFT_DTOA_FLOAT80_SUPPORT
// The power-of-10 tables do not directly store the associated binary
// exponent.  That's because the binary exponent is a simple linear
// function of the decimal power (and vice versa), so it's just as
// fast (and uses much less memory) to compute it:

// The binary exponent corresponding to a particular power of 10.
// This matches the power-of-10 tables across the full range of Float80.
static int binaryExponentFor10ToThe(int p) {
    return (int)(((((int64_t)p) * 55732705) >> 24) + 1);
}

// A decimal exponent that approximates a particular binary power.
static int decimalExponentFor2ToThe(int e) {
    return (int)(((int64_t)e * 20201781) >> 26);
}
#endif

#if SWIFT_DTOA_FLOAT_SUPPORT
// Given a power `p`, this returns three values:
// * 64-bit fractions `lower` and `upper`
// * integer `exponent`
//
// The returned values satisty the following:
// ```
//    lower * 2^exponent <= 10^p <= upper * 2^exponent
// ```
//
// Note: Max(*upper - *lower) = 3
static void intervalContainingPowerOf10_Float(int p, uint64_t *lower, uint64_t *upper, int *exponent) {
    if (p < 0) {
        uint64_t base = powersOf10_Float[p + 40];
        int baseExponent = binaryExponentFor10ToThe(p + 40);
        uint64_t tenToTheMinus40 = 0x8b61313bbabce2c6; // x 2^-132 ~= 10^-40
        *upper = multiply64x64RoundingUp(base + 1, tenToTheMinus40 + 1);
        *lower = multiply64x64RoundingDown(base, tenToTheMinus40);
        *exponent = baseExponent - 132;
    } else {
        uint64_t base = powersOf10_Float[p];
        *upper = base + 1;
        *lower = base;
        *exponent = binaryExponentFor10ToThe(p);
    }
}
#endif

#if SWIFT_DTOA_DOUBLE_SUPPORT
// As above, but returning 128-bit fractions suitable for
// converting doubles.
static void intervalContainingPowerOf10_Double(int p, swift_uint128_t *lower, swift_uint128_t *upper, int *exponent) {
    if (p >= 0 && p <= 54) {
        if (p <= 27) {
            // Use one 64-bit exact value
            swift_uint128_t exact;
            initialize128WithHigh64(exact, powersOf10_Float[p]);
            *upper = exact;
            *lower = exact;
            *exponent = binaryExponentFor10ToThe(p);
            return;
        } else {
            // Multiply two 64-bit exact values to get a 128-bit exact value
            swift_uint128_t base;
            initialize128WithHigh64(base, powersOf10_Float[p - 27]);
            int baseExponent = binaryExponentFor10ToThe(p - 27);
            uint64_t extra = powersOf10_Float[27];
            int extraExponent = binaryExponentFor10ToThe(27);
            swift_uint128_t exact = multiply128x64RoundingDown(base, extra);
            *upper = exact;
            *lower = exact;
            *exponent = baseExponent + extraExponent;
            return;
        }
    }

    // Multiply a 128-bit approximate value with a 64-bit exact value
    int index = p + 400;
    // Copy a pair of uint64_t into a swift_uint128_t
    const uint64_t *base_p = powersOf10_Double + (index / 28) * 2;
    swift_uint128_t base;
    initialize128WithHighLow64(base, base_p[1], base_p[0]);
    int extraPower = index % 28;
    int baseExponent = binaryExponentFor10ToThe(p - extraPower);

    int e = baseExponent;
    if (extraPower > 0) {
        int64_t extra = powersOf10_Float[extraPower];
        e += binaryExponentFor10ToThe(extraPower);
        *lower = multiply128x64RoundingDown(base, extra);
        increment128(base);
        *upper = multiply128x64RoundingUp(base, extra);
    } else {
        *lower = base;
        increment128(base);
        *upper = base;
    }
    *exponent = e;
}
#endif

#if SWIFT_DTOA_FLOAT80_SUPPORT
// As above, but returning 192-bit fractions suitable for
// converting float80.
static void intervalContainingPowerOf10_Float80(int p, swift_uint192_t *lower, swift_uint192_t *upper, int *exponent) {
    if (p >= 0 && p <= 27) {
        // We have an exact form, return a zero-width interval.
        uint64_t exact = powersOf10_Float[p];
        initialize192WithHighMidLow64(*upper, exact, 0, 0);
        initialize192WithHighMidLow64(*lower, exact, 0, 0);
        *exponent = binaryExponentFor10ToThe(p);
        return;
    }

    int index = p + 5063;
    const uint64_t *base_p = powersOf10_Float80 + (index / 83) * 3;
    // Note: The low-order value in the Float80 table above
    // is never UINT64_MAX, so there's never a carry from
    // the increment here.
    initialize192WithHighMidLow64(*upper, base_p[2], base_p[1], base_p[0] + 1);
    initialize192WithHighMidLow64(*lower, base_p[2], base_p[1], base_p[0]);
    int extraPower = index % 83;
    int e = binaryExponentFor10ToThe(p - extraPower);

    while (extraPower > 27) {
        uint64_t power27 = powersOf10_Float[27];
        multiply192x64RoundingDown(lower, power27);
        multiply192x64RoundingUp(upper, power27);
        e += binaryExponentFor10ToThe(27);
        extraPower -= 27;
    }
    if (extraPower > 0) {
        uint64_t extra = powersOf10_Float[extraPower];
        multiply192x64RoundingDown(lower, extra);
        multiply192x64RoundingUp(upper, extra);
        e += binaryExponentFor10ToThe(extraPower);
    }
    *exponent = e;
}
#endif