![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to floatconv.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Apple XNU] / GNUtools / libg++ / libio|
Return to floatconv.c CVS log | Up to [Apple XNU] / GNUtools / libg++ / libio |
1.1 root 1: /*
2: Copyright (C) 1993 Free Software Foundation
3:
4: This file is part of the GNU IO Library. This library is free
5: software; you can redistribute it and/or modify it under the
6: terms of the GNU General Public License as published by the
7: Free Software Foundation; either version 2, or (at your option)
8: any later version.
9:
10: This library is distributed in the hope that it will be useful,
11: but WITHOUT ANY WARRANTY; without even the implied warranty of
12: MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13: GNU General Public License for more details.
14:
15: You should have received a copy of the GNU General Public License
16: along with GNU CC; see the file COPYING. If not, write to
17: the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
18:
19: As a special exception, if you link this library with files
20: compiled with a GNU compiler to produce an executable, this does not cause
21: the resulting executable to be covered by the GNU General Public License.
22: This exception does not however invalidate any other reasons why
23: the executable file might be covered by the GNU General Public License. */
24:
25: #include <libioP.h>
26: #ifdef USE_DTOA
27: /****************************************************************
28: *
29: * The author of this software is David M. Gay.
30: *
31: * Copyright (c) 1991 by AT&T.
32: *
33: * Permission to use, copy, modify, and distribute this software for any
34: * purpose without fee is hereby granted, provided that this entire notice
35: * is included in all copies of any software which is or includes a copy
36: * or modification of this software and in all copies of the supporting
37: * documentation for such software.
38: *
39: * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
40: * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
41: * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
42: * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
43: *
44: ***************************************************************/
45:
46: /* Some cleaning up by Per Bothner, [email protected], 1992, 1993.
47: Re-written to not need static variables
48: (except result, result_k, HIWORD, LOWORD). */
49:
50: /* Please send bug reports to
51: David M. Gay
52: AT&T Bell Laboratories, Room 2C-463
53: 600 Mountain Avenue
54: Murray Hill, NJ 07974-2070
55: U.S.A.
56: [email protected] or research!dmg
57: */
58:
59: /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
60: *
61: * This strtod returns a nearest machine number to the input decimal
62: * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
63: * broken by the IEEE round-even rule. Otherwise ties are broken by
64: * biased rounding (add half and chop).
65: *
66: * Inspired loosely by William D. Clinger's paper "How to Read Floating
67: * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
68: *
69: * Modifications:
70: *
71: * 1. We only require IEEE, IBM, or VAX double-precision
72: * arithmetic (not IEEE double-extended).
73: * 2. We get by with floating-point arithmetic in a case that
74: * Clinger missed -- when we're computing d * 10^n
75: * for a small integer d and the integer n is not too
76: * much larger than 22 (the maximum integer k for which
77: * we can represent 10^k exactly), we may be able to
78: * compute (d*10^k) * 10^(e-k) with just one roundoff.
79: * 3. Rather than a bit-at-a-time adjustment of the binary
80: * result in the hard case, we use floating-point
81: * arithmetic to determine the adjustment to within
82: * one bit; only in really hard cases do we need to
83: * compute a second residual.
84: * 4. Because of 3., we don't need a large table of powers of 10
85: * for ten-to-e (just some small tables, e.g. of 10^k
86: * for 0 <= k <= 22).
87: */
88:
89: /*
90: * #define IEEE_8087 for IEEE-arithmetic machines where the least
91: * significant byte has the lowest address.
92: * #define IEEE_MC68k for IEEE-arithmetic machines where the most
93: * significant byte has the lowest address.
94: * #define Sudden_Underflow for IEEE-format machines without gradual
95: * underflow (i.e., that flush to zero on underflow).
96: * #define IBM for IBM mainframe-style floating-point arithmetic.
97: * #define VAX for VAX-style floating-point arithmetic.
98: * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
99: * #define No_leftright to omit left-right logic in fast floating-point
100: * computation of dtoa.
101: * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
102: * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
103: * that use extended-precision instructions to compute rounded
104: * products and quotients) with IBM.
105: * #define ROUND_BIASED for IEEE-format with biased rounding.
106: * #define Inaccurate_Divide for IEEE-format with correctly rounded
107: * products but inaccurate quotients, e.g., for Intel i860.
108: * #define KR_headers for old-style C function headers.
109: */
110:
111: #ifdef DEBUG
112: #include <stdio.h>
113: #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
114: #endif
115:
116: #ifdef __STDC__
117: #include <stdlib.h>
118: #include <string.h>
119: #include <float.h>
120: #define CONST const
121: #else
122: #define CONST
123: #define KR_headers
124:
125: /* In this case, we assume IEEE floats. */
126: #define FLT_ROUNDS 1
127: #define FLT_RADIX 2
128: #define DBL_MANT_DIG 53
129: #define DBL_DIG 15
130: #define DBL_MAX_10_EXP 308
131: #define DBL_MAX_EXP 1024
132: #endif
133:
134: #include <errno.h>
135: #ifndef __MATH_H__
136: #include <math.h>
137: #endif
138:
139: #ifdef Unsigned_Shifts
140: #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
141: #else
142: #define Sign_Extend(a,b) /*no-op*/
143: #endif
144:
145: #if defined(__i386__) || defined(__i860__) || defined(clipper)
146: #define IEEE_8087
147: #endif
148: #if defined(MIPSEL) || defined(__alpha__)
149: #define IEEE_8087
150: #endif
151: #if defined(__sparc__) || defined(sparc) || defined(MIPSEB)
152: #define IEEE_MC68k
153: #endif
154:
155: #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
156:
157: #if FLT_RADIX==16
158: #define IBM
159: #else
160: #if DBL_MANT_DIG==56
161: #define VAX
162: #else
163: #if DBL_MANT_DIG==53 && DBL_MAX_10_EXP==308
164: #define IEEE_Unknown
165: #else
166: Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
167: #endif
168: #endif
169: #endif
170: #endif
171:
172: #ifdef IEEE_8087
173: #define HIWORD 1
174: #define LOWORD 0
175: #define TEST_ENDIANNESS /* nothing */
176: #else
177: #if defined(IEEE_MC68k)
178: #define HIWORD 0
179: #define LOWORD 1
180: #define TEST_ENDIANNESS /* nothing */
181: #else
182: static int HIWORD = -1, LOWORD;
183: static void test_endianness()
184: {
185: union doubleword {
186: double d;
187: unsigned long u[2];
188: } dw;
189: dw.d = 10;
190: if (dw.u[0] != 0) /* big-endian */
191: HIWORD=0, LOWORD=1;
192: else
193: HIWORD=1, LOWORD=0;
194: }
195: #define TEST_ENDIANNESS if (HIWORD<0) test_endianness();
196: #endif
197: #endif
198:
199: #if 0
200: union {
201: double d;
202: unsigned long x[2];
203: } _temp;
204: #endif
205: #define word0(x) ((unsigned long *)&x)[HIWORD]
206: #if 0
207: #define word0(X) (_temp.d = X, _temp.x[HIWORD])
208: #define setword0(D,X) (_temp.d = D, _temp.x[HIWORD] = X, D = _temp.d)
209: #endif
210: #define word1(x) ((unsigned long *)&x)[LOWORD]
211:
212: /* The following definition of Storeinc is appropriate for MIPS processors. */
213: #if defined(IEEE_8087) + defined(VAX)
214: #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
215: ((unsigned short *)a)[0] = (unsigned short)c, a++)
216: #else
217: #if defined(IEEE_MC68k)
218: #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
219: ((unsigned short *)a)[1] = (unsigned short)c, a++)
220: #else
221: #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
222: #endif
223: #endif
224:
225: /* #define P DBL_MANT_DIG */
226: /* Ten_pmax = floor(P*log(2)/log(5)) */
227: /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
228: /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
229: /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
230:
231: #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(IEEE_Unknown)
232: #define Exp_shift 20
233: #define Exp_shift1 20
234: #define Exp_msk1 0x100000
235: #define Exp_msk11 0x100000
236: #define Exp_mask 0x7ff00000
237: #define P 53
238: #define Bias 1023
239: #define IEEE_Arith
240: #define Emin (-1022)
241: #define Exp_1 0x3ff00000
242: #define Exp_11 0x3ff00000
243: #define Ebits 11
244: #define Frac_mask 0xfffff
245: #define Frac_mask1 0xfffff
246: #define Ten_pmax 22
247: #define Bletch 0x10
248: #define Bndry_mask 0xfffff
249: #define Bndry_mask1 0xfffff
250: #define LSB 1
251: #define Sign_bit 0x80000000
252: #define Log2P 1
253: #define Tiny0 0
254: #define Tiny1 1
255: #define Quick_max 14
256: #define Int_max 14
257: #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
258: #else
259: #undef Sudden_Underflow
260: #define Sudden_Underflow
261: #ifdef IBM
262: #define Exp_shift 24
263: #define Exp_shift1 24
264: #define Exp_msk1 0x1000000
265: #define Exp_msk11 0x1000000
266: #define Exp_mask 0x7f000000
267: #define P 14
268: #define Bias 65
269: #define Exp_1 0x41000000
270: #define Exp_11 0x41000000
271: #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
272: #define Frac_mask 0xffffff
273: #define Frac_mask1 0xffffff
274: #define Bletch 4
275: #define Ten_pmax 22
276: #define Bndry_mask 0xefffff
277: #define Bndry_mask1 0xffffff
278: #define LSB 1
279: #define Sign_bit 0x80000000
280: #define Log2P 4
281: #define Tiny0 0x100000
282: #define Tiny1 0
283: #define Quick_max 14
284: #define Int_max 15
285: #else /* VAX */
286: #define Exp_shift 23
287: #define Exp_shift1 7
288: #define Exp_msk1 0x80
289: #define Exp_msk11 0x800000
290: #define Exp_mask 0x7f80
291: #define P 56
292: #define Bias 129
293: #define Exp_1 0x40800000
294: #define Exp_11 0x4080
295: #define Ebits 8
296: #define Frac_mask 0x7fffff
297: #define Frac_mask1 0xffff007f
298: #define Ten_pmax 24
299: #define Bletch 2
300: #define Bndry_mask 0xffff007f
301: #define Bndry_mask1 0xffff007f
302: #define LSB 0x10000
303: #define Sign_bit 0x8000
304: #define Log2P 1
305: #define Tiny0 0x80
306: #define Tiny1 0
307: #define Quick_max 15
308: #define Int_max 15
309: #endif
310: #endif
311:
312: #ifndef IEEE_Arith
313: #define ROUND_BIASED
314: #endif
315:
316: #ifdef RND_PRODQUOT
317: #define rounded_product(a,b) a = rnd_prod(a, b)
318: #define rounded_quotient(a,b) a = rnd_quot(a, b)
319: extern double rnd_prod(double, double), rnd_quot(double, double);
320: #else
321: #define rounded_product(a,b) a *= b
322: #define rounded_quotient(a,b) a /= b
323: #endif
324:
325: #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
326: #define Big1 0xffffffff
327:
328: #define Kmax 15
329:
330: /* (1<<BIGINT_MINIMUM_K) is the minimum number of words to allocate
331: in a Bigint. dtoa usually manages with 1<<2, and has not been
332: known to need more than 1<<3. */
333:
334: #define BIGINT_MINIMUM_K 3
335:
336: struct Bigint {
337: struct Bigint *next;
338: int k; /* Parameter given to Balloc(k) */
339: int maxwds; /* Allocated space: equals 1<<k. */
340: short on_stack; /* 1 if stack-allocated. */
341: short sign; /* 0 if value is positive or zero; 1 if negative. */
342: int wds; /* Current length. */
343: unsigned long x[1<<BIGINT_MINIMUM_K]; /* Actually: x[maxwds] */
344: };
345:
346: #define BIGINT_HEADER_SIZE \
347: (sizeof(Bigint) - (1<<BIGINT_MINIMUM_K) * sizeof(unsigned long))
348:
349: typedef struct Bigint Bigint;
350:
351: /* Initialize a stack-allocated Bigint. */
352:
353: static Bigint *
354: Binit
355: #ifdef KR_headers
356: (v) Bigint *v;
357: #else
358: (Bigint *v)
359: #endif
360: {
361: v->on_stack = 1;
362: v->k = BIGINT_MINIMUM_K;
363: v->maxwds = 1 << BIGINT_MINIMUM_K;
364: v->sign = v->wds = 0;
365: return v;
366: }
367:
368: /* Allocate a Bigint with '1<<k' big digits. */
369:
370: static Bigint *
371: Balloc
372: #ifdef KR_headers
373: (k) int k;
374: #else
375: (int k)
376: #endif
377: {
378: int x;
379: Bigint *rv;
380:
381: if (k < BIGINT_MINIMUM_K)
382: k = BIGINT_MINIMUM_K;
383:
384: x = 1 << k;
385: rv = (Bigint *)
386: malloc(BIGINT_HEADER_SIZE + x * sizeof(unsigned long));
387: rv->k = k;
388: rv->maxwds = x;
389: rv->sign = rv->wds = 0;
390: rv->on_stack = 0;
391: return rv;
392: }
393:
394: static void
395: Bfree
396: #ifdef KR_headers
397: (v) Bigint *v;
398: #else
399: (Bigint *v)
400: #endif
401: {
402: if (v && !v->on_stack)
403: free (v);
404: }
405:
406: static void
407: Bcopy
408: #ifdef KR_headers
409: (x, y) Bigint *x, *y;
410: #else
411: (Bigint *x, Bigint *y)
412: #endif
413: {
414: register unsigned long *xp, *yp;
415: register int i = y->wds;
416: x->sign = y->sign;
417: x->wds = i;
418: for (xp = x->x, yp = y->x; --i >= 0; )
419: *xp++ = *yp++;
420: }
421:
422: /* Make sure b has room for at least 1<<k big digits. */
423:
424: static Bigint *
425: Brealloc
426: #ifdef KR_headers
427: (b, k) Bigint *b; int k;
428: #else
429: (Bigint * b, int k)
430: #endif
431: {
432: if (b == NULL)
433: return Balloc(k);
434: if (b->k >= k)
435: return b;
436: else
437: {
438: Bigint *rv = Balloc (k);
439: Bcopy(rv, b);
440: Bfree(b);
441: return rv;
442: }
443: }
444:
445: /* Return b*m+a. b is modified.
446: Assumption: 0xFFFF*m+a fits in 32 bits. */
447:
448: static Bigint *
449: multadd
450: #ifdef KR_headers
451: (b, m, a) Bigint *b; int m, a;
452: #else
453: (Bigint *b, int m, int a)
454: #endif
455: {
456: int i, wds;
457: unsigned long *x, y;
458: unsigned long xi, z;
459:
460: wds = b->wds;
461: x = b->x;
462: i = 0;
463: do {
464: xi = *x;
465: y = (xi & 0xffff) * m + a;
466: z = (xi >> 16) * m + (y >> 16);
467: a = (int)(z >> 16);
468: *x++ = (z << 16) + (y & 0xffff);
469: }
470: while(++i < wds);
471: if (a) {
472: if (wds >= b->maxwds)
473: b = Brealloc(b, b->k+1);
474: b->x[wds++] = a;
475: b->wds = wds;
476: }
477: return b;
478: }
479:
480: static Bigint *
481: s2b
482: #ifdef KR_headers
483: (result, s, nd0, nd, y9)
484: Bigint *result; CONST char *s; int nd0, nd; unsigned long y9;
485: #else
486: (Bigint *result, CONST char *s, int nd0, int nd, unsigned long y9)
487: #endif
488: {
489: int i, k;
490: long x, y;
491:
492: x = (nd + 8) / 9;
493: for(k = 0, y = 1; x > y; y <<= 1, k++) ;
494: result = Brealloc(result, k);
495: result->x[0] = y9;
496: result->wds = 1;
497:
498: i = 9;
499: if (9 < nd0)
500: {
501: s += 9;
502: do
503: result = multadd(result, 10, *s++ - '0');
504: while (++i < nd0);
505: s++;
506: }
507: else
508: s += 10;
509: for(; i < nd; i++)
510: result = multadd(result, 10, *s++ - '0');
511: return result;
512: }
513:
514: static int
515: hi0bits
516: #ifdef KR_headers
517: (x) register unsigned long x;
518: #else
519: (register unsigned long x)
520: #endif
521: {
522: register int k = 0;
523:
524: if (!(x & 0xffff0000)) {
525: k = 16;
526: x <<= 16;
527: }
528: if (!(x & 0xff000000)) {
529: k += 8;
530: x <<= 8;
531: }
532: if (!(x & 0xf0000000)) {
533: k += 4;
534: x <<= 4;
535: }
536: if (!(x & 0xc0000000)) {
537: k += 2;
538: x <<= 2;
539: }
540: if (!(x & 0x80000000)) {
541: k++;
542: if (!(x & 0x40000000))
543: return 32;
544: }
545: return k;
546: }
547:
548: static int
549: lo0bits
550: #ifdef KR_headers
551: (y) unsigned long *y;
552: #else
553: (unsigned long *y)
554: #endif
555: {
556: register int k;
557: register unsigned long x = *y;
558:
559: if (x & 7) {
560: if (x & 1)
561: return 0;
562: if (x & 2) {
563: *y = x >> 1;
564: return 1;
565: }
566: *y = x >> 2;
567: return 2;
568: }
569: k = 0;
570: if (!(x & 0xffff)) {
571: k = 16;
572: x >>= 16;
573: }
574: if (!(x & 0xff)) {
575: k += 8;
576: x >>= 8;
577: }
578: if (!(x & 0xf)) {
579: k += 4;
580: x >>= 4;
581: }
582: if (!(x & 0x3)) {
583: k += 2;
584: x >>= 2;
585: }
586: if (!(x & 1)) {
587: k++;
588: x >>= 1;
589: if (!x & 1)
590: return 32;
591: }
592: *y = x;
593: return k;
594: }
595:
596: static Bigint *
597: i2b
598: #ifdef KR_headers
599: (result, i) Bigint *result; int i;
600: #else
601: (Bigint* result, int i)
602: #endif
603: {
604: result = Brealloc(result, 1);
605: result->x[0] = i;
606: result->wds = 1;
607: return result;
608: }
609:
610: /* Do: c = a * b. */
611:
612: static Bigint *
613: mult
614: #ifdef KR_headers
615: (c, a, b) Bigint *a, *b, *c;
616: #else
617: (Bigint *c, Bigint *a, Bigint *b)
618: #endif
619: {
620: int k, wa, wb, wc;
621: unsigned long carry, y, z;
622: unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
623: unsigned long z2;
624: if (a->wds < b->wds) {
625: Bigint *tmp = a;
626: a = b;
627: b = tmp;
628: }
629: k = a->k;
630: wa = a->wds;
631: wb = b->wds;
632: wc = wa + wb;
633: if (wc > a->maxwds)
634: k++;
635: c = Brealloc(c, k);
636: for(x = c->x, xa = x + wc; x < xa; x++)
637: *x = 0;
638: xa = a->x;
639: xae = xa + wa;
640: xb = b->x;
641: xbe = xb + wb;
642: xc0 = c->x;
643: for(; xb < xbe; xb++, xc0++) {
644: if (y = *xb & 0xffff) {
645: x = xa;
646: xc = xc0;
647: carry = 0;
648: do {
649: z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
650: carry = z >> 16;
651: z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
652: carry = z2 >> 16;
653: Storeinc(xc, z2, z);
654: }
655: while(x < xae);
656: *xc = carry;
657: }
658: if (y = *xb >> 16) {
659: x = xa;
660: xc = xc0;
661: carry = 0;
662: z2 = *xc;
663: do {
664: z = (*x & 0xffff) * y + (*xc >> 16) + carry;
665: carry = z >> 16;
666: Storeinc(xc, z, z2);
667: z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
668: carry = z2 >> 16;
669: }
670: while(x < xae);
671: *xc = z2;
672: }
673: }
674: for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
675: c->wds = wc;
676: return c;
677: }
678:
679: /* Returns b*(5**k). b is modified. */
680: /* Re-written by Per Bothner to not need a static list. */
681:
682: static Bigint *
683: pow5mult
684: #ifdef KR_headers
685: (b, k) Bigint *b; int k;
686: #else
687: (Bigint *b, int k)
688: #endif
689: {
690: static int p05[6] = { 5, 25, 125, 625, 3125, 15625 };
691:
692: for (; k > 6; k -= 6)
693: b = multadd(b, 15625, 0); /* b *= 5**6 */
694: if (k == 0)
695: return b;
696: else
697: return multadd(b, p05[k-1], 0);
698: }
699:
700: /* Re-written by Per Bothner so shift can be in place. */
701:
702: static Bigint *
703: lshift
704: #ifdef KR_headers
705: (b, k) Bigint *b; int k;
706: #else
707: (Bigint *b, int k)
708: #endif
709: {
710: int i;
711: unsigned long *x, *x1, *xe;
712: int old_wds = b->wds;
713: int n = k >> 5;
714: int k1 = b->k;
715: int n1 = n + old_wds + 1;
716:
717: if (k == 0)
718: return b;
719:
720: for(i = b->maxwds; n1 > i; i <<= 1)
721: k1++;
722: b = Brealloc(b, k1);
723:
724: xe = b->x; /* Source limit */
725: x = xe + old_wds; /* Source pointer */
726: x1 = x + n; /* Destination pointer */
727: if (k &= 0x1f) {
728: int k1 = 32 - k;
729: unsigned long z = *--x;
730: if ((*x1 = (z >> k1)) != 0) {
731: ++n1;
732: }
733: while (x > xe) {
734: unsigned long w = *--x;
735: *--x1 = (z << k) | (w >> k1);
736: z = w;
737: }
738: *--x1 = z << k;
739: }
740: else
741: do {
742: *--x1 = *--x;
743: } while(x > xe);
744: while (x1 > xe)
745: *--x1 = 0;
746: b->wds = n1 - 1;
747: return b;
748: }
749:
750: static int
751: cmp
752: #ifdef KR_headers
753: (a, b) Bigint *a, *b;
754: #else
755: (Bigint *a, Bigint *b)
756: #endif
757: {
758: unsigned long *xa, *xa0, *xb, *xb0;
759: int i, j;
760:
761: i = a->wds;
762: j = b->wds;
763: #ifdef DEBUG
764: if (i > 1 && !a->x[i-1])
765: Bug("cmp called with a->x[a->wds-1] == 0");
766: if (j > 1 && !b->x[j-1])
767: Bug("cmp called with b->x[b->wds-1] == 0");
768: #endif
769: if (i -= j)
770: return i;
771: xa0 = a->x;
772: xa = xa0 + j;
773: xb0 = b->x;
774: xb = xb0 + j;
775: for(;;) {
776: if (*--xa != *--xb)
777: return *xa < *xb ? -1 : 1;
778: if (xa <= xa0)
779: break;
780: }
781: return 0;
782: }
783:
784: /* Do: c = a-b. */
785:
786: static Bigint *
787: diff
788: #ifdef KR_headers
789: (c, a, b) Bigint *c, *a, *b;
790: #else
791: (Bigint *c, Bigint *a, Bigint *b)
792: #endif
793: {
794: int i, wa, wb;
795: long borrow, y; /* We need signed shifts here. */
796: unsigned long *xa, *xae, *xb, *xbe, *xc;
797: long z;
798:
799: i = cmp(a,b);
800: if (!i) {
801: c = Brealloc(c, 0);
802: c->wds = 1;
803: c->x[0] = 0;
804: return c;
805: }
806: if (i < 0) {
807: Bigint *tmp = a;
808: a = b;
809: b = tmp;
810: i = 1;
811: }
812: else
813: i = 0;
814: c = Brealloc(c, a->k);
815: c->sign = i;
816: wa = a->wds;
817: xa = a->x;
818: xae = xa + wa;
819: wb = b->wds;
820: xb = b->x;
821: xbe = xb + wb;
822: xc = c->x;
823: borrow = 0;
824: do {
825: y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
826: borrow = y >> 16;
827: Sign_Extend(borrow, y);
828: z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
829: borrow = z >> 16;
830: Sign_Extend(borrow, z);
831: Storeinc(xc, z, y);
832: }
833: while(xb < xbe);
834: while(xa < xae) {
835: y = (*xa & 0xffff) + borrow;
836: borrow = y >> 16;
837: Sign_Extend(borrow, y);
838: z = (*xa++ >> 16) + borrow;
839: borrow = z >> 16;
840: Sign_Extend(borrow, z);
841: Storeinc(xc, z, y);
842: }
843: while(!*--xc)
844: wa--;
845: c->wds = wa;
846: return c;
847: }
848:
849: static double
850: ulp
851: #ifdef KR_headers
852: (x) double x;
853: #else
854: (double x)
855: #endif
856: {
857: register long L;
858: double a;
859:
860: L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
861: #ifndef Sudden_Underflow
862: if (L > 0) {
863: #endif
864: #ifdef IBM
865: L |= Exp_msk1 >> 4;
866: #endif
867: word0(a) = L;
868: word1(a) = 0;
869: #ifndef Sudden_Underflow
870: }
871: else {
872: L = -L >> Exp_shift;
873: if (L < Exp_shift) {
874: word0(a) = 0x80000 >> L;
875: word1(a) = 0;
876: }
877: else {
878: word0(a) = 0;
879: L -= Exp_shift;
880: word1(a) = L >= 31 ? 1 : 1 << 31 - L;
881: }
882: }
883: #endif
884: return a;
885: }
886:
887: static double
888: b2d
889: #ifdef KR_headers
890: (a, e) Bigint *a; int *e;
891: #else
892: (Bigint *a, int *e)
893: #endif
894: {
895: unsigned long *xa, *xa0, w, y, z;
896: int k;
897: double d;
898: #ifdef VAX
899: unsigned long d0, d1;
900: #else
901: #define d0 word0(d)
902: #define d1 word1(d)
903: #endif
904:
905: xa0 = a->x;
906: xa = xa0 + a->wds;
907: y = *--xa;
908: #ifdef DEBUG
909: if (!y) Bug("zero y in b2d");
910: #endif
911: k = hi0bits(y);
912: *e = 32 - k;
913: if (k < Ebits) {
914: d0 = Exp_1 | y >> Ebits - k;
915: w = xa > xa0 ? *--xa : 0;
916: d1 = y << (32-Ebits) + k | w >> Ebits - k;
917: goto ret_d;
918: }
919: z = xa > xa0 ? *--xa : 0;
920: if (k -= Ebits) {
921: d0 = Exp_1 | y << k | z >> 32 - k;
922: y = xa > xa0 ? *--xa : 0;
923: d1 = z << k | y >> 32 - k;
924: }
925: else {
926: d0 = Exp_1 | y;
927: d1 = z;
928: }
929: ret_d:
930: #ifdef VAX
931: word0(d) = d0 >> 16 | d0 << 16;
932: word1(d) = d1 >> 16 | d1 << 16;
933: #else
934: #undef d0
935: #undef d1
936: #endif
937: return d;
938: }
939:
940: static Bigint *
941: d2b
942: #ifdef KR_headers
943: (result, d, e, bits) Bigint *result; double d; int *e, *bits;
944: #else
945: (Bigint *result, double d, int *e, int *bits)
946: #endif
947: {
948: int de, i, k;
949: unsigned long *x, y, z;
950: #ifdef VAX
951: unsigned long d0, d1;
952: d0 = word0(d) >> 16 | word0(d) << 16;
953: d1 = word1(d) >> 16 | word1(d) << 16;
954: #else
955: #define d0 word0(d)
956: #define d1 word1(d)
957: #endif
958:
959: result = Brealloc(result, 1);
960: x = result->x;
961:
962: z = d0 & Frac_mask;
963: d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
964:
965: de = (int)(d0 >> Exp_shift); /* The exponent part of d. */
966:
967: /* Put back the suppressed high-order bit, if normalized. */
968: #ifndef IBM
969: #ifndef Sudden_Underflow
970: if (de)
971: #endif
972: z |= Exp_msk11;
973: #endif
974:
975: if (y = d1) {
976: if (k = lo0bits(&y)) {
977: x[0] = y | z << 32 - k;
978: z >>= k;
979: }
980: else
981: x[0] = y;
982: i = result->wds = (x[1] = z) ? 2 : 1;
983: }
984: else {
985: #ifdef DEBUG
986: if (!z)
987: Bug("Zero passed to d2b");
988: #endif
989: k = lo0bits(&z);
990: x[0] = z;
991: i = result->wds = 1;
992: k += 32;
993: }
994: #ifndef Sudden_Underflow
995: if (de) {
996: #endif
997: #ifdef IBM
998: *e = (de - Bias - (P-1) << 2) + k;
999: *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1000: #else
1001: *e = de - Bias - (P-1) + k;
1002: *bits = P - k;
1003: #endif
1004: #ifndef Sudden_Underflow
1005: }
1006: else {
1007: *e = de - Bias - (P-1) + 1 + k;
1008: *bits = 32*i - hi0bits(x[i-1]);
1009: }
1010: #endif
1011: return result;
1012: }
1013: #undef d0
1014: #undef d1
1015:
1016: static double
1017: ratio
1018: #ifdef KR_headers
1019: (a, b) Bigint *a, *b;
1020: #else
1021: (Bigint *a, Bigint *b)
1022: #endif
1023: {
1024: double da, db;
1025: int k, ka, kb;
1026:
1027: da = b2d(a, &ka);
1028: db = b2d(b, &kb);
1029: k = ka - kb + 32*(a->wds - b->wds);
1030: #ifdef IBM
1031: if (k > 0) {
1032: word0(da) += (k >> 2)*Exp_msk1;
1033: if (k &= 3)
1034: da *= 1 << k;
1035: }
1036: else {
1037: k = -k;
1038: word0(db) += (k >> 2)*Exp_msk1;
1039: if (k &= 3)
1040: db *= 1 << k;
1041: }
1042: #else
1043: if (k > 0)
1044: word0(da) += k*Exp_msk1;
1045: else {
1046: k = -k;
1047: word0(db) += k*Exp_msk1;
1048: }
1049: #endif
1050: return da / db;
1051: }
1052:
1053: static double
1054: tens[] = {
1055: 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1056: 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1057: 1e20, 1e21, 1e22
1058: #ifdef VAX
1059: , 1e23, 1e24
1060: #endif
1061: };
1062:
1063: static double
1064: #ifdef IEEE_Arith
1065: bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1066: static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1067: #define n_bigtens 5
1068: #else
1069: #ifdef IBM
1070: bigtens[] = { 1e16, 1e32, 1e64 };
1071: static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1072: #define n_bigtens 3
1073: #else
1074: bigtens[] = { 1e16, 1e32 };
1075: static double tinytens[] = { 1e-16, 1e-32 };
1076: #define n_bigtens 2
1077: #endif
1078: #endif
1079:
1080: double
1081: _IO_strtod
1082: #ifdef KR_headers
1083: (s00, se) CONST char *s00; char **se;
1084: #else
1085: (CONST char *s00, char **se)
1086: #endif
1087: {
1088: int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1089: e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1090: CONST char *s, *s0, *s1;
1091: double aadj, aadj1, adj, rv, rv0;
1092: long L;
1093: unsigned long y, z;
1094: Bigint _bb, _b_avail, _bd, _bd0, _bs, _delta;
1095: Bigint *bb = Binit(&_bb);
1096: Bigint *bd = Binit(&_bd);
1097: Bigint *bd0 = Binit(&_bd0);
1098: Bigint *bs = Binit(&_bs);
1099: Bigint *b_avail = Binit(&_b_avail);
1100: Bigint *delta = Binit(&_delta);
1101:
1102: TEST_ENDIANNESS;
1103: sign = nz0 = nz = 0;
1104: rv = 0.;
1105: for(s = s00;;s++) switch(*s) {
1106: case '-':
1107: sign = 1;
1108: /* no break */
1109: case '+':
1110: if (*++s)
1111: goto break2;
1112: /* no break */
1113: case 0:
1114: s = s00;
1115: goto ret;
1116: case '\t':
1117: case '\n':
1118: case '\v':
1119: case '\f':
1120: case '\r':
1121: case ' ':
1122: continue;
1123: default:
1124: goto break2;
1125: }
1126: break2:
1127: if (*s == '0') {
1128: nz0 = 1;
1129: while(*++s == '0') ;
1130: if (!*s)
1131: goto ret;
1132: }
1133: s0 = s;
1134: y = z = 0;
1135: for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1136: if (nd < 9)
1137: y = 10*y + c - '0';
1138: else if (nd < 16)
1139: z = 10*z + c - '0';
1140: nd0 = nd;
1141: if (c == '.') {
1142: c = *++s;
1143: if (!nd) {
1144: for(; c == '0'; c = *++s)
1145: nz++;
1146: if (c > '0' && c <= '9') {
1147: s0 = s;
1148: nf += nz;
1149: nz = 0;
1150: goto have_dig;
1151: }
1152: goto dig_done;
1153: }
1154: for(; c >= '0' && c <= '9'; c = *++s) {
1155: have_dig:
1156: nz++;
1157: if (c -= '0') {
1158: nf += nz;
1159: for(i = 1; i < nz; i++)
1160: if (nd++ < 9)
1161: y *= 10;
1162: else if (nd <= DBL_DIG + 1)
1163: z *= 10;
1164: if (nd++ < 9)
1165: y = 10*y + c;
1166: else if (nd <= DBL_DIG + 1)
1167: z = 10*z + c;
1168: nz = 0;
1169: }
1170: }
1171: }
1172: dig_done:
1173: e = 0;
1174: if (c == 'e' || c == 'E') {
1175: if (!nd && !nz && !nz0) {
1176: s = s00;
1177: goto ret;
1178: }
1179: s00 = s;
1180: esign = 0;
1181: switch(c = *++s) {
1182: case '-':
1183: esign = 1;
1184: case '+':
1185: c = *++s;
1186: }
1187: if (c >= '0' && c <= '9') {
1188: while(c == '0')
1189: c = *++s;
1190: if (c > '0' && c <= '9') {
1191: e = c - '0';
1192: s1 = s;
1193: while((c = *++s) >= '0' && c <= '9')
1194: e = 10*e + c - '0';
1195: if (s - s1 > 8)
1196: /* Avoid confusion from exponents
1197: * so large that e might overflow.
1198: */
1199: e = 9999999;
1200: if (esign)
1201: e = -e;
1202: }
1203: else
1204: e = 0;
1205: }
1206: else
1207: s = s00;
1208: }
1209: if (!nd) {
1210: if (!nz && !nz0)
1211: s = s00;
1212: goto ret;
1213: }
1214: e1 = e -= nf;
1215:
1216: /* Now we have nd0 digits, starting at s0, followed by a
1217: * decimal point, followed by nd-nd0 digits. The number we're
1218: * after is the integer represented by those digits times
1219: * 10**e */
1220:
1221: if (!nd0)
1222: nd0 = nd;
1223: k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1224: rv = y;
1225: if (k > 9)
1226: rv = tens[k - 9] * rv + z;
1227: if (nd <= DBL_DIG
1228: #ifndef RND_PRODQUOT
1229: && FLT_ROUNDS == 1
1230: #endif
1231: ) {
1232: if (!e)
1233: goto ret;
1234: if (e > 0) {
1235: if (e <= Ten_pmax) {
1236: #ifdef VAX
1237: goto vax_ovfl_check;
1238: #else
1239: /* rv = */ rounded_product(rv, tens[e]);
1240: goto ret;
1241: #endif
1242: }
1243: i = DBL_DIG - nd;
1244: if (e <= Ten_pmax + i) {
1245: /* A fancier test would sometimes let us do
1246: * this for larger i values.
1247: */
1248: e -= i;
1249: rv *= tens[i];
1250: #ifdef VAX
1251: /* VAX exponent range is so narrow we must
1252: * worry about overflow here...
1253: */
1254: vax_ovfl_check:
1255: word0(rv) -= P*Exp_msk1;
1256: /* rv = */ rounded_product(rv, tens[e]);
1257: if ((word0(rv) & Exp_mask)
1258: > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1259: goto ovfl;
1260: word0(rv) += P*Exp_msk1;
1261: #else
1262: /* rv = */ rounded_product(rv, tens[e]);
1263: #endif
1264: goto ret;
1265: }
1266: }
1267: #ifndef Inaccurate_Divide
1268: else if (e >= -Ten_pmax) {
1269: /* rv = */ rounded_quotient(rv, tens[-e]);
1270: goto ret;
1271: }
1272: #endif
1273: }
1274: e1 += nd - k;
1275:
1276: /* Get starting approximation = rv * 10**e1 */
1277:
1278: if (e1 > 0) {
1279: if (i = e1 & 15)
1280: rv *= tens[i];
1281: if (e1 &= ~15) {
1282: if (e1 > DBL_MAX_10_EXP) {
1283: ovfl:
1284: errno = ERANGE;
1285: #if defined(sun) && !defined(__svr4__)
1286: /* SunOS defines HUGE_VAL as __infinity(), which is in libm. */
1287: #undef HUGE_VAL
1288: #endif
1289: #ifndef HUGE_VAL
1290: #define HUGE_VAL 1.7976931348623157E+308
1291: #endif
1292: rv = HUGE_VAL;
1293: goto ret;
1294: }
1295: if (e1 >>= 4) {
1296: for(j = 0; e1 > 1; j++, e1 >>= 1)
1297: if (e1 & 1)
1298: rv *= bigtens[j];
1299: /* The last multiplication could overflow. */
1300: word0(rv) -= P*Exp_msk1;
1301: rv *= bigtens[j];
1302: if ((z = word0(rv) & Exp_mask)
1303: > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1304: goto ovfl;
1305: if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1306: /* set to largest number */
1307: /* (Can't trust DBL_MAX) */
1308: word0(rv) = Big0;
1309: word1(rv) = Big1;
1310: }
1311: else
1312: word0(rv) += P*Exp_msk1;
1313: }
1314:
1315: }
1316: }
1317: else if (e1 < 0) {
1318: e1 = -e1;
1319: if (i = e1 & 15)
1320: rv /= tens[i];
1321: if (e1 &= ~15) {
1322: e1 >>= 4;
1323: for(j = 0; e1 > 1; j++, e1 >>= 1)
1324: if (e1 & 1)
1325: rv *= tinytens[j];
1326: /* The last multiplication could underflow. */
1327: rv0 = rv;
1328: rv *= tinytens[j];
1329: if (!rv) {
1330: rv = 2.*rv0;
1331: rv *= tinytens[j];
1332: if (!rv) {
1333: undfl:
1334: rv = 0.;
1335: errno = ERANGE;
1336: goto ret;
1337: }
1338: word0(rv) = Tiny0;
1339: word1(rv) = Tiny1;
1340: /* The refinement below will clean
1341: * this approximation up.
1342: */
1343: }
1344: }
1345: }
1346:
1347: /* Now the hard part -- adjusting rv to the correct value.*/
1348:
1349: /* Put digits into bd: true value = bd * 10^e */
1350:
1351: bd0 = s2b(bd0, s0, nd0, nd, y);
1352: bd = Brealloc(bd, bd0->k);
1353:
1354: for(;;) {
1355: Bcopy(bd, bd0);
1356: bb = d2b(bb, rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1357: bs = i2b(bs, 1);
1358:
1359: if (e >= 0) {
1360: bb2 = bb5 = 0;
1361: bd2 = bd5 = e;
1362: }
1363: else {
1364: bb2 = bb5 = -e;
1365: bd2 = bd5 = 0;
1366: }
1367: if (bbe >= 0)
1368: bb2 += bbe;
1369: else
1370: bd2 -= bbe;
1371: bs2 = bb2;
1372: #ifdef Sudden_Underflow
1373: #ifdef IBM
1374: j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1375: #else
1376: j = P + 1 - bbbits;
1377: #endif
1378: #else
1379: i = bbe + bbbits - 1; /* logb(rv) */
1380: if (i < Emin) /* denormal */
1381: j = bbe + (P-Emin);
1382: else
1383: j = P + 1 - bbbits;
1384: #endif
1385: bb2 += j;
1386: bd2 += j;
1387: i = bb2 < bd2 ? bb2 : bd2;
1388: if (i > bs2)
1389: i = bs2;
1390: if (i > 0) {
1391: bb2 -= i;
1392: bd2 -= i;
1393: bs2 -= i;
1394: }
1395: if (bb5 > 0) {
1396: Bigint *b_tmp;
1397: bs = pow5mult(bs, bb5);
1398: b_tmp = mult(b_avail, bs, bb);
1399: b_avail = bb;
1400: bb = b_tmp;
1401: }
1402: if (bb2 > 0)
1403: bb = lshift(bb, bb2);
1404: if (bd5 > 0)
1405: bd = pow5mult(bd, bd5);
1406: if (bd2 > 0)
1407: bd = lshift(bd, bd2);
1408: if (bs2 > 0)
1409: bs = lshift(bs, bs2);
1410: delta = diff(delta, bb, bd);
1411: dsign = delta->sign;
1412: delta->sign = 0;
1413: i = cmp(delta, bs);
1414: if (i < 0) {
1415: /* Error is less than half an ulp -- check for
1416: * special case of mantissa a power of two.
1417: */
1418: if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1419: break;
1420: delta = lshift(delta,Log2P);
1421: if (cmp(delta, bs) > 0)
1422: goto drop_down;
1423: break;
1424: }
1425: if (i == 0) {
1426: /* exactly half-way between */
1427: if (dsign) {
1428: if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1429: && word1(rv) == 0xffffffff) {
1430: /*boundary case -- increment exponent*/
1431: word0(rv) = (word0(rv) & Exp_mask)
1432: + Exp_msk1
1433: #ifdef IBM
1434: | Exp_msk1 >> 4
1435: #endif
1436: ;
1437: word1(rv) = 0;
1438: break;
1439: }
1440: }
1441: else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1442: drop_down:
1443: /* boundary case -- decrement exponent */
1444: #ifdef Sudden_Underflow
1445: L = word0(rv) & Exp_mask;
1446: #ifdef IBM
1447: if (L < Exp_msk1)
1448: #else
1449: if (L <= Exp_msk1)
1450: #endif
1451: goto undfl;
1452: L -= Exp_msk1;
1453: #else
1454: L = (word0(rv) & Exp_mask) - Exp_msk1;
1455: #endif
1456: word0(rv) = L | Bndry_mask1;
1457: word1(rv) = 0xffffffff;
1458: #ifdef IBM
1459: continue;
1460: #else
1461: break;
1462: #endif
1463: }
1464: #ifndef ROUND_BIASED
1465: if (!(word1(rv) & LSB))
1466: break;
1467: #endif
1468: if (dsign)
1469: rv += ulp(rv);
1470: #ifndef ROUND_BIASED
1471: else {
1472: rv -= ulp(rv);
1473: #ifndef Sudden_Underflow
1474: if (!rv)
1475: goto undfl;
1476: #endif
1477: }
1478: #endif
1479: break;
1480: }
1481: if ((aadj = ratio(delta, bs)) <= 2.) {
1482: if (dsign)
1483: aadj = aadj1 = 1.;
1484: else if (word1(rv) || word0(rv) & Bndry_mask) {
1485: #ifndef Sudden_Underflow
1486: if (word1(rv) == Tiny1 && !word0(rv))
1487: goto undfl;
1488: #endif
1489: aadj = 1.;
1490: aadj1 = -1.;
1491: }
1492: else {
1493: /* special case -- power of FLT_RADIX to be */
1494: /* rounded down... */
1495:
1496: if (aadj < 2./FLT_RADIX)
1497: aadj = 1./FLT_RADIX;
1498: else
1499: aadj *= 0.5;
1500: aadj1 = -aadj;
1501: }
1502: }
1503: else {
1504: aadj *= 0.5;
1505: aadj1 = dsign ? aadj : -aadj;
1506: #ifdef Check_FLT_ROUNDS
1507: switch(FLT_ROUNDS) {
1508: case 2: /* towards +infinity */
1509: aadj1 -= 0.5;
1510: break;
1511: case 0: /* towards 0 */
1512: case 3: /* towards -infinity */
1513: aadj1 += 0.5;
1514: }
1515: #else
1516: if (FLT_ROUNDS == 0)
1517: aadj1 += 0.5;
1518: #endif
1519: }
1520: y = word0(rv) & Exp_mask;
1521:
1522: /* Check for overflow */
1523:
1524: if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1525: rv0 = rv;
1526: word0(rv) -= P*Exp_msk1;
1527: adj = aadj1 * ulp(rv);
1528: rv += adj;
1529: if ((word0(rv) & Exp_mask) >=
1530: Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1531: if (word0(rv0) == Big0 && word1(rv0) == Big1)
1532: goto ovfl;
1533: word0(rv) = Big0;
1534: word1(rv) = Big1;
1535: continue;
1536: }
1537: else
1538: word0(rv) += P*Exp_msk1;
1539: }
1540: else {
1541: #ifdef Sudden_Underflow
1542: if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1543: rv0 = rv;
1544: word0(rv) += P*Exp_msk1;
1545: adj = aadj1 * ulp(rv);
1546: rv += adj;
1547: #ifdef IBM
1548: if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1549: #else
1550: if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1551: #endif
1552: {
1553: if (word0(rv0) == Tiny0
1554: && word1(rv0) == Tiny1)
1555: goto undfl;
1556: word0(rv) = Tiny0;
1557: word1(rv) = Tiny1;
1558: continue;
1559: }
1560: else
1561: word0(rv) -= P*Exp_msk1;
1562: }
1563: else {
1564: adj = aadj1 * ulp(rv);
1565: rv += adj;
1566: }
1567: #else
1568: /* Compute adj so that the IEEE rounding rules will
1569: * correctly round rv + adj in some half-way cases.
1570: * If rv * ulp(rv) is denormalized (i.e.,
1571: * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1572: * trouble from bits lost to denormalization;
1573: * example: 1.2e-307 .
1574: */
1575: if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1576: aadj1 = (double)(int)(aadj + 0.5);
1577: if (!dsign)
1578: aadj1 = -aadj1;
1579: }
1580: adj = aadj1 * ulp(rv);
1581: rv += adj;
1582: #endif
1583: }
1584: z = word0(rv) & Exp_mask;
1585: if (y == z) {
1586: /* Can we stop now? */
1587: L = (long)aadj;
1588: aadj -= L;
1589: /* The tolerances below are conservative. */
1590: if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1591: if (aadj < .4999999 || aadj > .5000001)
1592: break;
1593: }
1594: else if (aadj < .4999999/FLT_RADIX)
1595: break;
1596: }
1597: }
1598: Bfree(bb);
1599: Bfree(bd);
1600: Bfree(bs);
1601: Bfree(bd0);
1602: Bfree(delta);
1603: Bfree(b_avail);
1604: ret:
1605: if (se)
1606: *se = (char *)s;
1607: return sign ? -rv : rv;
1608: }
1609:
1610: static int
1611: quorem
1612: #ifdef KR_headers
1613: (b, S) Bigint *b, *S;
1614: #else
1615: (Bigint *b, Bigint *S)
1616: #endif
1617: {
1618: int n;
1619: long borrow, y;
1620: unsigned long carry, q, ys;
1621: unsigned long *bx, *bxe, *sx, *sxe;
1622: long z;
1623: unsigned long si, zs;
1624:
1625: n = S->wds;
1626: #ifdef DEBUG
1627: /*debug*/ if (b->wds > n)
1628: /*debug*/ Bug("oversize b in quorem");
1629: #endif
1630: if (b->wds < n)
1631: return 0;
1632: sx = S->x;
1633: sxe = sx + --n;
1634: bx = b->x;
1635: bxe = bx + n;
1636: q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1637: #ifdef DEBUG
1638: /*debug*/ if (q > 9)
1639: /*debug*/ Bug("oversized quotient in quorem");
1640: #endif
1641: if (q) {
1642: borrow = 0;
1643: carry = 0;
1644: do {
1645: si = *sx++;
1646: ys = (si & 0xffff) * q + carry;
1647: zs = (si >> 16) * q + (ys >> 16);
1648: carry = zs >> 16;
1649: y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1650: borrow = y >> 16;
1651: Sign_Extend(borrow, y);
1652: z = (*bx >> 16) - (zs & 0xffff) + borrow;
1653: borrow = z >> 16;
1654: Sign_Extend(borrow, z);
1655: Storeinc(bx, z, y);
1656: }
1657: while(sx <= sxe);
1658: if (!*bxe) {
1659: bx = b->x;
1660: while(--bxe > bx && !*bxe)
1661: --n;
1662: b->wds = n;
1663: }
1664: }
1665: if (cmp(b, S) >= 0) {
1666: q++;
1667: borrow = 0;
1668: carry = 0;
1669: bx = b->x;
1670: sx = S->x;
1671: do {
1672: si = *sx++;
1673: ys = (si & 0xffff) + carry;
1674: zs = (si >> 16) + (ys >> 16);
1675: carry = zs >> 16;
1676: y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1677: borrow = y >> 16;
1678: Sign_Extend(borrow, y);
1679: z = (*bx >> 16) - (zs & 0xffff) + borrow;
1680: borrow = z >> 16;
1681: Sign_Extend(borrow, z);
1682: Storeinc(bx, z, y);
1683: }
1684: while(sx <= sxe);
1685: bx = b->x;
1686: bxe = bx + n;
1687: if (!*bxe) {
1688: while(--bxe > bx && !*bxe)
1689: --n;
1690: b->wds = n;
1691: }
1692: }
1693: return q;
1694: }
1695:
1696: /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1697: *
1698: * Inspired by "How to Print Floating-Point Numbers Accurately" by
1699: * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1700: *
1701: * Modifications:
1702: * 1. Rather than iterating, we use a simple numeric overestimate
1703: * to determine k = floor(log10(d)). We scale relevant
1704: * quantities using O(log2(k)) rather than O(k) multiplications.
1705: * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1706: * try to generate digits strictly left to right. Instead, we
1707: * compute with fewer bits and propagate the carry if necessary
1708: * when rounding the final digit up. This is often faster.
1709: * 3. Under the assumption that input will be rounded nearest,
1710: * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1711: * That is, we allow equality in stopping tests when the
1712: * round-nearest rule will give the same floating-point value
1713: * as would satisfaction of the stopping test with strict
1714: * inequality.
1715: * 4. We remove common factors of powers of 2 from relevant
1716: * quantities.
1717: * 5. When converting floating-point integers less than 1e16,
1718: * we use floating-point arithmetic rather than resorting
1719: * to multiple-precision integers.
1720: * 6. When asked to produce fewer than 15 digits, we first try
1721: * to get by with floating-point arithmetic; we resort to
1722: * multiple-precision integer arithmetic only if we cannot
1723: * guarantee that the floating-point calculation has given
1724: * the correctly rounded result. For k requested digits and
1725: * "uniformly" distributed input, the probability is
1726: * something like 10^(k-15) that we must resort to the long
1727: * calculation.
1728: */
1729:
1730: char *
1731: _IO_dtoa
1732: #ifdef KR_headers
1733: (d, mode, ndigits, decpt, sign, rve)
1734: double d; int mode, ndigits, *decpt, *sign; char **rve;
1735: #else
1736: (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
1737: #endif
1738: {
1739: /* Arguments ndigits, decpt, sign are similar to those
1740: of ecvt and fcvt; trailing zeros are suppressed from
1741: the returned string. If not null, *rve is set to point
1742: to the end of the return value. If d is +-Infinity or NaN,
1743: then *decpt is set to 9999.
1744:
1745: mode:
1746: 0 ==> shortest string that yields d when read in
1747: and rounded to nearest.
1748: 1 ==> like 0, but with Steele & White stopping rule;
1749: e.g. with IEEE P754 arithmetic , mode 0 gives
1750: 1e23 whereas mode 1 gives 9.999999999999999e22.
1751: 2 ==> max(1,ndigits) significant digits. This gives a
1752: return value similar to that of ecvt, except
1753: that trailing zeros are suppressed.
1754: 3 ==> through ndigits past the decimal point. This
1755: gives a return value similar to that from fcvt,
1756: except that trailing zeros are suppressed, and
1757: ndigits can be negative.
1758: 4-9 should give the same return values as 2-3, i.e.,
1759: 4 <= mode <= 9 ==> same return as mode
1760: 2 + (mode & 1). These modes are mainly for
1761: debugging; often they run slower but sometimes
1762: faster than modes 2-3.
1763: 4,5,8,9 ==> left-to-right digit generation.
1764: 6-9 ==> don't try fast floating-point estimate
1765: (if applicable).
1766:
1767: Values of mode other than 0-9 are treated as mode 0.
1768:
1769: Sufficient space is allocated to the return value
1770: to hold the suppressed trailing zeros.
1771: */
1772:
1773: int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1774: j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1775: spec_case, try_quick;
1776: long L;
1777: #ifndef Sudden_Underflow
1778: int denorm;
1779: #endif
1780: Bigint _b_avail, _b, _mhi, _mlo, _S;
1781: Bigint *b_avail = Binit(&_b_avail);
1782: Bigint *b = Binit(&_b);
1783: Bigint *S = Binit(&_S);
1784: /* mhi and mlo are only set and used if leftright. */
1785: Bigint *mhi = NULL, *mlo = NULL;
1786: double d2, ds, eps;
1787: char *s, *s0;
1788: static Bigint *result = NULL;
1789: static int result_k;
1790:
1791: TEST_ENDIANNESS;
1792: if (result) {
1793: result->k = result_k;
1794: result->maxwds = 1 << result_k;
1795: }
1796:
1797: if (word0(d) & Sign_bit) {
1798: /* set sign for everything, including 0's and NaNs */
1799: *sign = 1;
1800: word0(d) &= ~Sign_bit; /* clear sign bit */
1801: }
1802: else
1803: *sign = 0;
1804:
1805: #if defined(IEEE_Arith) + defined(VAX)
1806: #ifdef IEEE_Arith
1807: if ((word0(d) & Exp_mask) == Exp_mask)
1808: #else
1809: if (word0(d) == 0x8000)
1810: #endif
1811: {
1812: /* Infinity or NaN */
1813: *decpt = 9999;
1814: #ifdef IEEE_Arith
1815: if (!word1(d) && !(word0(d) & 0xfffff))
1816: {
1817: s = "Infinity";
1818: if (rve)
1819: *rve = s + 8;
1820: }
1821: else
1822: #endif
1823: {
1824: s = "NaN";
1825: if (rve)
1826: *rve = s +3;
1827: }
1828: return s;
1829: }
1830: #endif
1831: #ifdef IBM
1832: d += 0; /* normalize */
1833: #endif
1834: if (!d) {
1835: *decpt = 1;
1836: s = "0";
1837: if (rve)
1838: *rve = s + 1;
1839: return s;
1840: }
1841:
1842: b = d2b(b, d, &be, &bbits);
1843: i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1844: #ifndef Sudden_Underflow
1845: if (i) {
1846: #endif
1847: d2 = d;
1848: word0(d2) &= Frac_mask1;
1849: word0(d2) |= Exp_11;
1850: #ifdef IBM
1851: if (j = 11 - hi0bits(word0(d2) & Frac_mask))
1852: d2 /= 1 << j;
1853: #endif
1854:
1855: i -= Bias;
1856: #ifdef IBM
1857: i <<= 2;
1858: i += j;
1859: #endif
1860: #ifndef Sudden_Underflow
1861: denorm = 0;
1862: }
1863: else {
1864: /* d is denormalized */
1865: unsigned long x;
1866:
1867: i = bbits + be + (Bias + (P-1) - 1);
1868: x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32
1869: : word1(d) << 32 - i;
1870: d2 = x;
1871: word0(d2) -= 31*Exp_msk1; /* adjust exponent */
1872: i -= (Bias + (P-1) - 1) + 1;
1873: denorm = 1;
1874: }
1875: #endif
1876:
1877: /* Now i is the unbiased base-2 exponent. */
1878:
1879: /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1880: * log10(x) = log(x) / log(10)
1881: * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1882: * log10(d) = i*log(2)/log(10) + log10(d2)
1883: *
1884: * This suggests computing an approximation k to log10(d) by
1885: *
1886: * k = i*0.301029995663981
1887: * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1888: *
1889: * We want k to be too large rather than too small.
1890: * The error in the first-order Taylor series approximation
1891: * is in our favor, so we just round up the constant enough
1892: * to compensate for any error in the multiplication of
1893: * (i) by 0.301029995663981; since |i| <= 1077,
1894: * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1895: * adding 1e-13 to the constant term more than suffices.
1896: * Hence we adjust the constant term to 0.1760912590558.
1897: * (We could get a more accurate k by invoking log10,
1898: * but this is probably not worthwhile.)
1899: */
1900:
1901: ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1902: k = (int)ds;
1903: if (ds < 0. && ds != k)
1904: k--; /* want k = floor(ds) */
1905: k_check = 1;
1906: if (k >= 0 && k <= Ten_pmax) {
1907: if (d < tens[k])
1908: k--;
1909: k_check = 0;
1910: }
1911: j = bbits - i - 1;
1912: if (j >= 0) {
1913: b2 = 0;
1914: s2 = j;
1915: }
1916: else {
1917: b2 = -j;
1918: s2 = 0;
1919: }
1920: if (k >= 0) {
1921: b5 = 0;
1922: s5 = k;
1923: s2 += k;
1924: }
1925: else {
1926: b2 -= k;
1927: b5 = -k;
1928: s5 = 0;
1929: }
1930: if (mode < 0 || mode > 9)
1931: mode = 0;
1932: try_quick = 1;
1933: if (mode > 5) {
1934: mode -= 4;
1935: try_quick = 0;
1936: }
1937: leftright = 1;
1938: switch(mode) {
1939: case 0:
1940: case 1:
1941: ilim = ilim1 = -1;
1942: i = 18;
1943: ndigits = 0;
1944: break;
1945: case 2:
1946: leftright = 0;
1947: /* no break */
1948: case 4:
1949: if (ndigits <= 0)
1950: ndigits = 1;
1951: ilim = ilim1 = i = ndigits;
1952: break;
1953: case 3:
1954: leftright = 0;
1955: /* no break */
1956: case 5:
1957: i = ndigits + k + 1;
1958: ilim = i;
1959: ilim1 = i - 1;
1960: if (i <= 0)
1961: i = 1;
1962: }
1963: /* i is now an upper bound of the number of digits to generate. */
1964: j = sizeof(unsigned long) * (1<<BIGINT_MINIMUM_K);
1965: /* The test is <= so as to allow room for the final '\0'. */
1966: for(result_k = BIGINT_MINIMUM_K; BIGINT_HEADER_SIZE + j <= i;
1967: j <<= 1) result_k++;
1968: result = Brealloc(result, result_k);
1969: s = s0 = (char *)result;
1970:
1971: if (ilim >= 0 && ilim <= Quick_max && try_quick) {
1972:
1973: /* Try to get by with floating-point arithmetic. */
1974:
1975: i = 0;
1976: d2 = d;
1977: k0 = k;
1978: ilim0 = ilim;
1979: ieps = 2; /* conservative */
1980: if (k > 0) {
1981: ds = tens[k&0xf];
1982: j = k >> 4;
1983: if (j & Bletch) {
1984: /* prevent overflows */
1985: j &= Bletch - 1;
1986: d /= bigtens[n_bigtens-1];
1987: ieps++;
1988: }
1989: for(; j; j >>= 1, i++)
1990: if (j & 1) {
1991: ieps++;
1992: ds *= bigtens[i];
1993: }
1994: d /= ds;
1995: }
1996: else if (j1 = -k) {
1997: d *= tens[j1 & 0xf];
1998: for(j = j1 >> 4; j; j >>= 1, i++)
1999: if (j & 1) {
2000: ieps++;
2001: d *= bigtens[i];
2002: }
2003: }
2004: if (k_check && d < 1. && ilim > 0) {
2005: if (ilim1 <= 0)
2006: goto fast_failed;
2007: ilim = ilim1;
2008: k--;
2009: d *= 10.;
2010: ieps++;
2011: }
2012: eps = ieps*d + 7.;
2013: word0(eps) -= (P-1)*Exp_msk1;
2014: if (ilim == 0) {
2015: d -= 5.;
2016: if (d > eps)
2017: goto one_digit;
2018: if (d < -eps)
2019: goto no_digits;
2020: goto fast_failed;
2021: }
2022: #ifndef No_leftright
2023: if (leftright) {
2024: /* Use Steele & White method of only
2025: * generating digits needed.
2026: */
2027: eps = 0.5/tens[ilim-1] - eps;
2028: for(i = 0;;) {
2029: L = (long)d;
2030: d -= L;
2031: *s++ = '0' + (int)L;
2032: if (d < eps)
2033: goto ret1;
2034: if (1. - d < eps)
2035: goto bump_up;
2036: if (++i >= ilim)
2037: break;
2038: eps *= 10.;
2039: d *= 10.;
2040: }
2041: }
2042: else {
2043: #endif
2044: /* Generate ilim digits, then fix them up. */
2045: eps *= tens[ilim-1];
2046: for(i = 1;; i++, d *= 10.) {
2047: L = (long)d;
2048: d -= L;
2049: *s++ = '0' + (int)L;
2050: if (i == ilim) {
2051: if (d > 0.5 + eps)
2052: goto bump_up;
2053: else if (d < 0.5 - eps) {
2054: while(*--s == '0');
2055: s++;
2056: goto ret1;
2057: }
2058: break;
2059: }
2060: }
2061: #ifndef No_leftright
2062: }
2063: #endif
2064: fast_failed:
2065: s = s0;
2066: d = d2;
2067: k = k0;
2068: ilim = ilim0;
2069: }
2070:
2071: /* Do we have a "small" integer? */
2072:
2073: if (be >= 0 && k <= Int_max) {
2074: /* Yes. */
2075: ds = tens[k];
2076: if (ndigits < 0 && ilim <= 0) {
2077: if (ilim < 0 || d <= 5*ds)
2078: goto no_digits;
2079: goto one_digit;
2080: }
2081: for(i = 1;; i++) {
2082: L = (long)(d / ds);
2083: d -= L*ds;
2084: #ifdef Check_FLT_ROUNDS
2085: /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2086: if (d < 0) {
2087: L--;
2088: d += ds;
2089: }
2090: #endif
2091: *s++ = '0' + (int)L;
2092: if (i == ilim) {
2093: d += d;
2094: if (d > ds || d == ds && L & 1) {
2095: bump_up:
2096: while(*--s == '9')
2097: if (s == s0) {
2098: k++;
2099: *s = '0';
2100: break;
2101: }
2102: ++*s++;
2103: }
2104: break;
2105: }
2106: if (!(d *= 10.))
2107: break;
2108: }
2109: goto ret1;
2110: }
2111:
2112: m2 = b2;
2113: m5 = b5;
2114: if (leftright) {
2115: if (mode < 2) {
2116: i =
2117: #ifndef Sudden_Underflow
2118: denorm ? be + (Bias + (P-1) - 1 + 1) :
2119: #endif
2120: #ifdef IBM
2121: 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2122: #else
2123: 1 + P - bbits;
2124: #endif
2125: }
2126: else {
2127: j = ilim - 1;
2128: if (m5 >= j)
2129: m5 -= j;
2130: else {
2131: s5 += j -= m5;
2132: b5 += j;
2133: m5 = 0;
2134: }
2135: if ((i = ilim) < 0) {
2136: m2 -= i;
2137: i = 0;
2138: }
2139: }
2140: b2 += i;
2141: s2 += i;
2142: mhi = i2b(Binit(&_mhi), 1);
2143: }
2144: if (m2 > 0 && s2 > 0) {
2145: i = m2 < s2 ? m2 : s2;
2146: b2 -= i;
2147: m2 -= i;
2148: s2 -= i;
2149: }
2150: if (b5 > 0) {
2151: if (leftright) {
2152: if (m5 > 0) {
2153: Bigint *b_tmp;
2154: mhi = pow5mult(mhi, m5);
2155: b_tmp = mult(b_avail, mhi, b);
2156: b_avail = b;
2157: b = b_tmp;
2158: }
2159: if (j = b5 - m5)
2160: b = pow5mult(b, j);
2161: }
2162: else
2163: b = pow5mult(b, b5);
2164: }
2165: S = i2b(S, 1);
2166: if (s5 > 0)
2167: S = pow5mult(S, s5);
2168:
2169: /* Check for special case that d is a normalized power of 2. */
2170:
2171: if (mode < 2) {
2172: if (!word1(d) && !(word0(d) & Bndry_mask)
2173: #ifndef Sudden_Underflow
2174: && word0(d) & Exp_mask
2175: #endif
2176: ) {
2177: /* The special case */
2178: b2 += Log2P;
2179: s2 += Log2P;
2180: spec_case = 1;
2181: }
2182: else
2183: spec_case = 0;
2184: }
2185:
2186: /* Arrange for convenient computation of quotients:
2187: * shift left if necessary so divisor has 4 leading 0 bits.
2188: *
2189: * Perhaps we should just compute leading 28 bits of S once
2190: * and for all and pass them and a shift to quorem, so it
2191: * can do shifts and ors to compute the numerator for q.
2192: */
2193: if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
2194: i = 32 - i;
2195: if (i > 4) {
2196: i -= 4;
2197: b2 += i;
2198: m2 += i;
2199: s2 += i;
2200: }
2201: else if (i < 4) {
2202: i += 28;
2203: b2 += i;
2204: m2 += i;
2205: s2 += i;
2206: }
2207: if (b2 > 0)
2208: b = lshift(b, b2);
2209: if (s2 > 0)
2210: S = lshift(S, s2);
2211: if (k_check) {
2212: if (cmp(b,S) < 0) {
2213: k--;
2214: b = multadd(b, 10, 0); /* we botched the k estimate */
2215: if (leftright)
2216: mhi = multadd(mhi, 10, 0);
2217: ilim = ilim1;
2218: }
2219: }
2220: if (ilim <= 0 && mode > 2) {
2221: if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2222: /* no digits, fcvt style */
2223: no_digits:
2224: k = -1 - ndigits;
2225: goto ret;
2226: }
2227: one_digit:
2228: *s++ = '1';
2229: k++;
2230: goto ret;
2231: }
2232: if (leftright) {
2233: if (m2 > 0)
2234: mhi = lshift(mhi, m2);
2235:
2236: /* Compute mlo -- check for special case
2237: * that d is a normalized power of 2.
2238: */
2239:
2240: if (spec_case) {
2241: mlo = Brealloc(Binit(&_mlo), mhi->k);
2242: Bcopy(mlo, mhi);
2243: mhi = lshift(mhi, Log2P);
2244: }
2245: else
2246: mlo = mhi;
2247:
2248: for(i = 1;;i++) {
2249: dig = quorem(b,S) + '0';
2250: /* Do we yet have the shortest decimal string
2251: * that will round to d?
2252: */
2253: j = cmp(b, mlo);
2254: b_avail = diff(b_avail, S, mhi); /* b_avail = S - mi */
2255: j1 = b_avail->sign ? 1 : cmp(b, b_avail);
2256: #ifndef ROUND_BIASED
2257: if (j1 == 0 && !mode && !(word1(d) & 1)) {
2258: if (dig == '9')
2259: goto round_9_up;
2260: if (j > 0)
2261: dig++;
2262: *s++ = dig;
2263: goto ret;
2264: }
2265: #endif
2266: if (j < 0 || j == 0 && !mode
2267: #ifndef ROUND_BIASED
2268: && !(word1(d) & 1)
2269: #endif
2270: ) {
2271: if (j1 > 0) {
2272: b = lshift(b, 1);
2273: j1 = cmp(b, S);
2274: if ((j1 > 0 || j1 == 0 && dig & 1)
2275: && dig++ == '9')
2276: goto round_9_up;
2277: }
2278: *s++ = dig;
2279: goto ret;
2280: }
2281: if (j1 > 0) {
2282: if (dig == '9') { /* possible if i == 1 */
2283: round_9_up:
2284: *s++ = '9';
2285: goto roundoff;
2286: }
2287: *s++ = dig + 1;
2288: goto ret;
2289: }
2290: *s++ = dig;
2291: if (i == ilim)
2292: break;
2293: b = multadd(b, 10, 0);
2294: if (mlo == mhi)
2295: mlo = mhi = multadd(mhi, 10, 0);
2296: else {
2297: mlo = multadd(mlo, 10, 0);
2298: mhi = multadd(mhi, 10, 0);
2299: }
2300: }
2301: }
2302: else
2303: for(i = 1;; i++) {
2304: *s++ = dig = quorem(b,S) + '0';
2305: if (i >= ilim)
2306: break;
2307: b = multadd(b, 10, 0);
2308: }
2309:
2310: /* Round off last digit */
2311:
2312: b = lshift(b, 1);
2313: j = cmp(b, S);
2314: if (j > 0 || j == 0 && dig & 1) {
2315: roundoff:
2316: while(*--s == '9')
2317: if (s == s0) {
2318: k++;
2319: *s++ = '1';
2320: goto ret;
2321: }
2322: ++*s++;
2323: }
2324: else {
2325: while(*--s == '0');
2326: s++;
2327: }
2328: ret:
2329: Bfree(b_avail);
2330: Bfree(S);
2331: if (mhi) {
2332: if (mlo && mlo != mhi)
2333: Bfree(mlo);
2334: Bfree(mhi);
2335: }
2336: ret1:
2337: Bfree(b);
2338: *s = 0;
2339: *decpt = k + 1;
2340: if (rve)
2341: *rve = s;
2342: return s0;
2343: }
2344: #endif /* USE_DTOA */
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.