1 | /*--------------------------------------------------------------------------- |
2 | * |
3 | * Ryu floating-point output for single precision. |
4 | * |
5 | * Portions Copyright (c) 2018-2019, PostgreSQL Global Development Group |
6 | * |
7 | * IDENTIFICATION |
8 | * src/common/f2s.c |
9 | * |
10 | * This is a modification of code taken from github.com/ulfjack/ryu under the |
11 | * terms of the Boost license (not the Apache license). The original copyright |
12 | * notice follows: |
13 | * |
14 | * Copyright 2018 Ulf Adams |
15 | * |
16 | * The contents of this file may be used under the terms of the Apache |
17 | * License, Version 2.0. |
18 | * |
19 | * (See accompanying file LICENSE-Apache or copy at |
20 | * http://www.apache.org/licenses/LICENSE-2.0) |
21 | * |
22 | * Alternatively, the contents of this file may be used under the terms of the |
23 | * Boost Software License, Version 1.0. |
24 | * |
25 | * (See accompanying file LICENSE-Boost or copy at |
26 | * https://www.boost.org/LICENSE_1_0.txt) |
27 | * |
28 | * Unless required by applicable law or agreed to in writing, this software is |
29 | * distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY |
30 | * KIND, either express or implied. |
31 | * |
32 | *--------------------------------------------------------------------------- |
33 | */ |
34 | |
35 | #ifndef FRONTEND |
36 | #include "postgres.h" |
37 | #else |
38 | #include "postgres_fe.h" |
39 | #endif |
40 | |
41 | #include "common/shortest_dec.h" |
42 | |
43 | #include "ryu_common.h" |
44 | #include "digit_table.h" |
45 | |
46 | #define FLOAT_MANTISSA_BITS 23 |
47 | #define FLOAT_EXPONENT_BITS 8 |
48 | #define FLOAT_BIAS 127 |
49 | |
50 | /* |
51 | * This table is generated (by the upstream) by PrintFloatLookupTable, |
52 | * and modified (by us) to add UINT64CONST. |
53 | */ |
54 | #define FLOAT_POW5_INV_BITCOUNT 59 |
55 | static const uint64 FLOAT_POW5_INV_SPLIT[31] = { |
56 | UINT64CONST(576460752303423489), UINT64CONST(461168601842738791), UINT64CONST(368934881474191033), UINT64CONST(295147905179352826), |
57 | UINT64CONST(472236648286964522), UINT64CONST(377789318629571618), UINT64CONST(302231454903657294), UINT64CONST(483570327845851670), |
58 | UINT64CONST(386856262276681336), UINT64CONST(309485009821345069), UINT64CONST(495176015714152110), UINT64CONST(396140812571321688), |
59 | UINT64CONST(316912650057057351), UINT64CONST(507060240091291761), UINT64CONST(405648192073033409), UINT64CONST(324518553658426727), |
60 | UINT64CONST(519229685853482763), UINT64CONST(415383748682786211), UINT64CONST(332306998946228969), UINT64CONST(531691198313966350), |
61 | UINT64CONST(425352958651173080), UINT64CONST(340282366920938464), UINT64CONST(544451787073501542), UINT64CONST(435561429658801234), |
62 | UINT64CONST(348449143727040987), UINT64CONST(557518629963265579), UINT64CONST(446014903970612463), UINT64CONST(356811923176489971), |
63 | UINT64CONST(570899077082383953), UINT64CONST(456719261665907162), UINT64CONST(365375409332725730) |
64 | }; |
65 | #define FLOAT_POW5_BITCOUNT 61 |
66 | static const uint64 FLOAT_POW5_SPLIT[47] = { |
67 | UINT64CONST(1152921504606846976), UINT64CONST(1441151880758558720), UINT64CONST(1801439850948198400), UINT64CONST(2251799813685248000), |
68 | UINT64CONST(1407374883553280000), UINT64CONST(1759218604441600000), UINT64CONST(2199023255552000000), UINT64CONST(1374389534720000000), |
69 | UINT64CONST(1717986918400000000), UINT64CONST(2147483648000000000), UINT64CONST(1342177280000000000), UINT64CONST(1677721600000000000), |
70 | UINT64CONST(2097152000000000000), UINT64CONST(1310720000000000000), UINT64CONST(1638400000000000000), UINT64CONST(2048000000000000000), |
71 | UINT64CONST(1280000000000000000), UINT64CONST(1600000000000000000), UINT64CONST(2000000000000000000), UINT64CONST(1250000000000000000), |
72 | UINT64CONST(1562500000000000000), UINT64CONST(1953125000000000000), UINT64CONST(1220703125000000000), UINT64CONST(1525878906250000000), |
73 | UINT64CONST(1907348632812500000), UINT64CONST(1192092895507812500), UINT64CONST(1490116119384765625), UINT64CONST(1862645149230957031), |
74 | UINT64CONST(1164153218269348144), UINT64CONST(1455191522836685180), UINT64CONST(1818989403545856475), UINT64CONST(2273736754432320594), |
75 | UINT64CONST(1421085471520200371), UINT64CONST(1776356839400250464), UINT64CONST(2220446049250313080), UINT64CONST(1387778780781445675), |
76 | UINT64CONST(1734723475976807094), UINT64CONST(2168404344971008868), UINT64CONST(1355252715606880542), UINT64CONST(1694065894508600678), |
77 | UINT64CONST(2117582368135750847), UINT64CONST(1323488980084844279), UINT64CONST(1654361225106055349), UINT64CONST(2067951531382569187), |
78 | UINT64CONST(1292469707114105741), UINT64CONST(1615587133892632177), UINT64CONST(2019483917365790221) |
79 | }; |
80 | |
81 | static inline uint32 |
82 | pow5Factor(uint32 value) |
83 | { |
84 | uint32 count = 0; |
85 | |
86 | for (;;) |
87 | { |
88 | Assert(value != 0); |
89 | const uint32 q = value / 5; |
90 | const uint32 r = value % 5; |
91 | |
92 | if (r != 0) |
93 | break; |
94 | |
95 | value = q; |
96 | ++count; |
97 | } |
98 | return count; |
99 | } |
100 | |
101 | /* Returns true if value is divisible by 5^p. */ |
102 | static inline bool |
103 | multipleOfPowerOf5(const uint32 value, const uint32 p) |
104 | { |
105 | return pow5Factor(value) >= p; |
106 | } |
107 | |
108 | /* Returns true if value is divisible by 2^p. */ |
109 | static inline bool |
110 | multipleOfPowerOf2(const uint32 value, const uint32 p) |
111 | { |
112 | /* return __builtin_ctz(value) >= p; */ |
113 | return (value & ((1u << p) - 1)) == 0; |
114 | } |
115 | |
116 | /* |
117 | * It seems to be slightly faster to avoid uint128_t here, although the |
118 | * generated code for uint128_t looks slightly nicer. |
119 | */ |
120 | static inline uint32 |
121 | mulShift(const uint32 m, const uint64 factor, const int32 shift) |
122 | { |
123 | /* |
124 | * The casts here help MSVC to avoid calls to the __allmul library |
125 | * function. |
126 | */ |
127 | const uint32 factorLo = (uint32) (factor); |
128 | const uint32 factorHi = (uint32) (factor >> 32); |
129 | const uint64 bits0 = (uint64) m * factorLo; |
130 | const uint64 bits1 = (uint64) m * factorHi; |
131 | |
132 | Assert(shift > 32); |
133 | |
134 | #ifdef RYU_32_BIT_PLATFORM |
135 | |
136 | /* |
137 | * On 32-bit platforms we can avoid a 64-bit shift-right since we only |
138 | * need the upper 32 bits of the result and the shift value is > 32. |
139 | */ |
140 | const uint32 bits0Hi = (uint32) (bits0 >> 32); |
141 | uint32 bits1Lo = (uint32) (bits1); |
142 | uint32 bits1Hi = (uint32) (bits1 >> 32); |
143 | |
144 | bits1Lo += bits0Hi; |
145 | bits1Hi += (bits1Lo < bits0Hi); |
146 | |
147 | const int32 s = shift - 32; |
148 | |
149 | return (bits1Hi << (32 - s)) | (bits1Lo >> s); |
150 | |
151 | #else /* RYU_32_BIT_PLATFORM */ |
152 | |
153 | const uint64 sum = (bits0 >> 32) + bits1; |
154 | const uint64 shiftedSum = sum >> (shift - 32); |
155 | |
156 | Assert(shiftedSum <= PG_UINT32_MAX); |
157 | return (uint32) shiftedSum; |
158 | |
159 | #endif /* RYU_32_BIT_PLATFORM */ |
160 | } |
161 | |
162 | static inline uint32 |
163 | mulPow5InvDivPow2(const uint32 m, const uint32 q, const int32 j) |
164 | { |
165 | return mulShift(m, FLOAT_POW5_INV_SPLIT[q], j); |
166 | } |
167 | |
168 | static inline uint32 |
169 | mulPow5divPow2(const uint32 m, const uint32 i, const int32 j) |
170 | { |
171 | return mulShift(m, FLOAT_POW5_SPLIT[i], j); |
172 | } |
173 | |
174 | static inline uint32 |
175 | decimalLength(const uint32 v) |
176 | { |
177 | /* Function precondition: v is not a 10-digit number. */ |
178 | /* (9 digits are sufficient for round-tripping.) */ |
179 | Assert(v < 1000000000); |
180 | if (v >= 100000000) |
181 | { |
182 | return 9; |
183 | } |
184 | if (v >= 10000000) |
185 | { |
186 | return 8; |
187 | } |
188 | if (v >= 1000000) |
189 | { |
190 | return 7; |
191 | } |
192 | if (v >= 100000) |
193 | { |
194 | return 6; |
195 | } |
196 | if (v >= 10000) |
197 | { |
198 | return 5; |
199 | } |
200 | if (v >= 1000) |
201 | { |
202 | return 4; |
203 | } |
204 | if (v >= 100) |
205 | { |
206 | return 3; |
207 | } |
208 | if (v >= 10) |
209 | { |
210 | return 2; |
211 | } |
212 | return 1; |
213 | } |
214 | |
215 | /* A floating decimal representing m * 10^e. */ |
216 | typedef struct floating_decimal_32 |
217 | { |
218 | uint32 mantissa; |
219 | int32 exponent; |
220 | } floating_decimal_32; |
221 | |
222 | static inline floating_decimal_32 |
223 | f2d(const uint32 ieeeMantissa, const uint32 ieeeExponent) |
224 | { |
225 | int32 e2; |
226 | uint32 m2; |
227 | |
228 | if (ieeeExponent == 0) |
229 | { |
230 | /* We subtract 2 so that the bounds computation has 2 additional bits. */ |
231 | e2 = 1 - FLOAT_BIAS - FLOAT_MANTISSA_BITS - 2; |
232 | m2 = ieeeMantissa; |
233 | } |
234 | else |
235 | { |
236 | e2 = ieeeExponent - FLOAT_BIAS - FLOAT_MANTISSA_BITS - 2; |
237 | m2 = (1u << FLOAT_MANTISSA_BITS) | ieeeMantissa; |
238 | } |
239 | |
240 | #if STRICTLY_SHORTEST |
241 | const bool even = (m2 & 1) == 0; |
242 | const bool acceptBounds = even; |
243 | #else |
244 | const bool acceptBounds = false; |
245 | #endif |
246 | |
247 | /* Step 2: Determine the interval of legal decimal representations. */ |
248 | const uint32 mv = 4 * m2; |
249 | const uint32 mp = 4 * m2 + 2; |
250 | |
251 | /* Implicit bool -> int conversion. True is 1, false is 0. */ |
252 | const uint32 mmShift = ieeeMantissa != 0 || ieeeExponent <= 1; |
253 | const uint32 mm = 4 * m2 - 1 - mmShift; |
254 | |
255 | /* Step 3: Convert to a decimal power base using 64-bit arithmetic. */ |
256 | uint32 vr, |
257 | vp, |
258 | vm; |
259 | int32 e10; |
260 | bool vmIsTrailingZeros = false; |
261 | bool vrIsTrailingZeros = false; |
262 | uint8 lastRemovedDigit = 0; |
263 | |
264 | if (e2 >= 0) |
265 | { |
266 | const uint32 q = log10Pow2(e2); |
267 | |
268 | e10 = q; |
269 | |
270 | const int32 k = FLOAT_POW5_INV_BITCOUNT + pow5bits(q) - 1; |
271 | const int32 i = -e2 + q + k; |
272 | |
273 | vr = mulPow5InvDivPow2(mv, q, i); |
274 | vp = mulPow5InvDivPow2(mp, q, i); |
275 | vm = mulPow5InvDivPow2(mm, q, i); |
276 | |
277 | if (q != 0 && (vp - 1) / 10 <= vm / 10) |
278 | { |
279 | /* |
280 | * We need to know one removed digit even if we are not going to |
281 | * loop below. We could use q = X - 1 above, except that would |
282 | * require 33 bits for the result, and we've found that 32-bit |
283 | * arithmetic is faster even on 64-bit machines. |
284 | */ |
285 | const int32 l = FLOAT_POW5_INV_BITCOUNT + pow5bits(q - 1) - 1; |
286 | |
287 | lastRemovedDigit = (uint8) (mulPow5InvDivPow2(mv, q - 1, -e2 + q - 1 + l) % 10); |
288 | } |
289 | if (q <= 9) |
290 | { |
291 | /* |
292 | * The largest power of 5 that fits in 24 bits is 5^10, but q <= 9 |
293 | * seems to be safe as well. |
294 | * |
295 | * Only one of mp, mv, and mm can be a multiple of 5, if any. |
296 | */ |
297 | if (mv % 5 == 0) |
298 | { |
299 | vrIsTrailingZeros = multipleOfPowerOf5(mv, q); |
300 | } |
301 | else if (acceptBounds) |
302 | { |
303 | vmIsTrailingZeros = multipleOfPowerOf5(mm, q); |
304 | } |
305 | else |
306 | { |
307 | vp -= multipleOfPowerOf5(mp, q); |
308 | } |
309 | } |
310 | } |
311 | else |
312 | { |
313 | const uint32 q = log10Pow5(-e2); |
314 | |
315 | e10 = q + e2; |
316 | |
317 | const int32 i = -e2 - q; |
318 | const int32 k = pow5bits(i) - FLOAT_POW5_BITCOUNT; |
319 | int32 j = q - k; |
320 | |
321 | vr = mulPow5divPow2(mv, i, j); |
322 | vp = mulPow5divPow2(mp, i, j); |
323 | vm = mulPow5divPow2(mm, i, j); |
324 | |
325 | if (q != 0 && (vp - 1) / 10 <= vm / 10) |
326 | { |
327 | j = q - 1 - (pow5bits(i + 1) - FLOAT_POW5_BITCOUNT); |
328 | lastRemovedDigit = (uint8) (mulPow5divPow2(mv, i + 1, j) % 10); |
329 | } |
330 | if (q <= 1) |
331 | { |
332 | /* |
333 | * {vr,vp,vm} is trailing zeros if {mv,mp,mm} has at least q |
334 | * trailing 0 bits. |
335 | */ |
336 | /* mv = 4 * m2, so it always has at least two trailing 0 bits. */ |
337 | vrIsTrailingZeros = true; |
338 | if (acceptBounds) |
339 | { |
340 | /* |
341 | * mm = mv - 1 - mmShift, so it has 1 trailing 0 bit iff |
342 | * mmShift == 1. |
343 | */ |
344 | vmIsTrailingZeros = mmShift == 1; |
345 | } |
346 | else |
347 | { |
348 | /* |
349 | * mp = mv + 2, so it always has at least one trailing 0 bit. |
350 | */ |
351 | --vp; |
352 | } |
353 | } |
354 | else if (q < 31) |
355 | { |
356 | /* TODO(ulfjack):Use a tighter bound here. */ |
357 | vrIsTrailingZeros = multipleOfPowerOf2(mv, q - 1); |
358 | } |
359 | } |
360 | |
361 | /* |
362 | * Step 4: Find the shortest decimal representation in the interval of |
363 | * legal representations. |
364 | */ |
365 | uint32 removed = 0; |
366 | uint32 output; |
367 | |
368 | if (vmIsTrailingZeros || vrIsTrailingZeros) |
369 | { |
370 | /* General case, which happens rarely (~4.0%). */ |
371 | while (vp / 10 > vm / 10) |
372 | { |
373 | vmIsTrailingZeros &= vm - (vm / 10) * 10 == 0; |
374 | vrIsTrailingZeros &= lastRemovedDigit == 0; |
375 | lastRemovedDigit = (uint8) (vr % 10); |
376 | vr /= 10; |
377 | vp /= 10; |
378 | vm /= 10; |
379 | ++removed; |
380 | } |
381 | if (vmIsTrailingZeros) |
382 | { |
383 | while (vm % 10 == 0) |
384 | { |
385 | vrIsTrailingZeros &= lastRemovedDigit == 0; |
386 | lastRemovedDigit = (uint8) (vr % 10); |
387 | vr /= 10; |
388 | vp /= 10; |
389 | vm /= 10; |
390 | ++removed; |
391 | } |
392 | } |
393 | |
394 | if (vrIsTrailingZeros && lastRemovedDigit == 5 && vr % 2 == 0) |
395 | { |
396 | /* Round even if the exact number is .....50..0. */ |
397 | lastRemovedDigit = 4; |
398 | } |
399 | |
400 | /* |
401 | * We need to take vr + 1 if vr is outside bounds or we need to round |
402 | * up. |
403 | */ |
404 | output = vr + ((vr == vm && (!acceptBounds || !vmIsTrailingZeros)) || lastRemovedDigit >= 5); |
405 | } |
406 | else |
407 | { |
408 | /* |
409 | * Specialized for the common case (~96.0%). Percentages below are |
410 | * relative to this. |
411 | * |
412 | * Loop iterations below (approximately): 0: 13.6%, 1: 70.7%, 2: |
413 | * 14.1%, 3: 1.39%, 4: 0.14%, 5+: 0.01% |
414 | */ |
415 | while (vp / 10 > vm / 10) |
416 | { |
417 | lastRemovedDigit = (uint8) (vr % 10); |
418 | vr /= 10; |
419 | vp /= 10; |
420 | vm /= 10; |
421 | ++removed; |
422 | } |
423 | |
424 | /* |
425 | * We need to take vr + 1 if vr is outside bounds or we need to round |
426 | * up. |
427 | */ |
428 | output = vr + (vr == vm || lastRemovedDigit >= 5); |
429 | } |
430 | |
431 | const int32 exp = e10 + removed; |
432 | |
433 | floating_decimal_32 fd; |
434 | |
435 | fd.exponent = exp; |
436 | fd.mantissa = output; |
437 | return fd; |
438 | } |
439 | |
440 | static inline int |
441 | to_chars_f(const floating_decimal_32 v, const uint32 olength, char *const result) |
442 | { |
443 | /* Step 5: Print the decimal representation. */ |
444 | int index = 0; |
445 | |
446 | uint32 output = v.mantissa; |
447 | int32 exp = v.exponent; |
448 | |
449 | /*---- |
450 | * On entry, mantissa * 10^exp is the result to be output. |
451 | * Caller has already done the - sign if needed. |
452 | * |
453 | * We want to insert the point somewhere depending on the output length |
454 | * and exponent, which might mean adding zeros: |
455 | * |
456 | * exp | format |
457 | * 1+ | ddddddddd000000 |
458 | * 0 | ddddddddd |
459 | * -1 .. -len+1 | dddddddd.d to d.ddddddddd |
460 | * -len ... | 0.ddddddddd to 0.000dddddd |
461 | */ |
462 | uint32 i = 0; |
463 | int32 nexp = exp + olength; |
464 | |
465 | if (nexp <= 0) |
466 | { |
467 | /* -nexp is number of 0s to add after '.' */ |
468 | Assert(nexp >= -3); |
469 | /* 0.000ddddd */ |
470 | index = 2 - nexp; |
471 | /* copy 8 bytes rather than 5 to let compiler optimize */ |
472 | memcpy(result, "0.000000" , 8); |
473 | } |
474 | else if (exp < 0) |
475 | { |
476 | /* |
477 | * dddd.dddd; leave space at the start and move the '.' in after |
478 | */ |
479 | index = 1; |
480 | } |
481 | else |
482 | { |
483 | /* |
484 | * We can save some code later by pre-filling with zeros. We know that |
485 | * there can be no more than 6 output digits in this form, otherwise |
486 | * we would not choose fixed-point output. memset 8 rather than 6 |
487 | * bytes to let the compiler optimize it. |
488 | */ |
489 | Assert(exp < 6 && exp + olength <= 6); |
490 | memset(result, '0', 8); |
491 | } |
492 | |
493 | while (output >= 10000) |
494 | { |
495 | const uint32 c = output - 10000 * (output / 10000); |
496 | const uint32 c0 = (c % 100) << 1; |
497 | const uint32 c1 = (c / 100) << 1; |
498 | |
499 | output /= 10000; |
500 | |
501 | memcpy(result + index + olength - i - 2, DIGIT_TABLE + c0, 2); |
502 | memcpy(result + index + olength - i - 4, DIGIT_TABLE + c1, 2); |
503 | i += 4; |
504 | } |
505 | if (output >= 100) |
506 | { |
507 | const uint32 c = (output % 100) << 1; |
508 | |
509 | output /= 100; |
510 | memcpy(result + index + olength - i - 2, DIGIT_TABLE + c, 2); |
511 | i += 2; |
512 | } |
513 | if (output >= 10) |
514 | { |
515 | const uint32 c = output << 1; |
516 | |
517 | memcpy(result + index + olength - i - 2, DIGIT_TABLE + c, 2); |
518 | } |
519 | else |
520 | { |
521 | result[index] = (char) ('0' + output); |
522 | } |
523 | |
524 | if (index == 1) |
525 | { |
526 | /* |
527 | * nexp is 1..6 here, representing the number of digits before the |
528 | * point. A value of 7+ is not possible because we switch to |
529 | * scientific notation when the display exponent reaches 6. |
530 | */ |
531 | Assert(nexp < 7); |
532 | /* gcc only seems to want to optimize memmove for small 2^n */ |
533 | if (nexp & 4) |
534 | { |
535 | memmove(result + index - 1, result + index, 4); |
536 | index += 4; |
537 | } |
538 | if (nexp & 2) |
539 | { |
540 | memmove(result + index - 1, result + index, 2); |
541 | index += 2; |
542 | } |
543 | if (nexp & 1) |
544 | { |
545 | result[index - 1] = result[index]; |
546 | } |
547 | result[nexp] = '.'; |
548 | index = olength + 1; |
549 | } |
550 | else if (exp >= 0) |
551 | { |
552 | /* we supplied the trailing zeros earlier, now just set the length. */ |
553 | index = olength + exp; |
554 | } |
555 | else |
556 | { |
557 | index = olength + (2 - nexp); |
558 | } |
559 | |
560 | return index; |
561 | } |
562 | |
563 | static inline int |
564 | to_chars(const floating_decimal_32 v, const bool sign, char *const result) |
565 | { |
566 | /* Step 5: Print the decimal representation. */ |
567 | int index = 0; |
568 | |
569 | uint32 output = v.mantissa; |
570 | uint32 olength = decimalLength(output); |
571 | int32 exp = v.exponent + olength - 1; |
572 | |
573 | if (sign) |
574 | result[index++] = '-'; |
575 | |
576 | /* |
577 | * The thresholds for fixed-point output are chosen to match printf |
578 | * defaults. Beware that both the code of to_chars_f and the value of |
579 | * FLOAT_SHORTEST_DECIMAL_LEN are sensitive to these thresholds. |
580 | */ |
581 | if (exp >= -4 && exp < 6) |
582 | return to_chars_f(v, olength, result + index) + sign; |
583 | |
584 | /* |
585 | * If v.exponent is exactly 0, we might have reached here via the small |
586 | * integer fast path, in which case v.mantissa might contain trailing |
587 | * (decimal) zeros. For scientific notation we need to move these zeros |
588 | * into the exponent. (For fixed point this doesn't matter, which is why |
589 | * we do this here rather than above.) |
590 | * |
591 | * Since we already calculated the display exponent (exp) above based on |
592 | * the old decimal length, that value does not change here. Instead, we |
593 | * just reduce the display length for each digit removed. |
594 | * |
595 | * If we didn't get here via the fast path, the raw exponent will not |
596 | * usually be 0, and there will be no trailing zeros, so we pay no more |
597 | * than one div10/multiply extra cost. We claw back half of that by |
598 | * checking for divisibility by 2 before dividing by 10. |
599 | */ |
600 | if (v.exponent == 0) |
601 | { |
602 | while ((output & 1) == 0) |
603 | { |
604 | const uint32 q = output / 10; |
605 | const uint32 r = output - 10 * q; |
606 | |
607 | if (r != 0) |
608 | break; |
609 | output = q; |
610 | --olength; |
611 | } |
612 | } |
613 | |
614 | /*---- |
615 | * Print the decimal digits. |
616 | * The following code is equivalent to: |
617 | * |
618 | * for (uint32 i = 0; i < olength - 1; ++i) { |
619 | * const uint32 c = output % 10; output /= 10; |
620 | * result[index + olength - i] = (char) ('0' + c); |
621 | * } |
622 | * result[index] = '0' + output % 10; |
623 | */ |
624 | uint32 i = 0; |
625 | |
626 | while (output >= 10000) |
627 | { |
628 | const uint32 c = output - 10000 * (output / 10000); |
629 | const uint32 c0 = (c % 100) << 1; |
630 | const uint32 c1 = (c / 100) << 1; |
631 | |
632 | output /= 10000; |
633 | |
634 | memcpy(result + index + olength - i - 1, DIGIT_TABLE + c0, 2); |
635 | memcpy(result + index + olength - i - 3, DIGIT_TABLE + c1, 2); |
636 | i += 4; |
637 | } |
638 | if (output >= 100) |
639 | { |
640 | const uint32 c = (output % 100) << 1; |
641 | |
642 | output /= 100; |
643 | memcpy(result + index + olength - i - 1, DIGIT_TABLE + c, 2); |
644 | i += 2; |
645 | } |
646 | if (output >= 10) |
647 | { |
648 | const uint32 c = output << 1; |
649 | |
650 | /* |
651 | * We can't use memcpy here: the decimal dot goes between these two |
652 | * digits. |
653 | */ |
654 | result[index + olength - i] = DIGIT_TABLE[c + 1]; |
655 | result[index] = DIGIT_TABLE[c]; |
656 | } |
657 | else |
658 | { |
659 | result[index] = (char) ('0' + output); |
660 | } |
661 | |
662 | /* Print decimal point if needed. */ |
663 | if (olength > 1) |
664 | { |
665 | result[index + 1] = '.'; |
666 | index += olength + 1; |
667 | } |
668 | else |
669 | { |
670 | ++index; |
671 | } |
672 | |
673 | /* Print the exponent. */ |
674 | result[index++] = 'e'; |
675 | if (exp < 0) |
676 | { |
677 | result[index++] = '-'; |
678 | exp = -exp; |
679 | } |
680 | else |
681 | result[index++] = '+'; |
682 | |
683 | memcpy(result + index, DIGIT_TABLE + 2 * exp, 2); |
684 | index += 2; |
685 | |
686 | return index; |
687 | } |
688 | |
689 | static inline bool |
690 | f2d_small_int(const uint32 ieeeMantissa, |
691 | const uint32 ieeeExponent, |
692 | floating_decimal_32 *v) |
693 | { |
694 | const int32 e2 = (int32) ieeeExponent - FLOAT_BIAS - FLOAT_MANTISSA_BITS; |
695 | |
696 | /* |
697 | * Avoid using multiple "return false;" here since it tends to provoke the |
698 | * compiler into inlining multiple copies of f2d, which is undesirable. |
699 | */ |
700 | |
701 | if (e2 >= -FLOAT_MANTISSA_BITS && e2 <= 0) |
702 | { |
703 | /*---- |
704 | * Since 2^23 <= m2 < 2^24 and 0 <= -e2 <= 23: |
705 | * 1 <= f = m2 / 2^-e2 < 2^24. |
706 | * |
707 | * Test if the lower -e2 bits of the significand are 0, i.e. whether |
708 | * the fraction is 0. We can use ieeeMantissa here, since the implied |
709 | * 1 bit can never be tested by this; the implied 1 can only be part |
710 | * of a fraction if e2 < -FLOAT_MANTISSA_BITS which we already |
711 | * checked. (e.g. 0.5 gives ieeeMantissa == 0 and e2 == -24) |
712 | */ |
713 | const uint32 mask = (1U << -e2) - 1; |
714 | const uint32 fraction = ieeeMantissa & mask; |
715 | |
716 | if (fraction == 0) |
717 | { |
718 | /*---- |
719 | * f is an integer in the range [1, 2^24). |
720 | * Note: mantissa might contain trailing (decimal) 0's. |
721 | * Note: since 2^24 < 10^9, there is no need to adjust |
722 | * decimalLength(). |
723 | */ |
724 | const uint32 m2 = (1U << FLOAT_MANTISSA_BITS) | ieeeMantissa; |
725 | |
726 | v->mantissa = m2 >> -e2; |
727 | v->exponent = 0; |
728 | return true; |
729 | } |
730 | } |
731 | |
732 | return false; |
733 | } |
734 | |
735 | /* |
736 | * Store the shortest decimal representation of the given float as an |
737 | * UNTERMINATED string in the caller's supplied buffer (which must be at least |
738 | * FLOAT_SHORTEST_DECIMAL_LEN-1 bytes long). |
739 | * |
740 | * Returns the number of bytes stored. |
741 | */ |
742 | int |
743 | float_to_shortest_decimal_bufn(float f, char *result) |
744 | { |
745 | /* |
746 | * Step 1: Decode the floating-point number, and unify normalized and |
747 | * subnormal cases. |
748 | */ |
749 | const uint32 bits = float_to_bits(f); |
750 | |
751 | /* Decode bits into sign, mantissa, and exponent. */ |
752 | const bool ieeeSign = ((bits >> (FLOAT_MANTISSA_BITS + FLOAT_EXPONENT_BITS)) & 1) != 0; |
753 | const uint32 ieeeMantissa = bits & ((1u << FLOAT_MANTISSA_BITS) - 1); |
754 | const uint32 ieeeExponent = (bits >> FLOAT_MANTISSA_BITS) & ((1u << FLOAT_EXPONENT_BITS) - 1); |
755 | |
756 | /* Case distinction; exit early for the easy cases. */ |
757 | if (ieeeExponent == ((1u << FLOAT_EXPONENT_BITS) - 1u) || (ieeeExponent == 0 && ieeeMantissa == 0)) |
758 | { |
759 | return copy_special_str(result, ieeeSign, (ieeeExponent != 0), (ieeeMantissa != 0)); |
760 | } |
761 | |
762 | floating_decimal_32 v; |
763 | const bool isSmallInt = f2d_small_int(ieeeMantissa, ieeeExponent, &v); |
764 | |
765 | if (!isSmallInt) |
766 | { |
767 | v = f2d(ieeeMantissa, ieeeExponent); |
768 | } |
769 | |
770 | return to_chars(v, ieeeSign, result); |
771 | } |
772 | |
773 | /* |
774 | * Store the shortest decimal representation of the given float as a |
775 | * null-terminated string in the caller's supplied buffer (which must be at |
776 | * least FLOAT_SHORTEST_DECIMAL_LEN bytes long). |
777 | * |
778 | * Returns the string length. |
779 | */ |
780 | int |
781 | float_to_shortest_decimal_buf(float f, char *result) |
782 | { |
783 | const int index = float_to_shortest_decimal_bufn(f, result); |
784 | |
785 | /* Terminate the string. */ |
786 | Assert(index < FLOAT_SHORTEST_DECIMAL_LEN); |
787 | result[index] = '\0'; |
788 | return index; |
789 | } |
790 | |
791 | /* |
792 | * Return the shortest decimal representation as a null-terminated palloc'd |
793 | * string (outside the backend, uses malloc() instead). |
794 | * |
795 | * Caller is responsible for freeing the result. |
796 | */ |
797 | char * |
798 | float_to_shortest_decimal(float f) |
799 | { |
800 | char *const result = (char *) palloc(FLOAT_SHORTEST_DECIMAL_LEN); |
801 | |
802 | float_to_shortest_decimal_buf(f, result); |
803 | return result; |
804 | } |
805 | |