drm/i915: Update to Linux 4.6
[dragonfly.git] / contrib / gdtoa / gdtoa.c
blobe15bb2ad862ed0435961dd4f52cc4562e59f26b2
1 /****************************************************************
3 The author of this software is David M. Gay.
5 Copyright (C) 1998, 1999 by Lucent Technologies
6 All Rights Reserved
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.
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.
27 ****************************************************************/
29 /* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
32 #include "gdtoaimp.h"
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
41 int i, k;
42 Bigint *b;
43 ULong *be, *x, *x0;
45 i = ULbits;
46 k = 0;
47 while(i < nbits) {
48 i <<= 1;
49 k++;
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;
71 b->wds = i + 1;
72 *bbits = i*ULbits + 32 - hi0bits(b->x[i]);
73 ret:
74 return b;
77 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
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].
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.
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
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 be = exponent: value = (integer represented by bits) * (2 to the power of be).
128 mode:
129 0 ==> shortest string that yields d when read in
130 and rounded to nearest.
131 1 ==> like 0, but with Steele & White stopping rule;
132 e.g. with IEEE P754 arithmetic , mode 0 gives
133 1e23 whereas mode 1 gives 9.999999999999999e22.
134 2 ==> max(1,ndigits) significant digits. This gives a
135 return value similar to that of ecvt, except
136 that trailing zeros are suppressed.
137 3 ==> through ndigits past the decimal point. This
138 gives a return value similar to that from fcvt,
139 except that trailing zeros are suppressed, and
140 ndigits can be negative.
141 4-9 should give the same return values as 2-3, i.e.,
142 4 <= mode <= 9 ==> same return as mode
143 2 + (mode & 1). These modes are mainly for
144 debugging; often they run slower but sometimes
145 faster than modes 2-3.
146 4,5,8,9 ==> left-to-right digit generation.
147 6-9 ==> don't try fast floating-point estimate
148 (if applicable).
150 Values of mode other than 0-9 are treated as mode 0.
152 Sufficient space is allocated to the return value
153 to hold the suppressed trailing zeros.
156 int bbits, b2, b5, be0, dig, i, ieps, ilim, ilim0, ilim1, inex;
157 int j, j1, k, k0, k_check, kind, leftright, m2, m5, nbits;
158 int rdir, s2, s5, spec_case, try_quick;
159 Long L;
160 Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
161 double d2, ds;
162 char *s, *s0;
163 U d, eps;
165 #ifndef MULTIPLE_THREADS
166 if (dtoa_result) {
167 freedtoa(dtoa_result);
168 dtoa_result = 0;
170 #endif
171 inex = 0;
172 kind = *kindp &= ~STRTOG_Inexact;
173 switch(kind & STRTOG_Retmask) {
174 case STRTOG_Zero:
175 goto ret_zero;
176 case STRTOG_Normal:
177 case STRTOG_Denormal:
178 break;
179 case STRTOG_Infinite:
180 *decpt = -32768;
181 return nrv_alloc("Infinity", rve, 8);
182 case STRTOG_NaN:
183 *decpt = -32768;
184 return nrv_alloc("NaN", rve, 3);
185 default:
186 return 0;
188 b = bitstob(bits, nbits = fpi->nbits, &bbits);
189 be0 = be;
190 if ( (i = trailz(b)) !=0) {
191 rshift(b, i);
192 be += i;
193 bbits -= i;
195 if (!b->wds) {
196 Bfree(b);
197 ret_zero:
198 *decpt = 1;
199 return nrv_alloc("0", rve, 1);
202 dval(&d) = b2d(b, &i);
203 i = be + bbits - 1;
204 word0(&d) &= Frac_mask1;
205 word0(&d) |= Exp_11;
206 #ifdef IBM
207 if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
208 dval(&d) /= 1 << j;
209 #endif
211 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
212 * log10(x) = log(x) / log(10)
213 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
214 * log10(&d) = (i-Bias)*log(2)/log(10) + log10(d2)
216 * This suggests computing an approximation k to log10(&d) by
218 * k = (i - Bias)*0.301029995663981
219 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
221 * We want k to be too large rather than too small.
222 * The error in the first-order Taylor series approximation
223 * is in our favor, so we just round up the constant enough
224 * to compensate for any error in the multiplication of
225 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
226 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
227 * adding 1e-13 to the constant term more than suffices.
228 * Hence we adjust the constant term to 0.1760912590558.
229 * (We could get a more accurate k by invoking log10,
230 * but this is probably not worthwhile.)
232 #ifdef IBM
233 i <<= 2;
234 i += j;
235 #endif
236 ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
238 /* correct assumption about exponent range */
239 if ((j = i) < 0)
240 j = -j;
241 if ((j -= 1077) > 0)
242 ds += j * 7e-17;
244 k = (int)ds;
245 if (ds < 0. && ds != k)
246 k--; /* want k = floor(ds) */
247 k_check = 1;
248 #ifdef IBM
249 j = be + bbits - 1;
250 if ( (j1 = j & 3) !=0)
251 dval(&d) *= 1 << j1;
252 word0(&d) += j << Exp_shift - 2 & Exp_mask;
253 #else
254 word0(&d) += (be + bbits - 1) << Exp_shift;
255 #endif
256 if (k >= 0 && k <= Ten_pmax) {
257 if (dval(&d) < tens[k])
258 k--;
259 k_check = 0;
261 j = bbits - i - 1;
262 if (j >= 0) {
263 b2 = 0;
264 s2 = j;
266 else {
267 b2 = -j;
268 s2 = 0;
270 if (k >= 0) {
271 b5 = 0;
272 s5 = k;
273 s2 += k;
275 else {
276 b2 -= k;
277 b5 = -k;
278 s5 = 0;
280 if (mode < 0 || mode > 9)
281 mode = 0;
282 try_quick = 1;
283 if (mode > 5) {
284 mode -= 4;
285 try_quick = 0;
287 else if (i >= -4 - Emin || i < Emin)
288 try_quick = 0;
289 leftright = 1;
290 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
291 /* silence erroneous "gcc -Wall" warning. */
292 switch(mode) {
293 case 0:
294 case 1:
295 i = (int)(nbits * .30103) + 3;
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;
316 s = s0 = rv_alloc(i);
318 if ( (rdir = fpi->rounding - 1) !=0) {
319 if (rdir < 0)
320 rdir = 2;
321 if (kind & STRTOG_Neg)
322 rdir = 3 - rdir;
325 /* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
327 if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir
328 #ifndef IMPRECISE_INEXACT
329 && k == 0
330 #endif
333 /* Try to get by with floating-point arithmetic. */
335 i = 0;
336 d2 = dval(&d);
337 #ifdef IBM
338 if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
339 dval(&d) /= 1 << j;
340 #endif
341 k0 = k;
342 ilim0 = ilim;
343 ieps = 2; /* conservative */
344 if (k > 0) {
345 ds = tens[k&0xf];
346 j = k >> 4;
347 if (j & Bletch) {
348 /* prevent overflows */
349 j &= Bletch - 1;
350 dval(&d) /= bigtens[n_bigtens-1];
351 ieps++;
353 for(; j; j >>= 1, i++)
354 if (j & 1) {
355 ieps++;
356 ds *= bigtens[i];
359 else {
360 ds = 1.;
361 if ( (j1 = -k) !=0) {
362 dval(&d) *= tens[j1 & 0xf];
363 for(j = j1 >> 4; j; j >>= 1, i++)
364 if (j & 1) {
365 ieps++;
366 dval(&d) *= bigtens[i];
370 if (k_check && dval(&d) < 1. && ilim > 0) {
371 if (ilim1 <= 0)
372 goto fast_failed;
373 ilim = ilim1;
374 k--;
375 dval(&d) *= 10.;
376 ieps++;
378 dval(&eps) = ieps*dval(&d) + 7.;
379 word0(&eps) -= (P-1)*Exp_msk1;
380 if (ilim == 0) {
381 S = mhi = 0;
382 dval(&d) -= 5.;
383 if (dval(&d) > dval(&eps))
384 goto one_digit;
385 if (dval(&d) < -dval(&eps))
386 goto no_digits;
387 goto fast_failed;
389 #ifndef No_leftright
390 if (leftright) {
391 /* Use Steele & White method of only
392 * generating digits needed.
394 dval(&eps) = ds*0.5/tens[ilim-1] - dval(&eps);
395 for(i = 0;;) {
396 L = (Long)(dval(&d)/ds);
397 dval(&d) -= L*ds;
398 *s++ = '0' + (int)L;
399 if (dval(&d) < dval(&eps)) {
400 if (dval(&d))
401 inex = STRTOG_Inexlo;
402 goto ret1;
404 if (ds - dval(&d) < dval(&eps))
405 goto bump_up;
406 if (++i >= ilim)
407 break;
408 dval(&eps) *= 10.;
409 dval(&d) *= 10.;
412 else {
413 #endif
414 /* Generate ilim digits, then fix them up. */
415 dval(&eps) *= tens[ilim-1];
416 for(i = 1;; i++, dval(&d) *= 10.) {
417 if ( (L = (Long)(dval(&d)/ds)) !=0)
418 dval(&d) -= L*ds;
419 *s++ = '0' + (int)L;
420 if (i == ilim) {
421 ds *= 0.5;
422 if (dval(&d) > ds + dval(&eps))
423 goto bump_up;
424 else if (dval(&d) < ds - dval(&eps)) {
425 if (dval(&d))
426 inex = STRTOG_Inexlo;
427 goto clear_trailing0;
429 break;
432 #ifndef No_leftright
434 #endif
435 fast_failed:
436 s = s0;
437 dval(&d) = d2;
438 k = k0;
439 ilim = ilim0;
442 /* Do we have a "small" integer? */
444 if (be >= 0 && k <= Int_max) {
445 /* Yes. */
446 ds = tens[k];
447 if (ndigits < 0 && ilim <= 0) {
448 S = mhi = 0;
449 if (ilim < 0 || dval(&d) <= 5*ds)
450 goto no_digits;
451 goto one_digit;
453 for(i = 1;; i++, dval(&d) *= 10.) {
454 L = dval(&d) / ds;
455 dval(&d) -= L*ds;
456 #ifdef Check_FLT_ROUNDS
457 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
458 if (dval(&d) < 0) {
459 L--;
460 dval(&d) += ds;
462 #endif
463 *s++ = '0' + (int)L;
464 if (dval(&d) == 0.)
465 break;
466 if (i == ilim) {
467 if (rdir) {
468 if (rdir == 1)
469 goto bump_up;
470 inex = STRTOG_Inexlo;
471 goto ret1;
473 dval(&d) += dval(&d);
474 #ifdef ROUND_BIASED
475 if (dval(&d) >= ds)
476 #else
477 if (dval(&d) > ds || (dval(&d) == ds && L & 1))
478 #endif
480 bump_up:
481 inex = STRTOG_Inexhi;
482 while(*--s == '9')
483 if (s == s0) {
484 k++;
485 *s = '0';
486 break;
488 ++*s++;
490 else {
491 inex = STRTOG_Inexlo;
492 clear_trailing0:
493 while(*--s == '0'){}
494 ++s;
496 break;
499 goto ret1;
502 m2 = b2;
503 m5 = b5;
504 mhi = mlo = 0;
505 if (leftright) {
506 i = nbits - bbits;
507 if (be - i++ < fpi->emin && mode != 3 && mode != 5) {
508 /* denormal */
509 i = be - fpi->emin + 1;
510 if (mode >= 2 && ilim > 0 && ilim < i)
511 goto small_ilim;
513 else if (mode >= 2) {
514 small_ilim:
515 j = ilim - 1;
516 if (m5 >= j)
517 m5 -= j;
518 else {
519 s5 += j -= m5;
520 b5 += j;
521 m5 = 0;
523 if ((i = ilim) < 0) {
524 m2 -= i;
525 i = 0;
528 b2 += i;
529 s2 += i;
530 mhi = i2b(1);
532 if (m2 > 0 && s2 > 0) {
533 i = m2 < s2 ? m2 : s2;
534 b2 -= i;
535 m2 -= i;
536 s2 -= i;
538 if (b5 > 0) {
539 if (leftright) {
540 if (m5 > 0) {
541 mhi = pow5mult(mhi, m5);
542 b1 = mult(mhi, b);
543 Bfree(b);
544 b = b1;
546 if ( (j = b5 - m5) !=0)
547 b = pow5mult(b, j);
549 else
550 b = pow5mult(b, b5);
552 S = i2b(1);
553 if (s5 > 0)
554 S = pow5mult(S, s5);
556 /* Check for special case that d is a normalized power of 2. */
558 spec_case = 0;
559 if (mode < 2) {
560 if (bbits == 1 && be0 > fpi->emin + 1) {
561 /* The special case */
562 b2++;
563 s2++;
564 spec_case = 1;
568 /* Arrange for convenient computation of quotients:
569 * shift left if necessary so divisor has 4 leading 0 bits.
571 * Perhaps we should just compute leading 28 bits of S once
572 * and for all and pass them and a shift to quorem, so it
573 * can do shifts and ors to compute the numerator for q.
575 i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;
576 m2 += i;
577 if ((b2 += i) > 0)
578 b = lshift(b, b2);
579 if ((s2 += i) > 0)
580 S = lshift(S, s2);
581 if (k_check) {
582 if (cmp(b,S) < 0) {
583 k--;
584 b = multadd(b, 10, 0); /* we botched the k estimate */
585 if (leftright)
586 mhi = multadd(mhi, 10, 0);
587 ilim = ilim1;
590 if (ilim <= 0 && mode > 2) {
591 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
592 /* no digits, fcvt style */
593 no_digits:
594 k = -1 - ndigits;
595 inex = STRTOG_Inexlo;
596 goto ret;
598 one_digit:
599 inex = STRTOG_Inexhi;
600 *s++ = '1';
601 k++;
602 goto ret;
604 if (leftright) {
605 if (m2 > 0)
606 mhi = lshift(mhi, m2);
608 /* Compute mlo -- check for special case
609 * that d is a normalized power of 2.
612 mlo = mhi;
613 if (spec_case) {
614 mhi = Balloc(mhi->k);
615 Bcopy(mhi, mlo);
616 mhi = lshift(mhi, 1);
619 for(i = 1;;i++) {
620 dig = quorem(b,S) + '0';
621 /* Do we yet have the shortest decimal string
622 * that will round to d?
624 j = cmp(b, mlo);
625 delta = diff(S, mhi);
626 j1 = delta->sign ? 1 : cmp(b, delta);
627 Bfree(delta);
628 #ifndef ROUND_BIASED
629 if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
630 if (dig == '9')
631 goto round_9_up;
632 if (j <= 0) {
633 if (b->wds > 1 || b->x[0])
634 inex = STRTOG_Inexlo;
636 else {
637 dig++;
638 inex = STRTOG_Inexhi;
640 *s++ = dig;
641 goto ret;
643 #endif
644 if (j < 0 || (j == 0 && !mode
645 #ifndef ROUND_BIASED
646 && !(bits[0] & 1)
647 #endif
648 )) {
649 if (rdir && (b->wds > 1 || b->x[0])) {
650 if (rdir == 2) {
651 inex = STRTOG_Inexlo;
652 goto accept;
654 while (cmp(S,mhi) > 0) {
655 *s++ = dig;
656 mhi1 = multadd(mhi, 10, 0);
657 if (mlo == mhi)
658 mlo = mhi1;
659 mhi = mhi1;
660 b = multadd(b, 10, 0);
661 dig = quorem(b,S) + '0';
663 if (dig++ == '9')
664 goto round_9_up;
665 inex = STRTOG_Inexhi;
666 goto accept;
668 if (j1 > 0) {
669 b = lshift(b, 1);
670 j1 = cmp(b, S);
671 #ifdef ROUND_BIASED
672 if (j1 >= 0 /*)*/
673 #else
674 if ((j1 > 0 || (j1 == 0 && dig & 1))
675 #endif
676 && dig++ == '9')
677 goto round_9_up;
678 inex = STRTOG_Inexhi;
680 if (b->wds > 1 || b->x[0])
681 inex = STRTOG_Inexlo;
682 accept:
683 *s++ = dig;
684 goto ret;
686 if (j1 > 0 && rdir != 2) {
687 if (dig == '9') { /* possible if i == 1 */
688 round_9_up:
689 *s++ = '9';
690 inex = STRTOG_Inexhi;
691 goto roundoff;
693 inex = STRTOG_Inexhi;
694 *s++ = dig + 1;
695 goto ret;
697 *s++ = dig;
698 if (i == ilim)
699 break;
700 b = multadd(b, 10, 0);
701 if (mlo == mhi)
702 mlo = mhi = multadd(mhi, 10, 0);
703 else {
704 mlo = multadd(mlo, 10, 0);
705 mhi = multadd(mhi, 10, 0);
709 else
710 for(i = 1;; i++) {
711 *s++ = dig = quorem(b,S) + '0';
712 if (i >= ilim)
713 break;
714 b = multadd(b, 10, 0);
717 /* Round off last digit */
719 if (rdir) {
720 if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
721 goto chopzeros;
722 goto roundoff;
724 b = lshift(b, 1);
725 j = cmp(b, S);
726 #ifdef ROUND_BIASED
727 if (j >= 0)
728 #else
729 if (j > 0 || (j == 0 && dig & 1))
730 #endif
732 roundoff:
733 inex = STRTOG_Inexhi;
734 while(*--s == '9')
735 if (s == s0) {
736 k++;
737 *s++ = '1';
738 goto ret;
740 ++*s++;
742 else {
743 chopzeros:
744 if (b->wds > 1 || b->x[0])
745 inex = STRTOG_Inexlo;
746 while(*--s == '0'){}
747 ++s;
749 ret:
750 Bfree(S);
751 if (mhi) {
752 if (mlo && mlo != mhi)
753 Bfree(mlo);
754 Bfree(mhi);
756 ret1:
757 Bfree(b);
758 *s = 0;
759 *decpt = k + 1;
760 if (rve)
761 *rve = s;
762 *kindp |= inex;
763 return s0;