]>
Commit | Line | Data |
---|---|---|
2aa62f2b | 1 | /* $NetBSD: dtoa.c,v 1.3.4.1.4.1 2008/04/08 21:10:55 jdc Exp $ */\r |
2 | \r | |
3 | /****************************************************************\r | |
4 | \r | |
5 | The author of this software is David M. Gay.\r | |
6 | \r | |
7 | Copyright (C) 1998, 1999 by Lucent Technologies\r | |
8 | All Rights Reserved\r | |
9 | \r | |
10 | Permission to use, copy, modify, and distribute this software and\r | |
11 | its documentation for any purpose and without fee is hereby\r | |
12 | granted, provided that the above copyright notice appear in all\r | |
13 | copies and that both that the copyright notice and this\r | |
14 | permission notice and warranty disclaimer appear in supporting\r | |
15 | documentation, and that the name of Lucent or any of its entities\r | |
16 | not be used in advertising or publicity pertaining to\r | |
17 | distribution of the software without specific, written prior\r | |
18 | permission.\r | |
19 | \r | |
20 | LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,\r | |
21 | INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.\r | |
22 | IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY\r | |
23 | SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES\r | |
24 | WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER\r | |
25 | IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,\r | |
26 | ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF\r | |
27 | THIS SOFTWARE.\r | |
28 | \r | |
29 | ****************************************************************/\r | |
30 | \r | |
31 | /* Please send bug reports to David M. Gay (dmg at acm dot org,\r | |
32 | * with " at " changed at "@" and " dot " changed to "."). */\r | |
33 | #include <LibConfig.h>\r | |
34 | \r | |
35 | #include "gdtoaimp.h"\r | |
36 | \r | |
37 | /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.\r | |
38 | *\r | |
39 | * Inspired by "How to Print Floating-Point Numbers Accurately" by\r | |
40 | * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].\r | |
41 | *\r | |
42 | * Modifications:\r | |
43 | * 1. Rather than iterating, we use a simple numeric overestimate\r | |
44 | * to determine k = floor(log10(d)). We scale relevant\r | |
45 | * quantities using O(log2(k)) rather than O(k) multiplications.\r | |
46 | * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't\r | |
47 | * try to generate digits strictly left to right. Instead, we\r | |
48 | * compute with fewer bits and propagate the carry if necessary\r | |
49 | * when rounding the final digit up. This is often faster.\r | |
50 | * 3. Under the assumption that input will be rounded nearest,\r | |
51 | * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.\r | |
52 | * That is, we allow equality in stopping tests when the\r | |
53 | * round-nearest rule will give the same floating-point value\r | |
54 | * as would satisfaction of the stopping test with strict\r | |
55 | * inequality.\r | |
56 | * 4. We remove common factors of powers of 2 from relevant\r | |
57 | * quantities.\r | |
58 | * 5. When converting floating-point integers less than 1e16,\r | |
59 | * we use floating-point arithmetic rather than resorting\r | |
60 | * to multiple-precision integers.\r | |
61 | * 6. When asked to produce fewer than 15 digits, we first try\r | |
62 | * to get by with floating-point arithmetic; we resort to\r | |
63 | * multiple-precision integer arithmetic only if we cannot\r | |
64 | * guarantee that the floating-point calculation has given\r | |
65 | * the correctly rounded result. For k requested digits and\r | |
66 | * "uniformly" distributed input, the probability is\r | |
67 | * something like 10^(k-15) that we must resort to the Long\r | |
68 | * calculation.\r | |
69 | */\r | |
70 | \r | |
71 | #ifdef Honor_FLT_ROUNDS\r | |
72 | #define Rounding rounding\r | |
73 | #undef Check_FLT_ROUNDS\r | |
74 | #define Check_FLT_ROUNDS\r | |
75 | #else\r | |
76 | #define Rounding Flt_Rounds\r | |
77 | #endif\r | |
78 | \r | |
79 | #if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */\r | |
80 | // Disable: warning C4700: uninitialized local variable 'xx' used\r | |
81 | #pragma warning ( disable : 4700 )\r | |
82 | #endif /* defined(_MSC_VER) */\r | |
83 | \r | |
84 | char *\r | |
85 | dtoa\r | |
86 | #ifdef KR_headers\r | |
87 | (d, mode, ndigits, decpt, sign, rve)\r | |
88 | double d; int mode, ndigits, *decpt, *sign; char **rve;\r | |
89 | #else\r | |
90 | (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)\r | |
91 | #endif\r | |
92 | {\r | |
93 | /* Arguments ndigits, decpt, sign are similar to those\r | |
94 | of ecvt and fcvt; trailing zeros are suppressed from\r | |
95 | the returned string. If not null, *rve is set to point\r | |
96 | to the end of the return value. If d is +-Infinity or NaN,\r | |
97 | then *decpt is set to 9999.\r | |
98 | \r | |
99 | mode:\r | |
100 | 0 ==> shortest string that yields d when read in\r | |
101 | and rounded to nearest.\r | |
102 | 1 ==> like 0, but with Steele & White stopping rule;\r | |
103 | e.g. with IEEE P754 arithmetic , mode 0 gives\r | |
104 | 1e23 whereas mode 1 gives 9.999999999999999e22.\r | |
105 | 2 ==> max(1,ndigits) significant digits. This gives a\r | |
106 | return value similar to that of ecvt, except\r | |
107 | that trailing zeros are suppressed.\r | |
108 | 3 ==> through ndigits past the decimal point. This\r | |
109 | gives a return value similar to that from fcvt,\r | |
110 | except that trailing zeros are suppressed, and\r | |
111 | ndigits can be negative.\r | |
112 | 4,5 ==> similar to 2 and 3, respectively, but (in\r | |
113 | round-nearest mode) with the tests of mode 0 to\r | |
114 | possibly return a shorter string that rounds to d.\r | |
115 | With IEEE arithmetic and compilation with\r | |
116 | -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same\r | |
117 | as modes 2 and 3 when FLT_ROUNDS != 1.\r | |
118 | 6-9 ==> Debugging modes similar to mode - 4: don't try\r | |
119 | fast floating-point estimate (if applicable).\r | |
120 | \r | |
121 | Values of mode other than 0-9 are treated as mode 0.\r | |
122 | \r | |
123 | Sufficient space is allocated to the return value\r | |
124 | to hold the suppressed trailing zeros.\r | |
125 | */\r | |
126 | \r | |
127 | int bbits, b2, b5, be, dig, i, ieps, ilim0,\r | |
128 | j, jj1, k, k0, k_check, leftright, m2, m5, s2, s5,\r | |
129 | spec_case, try_quick;\r | |
130 | int ilim = 0, ilim1 = 0; /* pacify gcc */\r | |
131 | Long L;\r | |
132 | #ifndef Sudden_Underflow\r | |
133 | int denorm;\r | |
134 | ULong x;\r | |
135 | #endif\r | |
136 | Bigint *b, *b1, *delta, *mhi, *S;\r | |
137 | Bigint *mlo = NULL; /* pacify gcc */\r | |
138 | double d2, ds, eps;\r | |
139 | char *s, *s0;\r | |
140 | #ifdef Honor_FLT_ROUNDS\r | |
141 | int rounding;\r | |
142 | #endif\r | |
143 | #ifdef SET_INEXACT\r | |
144 | int inexact, oldinexact;\r | |
145 | #endif\r | |
146 | \r | |
147 | #ifndef MULTIPLE_THREADS\r | |
148 | if (dtoa_result) {\r | |
149 | freedtoa(dtoa_result);\r | |
150 | dtoa_result = 0;\r | |
151 | }\r | |
152 | #endif\r | |
153 | \r | |
154 | if (word0(d) & Sign_bit) {\r | |
155 | /* set sign for everything, including 0's and NaNs */\r | |
156 | *sign = 1;\r | |
157 | word0(d) &= ~Sign_bit; /* clear sign bit */\r | |
158 | }\r | |
159 | else\r | |
160 | *sign = 0;\r | |
161 | \r | |
162 | #if defined(IEEE_Arith) + defined(VAX)\r | |
163 | #ifdef IEEE_Arith\r | |
164 | if ((word0(d) & Exp_mask) == Exp_mask)\r | |
165 | #else\r | |
166 | if (word0(d) == 0x8000)\r | |
167 | #endif\r | |
168 | {\r | |
169 | /* Infinity or NaN */\r | |
170 | *decpt = 9999;\r | |
171 | #ifdef IEEE_Arith\r | |
172 | if (!word1(d) && !(word0(d) & 0xfffff))\r | |
173 | return nrv_alloc("Infinity", rve, 8);\r | |
174 | #endif\r | |
175 | return nrv_alloc("NaN", rve, 3);\r | |
176 | }\r | |
177 | #endif\r | |
178 | #ifdef IBM\r | |
179 | dval(d) += 0; /* normalize */\r | |
180 | #endif\r | |
181 | if (!dval(d)) {\r | |
182 | *decpt = 1;\r | |
183 | return nrv_alloc("0", rve, 1);\r | |
184 | }\r | |
185 | \r | |
186 | #ifdef SET_INEXACT\r | |
187 | try_quick = oldinexact = get_inexact();\r | |
188 | inexact = 1;\r | |
189 | #endif\r | |
190 | #ifdef Honor_FLT_ROUNDS\r | |
191 | if ((rounding = Flt_Rounds) >= 2) {\r | |
192 | if (*sign)\r | |
193 | rounding = rounding == 2 ? 0 : 2;\r | |
194 | else\r | |
195 | if (rounding != 2)\r | |
196 | rounding = 0;\r | |
197 | }\r | |
198 | #endif\r | |
199 | \r | |
200 | b = d2b(dval(d), &be, &bbits);\r | |
201 | if (b == NULL)\r | |
202 | return NULL;\r | |
203 | #ifdef Sudden_Underflow\r | |
204 | i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));\r | |
205 | #else\r | |
206 | if (( i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {\r | |
207 | #endif\r | |
208 | dval(d2) = dval(d);\r | |
209 | word0(d2) &= Frac_mask1;\r | |
210 | word0(d2) |= Exp_11;\r | |
211 | #ifdef IBM\r | |
212 | if (( j = 11 - hi0bits(word0(d2) & Frac_mask) )!=0)\r | |
213 | dval(d2) /= 1 << j;\r | |
214 | #endif\r | |
215 | \r | |
216 | /* log(x) ~=~ log(1.5) + (x-1.5)/1.5\r | |
217 | * log10(x) = log(x) / log(10)\r | |
218 | * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))\r | |
219 | * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)\r | |
220 | *\r | |
221 | * This suggests computing an approximation k to log10(d) by\r | |
222 | *\r | |
223 | * k = (i - Bias)*0.301029995663981\r | |
224 | * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );\r | |
225 | *\r | |
226 | * We want k to be too large rather than too small.\r | |
227 | * The error in the first-order Taylor series approximation\r | |
228 | * is in our favor, so we just round up the constant enough\r | |
229 | * to compensate for any error in the multiplication of\r | |
230 | * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,\r | |
231 | * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,\r | |
232 | * adding 1e-13 to the constant term more than suffices.\r | |
233 | * Hence we adjust the constant term to 0.1760912590558.\r | |
234 | * (We could get a more accurate k by invoking log10,\r | |
235 | * but this is probably not worthwhile.)\r | |
236 | */\r | |
237 | \r | |
238 | i -= Bias;\r | |
239 | #ifdef IBM\r | |
240 | i <<= 2;\r | |
241 | i += j;\r | |
242 | #endif\r | |
243 | #ifndef Sudden_Underflow\r | |
244 | denorm = 0;\r | |
245 | }\r | |
246 | else {\r | |
247 | /* d is denormalized */\r | |
248 | \r | |
249 | i = bbits + be + (Bias + (P-1) - 1);\r | |
250 | x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)\r | |
251 | : word1(d) << (32 - i);\r | |
252 | dval(d2) = (double)x;\r | |
253 | word0(d2) -= 31*Exp_msk1; /* adjust exponent */\r | |
254 | i -= (Bias + (P-1) - 1) + 1;\r | |
255 | denorm = 1;\r | |
256 | }\r | |
257 | #endif\r | |
258 | ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;\r | |
259 | k = (int)ds;\r | |
260 | if (ds < 0. && ds != k)\r | |
261 | k--; /* want k = floor(ds) */\r | |
262 | k_check = 1;\r | |
263 | if (k >= 0 && k <= Ten_pmax) {\r | |
264 | if (dval(d) < tens[k])\r | |
265 | k--;\r | |
266 | k_check = 0;\r | |
267 | }\r | |
268 | j = bbits - i - 1;\r | |
269 | if (j >= 0) {\r | |
270 | b2 = 0;\r | |
271 | s2 = j;\r | |
272 | }\r | |
273 | else {\r | |
274 | b2 = -j;\r | |
275 | s2 = 0;\r | |
276 | }\r | |
277 | if (k >= 0) {\r | |
278 | b5 = 0;\r | |
279 | s5 = k;\r | |
280 | s2 += k;\r | |
281 | }\r | |
282 | else {\r | |
283 | b2 -= k;\r | |
284 | b5 = -k;\r | |
285 | s5 = 0;\r | |
286 | }\r | |
287 | if (mode < 0 || mode > 9)\r | |
288 | mode = 0;\r | |
289 | \r | |
290 | #ifndef SET_INEXACT\r | |
291 | #ifdef Check_FLT_ROUNDS\r | |
292 | try_quick = Rounding == 1;\r | |
293 | #else\r | |
294 | try_quick = 1;\r | |
295 | #endif\r | |
296 | #endif /*SET_INEXACT*/\r | |
297 | \r | |
298 | if (mode > 5) {\r | |
299 | mode -= 4;\r | |
300 | try_quick = 0;\r | |
301 | }\r | |
302 | leftright = 1;\r | |
303 | switch(mode) {\r | |
304 | case 0:\r | |
305 | case 1:\r | |
306 | ilim = ilim1 = -1;\r | |
307 | i = 18;\r | |
308 | ndigits = 0;\r | |
309 | break;\r | |
310 | case 2:\r | |
311 | leftright = 0;\r | |
312 | /* FALLTHROUGH */\r | |
313 | case 4:\r | |
314 | if (ndigits <= 0)\r | |
315 | ndigits = 1;\r | |
316 | ilim = ilim1 = i = ndigits;\r | |
317 | break;\r | |
318 | case 3:\r | |
319 | leftright = 0;\r | |
320 | /* FALLTHROUGH */\r | |
321 | case 5:\r | |
322 | i = ndigits + k + 1;\r | |
323 | ilim = i;\r | |
324 | ilim1 = i - 1;\r | |
325 | if (i <= 0)\r | |
326 | i = 1;\r | |
327 | }\r | |
328 | s = s0 = rv_alloc((size_t)i);\r | |
329 | if (s == NULL)\r | |
330 | return NULL;\r | |
331 | \r | |
332 | #ifdef Honor_FLT_ROUNDS\r | |
333 | if (mode > 1 && rounding != 1)\r | |
334 | leftright = 0;\r | |
335 | #endif\r | |
336 | \r | |
337 | if (ilim >= 0 && ilim <= Quick_max && try_quick) {\r | |
338 | \r | |
339 | /* Try to get by with floating-point arithmetic. */\r | |
340 | \r | |
341 | i = 0;\r | |
342 | dval(d2) = dval(d);\r | |
343 | k0 = k;\r | |
344 | ilim0 = ilim;\r | |
345 | ieps = 2; /* conservative */\r | |
346 | if (k > 0) {\r | |
347 | ds = tens[k&0xf];\r | |
348 | j = (unsigned int)k >> 4;\r | |
349 | if (j & Bletch) {\r | |
350 | /* prevent overflows */\r | |
351 | j &= Bletch - 1;\r | |
352 | dval(d) /= bigtens[n_bigtens-1];\r | |
353 | ieps++;\r | |
354 | }\r | |
355 | for(; j; j = (unsigned int)j >> 1, i++)\r | |
356 | if (j & 1) {\r | |
357 | ieps++;\r | |
358 | ds *= bigtens[i];\r | |
359 | }\r | |
360 | dval(d) /= ds;\r | |
361 | }\r | |
362 | else if (( jj1 = -k )!=0) {\r | |
363 | dval(d) *= tens[jj1 & 0xf];\r | |
364 | for(j = jj1 >> 4; j; j >>= 1, i++)\r | |
365 | if (j & 1) {\r | |
366 | ieps++;\r | |
367 | dval(d) *= bigtens[i];\r | |
368 | }\r | |
369 | }\r | |
370 | if (k_check && dval(d) < 1. && ilim > 0) {\r | |
371 | if (ilim1 <= 0)\r | |
372 | goto fast_failed;\r | |
373 | ilim = ilim1;\r | |
374 | k--;\r | |
375 | dval(d) *= 10.;\r | |
376 | ieps++;\r | |
377 | }\r | |
378 | dval(eps) = ieps*dval(d) + 7.;\r | |
379 | word0(eps) -= (P-1)*Exp_msk1;\r | |
380 | if (ilim == 0) {\r | |
381 | S = mhi = 0;\r | |
382 | dval(d) -= 5.;\r | |
383 | if (dval(d) > dval(eps))\r | |
384 | goto one_digit;\r | |
385 | if (dval(d) < -dval(eps))\r | |
386 | goto no_digits;\r | |
387 | goto fast_failed;\r | |
388 | }\r | |
389 | #ifndef No_leftright\r | |
390 | if (leftright) {\r | |
391 | /* Use Steele & White method of only\r | |
392 | * generating digits needed.\r | |
393 | */\r | |
394 | dval(eps) = 0.5/tens[ilim-1] - dval(eps);\r | |
395 | for(i = 0;;) {\r | |
396 | L = (INT32)dval(d);\r | |
397 | dval(d) -= L;\r | |
398 | *s++ = (char)('0' + (int)L);\r | |
399 | if (dval(d) < dval(eps))\r | |
400 | goto ret1;\r | |
401 | if (1. - dval(d) < dval(eps))\r | |
402 | goto bump_up;\r | |
403 | if (++i >= ilim)\r | |
404 | break;\r | |
405 | dval(eps) *= 10.;\r | |
406 | dval(d) *= 10.;\r | |
407 | }\r | |
408 | }\r | |
409 | else {\r | |
410 | #endif\r | |
411 | /* Generate ilim digits, then fix them up. */\r | |
412 | dval(eps) *= tens[ilim-1];\r | |
413 | for(i = 1;; i++, dval(d) *= 10.) {\r | |
414 | L = (Long)(dval(d));\r | |
415 | if (!(dval(d) -= L))\r | |
416 | ilim = i;\r | |
417 | *s++ = (char)('0' + (int)L);\r | |
418 | if (i == ilim) {\r | |
419 | if (dval(d) > 0.5 + dval(eps))\r | |
420 | goto bump_up;\r | |
421 | else if (dval(d) < 0.5 - dval(eps)) {\r | |
422 | while(*--s == '0');\r | |
423 | s++;\r | |
424 | goto ret1;\r | |
425 | }\r | |
426 | break;\r | |
427 | }\r | |
428 | }\r | |
429 | #ifndef No_leftright\r | |
430 | }\r | |
431 | #endif\r | |
432 | fast_failed:\r | |
433 | s = s0;\r | |
434 | dval(d) = dval(d2);\r | |
435 | k = k0;\r | |
436 | ilim = ilim0;\r | |
437 | }\r | |
438 | \r | |
439 | /* Do we have a "small" integer? */\r | |
440 | \r | |
441 | if (be >= 0 && k <= Int_max) {\r | |
442 | /* Yes. */\r | |
443 | ds = tens[k];\r | |
444 | if (ndigits < 0 && ilim <= 0) {\r | |
445 | S = mhi = 0;\r | |
446 | if (ilim < 0 || dval(d) <= 5*ds)\r | |
447 | goto no_digits;\r | |
448 | goto one_digit;\r | |
449 | }\r | |
450 | for(i = 1;; i++, dval(d) *= 10.) {\r | |
451 | L = (Long)(dval(d) / ds);\r | |
452 | dval(d) -= L*ds;\r | |
453 | #ifdef Check_FLT_ROUNDS\r | |
454 | /* If FLT_ROUNDS == 2, L will usually be high by 1 */\r | |
455 | if (dval(d) < 0) {\r | |
456 | L--;\r | |
457 | dval(d) += ds;\r | |
458 | }\r | |
459 | #endif\r | |
460 | *s++ = (char)('0' + (int)L);\r | |
461 | if (!dval(d)) {\r | |
462 | #ifdef SET_INEXACT\r | |
463 | inexact = 0;\r | |
464 | #endif\r | |
465 | break;\r | |
466 | }\r | |
467 | if (i == ilim) {\r | |
468 | #ifdef Honor_FLT_ROUNDS\r | |
469 | if (mode > 1)\r | |
470 | switch(rounding) {\r | |
471 | case 0: goto ret1;\r | |
472 | case 2: goto bump_up;\r | |
473 | }\r | |
474 | #endif\r | |
475 | dval(d) += dval(d);\r | |
476 | if (dval(d) > ds || (dval(d) == ds && L & 1)) {\r | |
477 | bump_up:\r | |
478 | while(*--s == '9')\r | |
479 | if (s == s0) {\r | |
480 | k++;\r | |
481 | *s = '0';\r | |
482 | break;\r | |
483 | }\r | |
484 | ++*s++;\r | |
485 | }\r | |
486 | break;\r | |
487 | }\r | |
488 | }\r | |
489 | goto ret1;\r | |
490 | }\r | |
491 | \r | |
492 | m2 = b2;\r | |
493 | m5 = b5;\r | |
494 | mhi = mlo = 0;\r | |
495 | if (leftright) {\r | |
496 | i =\r | |
497 | #ifndef Sudden_Underflow\r | |
498 | denorm ? be + (Bias + (P-1) - 1 + 1) :\r | |
499 | #endif\r | |
500 | #ifdef IBM\r | |
501 | 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);\r | |
502 | #else\r | |
503 | 1 + P - bbits;\r | |
504 | #endif\r | |
505 | b2 += i;\r | |
506 | s2 += i;\r | |
507 | mhi = i2b(1);\r | |
508 | if (mhi == NULL)\r | |
509 | return NULL;\r | |
510 | }\r | |
511 | if (m2 > 0 && s2 > 0) {\r | |
512 | i = m2 < s2 ? m2 : s2;\r | |
513 | b2 -= i;\r | |
514 | m2 -= i;\r | |
515 | s2 -= i;\r | |
516 | }\r | |
517 | if (b5 > 0) {\r | |
518 | if (leftright) {\r | |
519 | if (m5 > 0) {\r | |
520 | mhi = pow5mult(mhi, m5);\r | |
521 | if (mhi == NULL)\r | |
522 | return NULL;\r | |
523 | b1 = mult(mhi, b);\r | |
524 | if (b1 == NULL)\r | |
525 | return NULL;\r | |
526 | Bfree(b);\r | |
527 | b = b1;\r | |
528 | }\r | |
65ed9d7f LL |
529 | if (( j = b5 - m5 )!=0)\r |
530 | b = pow5mult(b, j);\r | |
2aa62f2b | 531 | if (b == NULL)\r |
532 | return NULL;\r | |
533 | }\r | |
534 | else\r | |
535 | b = pow5mult(b, b5);\r | |
65ed9d7f LL |
536 | if (b == NULL)\r |
537 | return NULL;\r | |
2aa62f2b | 538 | }\r |
539 | S = i2b(1);\r | |
540 | if (S == NULL)\r | |
541 | return NULL;\r | |
542 | if (s5 > 0) {\r | |
543 | S = pow5mult(S, s5);\r | |
544 | if (S == NULL)\r | |
545 | return NULL;\r | |
546 | }\r | |
547 | \r | |
548 | /* Check for special case that d is a normalized power of 2. */\r | |
549 | \r | |
550 | spec_case = 0;\r | |
551 | if ((mode < 2 || leftright)\r | |
552 | #ifdef Honor_FLT_ROUNDS\r | |
553 | && rounding == 1\r | |
554 | #endif\r | |
555 | ) {\r | |
556 | if (!word1(d) && !(word0(d) & Bndry_mask)\r | |
557 | #ifndef Sudden_Underflow\r | |
558 | && word0(d) & (Exp_mask & ~Exp_msk1)\r | |
559 | #endif\r | |
560 | ) {\r | |
561 | /* The special case */\r | |
562 | b2 += Log2P;\r | |
563 | s2 += Log2P;\r | |
564 | spec_case = 1;\r | |
565 | }\r | |
566 | }\r | |
567 | \r | |
568 | /* Arrange for convenient computation of quotients:\r | |
569 | * shift left if necessary so divisor has 4 leading 0 bits.\r | |
570 | *\r | |
571 | * Perhaps we should just compute leading 28 bits of S once\r | |
572 | * and for all and pass them and a shift to quorem, so it\r | |
573 | * can do shifts and ors to compute the numerator for q.\r | |
574 | */\r | |
575 | #ifdef Pack_32\r | |
576 | if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f )!=0)\r | |
577 | i = 32 - i;\r | |
578 | #else\r | |
579 | if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf )!=0)\r | |
580 | i = 16 - i;\r | |
581 | #endif\r | |
582 | if (i > 4) {\r | |
583 | i -= 4;\r | |
584 | b2 += i;\r | |
585 | m2 += i;\r | |
586 | s2 += i;\r | |
587 | }\r | |
588 | else if (i < 4) {\r | |
589 | i += 28;\r | |
590 | b2 += i;\r | |
591 | m2 += i;\r | |
592 | s2 += i;\r | |
593 | }\r | |
594 | if (b2 > 0) {\r | |
595 | b = lshift(b, b2);\r | |
596 | if (b == NULL)\r | |
597 | return NULL;\r | |
598 | }\r | |
599 | if (s2 > 0) {\r | |
600 | S = lshift(S, s2);\r | |
601 | if (S == NULL)\r | |
602 | return NULL;\r | |
603 | }\r | |
604 | if (k_check) {\r | |
605 | if (cmp(b,S) < 0) {\r | |
606 | k--;\r | |
607 | b = multadd(b, 10, 0); /* we botched the k estimate */\r | |
608 | if (b == NULL)\r | |
609 | return NULL;\r | |
610 | if (leftright) {\r | |
611 | mhi = multadd(mhi, 10, 0);\r | |
612 | if (mhi == NULL)\r | |
613 | return NULL;\r | |
614 | }\r | |
615 | ilim = ilim1;\r | |
616 | }\r | |
617 | }\r | |
618 | if (ilim <= 0 && (mode == 3 || mode == 5)) {\r | |
619 | if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {\r | |
620 | /* no digits, fcvt style */\r | |
621 | no_digits:\r | |
622 | k = -1 - ndigits;\r | |
623 | goto ret;\r | |
624 | }\r | |
625 | one_digit:\r | |
626 | *s++ = '1';\r | |
627 | k++;\r | |
628 | goto ret;\r | |
629 | }\r | |
630 | if (leftright) {\r | |
631 | if (m2 > 0) {\r | |
632 | mhi = lshift(mhi, m2);\r | |
633 | if (mhi == NULL)\r | |
634 | return NULL;\r | |
635 | }\r | |
636 | \r | |
637 | /* Compute mlo -- check for special case\r | |
638 | * that d is a normalized power of 2.\r | |
639 | */\r | |
640 | \r | |
641 | mlo = mhi;\r | |
642 | if (spec_case) {\r | |
643 | mhi = Balloc(mhi->k);\r | |
644 | if (mhi == NULL)\r | |
645 | return NULL;\r | |
646 | Bcopy(mhi, mlo);\r | |
647 | mhi = lshift(mhi, Log2P);\r | |
648 | if (mhi == NULL)\r | |
649 | return NULL;\r | |
650 | }\r | |
651 | \r | |
652 | for(i = 1;;i++) {\r | |
653 | dig = quorem(b,S) + '0';\r | |
654 | /* Do we yet have the shortest decimal string\r | |
655 | * that will round to d?\r | |
656 | */\r | |
657 | j = cmp(b, mlo);\r | |
658 | delta = diff(S, mhi);\r | |
659 | if (delta == NULL)\r | |
660 | return NULL;\r | |
661 | jj1 = delta->sign ? 1 : cmp(b, delta);\r | |
662 | Bfree(delta);\r | |
663 | #ifndef ROUND_BIASED\r | |
664 | if (jj1 == 0 && mode != 1 && !(word1(d) & 1)\r | |
665 | #ifdef Honor_FLT_ROUNDS\r | |
666 | && rounding >= 1\r | |
667 | #endif\r | |
668 | ) {\r | |
669 | if (dig == '9')\r | |
670 | goto round_9_up;\r | |
671 | if (j > 0)\r | |
672 | dig++;\r | |
673 | #ifdef SET_INEXACT\r | |
674 | else if (!b->x[0] && b->wds <= 1)\r | |
675 | inexact = 0;\r | |
676 | #endif\r | |
677 | *s++ = (char)dig;\r | |
678 | goto ret;\r | |
679 | }\r | |
680 | #endif\r | |
681 | if (j < 0 || (j == 0 && mode != 1\r | |
682 | #ifndef ROUND_BIASED\r | |
683 | && !(word1(d) & 1)\r | |
684 | #endif\r | |
685 | )) {\r | |
686 | if (!b->x[0] && b->wds <= 1) {\r | |
687 | #ifdef SET_INEXACT\r | |
688 | inexact = 0;\r | |
689 | #endif\r | |
690 | goto accept_dig;\r | |
691 | }\r | |
692 | #ifdef Honor_FLT_ROUNDS\r | |
693 | if (mode > 1)\r | |
694 | switch(rounding) {\r | |
695 | case 0: goto accept_dig;\r | |
696 | case 2: goto keep_dig;\r | |
697 | }\r | |
698 | #endif /*Honor_FLT_ROUNDS*/\r | |
699 | if (jj1 > 0) {\r | |
700 | b = lshift(b, 1);\r | |
701 | if (b == NULL)\r | |
702 | return NULL;\r | |
703 | jj1 = cmp(b, S);\r | |
704 | if ((jj1 > 0 || (jj1 == 0 && dig & 1))\r | |
705 | && dig++ == '9')\r | |
706 | goto round_9_up;\r | |
707 | }\r | |
708 | accept_dig:\r | |
709 | *s++ = (char)dig;\r | |
710 | goto ret;\r | |
711 | }\r | |
712 | if (jj1 > 0) {\r | |
713 | #ifdef Honor_FLT_ROUNDS\r | |
714 | if (!rounding)\r | |
715 | goto accept_dig;\r | |
716 | #endif\r | |
717 | if (dig == '9') { /* possible if i == 1 */\r | |
718 | round_9_up:\r | |
719 | *s++ = '9';\r | |
720 | goto roundoff;\r | |
721 | }\r | |
722 | *s++ = (char)(dig + 1);\r | |
723 | goto ret;\r | |
724 | }\r | |
725 | #ifdef Honor_FLT_ROUNDS\r | |
726 | keep_dig:\r | |
727 | #endif\r | |
728 | *s++ = (char)dig;\r | |
729 | if (i == ilim)\r | |
730 | break;\r | |
731 | b = multadd(b, 10, 0);\r | |
732 | if (b == NULL)\r | |
733 | return NULL;\r | |
734 | if (mlo == mhi) {\r | |
735 | mlo = mhi = multadd(mhi, 10, 0);\r | |
736 | if (mlo == NULL)\r | |
737 | return NULL;\r | |
738 | }\r | |
739 | else {\r | |
740 | mlo = multadd(mlo, 10, 0);\r | |
741 | if (mlo == NULL)\r | |
742 | return NULL;\r | |
743 | mhi = multadd(mhi, 10, 0);\r | |
744 | if (mhi == NULL)\r | |
745 | return NULL;\r | |
746 | }\r | |
747 | }\r | |
748 | }\r | |
749 | else\r | |
750 | for(i = 1;; i++) {\r | |
751 | *s++ = (char)(dig = (int)(quorem(b,S) + '0'));\r | |
752 | if (!b->x[0] && b->wds <= 1) {\r | |
753 | #ifdef SET_INEXACT\r | |
754 | inexact = 0;\r | |
755 | #endif\r | |
756 | goto ret;\r | |
757 | }\r | |
758 | if (i >= ilim)\r | |
759 | break;\r | |
760 | b = multadd(b, 10, 0);\r | |
761 | if (b == NULL)\r | |
762 | return NULL;\r | |
763 | }\r | |
764 | \r | |
765 | /* Round off last digit */\r | |
766 | \r | |
767 | #ifdef Honor_FLT_ROUNDS\r | |
768 | switch(rounding) {\r | |
769 | case 0: goto trimzeros;\r | |
770 | case 2: goto roundoff;\r | |
771 | }\r | |
772 | #endif\r | |
773 | b = lshift(b, 1);\r | |
774 | j = cmp(b, S);\r | |
775 | if (j > 0 || (j == 0 && dig & 1)) {\r | |
776 | roundoff:\r | |
777 | while(*--s == '9')\r | |
778 | if (s == s0) {\r | |
779 | k++;\r | |
780 | *s++ = '1';\r | |
781 | goto ret;\r | |
782 | }\r | |
783 | ++*s++;\r | |
784 | }\r | |
785 | else {\r | |
786 | #ifdef Honor_FLT_ROUNDS\r | |
787 | trimzeros:\r | |
788 | #endif\r | |
789 | while(*--s == '0');\r | |
790 | s++;\r | |
791 | }\r | |
792 | ret:\r | |
793 | Bfree(S);\r | |
794 | if (mhi) {\r | |
795 | if (mlo && mlo != mhi)\r | |
796 | Bfree(mlo);\r | |
797 | Bfree(mhi);\r | |
798 | }\r | |
799 | ret1:\r | |
800 | #ifdef SET_INEXACT\r | |
801 | if (inexact) {\r | |
802 | if (!oldinexact) {\r | |
803 | word0(d) = Exp_1 + (70 << Exp_shift);\r | |
804 | word1(d) = 0;\r | |
805 | dval(d) += 1.;\r | |
806 | }\r | |
807 | }\r | |
808 | else if (!oldinexact)\r | |
809 | clear_inexact();\r | |
810 | #endif\r | |
811 | Bfree(b);\r | |
812 | if (s == s0) { /* don't return empty string */\r | |
813 | *s++ = '0';\r | |
814 | k = 0;\r | |
815 | }\r | |
816 | *s = 0;\r | |
817 | *decpt = k + 1;\r | |
818 | if (rve)\r | |
819 | *rve = s;\r | |
820 | return s0;\r | |
821 | }\r |