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