]>
Commit | Line | Data |
---|---|---|
59e0d9fe | 1 | /*- |
31185420 | 2 | * Copyright (c) 2004, 2005 David Schultz <das@FreeBSD.ORG> |
59e0d9fe A |
3 | * All rights reserved. |
4 | * | |
5 | * Redistribution and use in source and binary forms, with or without | |
6 | * modification, are permitted provided that the following conditions | |
7 | * are met: | |
8 | * 1. Redistributions of source code must retain the above copyright | |
9 | * notice, this list of conditions and the following disclaimer. | |
10 | * 2. Redistributions in binary form must reproduce the above copyright | |
11 | * notice, this list of conditions and the following disclaimer in the | |
12 | * documentation and/or other materials provided with the distribution. | |
13 | * | |
14 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND | |
15 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | |
16 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | |
17 | * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE | |
18 | * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL | |
19 | * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS | |
20 | * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) | |
21 | * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT | |
22 | * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY | |
23 | * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF | |
24 | * SUCH DAMAGE. | |
25 | */ | |
26 | ||
27 | #include <sys/cdefs.h> | |
31185420 | 28 | __FBSDID("$FreeBSD: src/lib/libc/gdtoa/_hdtoa.c,v 1.3 2005/01/18 18:44:07 das Exp $"); |
59e0d9fe A |
29 | |
30 | #include <float.h> | |
59e0d9fe A |
31 | #include <limits.h> |
32 | #include <math.h> | |
59e0d9fe A |
33 | #include "fpmath.h" |
34 | #include "gdtoaimp.h" | |
35 | ||
36 | /* Strings values used by dtoa() */ | |
37 | #define INFSTR "Infinity" | |
38 | #define NANSTR "NaN" | |
39 | ||
31185420 A |
40 | #define DBL_ADJ (DBL_MAX_EXP - 2 + ((DBL_MANT_DIG - 1) % 4)) |
41 | #define LDBL_ADJ (LDBL_MAX_EXP - 2 + ((LDBL_MANT_DIG - 1) % 4)) | |
59e0d9fe A |
42 | |
43 | /* | |
44 | * Round up the given digit string. If the digit string is fff...f, | |
45 | * this procedure sets it to 100...0 and returns 1 to indicate that | |
46 | * the exponent needs to be bumped. Otherwise, 0 is returned. | |
47 | */ | |
48 | static int | |
49 | roundup(char *s0, int ndigits) | |
50 | { | |
51 | char *s; | |
52 | ||
53 | for (s = s0 + ndigits - 1; *s == 0xf; s--) { | |
54 | if (s == s0) { | |
55 | *s = 1; | |
56 | return (1); | |
57 | } | |
58 | ++*s; | |
59 | } | |
60 | ++*s; | |
61 | return (0); | |
62 | } | |
63 | ||
64 | /* | |
65 | * Round the given digit string to ndigits digits according to the | |
66 | * current rounding mode. Note that this could produce a string whose | |
67 | * value is not representable in the corresponding floating-point | |
68 | * type. The exponent pointed to by decpt is adjusted if necessary. | |
69 | */ | |
70 | static void | |
71 | dorounding(char *s0, int ndigits, int sign, int *decpt) | |
72 | { | |
73 | int adjust = 0; /* do we need to adjust the exponent? */ | |
74 | ||
75 | switch (FLT_ROUNDS) { | |
76 | case 0: /* toward zero */ | |
77 | default: /* implementation-defined */ | |
78 | break; | |
79 | case 1: /* to nearest, halfway rounds to even */ | |
80 | if ((s0[ndigits] > 8) || | |
81 | (s0[ndigits] == 8 && s0[ndigits - 1] & 1)) | |
82 | adjust = roundup(s0, ndigits); | |
83 | break; | |
84 | case 2: /* toward +inf */ | |
85 | if (sign == 0) | |
86 | adjust = roundup(s0, ndigits); | |
87 | break; | |
88 | case 3: /* toward -inf */ | |
89 | if (sign != 0) | |
90 | adjust = roundup(s0, ndigits); | |
91 | break; | |
92 | } | |
93 | ||
94 | if (adjust) | |
95 | *decpt += 4; | |
96 | } | |
97 | ||
98 | /* | |
99 | * This procedure converts a double-precision number in IEEE format | |
100 | * into a string of hexadecimal digits and an exponent of 2. Its | |
101 | * behavior is bug-for-bug compatible with dtoa() in mode 2, with the | |
102 | * following exceptions: | |
103 | * | |
104 | * - An ndigits < 0 causes it to use as many digits as necessary to | |
105 | * represent the number exactly. | |
106 | * - The additional xdigs argument should point to either the string | |
107 | * "0123456789ABCDEF" or the string "0123456789abcdef", depending on | |
108 | * which case is desired. | |
109 | * - This routine does not repeat dtoa's mistake of setting decpt | |
110 | * to 9999 in the case of an infinity or NaN. INT_MAX is used | |
111 | * for this purpose instead. | |
112 | * | |
113 | * Note that the C99 standard does not specify what the leading digit | |
114 | * should be for non-zero numbers. For instance, 0x1.3p3 is the same | |
115 | * as 0x2.6p2 is the same as 0x4.cp3. This implementation chooses the | |
116 | * first digit so that subsequent digits are aligned on nibble | |
117 | * boundaries (before rounding). | |
118 | * | |
119 | * Inputs: d, xdigs, ndigits | |
120 | * Outputs: decpt, sign, rve | |
121 | */ | |
122 | char * | |
123 | __hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign, | |
124 | char **rve) | |
125 | { | |
31185420 | 126 | static const int sigfigs = (DBL_MANT_DIG + 3) / 4; |
59e0d9fe A |
127 | union IEEEd2bits u; |
128 | char *s, *s0; | |
129 | int bufsize; | |
59e0d9fe A |
130 | |
131 | u.d = d; | |
132 | *sign = u.bits.sign; | |
133 | ||
134 | switch (fpclassify(d)) { | |
135 | case FP_NORMAL: | |
31185420 | 136 | *decpt = u.bits.exp - DBL_ADJ; |
59e0d9fe A |
137 | break; |
138 | case FP_ZERO: | |
139 | *decpt = 1; | |
140 | return (nrv_alloc("0", rve, 1)); | |
141 | case FP_SUBNORMAL: | |
31185420 A |
142 | u.d *= 0x1p514; |
143 | *decpt = u.bits.exp - (514 + DBL_ADJ); | |
59e0d9fe A |
144 | break; |
145 | case FP_INFINITE: | |
146 | *decpt = INT_MAX; | |
147 | return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1)); | |
148 | case FP_NAN: | |
149 | *decpt = INT_MAX; | |
150 | return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1)); | |
151 | default: | |
152 | abort(); | |
153 | } | |
154 | ||
155 | /* FP_NORMAL or FP_SUBNORMAL */ | |
156 | ||
157 | if (ndigits == 0) /* dtoa() compatibility */ | |
158 | ndigits = 1; | |
159 | ||
160 | /* | |
161 | * For simplicity, we generate all the digits even if the | |
162 | * caller has requested fewer. | |
163 | */ | |
164 | bufsize = (sigfigs > ndigits) ? sigfigs : ndigits; | |
165 | s0 = rv_alloc(bufsize); | |
166 | ||
167 | /* | |
168 | * We work from right to left, first adding any requested zero | |
169 | * padding, then the least significant portion of the | |
170 | * mantissa, followed by the most significant. The buffer is | |
171 | * filled with the byte values 0x0 through 0xf, which are | |
172 | * converted to xdigs[0x0] through xdigs[0xf] after the | |
173 | * rounding phase. | |
174 | */ | |
175 | for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--) | |
176 | *s = 0; | |
177 | for (; s > s0 + sigfigs - (DBL_MANL_SIZE / 4) - 1 && s > s0; s--) { | |
178 | *s = u.bits.manl & 0xf; | |
179 | u.bits.manl >>= 4; | |
180 | } | |
181 | for (; s > s0; s--) { | |
182 | *s = u.bits.manh & 0xf; | |
183 | u.bits.manh >>= 4; | |
184 | } | |
185 | ||
186 | /* | |
187 | * At this point, we have snarfed all the bits in the | |
188 | * mantissa, with the possible exception of the highest-order | |
189 | * (partial) nibble, which is dealt with by the next | |
31185420 | 190 | * statement. We also tack on the implicit normalization bit. |
59e0d9fe | 191 | */ |
31185420 | 192 | *s = u.bits.manh | (1U << ((DBL_MANT_DIG - 1) % 4)); |
59e0d9fe A |
193 | |
194 | /* If ndigits < 0, we are expected to auto-size the precision. */ | |
195 | if (ndigits < 0) { | |
196 | for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--) | |
197 | ; | |
198 | } | |
199 | ||
200 | if (sigfigs > ndigits && s0[ndigits] != 0) | |
201 | dorounding(s0, ndigits, u.bits.sign, decpt); | |
202 | ||
203 | s = s0 + ndigits; | |
204 | if (rve != NULL) | |
205 | *rve = s; | |
206 | *s-- = '\0'; | |
207 | for (; s >= s0; s--) | |
208 | *s = xdigs[(unsigned int)*s]; | |
209 | ||
210 | return (s0); | |
211 | } | |
212 | ||
213 | #if (LDBL_MANT_DIG > DBL_MANT_DIG) | |
214 | ||
215 | /* | |
216 | * This is the long double version of __hdtoa(). | |
59e0d9fe A |
217 | */ |
218 | char * | |
219 | __hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign, | |
220 | char **rve) | |
221 | { | |
31185420 | 222 | static const int sigfigs = (LDBL_MANT_DIG + 3) / 4; |
59e0d9fe A |
223 | union IEEEl2bits u; |
224 | char *s, *s0; | |
225 | int bufsize; | |
59e0d9fe A |
226 | |
227 | u.e = e; | |
228 | *sign = u.bits.sign; | |
229 | ||
230 | switch (fpclassify(e)) { | |
231 | case FP_NORMAL: | |
31185420 | 232 | *decpt = u.bits.exp - LDBL_ADJ; |
59e0d9fe A |
233 | break; |
234 | case FP_ZERO: | |
235 | *decpt = 1; | |
236 | return (nrv_alloc("0", rve, 1)); | |
237 | case FP_SUBNORMAL: | |
31185420 A |
238 | u.e *= 0x1p514L; |
239 | *decpt = u.bits.exp - (514 + LDBL_ADJ); | |
59e0d9fe A |
240 | break; |
241 | case FP_INFINITE: | |
242 | *decpt = INT_MAX; | |
243 | return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1)); | |
244 | case FP_NAN: | |
245 | *decpt = INT_MAX; | |
246 | return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1)); | |
247 | default: | |
248 | abort(); | |
249 | } | |
250 | ||
251 | /* FP_NORMAL or FP_SUBNORMAL */ | |
252 | ||
253 | if (ndigits == 0) /* dtoa() compatibility */ | |
254 | ndigits = 1; | |
255 | ||
256 | /* | |
257 | * For simplicity, we generate all the digits even if the | |
258 | * caller has requested fewer. | |
259 | */ | |
260 | bufsize = (sigfigs > ndigits) ? sigfigs : ndigits; | |
261 | s0 = rv_alloc(bufsize); | |
262 | ||
263 | /* | |
264 | * We work from right to left, first adding any requested zero | |
265 | * padding, then the least significant portion of the | |
266 | * mantissa, followed by the most significant. The buffer is | |
267 | * filled with the byte values 0x0 through 0xf, which are | |
268 | * converted to xdigs[0x0] through xdigs[0xf] after the | |
269 | * rounding phase. | |
270 | */ | |
271 | for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--) | |
272 | *s = 0; | |
273 | for (; s > s0 + sigfigs - (LDBL_MANL_SIZE / 4) - 1 && s > s0; s--) { | |
274 | *s = u.bits.manl & 0xf; | |
275 | u.bits.manl >>= 4; | |
276 | } | |
277 | for (; s > s0; s--) { | |
278 | *s = u.bits.manh & 0xf; | |
279 | u.bits.manh >>= 4; | |
280 | } | |
281 | ||
282 | /* | |
283 | * At this point, we have snarfed all the bits in the | |
284 | * mantissa, with the possible exception of the highest-order | |
285 | * (partial) nibble, which is dealt with by the next | |
31185420 | 286 | * statement. We also tack on the implicit normalization bit. |
59e0d9fe | 287 | */ |
31185420 | 288 | *s = u.bits.manh | (1U << ((LDBL_MANT_DIG - 1) % 4)); |
59e0d9fe A |
289 | |
290 | /* If ndigits < 0, we are expected to auto-size the precision. */ | |
291 | if (ndigits < 0) { | |
292 | for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--) | |
293 | ; | |
294 | } | |
295 | ||
296 | if (sigfigs > ndigits && s0[ndigits] != 0) | |
297 | dorounding(s0, ndigits, u.bits.sign, decpt); | |
298 | ||
299 | s = s0 + ndigits; | |
300 | if (rve != NULL) | |
301 | *rve = s; | |
302 | *s-- = '\0'; | |
303 | for (; s >= s0; s--) | |
304 | *s = xdigs[(unsigned int)*s]; | |
305 | ||
306 | return (s0); | |
307 | } | |
308 | ||
309 | #else /* (LDBL_MANT_DIG == DBL_MANT_DIG) */ | |
310 | ||
311 | char * | |
312 | __hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign, | |
313 | char **rve) | |
314 | { | |
315 | ||
316 | return (__hdtoa((double)e, xdigs, ndigits, decpt, sign, rve)); | |
317 | } | |
318 | ||
319 | #endif /* (LDBL_MANT_DIG == DBL_MANT_DIG) */ |