1 //===----------------------------------------------------------------------===//
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
7 //===----------------------------------------------------------------------===//
9 // Copyright (c) Microsoft Corporation.
10 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
12 // Copyright 2018 Ulf Adams
13 // Copyright (c) Microsoft Corporation. All rights reserved.
15 // Boost Software License - Version 1.0 - August 17th, 2003
17 // Permission is hereby granted, free of charge, to any person or organization
18 // obtaining a copy of the software and accompanying documentation covered by
19 // this license (the "Software") to use, reproduce, display, distribute,
20 // execute, and transmit the Software, and to prepare derivative works of the
21 // Software, and to permit third-parties to whom the Software is furnished to
22 // do so, all subject to the following:
24 // The copyright notices in the Software and this entire statement, including
25 // the above license grant, this restriction and the following disclaimer,
26 // must be included in all copies of the Software, in whole or in part, and
27 // all derivative works of the Software, unless such copies or derivative
28 // works are solely in the form of machine-executable object code generated by
29 // a source language processor.
31 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
32 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
33 // FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
34 // SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
35 // FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
36 // ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
37 // DEALINGS IN THE SOFTWARE.
39 // Avoid formatting to keep the changes with the original code minimal.
46 #include "include/ryu/common.h"
47 #include "include/ryu/d2fixed.h"
48 #include "include/ryu/d2s.h"
49 #include "include/ryu/d2s_full_table.h"
50 #include "include/ryu/d2s_intrinsics.h"
51 #include "include/ryu/digit_table.h"
52 #include "include/ryu/ryu.h"
54 _LIBCPP_BEGIN_NAMESPACE_STD
56 // We need a 64x128-bit multiplication and a subsequent 128-bit shift.
58 // The 64-bit factor is variable and passed in, the 128-bit factor comes
59 // from a lookup table. We know that the 64-bit factor only has 55
60 // significant bits (i.e., the 9 topmost bits are zeros). The 128-bit
61 // factor only has 124 significant bits (i.e., the 4 topmost bits are
64 // In principle, the multiplication result requires 55 + 124 = 179 bits to
65 // represent. However, we then shift this value to the right by __j, which is
66 // at least __j >= 115, so the result is guaranteed to fit into 179 - 115 = 64
67 // bits. This means that we only need the topmost 64 significant bits of
68 // the 64x128-bit multiplication.
70 // There are several ways to do this:
71 // 1. Best case: the compiler exposes a 128-bit type.
72 // We perform two 64x64-bit multiplications, add the higher 64 bits of the
73 // lower result to the higher result, and shift by __j - 64 bits.
75 // We explicitly cast from 64-bit to 128-bit, so the compiler can tell
76 // that these are only 64-bit inputs, and can map these to the best
77 // possible sequence of assembly instructions.
78 // x64 machines happen to have matching assembly instructions for
79 // 64x64-bit multiplications and 128-bit shifts.
81 // 2. Second best case: the compiler exposes intrinsics for the x64 assembly
82 // instructions mentioned in 1.
84 // 3. We only have 64x64 bit instructions that return the lower 64 bits of
85 // the result, i.e., we have to use plain C.
86 // Our inputs are less than the full width, so we have three options:
87 // a. Ignore this fact and just implement the intrinsics manually.
88 // b. Split both into 31-bit pieces, which guarantees no internal overflow,
89 // but requires extra work upfront (unless we change the lookup table).
90 // c. Split only the first factor into 31-bit pieces, which also guarantees
91 // no internal overflow, but requires extra work since the intermediate
92 // results are not perfectly aligned.
93 #ifdef _LIBCPP_INTRINSIC128
95 [[nodiscard
]] _LIBCPP_HIDE_FROM_ABI
inline uint64_t __mulShift(const uint64_t __m
, const uint64_t* const __mul
, const int32_t __j
) {
96 // __m is maximum 55 bits
97 uint64_t __high1
; // 128
98 const uint64_t __low1
= __ryu_umul128(__m
, __mul
[1], &__high1
); // 64
99 uint64_t __high0
; // 64
100 (void) __ryu_umul128(__m
, __mul
[0], &__high0
); // 0
101 const uint64_t __sum
= __high0
+ __low1
;
102 if (__sum
< __high0
) {
103 ++__high1
; // overflow into __high1
105 return __ryu_shiftright128(__sum
, __high1
, static_cast<uint32_t>(__j
- 64));
108 [[nodiscard
]] _LIBCPP_HIDE_FROM_ABI
inline uint64_t __mulShiftAll(const uint64_t __m
, const uint64_t* const __mul
, const int32_t __j
,
109 uint64_t* const __vp
, uint64_t* const __vm
, const uint32_t __mmShift
) {
110 *__vp
= __mulShift(4 * __m
+ 2, __mul
, __j
);
111 *__vm
= __mulShift(4 * __m
- 1 - __mmShift
, __mul
, __j
);
112 return __mulShift(4 * __m
, __mul
, __j
);
115 #else // ^^^ intrinsics available ^^^ / vvv intrinsics unavailable vvv
117 [[nodiscard
]] _LIBCPP_HIDE_FROM_ABI
inline _LIBCPP_ALWAYS_INLINE
uint64_t __mulShiftAll(uint64_t __m
, const uint64_t* const __mul
, const int32_t __j
,
118 uint64_t* const __vp
, uint64_t* const __vm
, const uint32_t __mmShift
) { // TRANSITION, VSO-634761
120 // __m is maximum 55 bits
122 const uint64_t __lo
= __ryu_umul128(__m
, __mul
[0], &__tmp
);
124 const uint64_t __mid
= __tmp
+ __ryu_umul128(__m
, __mul
[1], &__hi
);
125 __hi
+= __mid
< __tmp
; // overflow into __hi
127 const uint64_t __lo2
= __lo
+ __mul
[0];
128 const uint64_t __mid2
= __mid
+ __mul
[1] + (__lo2
< __lo
);
129 const uint64_t __hi2
= __hi
+ (__mid2
< __mid
);
130 *__vp
= __ryu_shiftright128(__mid2
, __hi2
, static_cast<uint32_t>(__j
- 64 - 1));
132 if (__mmShift
== 1) {
133 const uint64_t __lo3
= __lo
- __mul
[0];
134 const uint64_t __mid3
= __mid
- __mul
[1] - (__lo3
> __lo
);
135 const uint64_t __hi3
= __hi
- (__mid3
> __mid
);
136 *__vm
= __ryu_shiftright128(__mid3
, __hi3
, static_cast<uint32_t>(__j
- 64 - 1));
138 const uint64_t __lo3
= __lo
+ __lo
;
139 const uint64_t __mid3
= __mid
+ __mid
+ (__lo3
< __lo
);
140 const uint64_t __hi3
= __hi
+ __hi
+ (__mid3
< __mid
);
141 const uint64_t __lo4
= __lo3
- __mul
[0];
142 const uint64_t __mid4
= __mid3
- __mul
[1] - (__lo4
> __lo3
);
143 const uint64_t __hi4
= __hi3
- (__mid4
> __mid3
);
144 *__vm
= __ryu_shiftright128(__mid4
, __hi4
, static_cast<uint32_t>(__j
- 64));
147 return __ryu_shiftright128(__mid
, __hi
, static_cast<uint32_t>(__j
- 64 - 1));
150 #endif // ^^^ intrinsics unavailable ^^^
152 [[nodiscard
]] _LIBCPP_HIDE_FROM_ABI
inline uint32_t __decimalLength17(const uint64_t __v
) {
153 // This is slightly faster than a loop.
154 // The average output length is 16.38 digits, so we check high-to-low.
155 // Function precondition: __v is not an 18, 19, or 20-digit number.
156 // (17 digits are sufficient for round-tripping.)
157 _LIBCPP_ASSERT_UNCATEGORIZED(__v
< 100000000000000000u, "");
158 if (__v
>= 10000000000000000u) { return 17; }
159 if (__v
>= 1000000000000000u) { return 16; }
160 if (__v
>= 100000000000000u) { return 15; }
161 if (__v
>= 10000000000000u) { return 14; }
162 if (__v
>= 1000000000000u) { return 13; }
163 if (__v
>= 100000000000u) { return 12; }
164 if (__v
>= 10000000000u) { return 11; }
165 if (__v
>= 1000000000u) { return 10; }
166 if (__v
>= 100000000u) { return 9; }
167 if (__v
>= 10000000u) { return 8; }
168 if (__v
>= 1000000u) { return 7; }
169 if (__v
>= 100000u) { return 6; }
170 if (__v
>= 10000u) { return 5; }
171 if (__v
>= 1000u) { return 4; }
172 if (__v
>= 100u) { return 3; }
173 if (__v
>= 10u) { return 2; }
177 // A floating decimal representing m * 10^e.
178 struct __floating_decimal_64
{
183 [[nodiscard
]] _LIBCPP_HIDE_FROM_ABI
inline __floating_decimal_64
__d2d(const uint64_t __ieeeMantissa
, const uint32_t __ieeeExponent
) {
186 if (__ieeeExponent
== 0) {
187 // We subtract 2 so that the bounds computation has 2 additional bits.
188 __e2
= 1 - __DOUBLE_BIAS
- __DOUBLE_MANTISSA_BITS
- 2;
189 __m2
= __ieeeMantissa
;
191 __e2
= static_cast<int32_t>(__ieeeExponent
) - __DOUBLE_BIAS
- __DOUBLE_MANTISSA_BITS
- 2;
192 __m2
= (1ull << __DOUBLE_MANTISSA_BITS
) | __ieeeMantissa
;
194 const bool __even
= (__m2
& 1) == 0;
195 const bool __acceptBounds
= __even
;
197 // Step 2: Determine the interval of valid decimal representations.
198 const uint64_t __mv
= 4 * __m2
;
199 // Implicit bool -> int conversion. True is 1, false is 0.
200 const uint32_t __mmShift
= __ieeeMantissa
!= 0 || __ieeeExponent
<= 1;
201 // We would compute __mp and __mm like this:
202 // uint64_t __mp = 4 * __m2 + 2;
203 // uint64_t __mm = __mv - 1 - __mmShift;
205 // Step 3: Convert to a decimal power base using 128-bit arithmetic.
206 uint64_t __vr
, __vp
, __vm
;
208 bool __vmIsTrailingZeros
= false;
209 bool __vrIsTrailingZeros
= false;
211 // I tried special-casing __q == 0, but there was no effect on performance.
212 // This expression is slightly faster than max(0, __log10Pow2(__e2) - 1).
213 const uint32_t __q
= __log10Pow2(__e2
) - (__e2
> 3);
214 __e10
= static_cast<int32_t>(__q
);
215 const int32_t __k
= __DOUBLE_POW5_INV_BITCOUNT
+ __pow5bits(static_cast<int32_t>(__q
)) - 1;
216 const int32_t __i
= -__e2
+ static_cast<int32_t>(__q
) + __k
;
217 __vr
= __mulShiftAll(__m2
, __DOUBLE_POW5_INV_SPLIT
[__q
], __i
, &__vp
, &__vm
, __mmShift
);
219 // This should use __q <= 22, but I think 21 is also safe. Smaller values
220 // may still be safe, but it's more difficult to reason about them.
221 // Only one of __mp, __mv, and __mm can be a multiple of 5, if any.
222 const uint32_t __mvMod5
= static_cast<uint32_t>(__mv
) - 5 * static_cast<uint32_t>(__div5(__mv
));
224 __vrIsTrailingZeros
= __multipleOfPowerOf5(__mv
, __q
);
225 } else if (__acceptBounds
) {
226 // Same as min(__e2 + (~__mm & 1), __pow5Factor(__mm)) >= __q
227 // <=> __e2 + (~__mm & 1) >= __q && __pow5Factor(__mm) >= __q
228 // <=> true && __pow5Factor(__mm) >= __q, since __e2 >= __q.
229 __vmIsTrailingZeros
= __multipleOfPowerOf5(__mv
- 1 - __mmShift
, __q
);
231 // Same as min(__e2 + 1, __pow5Factor(__mp)) >= __q.
232 __vp
-= __multipleOfPowerOf5(__mv
+ 2, __q
);
236 // This expression is slightly faster than max(0, __log10Pow5(-__e2) - 1).
237 const uint32_t __q
= __log10Pow5(-__e2
) - (-__e2
> 1);
238 __e10
= static_cast<int32_t>(__q
) + __e2
;
239 const int32_t __i
= -__e2
- static_cast<int32_t>(__q
);
240 const int32_t __k
= __pow5bits(__i
) - __DOUBLE_POW5_BITCOUNT
;
241 const int32_t __j
= static_cast<int32_t>(__q
) - __k
;
242 __vr
= __mulShiftAll(__m2
, __DOUBLE_POW5_SPLIT
[__i
], __j
, &__vp
, &__vm
, __mmShift
);
244 // {__vr,__vp,__vm} is trailing zeros if {__mv,__mp,__mm} has at least __q trailing 0 bits.
245 // __mv = 4 * __m2, so it always has at least two trailing 0 bits.
246 __vrIsTrailingZeros
= true;
247 if (__acceptBounds
) {
248 // __mm = __mv - 1 - __mmShift, so it has 1 trailing 0 bit iff __mmShift == 1.
249 __vmIsTrailingZeros
= __mmShift
== 1;
251 // __mp = __mv + 2, so it always has at least one trailing 0 bit.
254 } else if (__q
< 63) { // TRANSITION(ulfjack): Use a tighter bound here.
255 // We need to compute min(ntz(__mv), __pow5Factor(__mv) - __e2) >= __q - 1
256 // <=> ntz(__mv) >= __q - 1 && __pow5Factor(__mv) - __e2 >= __q - 1
257 // <=> ntz(__mv) >= __q - 1 (__e2 is negative and -__e2 >= __q)
258 // <=> (__mv & ((1 << (__q - 1)) - 1)) == 0
259 // We also need to make sure that the left shift does not overflow.
260 __vrIsTrailingZeros
= __multipleOfPowerOf2(__mv
, __q
- 1);
264 // Step 4: Find the shortest decimal representation in the interval of valid representations.
265 int32_t __removed
= 0;
266 uint8_t __lastRemovedDigit
= 0;
268 // On average, we remove ~2 digits.
269 if (__vmIsTrailingZeros
|| __vrIsTrailingZeros
) {
270 // General case, which happens rarely (~0.7%).
272 const uint64_t __vpDiv10
= __div10(__vp
);
273 const uint64_t __vmDiv10
= __div10(__vm
);
274 if (__vpDiv10
<= __vmDiv10
) {
277 const uint32_t __vmMod10
= static_cast<uint32_t>(__vm
) - 10 * static_cast<uint32_t>(__vmDiv10
);
278 const uint64_t __vrDiv10
= __div10(__vr
);
279 const uint32_t __vrMod10
= static_cast<uint32_t>(__vr
) - 10 * static_cast<uint32_t>(__vrDiv10
);
280 __vmIsTrailingZeros
&= __vmMod10
== 0;
281 __vrIsTrailingZeros
&= __lastRemovedDigit
== 0;
282 __lastRemovedDigit
= static_cast<uint8_t>(__vrMod10
);
288 if (__vmIsTrailingZeros
) {
290 const uint64_t __vmDiv10
= __div10(__vm
);
291 const uint32_t __vmMod10
= static_cast<uint32_t>(__vm
) - 10 * static_cast<uint32_t>(__vmDiv10
);
292 if (__vmMod10
!= 0) {
295 const uint64_t __vpDiv10
= __div10(__vp
);
296 const uint64_t __vrDiv10
= __div10(__vr
);
297 const uint32_t __vrMod10
= static_cast<uint32_t>(__vr
) - 10 * static_cast<uint32_t>(__vrDiv10
);
298 __vrIsTrailingZeros
&= __lastRemovedDigit
== 0;
299 __lastRemovedDigit
= static_cast<uint8_t>(__vrMod10
);
306 if (__vrIsTrailingZeros
&& __lastRemovedDigit
== 5 && __vr
% 2 == 0) {
307 // Round even if the exact number is .....50..0.
308 __lastRemovedDigit
= 4;
310 // We need to take __vr + 1 if __vr is outside bounds or we need to round up.
311 _Output
= __vr
+ ((__vr
== __vm
&& (!__acceptBounds
|| !__vmIsTrailingZeros
)) || __lastRemovedDigit
>= 5);
313 // Specialized for the common case (~99.3%). Percentages below are relative to this.
314 bool __roundUp
= false;
315 const uint64_t __vpDiv100
= __div100(__vp
);
316 const uint64_t __vmDiv100
= __div100(__vm
);
317 if (__vpDiv100
> __vmDiv100
) { // Optimization: remove two digits at a time (~86.2%).
318 const uint64_t __vrDiv100
= __div100(__vr
);
319 const uint32_t __vrMod100
= static_cast<uint32_t>(__vr
) - 100 * static_cast<uint32_t>(__vrDiv100
);
320 __roundUp
= __vrMod100
>= 50;
326 // Loop iterations below (approximately), without optimization above:
327 // 0: 0.03%, 1: 13.8%, 2: 70.6%, 3: 14.0%, 4: 1.40%, 5: 0.14%, 6+: 0.02%
328 // Loop iterations below (approximately), with optimization above:
329 // 0: 70.6%, 1: 27.8%, 2: 1.40%, 3: 0.14%, 4+: 0.02%
331 const uint64_t __vpDiv10
= __div10(__vp
);
332 const uint64_t __vmDiv10
= __div10(__vm
);
333 if (__vpDiv10
<= __vmDiv10
) {
336 const uint64_t __vrDiv10
= __div10(__vr
);
337 const uint32_t __vrMod10
= static_cast<uint32_t>(__vr
) - 10 * static_cast<uint32_t>(__vrDiv10
);
338 __roundUp
= __vrMod10
>= 5;
344 // We need to take __vr + 1 if __vr is outside bounds or we need to round up.
345 _Output
= __vr
+ (__vr
== __vm
|| __roundUp
);
347 const int32_t __exp
= __e10
+ __removed
;
349 __floating_decimal_64 __fd
;
350 __fd
.__exponent
= __exp
;
351 __fd
.__mantissa
= _Output
;
355 [[nodiscard
]] _LIBCPP_HIDE_FROM_ABI
inline to_chars_result
__to_chars(char* const _First
, char* const _Last
, const __floating_decimal_64 __v
,
356 chars_format _Fmt
, const double __f
) {
357 // Step 5: Print the decimal representation.
358 uint64_t _Output
= __v
.__mantissa
;
359 int32_t _Ryu_exponent
= __v
.__exponent
;
360 const uint32_t __olength
= __decimalLength17(_Output
);
361 int32_t _Scientific_exponent
= _Ryu_exponent
+ static_cast<int32_t>(__olength
) - 1;
363 if (_Fmt
== chars_format
{}) {
367 if (__olength
== 1) {
368 // Value | Fixed | Scientific
369 // 1e-3 | "0.001" | "1e-03"
370 // 1e4 | "10000" | "1e+04"
374 // Value | Fixed | Scientific
375 // 1234e-7 | "0.0001234" | "1.234e-04"
376 // 1234e5 | "123400000" | "1.234e+08"
377 _Lower
= -static_cast<int32_t>(__olength
+ 3);
381 if (_Lower
<= _Ryu_exponent
&& _Ryu_exponent
<= _Upper
) {
382 _Fmt
= chars_format::fixed
;
384 _Fmt
= chars_format::scientific
;
386 } else if (_Fmt
== chars_format::general
) {
387 // C11 7.21.6.1 "The fprintf function"/8:
388 // "Let P equal [...] 6 if the precision is omitted [...].
389 // Then, if a conversion with style E would have an exponent of X:
390 // - if P > X >= -4, the conversion is with style f [...].
391 // - otherwise, the conversion is with style e [...]."
392 if (-4 <= _Scientific_exponent
&& _Scientific_exponent
< 6) {
393 _Fmt
= chars_format::fixed
;
395 _Fmt
= chars_format::scientific
;
399 if (_Fmt
== chars_format::fixed
) {
400 // Example: _Output == 1729, __olength == 4
402 // _Ryu_exponent | Printed | _Whole_digits | _Total_fixed_length | Notes
403 // --------------|----------|---------------|----------------------|---------------------------------------
404 // 2 | 172900 | 6 | _Whole_digits | Ryu can't be used for printing
405 // 1 | 17290 | 5 | (sometimes adjusted) | when the trimmed digits are nonzero.
406 // --------------|----------|---------------|----------------------|---------------------------------------
407 // 0 | 1729 | 4 | _Whole_digits | Unified length cases.
408 // --------------|----------|---------------|----------------------|---------------------------------------
409 // -1 | 172.9 | 3 | __olength + 1 | This case can't happen for
410 // -2 | 17.29 | 2 | | __olength == 1, but no additional
411 // -3 | 1.729 | 1 | | code is needed to avoid it.
412 // --------------|----------|---------------|----------------------|---------------------------------------
413 // -4 | 0.1729 | 0 | 2 - _Ryu_exponent | C11 7.21.6.1 "The fprintf function"/8:
414 // -5 | 0.01729 | -1 | | "If a decimal-point character appears,
415 // -6 | 0.001729 | -2 | | at least one digit appears before it."
417 const int32_t _Whole_digits
= static_cast<int32_t>(__olength
) + _Ryu_exponent
;
419 uint32_t _Total_fixed_length
;
420 if (_Ryu_exponent
>= 0) { // cases "172900" and "1729"
421 _Total_fixed_length
= static_cast<uint32_t>(_Whole_digits
);
423 // Rounding can affect the number of digits.
424 // For example, 1e23 is exactly "99999999999999991611392" which is 23 digits instead of 24.
425 // We can use a lookup table to detect this and adjust the total length.
426 static constexpr uint8_t _Adjustment
[309] = {
427 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,1,1,0,1,0,1,1,1,0,1,1,1,0,0,0,0,0,
428 1,1,0,0,1,0,1,1,1,0,0,0,0,1,1,1,1,0,0,0,1,1,1,1,0,0,0,1,1,1,1,0,1,0,1,0,1,1,0,0,0,0,1,1,1,
429 1,0,0,0,0,0,0,0,1,1,0,1,1,0,0,1,0,1,0,1,0,1,1,0,0,0,0,0,1,1,1,0,0,1,1,1,1,1,0,1,0,1,1,0,1,
430 1,0,0,0,0,0,0,0,0,0,1,1,1,0,0,1,0,0,1,0,0,1,1,1,1,0,0,1,1,0,1,1,0,1,1,0,1,0,0,0,1,0,0,0,1,
431 0,1,0,1,0,1,1,1,0,0,0,0,0,0,1,1,1,1,0,0,1,0,1,1,1,0,0,0,1,0,1,1,1,1,1,1,0,1,0,1,1,0,0,0,1,
432 1,1,0,1,1,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,1,0,1,1,0,0,1,1,1,0,0,0,1,0,1,0,0,0,0,0,1,1,0,
433 0,1,0,1,1,1,0,0,1,0,0,0,0,1,0,1,0,0,0,0,0,1,0,1,0,1,1,0,1,0,0,0,0,0,1,1,0,1,0 };
434 _Total_fixed_length
-= _Adjustment
[_Ryu_exponent
];
435 // _Whole_digits doesn't need to be adjusted because these cases won't refer to it later.
437 } else if (_Whole_digits
> 0) { // case "17.29"
438 _Total_fixed_length
= __olength
+ 1;
439 } else { // case "0.001729"
440 _Total_fixed_length
= static_cast<uint32_t>(2 - _Ryu_exponent
);
443 if (_Last
- _First
< static_cast<ptrdiff_t>(_Total_fixed_length
)) {
444 return { _Last
, errc::value_too_large
};
448 if (_Ryu_exponent
> 0) { // case "172900"
451 if (_Ryu_exponent
> 22) { // 10^22 is the largest power of 10 that's exactly representable as a double.
452 _Can_use_ryu
= false;
454 // Ryu generated X: __v.__mantissa * 10^_Ryu_exponent
455 // __v.__mantissa == 2^_Trailing_zero_bits * (__v.__mantissa >> _Trailing_zero_bits)
456 // 10^_Ryu_exponent == 2^_Ryu_exponent * 5^_Ryu_exponent
458 // _Trailing_zero_bits is [0, 56] (aside: because 2^56 is the largest power of 2
459 // with 17 decimal digits, which is double's round-trip limit.)
460 // _Ryu_exponent is [1, 22].
461 // Normalization adds [2, 52] (aside: at least 2 because the pre-normalized mantissa is at least 5).
462 // This adds up to [3, 130], which is well below double's maximum binary exponent 1023.
464 // Therefore, we just need to consider (__v.__mantissa >> _Trailing_zero_bits) * 5^_Ryu_exponent.
466 // If that product would exceed 53 bits, then X can't be exactly represented as a double.
467 // (That's not a problem for round-tripping, because X is close enough to the original double,
468 // but X isn't mathematically equal to the original double.) This requires a high-precision fallback.
470 // If the product is 53 bits or smaller, then X can be exactly represented as a double (and we don't
471 // need to re-synthesize it; the original double must have been X, because Ryu wouldn't produce the
472 // same output for two different doubles X and Y). This allows Ryu's output to be used (zero-filled).
474 // (2^53 - 1) / 5^0 (for indexing), (2^53 - 1) / 5^1, ..., (2^53 - 1) / 5^22
475 static constexpr uint64_t _Max_shifted_mantissa
[23] = {
476 9007199254740991u, 1801439850948198u, 360287970189639u, 72057594037927u, 14411518807585u,
477 2882303761517u, 576460752303u, 115292150460u, 23058430092u, 4611686018u, 922337203u, 184467440u,
478 36893488u, 7378697u, 1475739u, 295147u, 59029u, 11805u, 2361u, 472u, 94u, 18u, 3u };
480 unsigned long _Trailing_zero_bits
;
481 #ifdef _LIBCPP_HAS_BITSCAN64
482 (void) _BitScanForward64(&_Trailing_zero_bits
, __v
.__mantissa
); // __v.__mantissa is guaranteed nonzero
483 #else // ^^^ 64-bit ^^^ / vvv 32-bit vvv
484 const uint32_t _Low_mantissa
= static_cast<uint32_t>(__v
.__mantissa
);
485 if (_Low_mantissa
!= 0) {
486 (void) _BitScanForward(&_Trailing_zero_bits
, _Low_mantissa
);
488 const uint32_t _High_mantissa
= static_cast<uint32_t>(__v
.__mantissa
>> 32); // nonzero here
489 (void) _BitScanForward(&_Trailing_zero_bits
, _High_mantissa
);
490 _Trailing_zero_bits
+= 32;
492 #endif // ^^^ 32-bit ^^^
493 const uint64_t _Shifted_mantissa
= __v
.__mantissa
>> _Trailing_zero_bits
;
494 _Can_use_ryu
= _Shifted_mantissa
<= _Max_shifted_mantissa
[_Ryu_exponent
];
498 // Print the integer exactly.
499 // Performance note: This will redundantly perform bounds checking.
500 // Performance note: This will redundantly decompose the IEEE representation.
501 return __d2fixed_buffered_n(_First
, _Last
, __f
, 0);
505 // Print the decimal digits, left-aligned within [_First, _First + _Total_fixed_length).
506 _Mid
= _First
+ __olength
;
507 } else { // cases "1729", "17.29", and "0.001729"
508 // Print the decimal digits, right-aligned within [_First, _First + _Total_fixed_length).
509 _Mid
= _First
+ _Total_fixed_length
;
512 // We prefer 32-bit operations, even on 64-bit platforms.
513 // We have at most 17 digits, and uint32_t can store 9 digits.
514 // If _Output doesn't fit into uint32_t, we cut off 8 digits,
515 // so the rest will fit into uint32_t.
516 if ((_Output
>> 32) != 0) {
517 // Expensive 64-bit division.
518 const uint64_t __q
= __div1e8(_Output
);
519 uint32_t __output2
= static_cast<uint32_t>(_Output
- 100000000 * __q
);
522 const uint32_t __c
= __output2
% 10000;
524 const uint32_t __d
= __output2
% 10000;
525 const uint32_t __c0
= (__c
% 100) << 1;
526 const uint32_t __c1
= (__c
/ 100) << 1;
527 const uint32_t __d0
= (__d
% 100) << 1;
528 const uint32_t __d1
= (__d
/ 100) << 1;
530 std::memcpy(_Mid
-= 2, __DIGIT_TABLE
+ __c0
, 2);
531 std::memcpy(_Mid
-= 2, __DIGIT_TABLE
+ __c1
, 2);
532 std::memcpy(_Mid
-= 2, __DIGIT_TABLE
+ __d0
, 2);
533 std::memcpy(_Mid
-= 2, __DIGIT_TABLE
+ __d1
, 2);
535 uint32_t __output2
= static_cast<uint32_t>(_Output
);
536 while (__output2
>= 10000) {
537 #ifdef __clang__ // TRANSITION, LLVM-38217
538 const uint32_t __c
= __output2
- 10000 * (__output2
/ 10000);
540 const uint32_t __c
= __output2
% 10000;
543 const uint32_t __c0
= (__c
% 100) << 1;
544 const uint32_t __c1
= (__c
/ 100) << 1;
545 std::memcpy(_Mid
-= 2, __DIGIT_TABLE
+ __c0
, 2);
546 std::memcpy(_Mid
-= 2, __DIGIT_TABLE
+ __c1
, 2);
548 if (__output2
>= 100) {
549 const uint32_t __c
= (__output2
% 100) << 1;
551 std::memcpy(_Mid
-= 2, __DIGIT_TABLE
+ __c
, 2);
553 if (__output2
>= 10) {
554 const uint32_t __c
= __output2
<< 1;
555 std::memcpy(_Mid
-= 2, __DIGIT_TABLE
+ __c
, 2);
557 *--_Mid
= static_cast<char>('0' + __output2
);
560 if (_Ryu_exponent
> 0) { // case "172900" with _Can_use_ryu
561 // Performance note: it might be more efficient to do this immediately after setting _Mid.
562 std::memset(_First
+ __olength
, '0', static_cast<size_t>(_Ryu_exponent
));
563 } else if (_Ryu_exponent
== 0) { // case "1729"
565 } else if (_Whole_digits
> 0) { // case "17.29"
566 // Performance note: moving digits might not be optimal.
567 std::memmove(_First
, _First
+ 1, static_cast<size_t>(_Whole_digits
));
568 _First
[_Whole_digits
] = '.';
569 } else { // case "0.001729"
570 // Performance note: a larger memset() followed by overwriting '.' might be more efficient.
573 std::memset(_First
+ 2, '0', static_cast<size_t>(-_Whole_digits
));
576 return { _First
+ _Total_fixed_length
, errc
{} };
579 const uint32_t _Total_scientific_length
= __olength
+ (__olength
> 1) // digits + possible decimal point
580 + (-100 < _Scientific_exponent
&& _Scientific_exponent
< 100 ? 4 : 5); // + scientific exponent
581 if (_Last
- _First
< static_cast<ptrdiff_t>(_Total_scientific_length
)) {
582 return { _Last
, errc::value_too_large
};
584 char* const __result
= _First
;
586 // Print the decimal digits.
588 // We prefer 32-bit operations, even on 64-bit platforms.
589 // We have at most 17 digits, and uint32_t can store 9 digits.
590 // If _Output doesn't fit into uint32_t, we cut off 8 digits,
591 // so the rest will fit into uint32_t.
592 if ((_Output
>> 32) != 0) {
593 // Expensive 64-bit division.
594 const uint64_t __q
= __div1e8(_Output
);
595 uint32_t __output2
= static_cast<uint32_t>(_Output
) - 100000000 * static_cast<uint32_t>(__q
);
598 const uint32_t __c
= __output2
% 10000;
600 const uint32_t __d
= __output2
% 10000;
601 const uint32_t __c0
= (__c
% 100) << 1;
602 const uint32_t __c1
= (__c
/ 100) << 1;
603 const uint32_t __d0
= (__d
% 100) << 1;
604 const uint32_t __d1
= (__d
/ 100) << 1;
605 std::memcpy(__result
+ __olength
- __i
- 1, __DIGIT_TABLE
+ __c0
, 2);
606 std::memcpy(__result
+ __olength
- __i
- 3, __DIGIT_TABLE
+ __c1
, 2);
607 std::memcpy(__result
+ __olength
- __i
- 5, __DIGIT_TABLE
+ __d0
, 2);
608 std::memcpy(__result
+ __olength
- __i
- 7, __DIGIT_TABLE
+ __d1
, 2);
611 uint32_t __output2
= static_cast<uint32_t>(_Output
);
612 while (__output2
>= 10000) {
613 #ifdef __clang__ // TRANSITION, LLVM-38217
614 const uint32_t __c
= __output2
- 10000 * (__output2
/ 10000);
616 const uint32_t __c
= __output2
% 10000;
619 const uint32_t __c0
= (__c
% 100) << 1;
620 const uint32_t __c1
= (__c
/ 100) << 1;
621 std::memcpy(__result
+ __olength
- __i
- 1, __DIGIT_TABLE
+ __c0
, 2);
622 std::memcpy(__result
+ __olength
- __i
- 3, __DIGIT_TABLE
+ __c1
, 2);
625 if (__output2
>= 100) {
626 const uint32_t __c
= (__output2
% 100) << 1;
628 std::memcpy(__result
+ __olength
- __i
- 1, __DIGIT_TABLE
+ __c
, 2);
631 if (__output2
>= 10) {
632 const uint32_t __c
= __output2
<< 1;
633 // We can't use memcpy here: the decimal dot goes between these two digits.
634 __result
[2] = __DIGIT_TABLE
[__c
+ 1];
635 __result
[0] = __DIGIT_TABLE
[__c
];
637 __result
[0] = static_cast<char>('0' + __output2
);
640 // Print decimal point if needed.
644 __index
= __olength
+ 1;
649 // Print the exponent.
650 __result
[__index
++] = 'e';
651 if (_Scientific_exponent
< 0) {
652 __result
[__index
++] = '-';
653 _Scientific_exponent
= -_Scientific_exponent
;
655 __result
[__index
++] = '+';
658 if (_Scientific_exponent
>= 100) {
659 const int32_t __c
= _Scientific_exponent
% 10;
660 std::memcpy(__result
+ __index
, __DIGIT_TABLE
+ 2 * (_Scientific_exponent
/ 10), 2);
661 __result
[__index
+ 2] = static_cast<char>('0' + __c
);
664 std::memcpy(__result
+ __index
, __DIGIT_TABLE
+ 2 * _Scientific_exponent
, 2);
668 return { _First
+ _Total_scientific_length
, errc
{} };
671 [[nodiscard
]] _LIBCPP_HIDE_FROM_ABI
inline bool __d2d_small_int(const uint64_t __ieeeMantissa
, const uint32_t __ieeeExponent
,
672 __floating_decimal_64
* const __v
) {
673 const uint64_t __m2
= (1ull << __DOUBLE_MANTISSA_BITS
) | __ieeeMantissa
;
674 const int32_t __e2
= static_cast<int32_t>(__ieeeExponent
) - __DOUBLE_BIAS
- __DOUBLE_MANTISSA_BITS
;
677 // f = __m2 * 2^__e2 >= 2^53 is an integer.
678 // Ignore this case for now.
687 // Since 2^52 <= __m2 < 2^53 and 0 <= -__e2 <= 52: 1 <= f = __m2 / 2^-__e2 < 2^53.
688 // Test if the lower -__e2 bits of the significand are 0, i.e. whether the fraction is 0.
689 const uint64_t __mask
= (1ull << -__e2
) - 1;
690 const uint64_t __fraction
= __m2
& __mask
;
691 if (__fraction
!= 0) {
695 // f is an integer in the range [1, 2^53).
696 // Note: __mantissa might contain trailing (decimal) 0's.
697 // Note: since 2^53 < 10^16, there is no need to adjust __decimalLength17().
698 __v
->__mantissa
= __m2
>> -__e2
;
703 [[nodiscard
]] to_chars_result
__d2s_buffered_n(char* const _First
, char* const _Last
, const double __f
,
704 const chars_format _Fmt
) {
706 // Step 1: Decode the floating-point number, and unify normalized and subnormal cases.
707 const uint64_t __bits
= __double_to_bits(__f
);
709 // Case distinction; exit early for the easy cases.
711 if (_Fmt
== chars_format::scientific
) {
712 if (_Last
- _First
< 5) {
713 return { _Last
, errc::value_too_large
};
716 std::memcpy(_First
, "0e+00", 5);
718 return { _First
+ 5, errc
{} };
721 // Print "0" for chars_format::fixed, chars_format::general, and chars_format{}.
722 if (_First
== _Last
) {
723 return { _Last
, errc::value_too_large
};
728 return { _First
+ 1, errc
{} };
731 // Decode __bits into mantissa and exponent.
732 const uint64_t __ieeeMantissa
= __bits
& ((1ull << __DOUBLE_MANTISSA_BITS
) - 1);
733 const uint32_t __ieeeExponent
= static_cast<uint32_t>(__bits
>> __DOUBLE_MANTISSA_BITS
);
735 if (_Fmt
== chars_format::fixed
) {
736 // const uint64_t _Mantissa2 = __ieeeMantissa | (1ull << __DOUBLE_MANTISSA_BITS); // restore implicit bit
737 const int32_t _Exponent2
= static_cast<int32_t>(__ieeeExponent
)
738 - __DOUBLE_BIAS
- __DOUBLE_MANTISSA_BITS
; // bias and normalization
740 // Normal values are equal to _Mantissa2 * 2^_Exponent2.
741 // (Subnormals are different, but they'll be rejected by the _Exponent2 test here, so they can be ignored.)
743 // For nonzero integers, _Exponent2 >= -52. (The minimum value occurs when _Mantissa2 * 2^_Exponent2 is 1.
744 // In that case, _Mantissa2 is the implicit 1 bit followed by 52 zeros, so _Exponent2 is -52 to shift away
745 // the zeros.) The dense range of exactly representable integers has negative or zero exponents
746 // (as positive exponents make the range non-dense). For that dense range, Ryu will always be used:
747 // every digit is necessary to uniquely identify the value, so Ryu must print them all.
749 // Positive exponents are the non-dense range of exactly representable integers. This contains all of the values
750 // for which Ryu can't be used (and a few Ryu-friendly values). We can save time by detecting positive
751 // exponents here and skipping Ryu. Calling __d2fixed_buffered_n() with precision 0 is valid for all integers
752 // (so it's okay if we call it with a Ryu-friendly value).
753 if (_Exponent2
> 0) {
754 return __d2fixed_buffered_n(_First
, _Last
, __f
, 0);
758 __floating_decimal_64 __v
;
759 const bool __isSmallInt
= __d2d_small_int(__ieeeMantissa
, __ieeeExponent
, &__v
);
761 // For small integers in the range [1, 2^53), __v.__mantissa might contain trailing (decimal) zeros.
762 // For scientific notation we need to move these zeros into the exponent.
763 // (This is not needed for fixed-point notation, so it might be beneficial to trim
764 // trailing zeros in __to_chars only if needed - once fixed-point notation output is implemented.)
766 const uint64_t __q
= __div10(__v
.__mantissa
);
767 const uint32_t __r
= static_cast<uint32_t>(__v
.__mantissa
) - 10 * static_cast<uint32_t>(__q
);
771 __v
.__mantissa
= __q
;
775 __v
= __d2d(__ieeeMantissa
, __ieeeExponent
);
778 return __to_chars(_First
, _Last
, __v
, _Fmt
, __f
);
781 _LIBCPP_END_NAMESPACE_STD