static_assert(sizeof(Target) == sizeof(Source), "size mismatch");

Target target;

std::memcpy(&target, &source, sizeof(Source));

return target;

{

static constexpr int kPrecision = 64; // = q

std::uint64_t f = 0;

int e = 0;

constexpr diyfp(std::uint64_t f_, int e_) noexcept : f(f_), e(e_) {}

/*!

static diyfp sub(const diyfp &x, const diyfp &y) noexcept {

return {x.f - y.f, x.e};

}

/*!

static diyfp mul(const diyfp &x, const diyfp &y) noexcept {

static_assert(kPrecision == 64, "internal error");

// Computes:

// f = round((x.f * y.f) / 2^q)

// e = x.e + y.e + q

// Emulate the 64-bit * 64-bit multiplication:

//

// p = u * v

// = (u_lo + 2^32 u_hi) (v_lo + 2^32 v_hi)

// = (u_lo v_lo ) + 2^32 ((u_lo v_hi ) + (u_hi v_lo )) +

// 2^64 (u_hi v_hi ) = (p0 ) + 2^32 ((p1 ) + (p2 ))

// + 2^64 (p3 ) = (p0_lo + 2^32 p0_hi) + 2^32 ((p1_lo +

// 2^32 p1_hi) + (p2_lo + 2^32 p2_hi)) + 2^64 (p3 ) =

// (p0_lo ) + 2^32 (p0_hi + p1_lo + p2_lo ) + 2^64 (p1_hi +

// p2_hi + p3) = (p0_lo ) + 2^32 (Q ) + 2^64 (H ) = (p0_lo ) +

// 2^32 (Q_lo + 2^32 Q_hi ) + 2^64 (H )

//

// (Since Q might be larger than 2^32 - 1)

//

// = (p0_lo + 2^32 Q_lo) + 2^64 (Q_hi + H)

//

// (Q_hi + H does not overflow a 64-bit int)

//

// = p_lo + 2^64 p_hi

const std::uint64_t u_lo = x.f & 0xFFFFFFFFu;

const std::uint64_t u_hi = x.f >> 32u;

const std::uint64_t v_lo = y.f & 0xFFFFFFFFu;

const std::uint64_t v_hi = y.f >> 32u;

const std::uint64_t p0 = u_lo * v_lo;

const std::uint64_t p1 = u_lo * v_hi;

const std::uint64_t p2 = u_hi * v_lo;

const std::uint64_t p3 = u_hi * v_hi;

const std::uint64_t p0_hi = p0 >> 32u;

const std::uint64_t p1_lo = p1 & 0xFFFFFFFFu;

const std::uint64_t p1_hi = p1 >> 32u;

const std::uint64_t p2_lo = p2 & 0xFFFFFFFFu;

const std::uint64_t p2_hi = p2 >> 32u;

std::uint64_t Q = p0_hi + p1_lo + p2_lo;

// The full product might now be computed as

//

// p_hi = p3 + p2_hi + p1_hi + (Q >> 32)

// p_lo = p0_lo + (Q << 32)

//

// But in this particular case here, the full p_lo is not required.

// Effectively we only need to add the highest bit in p_lo to p_hi (and

// Q_hi + 1 does not overflow).

Q += std::uint64_t{1} << (64u - 32u - 1u); // round, ties up

const std::uint64_t h = p3 + p2_hi + p1_hi + (Q >> 32u);

return {h, x.e + y.e + 64};

}

/*!

static diyfp normalize(diyfp x) noexcept {

while ((x.f >> 63u) == 0) {

x.f <<= 1u;

x.e--;

}

return x;

}

/*!

static diyfp normalize_to(const diyfp &x,

const int target_exponent) noexcept {

const int delta = x.e - target_exponent;

return {x.f << delta, target_exponent};

}

diyfp w;

diyfp minus;

diyfp plus;

// Convert the IEEE representation into a diyfp.

//

// If v is denormal:

// value = 0.F * 2^(1 - bias) = ( F) * 2^(1 - bias - (p-1))

// If v is normalized:

// value = 1.F * 2^(E - bias) = (2^(p-1) + F) * 2^(E - bias - (p-1))

static_assert(std::numeric_limits<FloatType>::is_iec559,

"internal error: dtoa_short requires an IEEE-754 "

"floating-point implementation");

constexpr int kPrecision =

std::numeric_limits<FloatType>::digits; // = p (includes the hidden bit)

constexpr int kBias =

std::numeric_limits<FloatType>::max_exponent - 1 + (kPrecision - 1);

constexpr int kMinExp = 1 - kBias;

constexpr std::uint64_t kHiddenBit = std::uint64_t{1}

<< (kPrecision - 1); // = 2^(p-1)

using bits_type = typename std::conditional<kPrecision == 24, std::uint32_t,

std::uint64_t>::type;

const std::uint64_t bits = reinterpret_bits<bits_type>(value);

const std::uint64_t E = bits >> (kPrecision - 1);

const std::uint64_t F = bits & (kHiddenBit - 1);

const bool is_denormal = E == 0;

const diyfp v = is_denormal

? diyfp(F, kMinExp)

: diyfp(F + kHiddenBit, static_cast<int>(E) - kBias);

// Compute the boundaries m- and m+ of the floating-point value

// v = f * 2^e.

//

// Determine v- and v+, the floating-point predecessor and successor if v,

// respectively.

//

// v- = v - 2^e if f != 2^(p-1) or e == e_min (A)

// = v - 2^(e-1) if f == 2^(p-1) and e > e_min (B)

//

// v+ = v + 2^e

//

// Let m- = (v- + v) / 2 and m+ = (v + v+) / 2. All real numbers _strictly_

// between m- and m+ round to v, regardless of how the input rounding

// algorithm breaks ties.

//

// ---+-------------+-------------+-------------+-------------+--- (A)

// v- m- v m+ v+

//

// -----------------+------+------+-------------+-------------+--- (B)

// v- m- v m+ v+

const bool lower_boundary_is_closer = F == 0 && E > 1;

const diyfp m_plus = diyfp(2 * v.f + 1, v.e - 1);

const diyfp m_minus = lower_boundary_is_closer

? diyfp(4 * v.f - 1, v.e - 2) // (B)

: diyfp(2 * v.f - 1, v.e - 1); // (A)

// Determine the normalized w+ = m+.

const diyfp w_plus = diyfp::normalize(m_plus);

// Determine w- = m- such that e_(w-) = e_(w+).

const diyfp w_minus = diyfp::normalize_to(m_minus, w_plus.e);

return {diyfp::normalize(v), w_minus, w_plus};

{

std::uint64_t f;

int e;

int k;

// Now

//

// alpha <= e_c + e + q <= gamma (1)

// ==> f_c * 2^alpha <= c * 2^e * 2^q

//

// and since the c's are normalized, 2^(q-1) <= f_c,

//

// ==> 2^(q - 1 + alpha) <= c * 2^(e + q)

// ==> 2^(alpha - e - 1) <= c

//

// If c were an exact power of ten, i.e. c = 10^k, one may determine k as

//

// k = ceil( log_10( 2^(alpha - e - 1) ) )

// = ceil( (alpha - e - 1) * log_10(2) )

//

// From the paper:

// "In theory the result of the procedure could be wrong since c is rounded,

// and the computation itself is approximated [...]. In practice, however,

// this simple function is sufficient."

//

// For IEEE double precision floating-point numbers converted into

// normalized diyfp's w = f * 2^e, with q = 64,

//

// e >= -1022 (min IEEE exponent)

// -52 (p - 1)

// -52 (p - 1, possibly normalize denormal IEEE numbers)

// -11 (normalize the diyfp)

// = -1137

//

// and

//

// e <= +1023 (max IEEE exponent)

// -52 (p - 1)

// -11 (normalize the diyfp)

// = 960

//

// This binary exponent range [-1137,960] results in a decimal exponent

// range [-307,324]. One does not need to store a cached power for each

// k in this range. For each such k it suffices to find a cached power

// such that the exponent of the product lies in [alpha,gamma].

// This implies that the difference of the decimal exponents of adjacent

// table entries must be less than or equal to

//

// floor( (gamma - alpha) * log_10(2) ) = 8.

//

// (A smaller distance gamma-alpha would require a larger table.)

// NB:

// Actually this function returns c, such that -60 <= e_c + e + 64 <= -34.

constexpr int kCachedPowersMinDecExp = -300;

constexpr int kCachedPowersDecStep = 8;

static constexpr std::array<cached_power, 79> kCachedPowers = {{

{0xAB70FE17C79AC6CA, -1060, -300}, {0xFF77B1FCBEBCDC4F, -1034, -292},

{0xBE5691EF416BD60C, -1007, -284}, {0x8DD01FAD907FFC3C, -980, -276},

{0xD3515C2831559A83, -954, -268}, {0x9D71AC8FADA6C9B5, -927, -260},

{0xEA9C227723EE8BCB, -901, -252}, {0xAECC49914078536D, -874, -244},

{0x823C12795DB6CE57, -847, -236}, {0xC21094364DFB5637, -821, -228},

{0x9096EA6F3848984F, -794, -220}, {0xD77485CB25823AC7, -768, -212},

{0xA086CFCD97BF97F4, -741, -204}, {0xEF340A98172AACE5, -715, -196},

{0xB23867FB2A35B28E, -688, -188}, {0x84C8D4DFD2C63F3B, -661, -180},

{0xC5DD44271AD3CDBA, -635, -172}, {0x936B9FCEBB25C996, -608, -164},

{0xDBAC6C247D62A584, -582, -156}, {0xA3AB66580D5FDAF6, -555, -148},

{0xF3E2F893DEC3F126, -529, -140}, {0xB5B5ADA8AAFF80B8, -502, -132},

{0x87625F056C7C4A8B, -475, -124}, {0xC9BCFF6034C13053, -449, -116},

{0x964E858C91BA2655, -422, -108}, {0xDFF9772470297EBD, -396, -100},

{0xA6DFBD9FB8E5B88F, -369, -92}, {0xF8A95FCF88747D94, -343, -84},

{0xB94470938FA89BCF, -316, -76}, {0x8A08F0F8BF0F156B, -289, -68},

{0xCDB02555653131B6, -263, -60}, {0x993FE2C6D07B7FAC, -236, -52},

{0xE45C10C42A2B3B06, -210, -44}, {0xAA242499697392D3, -183, -36},

{0xFD87B5F28300CA0E, -157, -28}, {0xBCE5086492111AEB, -130, -20},

{0x8CBCCC096F5088CC, -103, -12}, {0xD1B71758E219652C, -77, -4},

{0x9C40000000000000, -50, 4}, {0xE8D4A51000000000, -24, 12},

{0xAD78EBC5AC620000, 3, 20}, {0x813F3978F8940984, 30, 28},

{0xC097CE7BC90715B3, 56, 36}, {0x8F7E32CE7BEA5C70, 83, 44},

{0xD5D238A4ABE98068, 109, 52}, {0x9F4F2726179A2245, 136, 60},

{0xED63A231D4C4FB27, 162, 68}, {0xB0DE65388CC8ADA8, 189, 76},

{0x83C7088E1AAB65DB, 216, 84}, {0xC45D1DF942711D9A, 242, 92},

{0x924D692CA61BE758, 269, 100}, {0xDA01EE641A708DEA, 295, 108},

{0xA26DA3999AEF774A, 322, 116}, {0xF209787BB47D6B85, 348, 124},

{0xB454E4A179DD1877, 375, 132}, {0x865B86925B9BC5C2, 402, 140},

{0xC83553C5C8965D3D, 428, 148}, {0x952AB45CFA97A0B3, 455, 156},

{0xDE469FBD99A05FE3, 481, 164}, {0xA59BC234DB398C25, 508, 172},

{0xF6C69A72A3989F5C, 534, 180}, {0xB7DCBF5354E9BECE, 561, 188},

{0x88FCF317F22241E2, 588, 196}, {0xCC20CE9BD35C78A5, 614, 204},

{0x98165AF37B2153DF, 641, 212}, {0xE2A0B5DC971F303A, 667, 220},

{0xA8D9D1535CE3B396, 694, 228}, {0xFB9B7CD9A4A7443C, 720, 236},

{0xBB764C4CA7A44410, 747, 244}, {0x8BAB8EEFB6409C1A, 774, 252},

{0xD01FEF10A657842C, 800, 260}, {0x9B10A4E5E9913129, 827, 268},

{0xE7109BFBA19C0C9D, 853, 276}, {0xAC2820D9623BF429, 880, 284},

{0x80444B5E7AA7CF85, 907, 292}, {0xBF21E44003ACDD2D, 933, 300},

{0x8E679C2F5E44FF8F, 960, 308}, {0xD433179D9C8CB841, 986, 316},

{0x9E19DB92B4E31BA9, 1013, 324},

}};

// This computation gives exactly the same results for k as

// k = ceil((kAlpha - e - 1) * 0.30102999566398114)

// for |e| <= 1500, but doesn't require floating-point operations.

// NB: log_10(2) ~= 78913 / 2^18

const int f = kAlpha - e - 1;

const int k = (f * 78913) / (1 << 18) + static_cast<int>(f > 0);

const int index = (-kCachedPowersMinDecExp + k + (kCachedPowersDecStep - 1)) /

kCachedPowersDecStep;

const cached_power cached = kCachedPowers[static_cast<std::size_t>(index)];

return cached;

// LCOV_EXCL_START

if (n >= 1000000000) {

pow10 = 1000000000;

return 10;

}

// LCOV_EXCL_STOP

else if (n >= 100000000) {

pow10 = 100000000;

return 9;

} else if (n >= 10000000) {

pow10 = 10000000;

return 8;

} else if (n >= 1000000) {

pow10 = 1000000;

return 7;

} else if (n >= 100000) {

pow10 = 100000;

return 6;

} else if (n >= 10000) {

pow10 = 10000;

return 5;

} else if (n >= 1000) {

pow10 = 1000;

return 4;

} else if (n >= 100) {

pow10 = 100;

return 3;

} else if (n >= 10) {

pow10 = 10;

return 2;

} else {

pow10 = 1;

return 1;

}

std::uint64_t delta, std::uint64_t rest,

std::uint64_t ten_k) {

// <--------------------------- delta ---->

// <---- dist --------->

// --------------[------------------+-------------------]--------------

// M- w M+

//

// ten_k

// <------>

// <---- rest ---->

// --------------[------------------+----+--------------]--------------

// w V

// = buf * 10^k

//

// ten_k represents a unit-in-the-last-place in the decimal representation

// stored in buf.

// Decrement buf by ten_k while this takes buf closer to w.

// The tests are written in this order to avoid overflow in unsigned

// integer arithmetic.

while (rest < dist && delta - rest >= ten_k &&

(rest + ten_k < dist || dist - rest > rest + ten_k - dist)) {

buf[len - 1]--;

rest += ten_k;

}

diyfp M_minus, diyfp w, diyfp M_plus) {

static_assert(kAlpha >= -60, "internal error");

static_assert(kGamma <= -32, "internal error");

// Generates the digits (and the exponent) of a decimal floating-point

// number V = buffer * 10^decimal_exponent in the range [M-, M+]. The diyfp's

// w, M- and M+ share the same exponent e, which satisfies alpha <= e <=

// gamma.

//

// <--------------------------- delta ---->

// <---- dist --------->

// --------------[------------------+-------------------]--------------

// M- w M+

//

// Grisu2 generates the digits of M+ from left to right and stops as soon as

// V is in [M-,M+].

std::uint64_t delta =

diyfp::sub(M_plus, M_minus)

.f; // (significand of (M+ - M-), implicit exponent is e)

std::uint64_t dist =

diyfp::sub(M_plus, w)

.f; // (significand of (M+ - w ), implicit exponent is e)

// Split M+ = f * 2^e into two parts p1 and p2 (note: e < 0):

//

// M+ = f * 2^e

// = ((f div 2^-e) * 2^-e + (f mod 2^-e)) * 2^e

// = ((p1 ) * 2^-e + (p2 )) * 2^e

// = p1 + p2 * 2^e

const diyfp one(std::uint64_t{1} << -M_plus.e, M_plus.e);

auto p1 = static_cast<std::uint32_t>(

M_plus.f >>

-one.e); // p1 = f div 2^-e (Since -e >= 32, p1 fits into a 32-bit int.)

std::uint64_t p2 = M_plus.f & (one.f - 1); // p2 = f mod 2^-e

// 1)

//

// Generate the digits of the integral part p1 = d[n-1]...d[1]d[0]

std::uint32_t pow10;

const int k = find_largest_pow10(p1, pow10);

// 10^(k-1) <= p1 < 10^k, pow10 = 10^(k-1)

//

// p1 = (p1 div 10^(k-1)) * 10^(k-1) + (p1 mod 10^(k-1))

// = (d[k-1] ) * 10^(k-1) + (p1 mod 10^(k-1))

//

// M+ = p1 + p2 * 2^e

// = d[k-1] * 10^(k-1) + (p1 mod 10^(k-1)) + p2 * 2^e

// = d[k-1] * 10^(k-1) + ((p1 mod 10^(k-1)) * 2^-e + p2) * 2^e

// = d[k-1] * 10^(k-1) + ( rest) * 2^e

//

// Now generate the digits d[n] of p1 from left to right (n = k-1,...,0)

//

// p1 = d[k-1]...d[n] * 10^n + d[n-1]...d[0]

//

// but stop as soon as

//

// rest * 2^e = (d[n-1]...d[0] * 2^-e + p2) * 2^e <= delta * 2^e

int n = k;

while (n > 0) {

// Invariants:

// M+ = buffer * 10^n + (p1 + p2 * 2^e) (buffer = 0 for n = k)

// pow10 = 10^(n-1) <= p1 < 10^n

//

const std::uint32_t d = p1 / pow10; // d = p1 div 10^(n-1)

const std::uint32_t r = p1 % pow10; // r = p1 mod 10^(n-1)

//

// M+ = buffer * 10^n + (d * 10^(n-1) + r) + p2 * 2^e

// = (buffer * 10 + d) * 10^(n-1) + (r + p2 * 2^e)

//

buffer[length++] = static_cast<char>('0' + d); // buffer := buffer * 10 + d

//

// M+ = buffer * 10^(n-1) + (r + p2 * 2^e)

//

p1 = r;

n--;

//

// M+ = buffer * 10^n + (p1 + p2 * 2^e)

// pow10 = 10^n

//

// Now check if enough digits have been generated.

// Compute

//

// p1 + p2 * 2^e = (p1 * 2^-e + p2) * 2^e = rest * 2^e

//

// Note:

// Since rest and delta share the same exponent e, it suffices to

// compare the significands.

const std::uint64_t rest = (std::uint64_t{p1} << -one.e) + p2;

if (rest <= delta) {

// V = buffer * 10^n, with M- <= V <= M+.

decimal_exponent += n;

// We may now just stop. But instead look if the buffer could be

// decremented to bring V closer to w.

//

// pow10 = 10^n is now 1 ulp in the decimal representation V.

// The rounding procedure works with diyfp's with an implicit

// exponent of e.

//

// 10^n = (10^n * 2^-e) * 2^e = ulp * 2^e

//

const std::uint64_t ten_n = std::uint64_t{pow10} << -one.e;

grisu2_round(buffer, length, dist, delta, rest, ten_n);

return;

}

pow10 /= 10;

//

// pow10 = 10^(n-1) <= p1 < 10^n

// Invariants restored.

}

// 2)

//

// The digits of the integral part have been generated:

//

// M+ = d[k-1]...d[1]d[0] + p2 * 2^e

// = buffer + p2 * 2^e

//

// Now generate the digits of the fractional part p2 * 2^e.

//

// Note:

// No decimal point is generated: the exponent is adjusted instead.

//

// p2 actually represents the fraction

//

// p2 * 2^e

// = p2 / 2^-e

// = d[-1] / 10^1 + d[-2] / 10^2 + ...

//

// Now generate the digits d[-m] of p1 from left to right (m = 1,2,...)

//

// p2 * 2^e = d[-1]d[-2]...d[-m] * 10^-m

// + 10^-m * (d[-m-1] / 10^1 + d[-m-2] / 10^2 + ...)

//

// using

//

// 10^m * p2 = ((10^m * p2) div 2^-e) * 2^-e + ((10^m * p2) mod 2^-e)

// = ( d) * 2^-e + ( r)

//

// or

// 10^m * p2 * 2^e = d + r * 2^e

//

// i.e.

//

// M+ = buffer + p2 * 2^e

// = buffer + 10^-m * (d + r * 2^e)

// = (buffer * 10^m + d) * 10^-m + 10^-m * r * 2^e

//

// and stop as soon as 10^-m * r * 2^e <= delta * 2^e

int m = 0;

for (;;) {

// Invariant:

// M+ = buffer * 10^-m + 10^-m * (d[-m-1] / 10 + d[-m-2] / 10^2 + ...)

// * 2^e

// = buffer * 10^-m + 10^-m * (p2 )

// * 2^e = buffer * 10^-m + 10^-m * (1/10 * (10 * p2) ) * 2^e =

// buffer * 10^-m + 10^-m * (1/10 * ((10*p2 div 2^-e) * 2^-e +

// (10*p2 mod 2^-e)) * 2^e

//

p2 *= 10;

const std::uint64_t d = p2 >> -one.e; // d = (10 * p2) div 2^-e

const std::uint64_t r = p2 & (one.f - 1); // r = (10 * p2) mod 2^-e

//

// M+ = buffer * 10^-m + 10^-m * (1/10 * (d * 2^-e + r) * 2^e

// = buffer * 10^-m + 10^-m * (1/10 * (d + r * 2^e))

// = (buffer * 10 + d) * 10^(-m-1) + 10^(-m-1) * r * 2^e

//

buffer[length++] = static_cast<char>('0' + d); // buffer := buffer * 10 + d

//

// M+ = buffer * 10^(-m-1) + 10^(-m-1) * r * 2^e

//

p2 = r;

m++;

//

// M+ = buffer * 10^-m + 10^-m * p2 * 2^e

// Invariant restored.

// Check if enough digits have been generated.

//

// 10^-m * p2 * 2^e <= delta * 2^e

// p2 * 2^e <= 10^m * delta * 2^e

// p2 <= 10^m * delta

delta *= 10;

dist *= 10;

if (p2 <= delta) {

break;

}

}

// V = buffer * 10^-m, with M- <= V <= M+.

decimal_exponent -= m;

// 1 ulp in the decimal representation is now 10^-m.

// Since delta and dist are now scaled by 10^m, we need to do the

// same with ulp in order to keep the units in sync.

//

// 10^m * 10^-m = 1 = 2^-e * 2^e = ten_m * 2^e

//

const std::uint64_t ten_m = one.f;

grisu2_round(buffer, length, dist, delta, p2, ten_m);

// By construction this algorithm generates the shortest possible decimal

// number (Loitsch, Theorem 6.2) which rounds back to w.

// For an input number of precision p, at least

//

// N = 1 + ceil(p * log_10(2))

//

// decimal digits are sufficient to identify all binary floating-point

// numbers (Matula, "In-and-Out conversions").

// This implies that the algorithm does not produce more than N decimal

// digits.

//

// N = 17 for p = 53 (IEEE double precision)

// N = 9 for p = 24 (IEEE single precision)

diyfp v, diyfp m_plus) {

// --------(-----------------------+-----------------------)-------- (A)

// m- v m+

//

// --------------------(-----------+-----------------------)-------- (B)

// m- v m+

//

// First scale v (and m- and m+) such that the exponent is in the range

// [alpha, gamma].

const cached_power cached = get_cached_power_for_binary_exponent(m_plus.e);

const diyfp c_minus_k(cached.f, cached.e); // = c ~= 10^-k

// The exponent of the products is = v.e + c_minus_k.e + q and is in the range

// [alpha,gamma]

const diyfp w = diyfp::mul(v, c_minus_k);

const diyfp w_minus = diyfp::mul(m_minus, c_minus_k);

const diyfp w_plus = diyfp::mul(m_plus, c_minus_k);

// ----(---+---)---------------(---+---)---------------(---+---)----

// w- w w+

// = c*m- = c*v = c*m+

//

// diyfp::mul rounds its result and c_minus_k is approximated too. w, w- and

// w+ are now off by a small amount.

// In fact:

//

// w - v * 10^k < 1 ulp

//

// To account for this inaccuracy, add resp. subtract 1 ulp.

//

// --------+---[---------------(---+---)---------------]---+--------

// w- M- w M+ w+

//

// Now any number in [M-, M+] (bounds included) will round to w when input,

// regardless of how the input rounding algorithm breaks ties.

//

// And digit_gen generates the shortest possible such number in [M-, M+].

// Note that this does not mean that Grisu2 always generates the shortest

// possible number in the interval (m-, m+).

const diyfp M_minus(w_minus.f + 1, w_minus.e);

const diyfp M_plus(w_plus.f - 1, w_plus.e);

decimal_exponent = -cached.k; // = -(-k) = k

grisu2_digit_gen(buf, len, decimal_exponent, M_minus, w, M_plus);

static_assert(diyfp::kPrecision >= std::numeric_limits<FloatType>::digits + 3,

"internal error: not enough precision");

// If the neighbors (and boundaries) of 'value' are always computed for

// double-precision numbers, all float's can be recovered using strtod (and

// strtof). However, the resulting decimal representations are not exactly

// "short".

//

// The documentation for 'std::to_chars'

// (https://en.cppreference.com/w/cpp/utility/to_chars) says "value is

// converted to a string as if by std::sprintf in the default ("C") locale"

// and since sprintf promotes float's to double's, I think this is exactly

// what 'std::to_chars' does. On the other hand, the documentation for

// 'std::to_chars' requires that "parsing the representation using the

// corresponding std::from_chars function recovers value exactly". That

// indicates that single precision floating-point numbers should be recovered

// using 'std::strtof'.

//

// NB: If the neighbors are computed for single-precision numbers, there is a

// single float

// (7.0385307e-26f) which can't be recovered using strtod. The resulting

// double precision value is off by 1 ulp.

#if 0

const boundaries w = compute_boundaries(static_cast<double>(value));

#else

const boundaries w = compute_boundaries(value);

#endif

grisu2(buf, len, decimal_exponent, w.minus, w.w, w.plus);