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