]> git.proxmox.com Git - mirror_edk2.git/blame - StdLib/LibC/gdtoa/dtoa.c
StdLib: GCC 6 build fixes
[mirror_edk2.git] / StdLib / LibC / gdtoa / dtoa.c
CommitLineData
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
5The author of this software is David M. Gay.\r
6\r
7Copyright (C) 1998, 1999 by Lucent Technologies\r
8All Rights Reserved\r
9\r
10Permission to use, copy, modify, and distribute this software and\r
11its documentation for any purpose and without fee is hereby\r
12granted, provided that the above copyright notice appear in all\r
13copies and that both that the copyright notice and this\r
14permission notice and warranty disclaimer appear in supporting\r
15documentation, and that the name of Lucent or any of its entities\r
16not be used in advertising or publicity pertaining to\r
17distribution of the software without specific, written prior\r
18permission.\r
19\r
20LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,\r
21INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.\r
22IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY\r
23SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES\r
24WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER\r
25IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,\r
26ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF\r
27THIS 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
85dtoa\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