]> git.saurik.com Git - apple/libc.git/blob - gdtoa/gdtoa-gdtoa-fbsd.c
projects / apple / libc.git / blob
? search:


4a488870dd7f3fa507d9bda82e4da313b81cc5bd
[apple/libc.git] / gdtoa / gdtoa-gdtoa-fbsd.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 while(*--s == '0'){}
421 s++;
422 if (dval(d))
423 inex = STRTOG_Inexlo;
424 goto ret1;
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 while(*--s == '0'){}
485 s++;
486 }
487 break;
488 }
489 }
490 goto ret1;
491 }
492
493 m2 = b2;
494 m5 = b5;
495 mhi = mlo = 0;
496 if (leftright) {
497 if (mode < 2) {
498 i = nbits - bbits;
499 if (be - i++ < fpi->emin)
500 /* denormal */
501 i = be - fpi->emin + 1;
502 }
503 else {
504 j = ilim - 1;
505 if (m5 >= j)
506 m5 -= j;
507 else {
508 s5 += j -= m5;
509 b5 += j;
510 m5 = 0;
511 }
512 if ((i = ilim) < 0) {
513 m2 -= i;
514 i = 0;
515 }
516 }
517 b2 += i;
518 s2 += i;
519 mhi = i2b(1);
520 }
521 if (m2 > 0 && s2 > 0) {
522 i = m2 < s2 ? m2 : s2;
523 b2 -= i;
524 m2 -= i;
525 s2 -= i;
526 }
527 if (b5 > 0) {
528 if (leftright) {
529 if (m5 > 0) {
530 mhi = pow5mult(mhi, m5);
531 b1 = mult(mhi, b);
532 Bfree(b);
533 b = b1;
534 }
535 if ( (j = b5 - m5) !=0)
536 b = pow5mult(b, j);
537 }
538 else
539 b = pow5mult(b, b5);
540 }
541 S = i2b(1);
542 if (s5 > 0)
543 S = pow5mult(S, s5);
544
545 /* Check for special case that d is a normalized power of 2. */
546
547 spec_case = 0;
548 if (mode < 2) {
549 if (bbits == 1 && be0 > fpi->emin + 1) {
550 /* The special case */
551 b2++;
552 s2++;
553 spec_case = 1;
554 }
555 }
556
557 /* Arrange for convenient computation of quotients:
558 * shift left if necessary so divisor has 4 leading 0 bits.
559 *
560 * Perhaps we should just compute leading 28 bits of S once
561 * and for all and pass them and a shift to quorem, so it
562 * can do shifts and ors to compute the numerator for q.
563 */
564 #ifdef Pack_32
565 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) !=0)
566 i = 32 - i;
567 #else
568 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) !=0)
569 i = 16 - i;
570 #endif
571 if (i > 4) {
572 i -= 4;
573 b2 += i;
574 m2 += i;
575 s2 += i;
576 }
577 else if (i < 4) {
578 i += 28;
579 b2 += i;
580 m2 += i;
581 s2 += i;
582 }
583 if (b2 > 0)
584 b = lshift(b, b2);
585 if (s2 > 0)
586 S = lshift(S, s2);
587 if (k_check) {
588 if (cmp(b,S) < 0) {
589 k--;
590 b = multadd(b, 10, 0); /* we botched the k estimate */
591 if (leftright)
592 mhi = multadd(mhi, 10, 0);
593 ilim = ilim1;
594 }
595 }
596 if (ilim <= 0 && mode > 2) {
597 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
598 /* no digits, fcvt style */
599 no_digits:
600 k = -1 - ndigits;
601 inex = STRTOG_Inexlo;
602 goto ret;
603 }
604 one_digit:
605 inex = STRTOG_Inexhi;
606 *s++ = '1';
607 k++;
608 goto ret;
609 }
610 if (leftright) {
611 if (m2 > 0)
612 mhi = lshift(mhi, m2);
613
614 /* Compute mlo -- check for special case
615 * that d is a normalized power of 2.
616 */
617
618 mlo = mhi;
619 if (spec_case) {
620 mhi = Balloc(mhi->k);
621 Bcopy(mhi, mlo);
622 mhi = lshift(mhi, 1);
623 }
624
625 for(i = 1;;i++) {
626 dig = quorem(b,S) + '0';
627 /* Do we yet have the shortest decimal string
628 * that will round to d?
629 */
630 j = cmp(b, mlo);
631 delta = diff(S, mhi);
632 j1 = delta->sign ? 1 : cmp(b, delta);
633 Bfree(delta);
634 #ifndef ROUND_BIASED
635 if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
636 if (dig == '9')
637 goto round_9_up;
638 if (j <= 0) {
639 if (b->wds > 1 || b->x[0])
640 inex = STRTOG_Inexlo;
641 }
642 else {
643 dig++;
644 inex = STRTOG_Inexhi;
645 }
646 *s++ = dig;
647 goto ret;
648 }
649 #endif
650 if (j < 0 || j == 0 && !mode
651 #ifndef ROUND_BIASED
652 && !(bits[0] & 1)
653 #endif
654 ) {
655 if (rdir && (b->wds > 1 || b->x[0])) {
656 if (rdir == 2) {
657 inex = STRTOG_Inexlo;
658 goto accept;
659 }
660 while (cmp(S,mhi) > 0) {
661 *s++ = dig;
662 mhi1 = multadd(mhi, 10, 0);
663 if (mlo == mhi)
664 mlo = mhi1;
665 mhi = mhi1;
666 b = multadd(b, 10, 0);
667 dig = quorem(b,S) + '0';
668 }
669 if (dig++ == '9')
670 goto round_9_up;
671 inex = STRTOG_Inexhi;
672 goto accept;
673 }
674 if (j1 > 0) {
675 b = lshift(b, 1);
676 j1 = cmp(b, S);
677 if ((j1 > 0 || j1 == 0 && dig & 1)
678 && dig++ == '9')
679 goto round_9_up;
680 inex = STRTOG_Inexhi;
681 }
682 if (b->wds > 1 || b->x[0])
683 inex = STRTOG_Inexlo;
684 accept:
685 *s++ = dig;
686 goto ret;
687 }
688 if (j1 > 0 && rdir != 2) {
689 if (dig == '9') { /* possible if i == 1 */
690 round_9_up:
691 *s++ = '9';
692 inex = STRTOG_Inexhi;
693 goto roundoff;
694 }
695 inex = STRTOG_Inexhi;
696 *s++ = dig + 1;
697 goto ret;
698 }
699 *s++ = dig;
700 if (i == ilim)
701 break;
702 b = multadd(b, 10, 0);
703 if (mlo == mhi)
704 mlo = mhi = multadd(mhi, 10, 0);
705 else {
706 mlo = multadd(mlo, 10, 0);
707 mhi = multadd(mhi, 10, 0);
708 }
709 }
710 }
711 else
712 for(i = 1;; i++) {
713 *s++ = dig = quorem(b,S) + '0';
714 if (i >= ilim)
715 break;
716 b = multadd(b, 10, 0);
717 }
718
719 /* Round off last digit */
720
721 if (rdir) {
722 if (rdir == 2 || b->wds <= 1 && !b->x[0])
723 goto chopzeros;
724 goto roundoff;
725 }
726 b = lshift(b, 1);
727 j = cmp(b, S);
728 if (j > 0 || j == 0 && dig & 1) {
729 roundoff:
730 inex = STRTOG_Inexhi;
731 while(*--s == '9')
732 if (s == s0) {
733 k++;
734 *s++ = '1';
735 goto ret;
736 }
737 ++*s++;
738 }
739 else {
740 chopzeros:
741 if (b->wds > 1 || b->x[0])
742 inex = STRTOG_Inexlo;
743 while(*--s == '0'){}
744 s++;
745 }
746 ret:
747 Bfree(S);
748 if (mhi) {
749 if (mlo && mlo != mhi)
750 Bfree(mlo);
751 Bfree(mhi);
752 }
753 ret1:
754 Bfree(b);
755 *s = 0;
756 *decpt = k + 1;
757 if (rve)
758 *rve = s;
759 *kindp |= inex;
760 return s0;
761 }
mirror of libc