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