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