]> git.saurik.com Git - apple/libc.git/blob - gdtoa/FreeBSD/gdtoa-gdtoa.c
Libc-594.9.5.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 d, d2, ds, eps;
161 char *s, *s0;
162
163 #ifndef MULTIPLE_THREADS
164 if (dtoa_result) {
165 freedtoa(dtoa_result);
166 dtoa_result = 0;
167 }
168 #endif
169 inex = 0;
170 kind = *kindp &= ~STRTOG_Inexact;
171 switch(kind & STRTOG_Retmask) {
172 case STRTOG_Zero:
173 goto ret_zero;
174 case STRTOG_Normal:
175 case STRTOG_Denormal:
176 break;
177 case STRTOG_Infinite:
178 *decpt = -32768;
179 return nrv_alloc("Infinity", rve, 8);
180 case STRTOG_NaN:
181 *decpt = -32768;
182 return nrv_alloc("NaN", rve, 3);
183 default:
184 return 0;
185 }
186 b = bitstob(bits, nbits = fpi->nbits, &bbits);
187 be0 = be;
188 if ( (i = trailz(b)) !=0) {
189 rshift(b, i);
190 be += i;
191 bbits -= i;
192 }
193 if (!b->wds) {
194 Bfree(b);
195 ret_zero:
196 *decpt = 1;
197 return nrv_alloc("0", rve, 1);
198 }
199
200 dval(d) = b2d(b, &i);
201 i = be + bbits - 1;
202 word0(d) &= Frac_mask1;
203 word0(d) |= Exp_11;
204 #ifdef IBM
205 if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
206 dval(d) /= 1 << j;
207 #endif
208
209 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
210 * log10(x) = log(x) / log(10)
211 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
212 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
213 *
214 * This suggests computing an approximation k to log10(d) by
215 *
216 * k = (i - Bias)*0.301029995663981
217 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
218 *
219 * We want k to be too large rather than too small.
220 * The error in the first-order Taylor series approximation
221 * is in our favor, so we just round up the constant enough
222 * to compensate for any error in the multiplication of
223 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
224 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
225 * adding 1e-13 to the constant term more than suffices.
226 * Hence we adjust the constant term to 0.1760912590558.
227 * (We could get a more accurate k by invoking log10,
228 * but this is probably not worthwhile.)
229 */
230 #ifdef IBM
231 i <<= 2;
232 i += j;
233 #endif
234 ds = (dval(d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
235
236 /* correct assumption about exponent range */
237 if ((j = i) < 0)
238 j = -j;
239 if ((j -= 1077) > 0)
240 ds += j * 7e-17;
241
242 k = (int)ds;
243 if (ds < 0. && ds != k)
244 k--; /* want k = floor(ds) */
245 k_check = 1;
246 #ifdef IBM
247 j = be + bbits - 1;
248 if ( (j1 = j & 3) !=0)
249 dval(d) *= 1 << j1;
250 word0(d) += j << Exp_shift - 2 & Exp_mask;
251 #else
252 word0(d) += (be + bbits - 1) << Exp_shift;
253 #endif
254 if (k >= 0 && k <= Ten_pmax) {
255 if (dval(d) < tens[k])
256 k--;
257 k_check = 0;
258 }
259 j = bbits - i - 1;
260 if (j >= 0) {
261 b2 = 0;
262 s2 = j;
263 }
264 else {
265 b2 = -j;
266 s2 = 0;
267 }
268 if (k >= 0) {
269 b5 = 0;
270 s5 = k;
271 s2 += k;
272 }
273 else {
274 b2 -= k;
275 b5 = -k;
276 s5 = 0;
277 }
278 if (mode < 0 || mode > 9)
279 mode = 0;
280 try_quick = 1;
281 if (mode > 5) {
282 mode -= 4;
283 try_quick = 0;
284 }
285 leftright = 1;
286 switch(mode) {
287 case 0:
288 case 1:
289 ilim = ilim1 = -1;
290 i = (int)(nbits * .30103) + 3;
291 ndigits = 0;
292 break;
293 case 2:
294 leftright = 0;
295 /* no break */
296 case 4:
297 if (ndigits <= 0)
298 ndigits = 1;
299 ilim = ilim1 = i = ndigits;
300 break;
301 case 3:
302 leftright = 0;
303 /* no break */
304 case 5:
305 i = ndigits + k + 1;
306 ilim = i;
307 ilim1 = i - 1;
308 if (i <= 0)
309 i = 1;
310 }
311 s = s0 = rv_alloc(i);
312
313 if ( (rdir = fpi->rounding - 1) !=0) {
314 if (rdir < 0)
315 rdir = 2;
316 if (kind & STRTOG_Neg)
317 rdir = 3 - rdir;
318 }
319
320 /* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
321
322 if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir
323 #ifndef IMPRECISE_INEXACT
324 && k == 0
325 #endif
326 ) {
327
328 /* Try to get by with floating-point arithmetic. */
329
330 i = 0;
331 d2 = dval(d);
332 #ifdef IBM
333 if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
334 dval(d) /= 1 << j;
335 #endif
336 k0 = k;
337 ilim0 = ilim;
338 ieps = 2; /* conservative */
339 if (k > 0) {
340 ds = tens[k&0xf];
341 j = k >> 4;
342 if (j & Bletch) {
343 /* prevent overflows */
344 j &= Bletch - 1;
345 dval(d) /= bigtens[n_bigtens-1];
346 ieps++;
347 }
348 for(; j; j >>= 1, i++)
349 if (j & 1) {
350 ieps++;
351 ds *= bigtens[i];
352 }
353 }
354 else {
355 ds = 1.;
356 if ( (j1 = -k) !=0) {
357 dval(d) *= tens[j1 & 0xf];
358 for(j = j1 >> 4; j; j >>= 1, i++)
359 if (j & 1) {
360 ieps++;
361 dval(d) *= bigtens[i];
362 }
363 }
364 }
365 if (k_check && dval(d) < 1. && ilim > 0) {
366 if (ilim1 <= 0)
367 goto fast_failed;
368 ilim = ilim1;
369 k--;
370 dval(d) *= 10.;
371 ieps++;
372 }
373 dval(eps) = ieps*dval(d) + 7.;
374 word0(eps) -= (P-1)*Exp_msk1;
375 if (ilim == 0) {
376 S = mhi = 0;
377 dval(d) -= 5.;
378 if (dval(d) > dval(eps))
379 goto one_digit;
380 if (dval(d) < -dval(eps))
381 goto no_digits;
382 goto fast_failed;
383 }
384 #ifndef No_leftright
385 if (leftright) {
386 /* Use Steele & White method of only
387 * generating digits needed.
388 */
389 dval(eps) = ds*0.5/tens[ilim-1] - dval(eps);
390 for(i = 0;;) {
391 L = (Long)(dval(d)/ds);
392 dval(d) -= L*ds;
393 *s++ = '0' + (int)L;
394 if (dval(d) < dval(eps)) {
395 if (dval(d))
396 inex = STRTOG_Inexlo;
397 goto ret1;
398 }
399 if (ds - dval(d) < dval(eps))
400 goto bump_up;
401 if (++i >= ilim)
402 break;
403 dval(eps) *= 10.;
404 dval(d) *= 10.;
405 }
406 }
407 else {
408 #endif
409 /* Generate ilim digits, then fix them up. */
410 dval(eps) *= tens[ilim-1];
411 for(i = 1;; i++, dval(d) *= 10.) {
412 if ( (L = (Long)(dval(d)/ds)) !=0)
413 dval(d) -= L*ds;
414 *s++ = '0' + (int)L;
415 if (i == ilim) {
416 ds *= 0.5;
417 if (dval(d) > ds + dval(eps))
418 goto bump_up;
419 else if (dval(d) < ds - dval(eps)) {
420 if (dval(d))
421 inex = STRTOG_Inexlo;
422 goto clear_trailing0;
423 }
424 break;
425 }
426 }
427 #ifndef No_leftright
428 }
429 #endif
430 fast_failed:
431 s = s0;
432 dval(d) = d2;
433 k = k0;
434 ilim = ilim0;
435 }
436
437 /* Do we have a "small" integer? */
438
439 if (be >= 0 && k <= Int_max) {
440 /* Yes. */
441 ds = tens[k];
442 if (ndigits < 0 && ilim <= 0) {
443 S = mhi = 0;
444 if (ilim < 0 || dval(d) <= 5*ds)
445 goto no_digits;
446 goto one_digit;
447 }
448 for(i = 1;; i++, dval(d) *= 10.) {
449 L = dval(d) / ds;
450 dval(d) -= L*ds;
451 #ifdef Check_FLT_ROUNDS
452 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
453 if (dval(d) < 0) {
454 L--;
455 dval(d) += ds;
456 }
457 #endif
458 *s++ = '0' + (int)L;
459 if (dval(d) == 0.)
460 break;
461 if (i == ilim) {
462 if (rdir) {
463 if (rdir == 1)
464 goto bump_up;
465 inex = STRTOG_Inexlo;
466 goto ret1;
467 }
468 dval(d) += dval(d);
469 if (dval(d) > ds || dval(d) == ds && L & 1) {
470 bump_up:
471 inex = STRTOG_Inexhi;
472 while(*--s == '9')
473 if (s == s0) {
474 k++;
475 *s = '0';
476 break;
477 }
478 ++*s++;
479 }
480 else {
481 inex = STRTOG_Inexlo;
482 clear_trailing0:
483 while(*--s == '0'){}
484 ++s;
485 }
486 break;
487 }
488 }
489 goto ret1;
490 }
491
492 m2 = b2;
493 m5 = b5;
494 mhi = mlo = 0;
495 if (leftright) {
496 if (mode < 2) {
497 i = nbits - bbits;
498 if (be - i++ < fpi->emin)
499 /* denormal */
500 i = be - fpi->emin + 1;
501 }
502 else {
503 j = ilim - 1;
504 if (m5 >= j)
505 m5 -= j;
506 else {
507 s5 += j -= m5;
508 b5 += j;
509 m5 = 0;
510 }
511 if ((i = ilim) < 0) {
512 m2 -= i;
513 i = 0;
514 }
515 }
516 b2 += i;
517 s2 += i;
518 mhi = i2b(1);
519 }
520 if (m2 > 0 && s2 > 0) {
521 i = m2 < s2 ? m2 : s2;
522 b2 -= i;
523 m2 -= i;
524 s2 -= i;
525 }
526 if (b5 > 0) {
527 if (leftright) {
528 if (m5 > 0) {
529 mhi = pow5mult(mhi, m5);
530 b1 = mult(mhi, b);
531 Bfree(b);
532 b = b1;
533 }
534 if ( (j = b5 - m5) !=0)
535 b = pow5mult(b, j);
536 }
537 else
538 b = pow5mult(b, b5);
539 }
540 S = i2b(1);
541 if (s5 > 0)
542 S = pow5mult(S, s5);
543
544 /* Check for special case that d is a normalized power of 2. */
545
546 spec_case = 0;
547 if (mode < 2) {
548 if (bbits == 1 && be0 > fpi->emin + 1) {
549 /* The special case */
550 b2++;
551 s2++;
552 spec_case = 1;
553 }
554 }
555
556 /* Arrange for convenient computation of quotients:
557 * shift left if necessary so divisor has 4 leading 0 bits.
558 *
559 * Perhaps we should just compute leading 28 bits of S once
560 * and for all and pass them and a shift to quorem, so it
561 * can do shifts and ors to compute the numerator for q.
562 */
563 #ifdef Pack_32
564 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) !=0)
565 i = 32 - i;
566 #else
567 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) !=0)
568 i = 16 - i;
569 #endif
570 if (i > 4) {
571 i -= 4;
572 b2 += i;
573 m2 += i;
574 s2 += i;
575 }
576 else if (i < 4) {
577 i += 28;
578 b2 += i;
579 m2 += i;
580 s2 += i;
581 }
582 if (b2 > 0)
583 b = lshift(b, b2);
584 if (s2 > 0)
585 S = lshift(S, s2);
586 if (k_check) {
587 if (cmp(b,S) < 0) {
588 k--;
589 b = multadd(b, 10, 0); /* we botched the k estimate */
590 if (leftright)
591 mhi = multadd(mhi, 10, 0);
592 ilim = ilim1;
593 }
594 }
595 if (ilim <= 0 && mode > 2) {
596 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
597 /* no digits, fcvt style */
598 no_digits:
599 k = -1 - ndigits;
600 inex = STRTOG_Inexlo;
601 goto ret;
602 }
603 one_digit:
604 inex = STRTOG_Inexhi;
605 *s++ = '1';
606 k++;
607 goto ret;
608 }
609 if (leftright) {
610 if (m2 > 0)
611 mhi = lshift(mhi, m2);
612
613 /* Compute mlo -- check for special case
614 * that d is a normalized power of 2.
615 */
616
617 mlo = mhi;
618 if (spec_case) {
619 mhi = Balloc(mhi->k);
620 Bcopy(mhi, mlo);
621 mhi = lshift(mhi, 1);
622 }
623
624 for(i = 1;;i++) {
625 dig = quorem(b,S) + '0';
626 /* Do we yet have the shortest decimal string
627 * that will round to d?
628 */
629 j = cmp(b, mlo);
630 delta = diff(S, mhi);
631 j1 = delta->sign ? 1 : cmp(b, delta);
632 Bfree(delta);
633 #ifndef ROUND_BIASED
634 if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
635 if (dig == '9')
636 goto round_9_up;
637 if (j <= 0) {
638 if (b->wds > 1 || b->x[0])
639 inex = STRTOG_Inexlo;
640 }
641 else {
642 dig++;
643 inex = STRTOG_Inexhi;
644 }
645 *s++ = dig;
646 goto ret;
647 }
648 #endif
649 if (j < 0 || j == 0 && !mode
650 #ifndef ROUND_BIASED
651 && !(bits[0] & 1)
652 #endif
653 ) {
654 if (rdir && (b->wds > 1 || b->x[0])) {
655 if (rdir == 2) {
656 inex = STRTOG_Inexlo;
657 goto accept;
658 }
659 while (cmp(S,mhi) > 0) {
660 *s++ = dig;
661 mhi1 = multadd(mhi, 10, 0);
662 if (mlo == mhi)
663 mlo = mhi1;
664 mhi = mhi1;
665 b = multadd(b, 10, 0);
666 dig = quorem(b,S) + '0';
667 }
668 if (dig++ == '9')
669 goto round_9_up;
670 inex = STRTOG_Inexhi;
671 goto accept;
672 }
673 if (j1 > 0) {
674 b = lshift(b, 1);
675 j1 = cmp(b, S);
676 if ((j1 > 0 || j1 == 0 && dig & 1)
677 && dig++ == '9')
678 goto round_9_up;
679 inex = STRTOG_Inexhi;
680 }
681 if (b->wds > 1 || b->x[0])
682 inex = STRTOG_Inexlo;
683 accept:
684 *s++ = dig;
685 goto ret;
686 }
687 if (j1 > 0 && rdir != 2) {
688 if (dig == '9') { /* possible if i == 1 */
689 round_9_up:
690 *s++ = '9';
691 inex = STRTOG_Inexhi;
692 goto roundoff;
693 }
694 inex = STRTOG_Inexhi;
695 *s++ = dig + 1;
696 goto ret;
697 }
698 *s++ = dig;
699 if (i == ilim)
700 break;
701 b = multadd(b, 10, 0);
702 if (mlo == mhi)
703 mlo = mhi = multadd(mhi, 10, 0);
704 else {
705 mlo = multadd(mlo, 10, 0);
706 mhi = multadd(mhi, 10, 0);
707 }
708 }
709 }
710 else
711 for(i = 1;; i++) {
712 *s++ = dig = quorem(b,S) + '0';
713 if (i >= ilim)
714 break;
715 b = multadd(b, 10, 0);
716 }
717
718 /* Round off last digit */
719
720 if (rdir) {
721 if (rdir == 2 || b->wds <= 1 && !b->x[0])
722 goto chopzeros;
723 goto roundoff;
724 }
725 b = lshift(b, 1);
726 j = cmp(b, S);
727 if (j > 0 || j == 0 && dig & 1) {
728 roundoff:
729 inex = STRTOG_Inexhi;
730 while(*--s == '9')
731 if (s == s0) {
732 k++;
733 *s++ = '1';
734 goto ret;
735 }
736 ++*s++;
737 }
738 else {
739 chopzeros:
740 if (b->wds > 1 || b->x[0])
741 inex = STRTOG_Inexlo;
742 while(*--s == '0'){}
743 ++s;
744 }
745 ret:
746 Bfree(S);
747 if (mhi) {
748 if (mlo && mlo != mhi)
749 Bfree(mlo);
750 Bfree(mhi);
751 }
752 ret1:
753 Bfree(b);
754 *s = 0;
755 *decpt = k + 1;
756 if (rve)
757 *rve = s;
758 *kindp |= inex;
759 return s0;
760 }