1 /* Copyright © 2018, Keith Packard
2 All rights reserved.
3
4 Redistribution and use in source and binary forms, with or without
5 modification, are permitted provided that the following conditions are met:
6
7 * Redistributions of source code must retain the above copyright
8 notice, this list of conditions and the following disclaimer.
9 * Redistributions in binary form must reproduce the above copyright
10 notice, this list of conditions and the following disclaimer in
11 the documentation and/or other materials provided with the
12 distribution.
13 * Neither the name of the copyright holders nor the names of
14 contributors may be used to endorse or promote products derived
15 from this software without specific prior written permission.
16
17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 POSSIBILITY OF SUCH DAMAGE. */
28
29 #include "stdio_private.h"
30
31 #ifdef FLOAT64
32
33 #define _NEED_IO_FLOAT64
34
35 #include "dtoa.h"
36
37 #define max(a, b) ({\
38 __typeof(a) _a = a;\
39 __typeof(b) _b = b;\
40 _a > _b ? _a : _b; })
41
42 #define min(a, b) ({\
43 __typeof(a) _a = a;\
44 __typeof(b) _b = b;\
45 _a < _b ? _a : _b; })
46
47 #define FRACTION_BITS (53 - 1)
48 #define FRACTION_MASK (((uint64_t) 1 << FRACTION_BITS) - 1)
49 #define EXPONENT_BITS (64 - FRACTION_BITS - 1)
50 #define EXPONENT_MASK (((uint64_t) 1 << EXPONENT_BITS) - 1)
51 #define SIGN_BIT ((uint64_t) 1 << (64 - 1))
52 #define BIT64(x) ((uint64_t) 1 << (x))
53
54 /*
55 * Check to see if the high bit is set by casting
56 * to signed and checking for negative
57 */
58 static inline bool
high_bit_set(uint64_t fract)59 high_bit_set(uint64_t fract)
60 {
61 return ((int64_t) fract < 0);
62 }
63
64 int
__dtoa_engine(FLOAT64 x,struct dtoa * dtoa,int max_digits,bool fmode,int max_decimals)65 __dtoa_engine(FLOAT64 x, struct dtoa *dtoa, int max_digits, bool fmode, int max_decimals)
66 {
67 uint64_t v = asuint64(x);
68 uint64_t fract = (v << (EXPONENT_BITS + 1)) >> 1;
69 int expo = (v << 1) >> (FRACTION_BITS + 1);
70 int decexp = 0;
71 int i;
72 uint8_t flags = 0;
73
74 if (high_bit_set(v))
75 flags |= DTOA_MINUS;
76 if (expo == EXPONENT_MASK) {
77 if (fract)
78 flags |= DTOA_NAN;
79 else
80 flags |= DTOA_INF;
81 } else {
82 if (expo == 0) {
83 if (fract == 0) {
84 flags |= DTOA_ZERO;
85 dtoa->digits[0] = '0';
86 max_digits = 1;
87 goto done;
88 }
89 /* normalize */
90 while ((int64_t) (fract <<= 1) >= 0)
91 expo--;
92 }
93
94 /* add implicit bit */
95 fract |= SIGN_BIT;
96 /* adjust exponent */
97 expo -= (DTOA_MAX_EXP - 2);
98 decexp = -1;
99
100 /*
101 * Let's consider:
102 *
103 * value = fract * 2^expo * 10^decexp
104 *
105 * Initially decexp = 0. The goal is to bring exp between
106 * 0 and -2 as the magnitude of a fractional decimal digit is 3 bits.
107 */
108
109 while (expo < -2) {
110 /*
111 * Make room to allow a multiplication by 5 without overflow.
112 * We test only the top part for faster code.
113 */
114 do {
115 /* Round this division to avoid accumulating errors */
116 fract = (fract >> 1) + (fract&1);
117 expo++;
118 } while ((uint32_t)(fract >> 32) >= (UINT32_MAX / 5U));
119
120 /* Perform fract * 5 * 2 / 10 */
121 fract *= 5U;
122 expo++;
123 decexp--;
124 }
125
126 while (expo > 0) {
127 /*
128 * Perform fract / 5 / 2 * 10.
129 * The +2 is there to do round the result of the division
130 * by 5 not to lose too much precision in extreme cases.
131 */
132 fract += 2;
133 fract /= 5U;
134 expo--;
135 decexp++;
136
137 /* Bring back our fractional number to full scale */
138 do {
139 fract <<= 1;
140 expo--;
141 } while (!high_bit_set(fract));
142 }
143
144
145 /*
146 * The binary fractional point is located somewhere above bit 63.
147 * Move it between bits 59 and 60 to give 4 bits of room to the
148 * integer part.
149 */
150 fract >>= (4 - expo);
151
152 /* If limiting decimals, then limit the max digits
153 * to no more than the number of digits left of the decimal
154 * plus the number of digits right of the decimal
155 *
156 * exp: exponent value. If negative, there are
157 * -exp - 1 zeros left of the first non-zero
158 * digit in 'f' format. If non-negative,
159 * there are exp digits to the left of
160 * the decimal point
161 *
162 * max_decimals: Only used in 'f' format. Round to no
163 * more than this many digits to the
164 * right of the decimal point (left if
165 * negative)
166 *
167 * max_digits: We can't convert more than this number
168 * of digits given the limits of the
169 * buffer
170 */
171
172 int save_max_digits = max_digits;
173 if (fmode) {
174 /*
175 * This covers two cases:
176 *
177 * When exp is < 0, there are -exp-1 zeros taking up
178 * space before we can display any of the real digits,
179 * so we have to subtract those off max_decimals before
180 * we round that (max_decimals - (-exp - 1)). This
181 * may end up less than zero, in which case we have
182 * no digits to display.
183 *
184 * When exp >= 0, there are exp + 1 digits left of the
185 * decimal point *plus* max_decimals right of the
186 * decimal point that need to be generated
187 *
188 * A single expression gives the right answer in both
189 * cases, which is kinda cool
190 */
191 /* max_decimals comes in biased by 1 to flag the 'f' case */
192 max_digits = min(max_digits, max(1, max_decimals + decexp + 1));
193 }
194
195 int decimals = max_digits;
196
197 /* Round the value to the last digit being printed. */
198 uint64_t round = BIT64(59); /* 0.5 */
199 while (decimals--) {
200 round /= 10U;
201 }
202 fract += round;
203 /* Make sure rounding didn't make fract >= 1.0 */
204 if (fract >= BIT64(60)) {
205 fract /= 10U;
206 decexp++;
207 max_digits = min(save_max_digits, max(1, max_decimals + decexp + 1));
208 }
209
210 /* Now convert mantissa to decimal. */
211
212 /* Compute digits */
213 for (i = 0; i < max_digits; i++) {
214 fract *= 10U;
215 dtoa->digits[i] = (fract >> 60) + '0';
216 fract &= BIT64(60) - 1;
217 }
218 }
219 done:
220 dtoa->flags = flags;
221 dtoa->exp = decexp;
222 return max_digits;
223 }
224
225 #endif
226
