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 "dtoa_engine.h"
30 #include <math.h>
31 #include "../../libm/common/math_config.h"
32
33 #define max(a, b) ({\
34 __typeof(a) _a = a;\
35 __typeof(b) _b = b;\
36 _a > _b ? _a : _b; })
37
38 #define min(a, b) ({\
39 __typeof(a) _a = a;\
40 __typeof(b) _b = b;\
41 _a < _b ? _a : _b; })
42
43 #define FRACTION_BITS (__DBL_MANT_DIG__ - 1)
44 #define FRACTION_MASK (((uint64_t) 1 << FRACTION_BITS) - 1)
45 #define EXPONENT_BITS (sizeof(double) * 8 - FRACTION_BITS - 1)
46 #define EXPONENT_MASK (((uint64_t) 1 << EXPONENT_BITS) - 1)
47 #define SIGN_BIT ((uint64_t) 1 << (sizeof(double) * 8 - 1))
48 #define BIT64(x) ((uint64_t) 1 << (x))
49
50 /*
51 * Check to see if the high bit is set by casting
52 * to signed and checking for negative
53 */
54 static inline bool
high_bit_set(uint64_t fract)55 high_bit_set(uint64_t fract)
56 {
57 return ((int64_t) fract < 0);
58 }
59
60 int
__dtoa_engine(double x,struct dtoa * dtoa,int max_digits,bool fmode,int max_decimals)61 __dtoa_engine(double x, struct dtoa *dtoa, int max_digits, bool fmode, int max_decimals)
62 {
63 uint64_t v = asuint64(x);
64 uint64_t fract = (v << (EXPONENT_BITS + 1)) >> 1;
65 int expo = (v << 1) >> (FRACTION_BITS + 1);
66 int decexp = 0;
67 int i;
68 uint8_t flags = 0;
69
70 if (high_bit_set(v))
71 flags |= DTOA_MINUS;
72 if (expo == EXPONENT_MASK) {
73 if (fract)
74 flags |= DTOA_NAN;
75 else
76 flags |= DTOA_INF;
77 } else {
78 if (expo == 0) {
79 if (fract == 0) {
80 flags |= DTOA_ZERO;
81 dtoa->digits[0] = '0';
82 max_digits = 1;
83 goto done;
84 }
85 /* normalize */
86 while ((int64_t) (fract <<= 1) >= 0)
87 expo--;
88 }
89
90 /* add implicit bit */
91 fract |= SIGN_BIT;
92 /* adjust exponent */
93 expo -= (__DBL_MAX_EXP__ - 2);
94 decexp = -1;
95
96 /*
97 * Let's consider:
98 *
99 * value = fract * 2^expo * 10^decexp
100 *
101 * Initially decexp = 0. The goal is to bring exp between
102 * 0 and -2 as the magnitude of a fractional decimal digit is 3 bits.
103 */
104
105 while (expo < -2) {
106 /*
107 * Make room to allow a multiplication by 5 without overflow.
108 * We test only the top part for faster code.
109 */
110 do {
111 /* Round this division to avoid accumulating errors */
112 fract = (fract >> 1) + (fract&1);
113 expo++;
114 } while ((uint32_t)(fract >> 32) >= (UINT32_MAX / 5U));
115
116 /* Perform fract * 5 * 2 / 10 */
117 fract *= 5U;
118 expo++;
119 decexp--;
120 }
121
122 while (expo > 0) {
123 /*
124 * Perform fract / 5 / 2 * 10.
125 * The +2 is there to do round the result of the division
126 * by 5 not to lose too much precision in extreme cases.
127 */
128 fract += 2;
129 fract /= 5U;
130 expo--;
131 decexp++;
132
133 /* Bring back our fractional number to full scale */
134 do {
135 fract <<= 1;
136 expo--;
137 } while (!high_bit_set(fract));
138 }
139
140
141 /*
142 * The binary fractional point is located somewhere above bit 63.
143 * Move it between bits 59 and 60 to give 4 bits of room to the
144 * integer part.
145 */
146 fract >>= (4 - expo);
147
148 /* If limiting decimals, then limit the max digits
149 * to no more than the number of digits left of the decimal
150 * plus the number of digits right of the decimal
151 *
152 * exp: exponent value. If negative, there are
153 * -exp - 1 zeros left of the first non-zero
154 * digit in 'f' format. If non-negative,
155 * there are exp digits to the left of
156 * the decimal point
157 *
158 * max_decimals: Only used in 'f' format. Round to no
159 * more than this many digits to the
160 * right of the decimal point (left if
161 * negative)
162 *
163 * max_digits: We can't convert more than this number
164 * of digits given the limits of the
165 * buffer
166 */
167
168 int save_max_digits = max_digits;
169 if (fmode) {
170 /*
171 * This covers two cases:
172 *
173 * When exp is < 0, there are -exp-1 zeros taking up
174 * space before we can display any of the real digits,
175 * so we have to subtract those off max_decimals before
176 * we round that (max_decimals - (-exp - 1)). This
177 * may end up less than zero, in which case we have
178 * no digits to display.
179 *
180 * When exp >= 0, there are exp + 1 digits left of the
181 * decimal point *plus* max_decimals right of the
182 * decimal point that need to be generated
183 *
184 * A single expression gives the right answer in both
185 * cases, which is kinda cool
186 */
187 /* max_decimals comes in biased by 1 to flag the 'f' case */
188 max_digits = min(max_digits, max(1, max_decimals + decexp + 1));
189 }
190
191 int decimals = max_digits;
192
193 /* Round the value to the last digit being printed. */
194 uint64_t round = BIT64(59); /* 0.5 */
195 while (decimals--) {
196 round /= 10U;
197 }
198 fract += round;
199 /* Make sure rounding didn't make fract >= 1.0 */
200 if (fract >= BIT64(60)) {
201 fract /= 10U;
202 decexp++;
203 max_digits = min(save_max_digits, max(1, max_decimals + decexp + 1));
204 }
205
206 /* Now convert mantissa to decimal. */
207
208 /* Compute digits */
209 for (i = 0; i < max_digits; i++) {
210 fract *= 10U;
211 dtoa->digits[i] = (fract >> 60) + '0';
212 fract &= BIT64(60) - 1;
213 }
214 }
215 done:
216 dtoa->digits[max_digits] = '\0';
217 dtoa->flags = flags;
218 dtoa->exp = decexp;
219 return max_digits;
220 }
221
