1 /****************************************************************
3 * The author of this software is David M. Gay.
5 * Copyright (c) 1991 by AT&T.
7 * Permission to use, copy, modify, and distribute this software for any
8 * purpose without fee is hereby granted, provided that this entire notice
9 * is included in all copies of any software which is or includes a copy
10 * or modification of this software and in all copies of the supporting
11 * documentation for such software.
13 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
15 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
18 ***************************************************************/
20 /* Please send bug reports to
22 AT&T Bell Laboratories, Room 2C-463
24 Murray Hill, NJ 07974-2070
26 dmg@research.att.com or research!dmg
29 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
31 * This strtod returns a nearest machine number to the input decimal
32 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
33 * broken by the IEEE round-even rule. Otherwise ties are broken by
34 * biased rounding (add half and chop).
36 * Inspired loosely by William D. Clinger's paper "How to Read Floating
37 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
41 * 1. We only require IEEE, IBM, or VAX double-precision
42 * arithmetic (not IEEE double-extended).
43 * 2. We get by with floating-point arithmetic in a case that
44 * Clinger missed -- when we're computing d * 10^n
45 * for a small integer d and the integer n is not too
46 * much larger than 22 (the maximum integer k for which
47 * we can represent 10^k exactly), we may be able to
48 * compute (d*10^k) * 10^(e-k) with just one roundoff.
49 * 3. Rather than a bit-at-a-time adjustment of the binary
50 * result in the hard case, we use floating-point
51 * arithmetic to determine the adjustment to within
52 * one bit; only in really hard cases do we need to
53 * compute a second residual.
54 * 4. Because of 3., we don't need a large table of powers of 10
55 * for ten-to-e (just some small tables, e.g. of 10^k
60 * #define IEEE_8087 for IEEE-arithmetic machines where the least
61 * significant byte has the lowest address.
62 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
63 * significant byte has the lowest address.
64 * #define Sudden_Underflow for IEEE-format machines without gradual
65 * underflow (i.e., that flush to zero on underflow).
66 * #define IBM for IBM mainframe-style floating-point arithmetic.
67 * #define VAX for VAX-style floating-point arithmetic.
68 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
69 * #define No_leftright to omit left-right logic in fast floating-point
70 * computation of dtoa.
71 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
72 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
73 * that use extended-precision instructions to compute rounded
74 * products and quotients) with IBM.
75 * #define ROUND_BIASED for IEEE-format with biased rounding.
76 * #define Inaccurate_Divide for IEEE-format with correctly rounded
77 * products but inaccurate quotients, e.g., for Intel i860.
78 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
79 * integer arithmetic. Whether this speeds things up or slows things
80 * down depends on the machine and the number being converted.
83 /*#include <_ansi.h>*/
87 /* #include <reent.h> */
90 /* reent.c knows this value */
91 /* #define _Kmax 15 */
93 #define _reent _Jv_reent
94 #define _Bigint _Jv_Bigint
96 #define _REENT_CHECK_MP(x)
97 #define _REENT_MP_FREELIST(x) ((x)->_freelist)
98 #define _REENT_MP_P5S(x) ((x)->_p5s)
100 typedef unsigned long __ULong;
104 mprec_calloc (void *ignore, size_t x1, size_t x2)
106 char *result = (char *) malloc (x1 * x2);
107 memset (result, 0, x1 * x2);
112 _DEFUN (Balloc, (ptr, k), struct _reent *ptr _AND int k)
118 _REENT_CHECK_MP(ptr);
119 if (_REENT_MP_FREELIST(ptr) == NULL)
121 /* Allocate a list of pointers to the mprec objects */
122 _REENT_MP_FREELIST(ptr) = (struct _Bigint **) mprec_calloc (ptr,
123 sizeof (struct _Bigint *),
125 if (_REENT_MP_FREELIST(ptr) == NULL)
131 else if (new_k > ptr->_max_k)
133 struct _Bigint **new_list
134 = (struct _Bigint **) realloc (ptr->_freelist,
135 new_k * sizeof (struct _Bigint *));
136 memset (&new_list[ptr->_max_k], 0,
137 (new_k - ptr->_max_k) * sizeof (struct _Bigint *));
138 ptr->_freelist = new_list;
143 assert (k <= ptr->_max_k);
145 if ((rv = _REENT_MP_FREELIST(ptr)[k]) != 0)
147 _REENT_MP_FREELIST(ptr)[k] = rv->_next;
152 /* Allocate an mprec Bigint and stick in in the freelist */
153 rv = (_Bigint *) mprec_calloc (ptr,
156 (x-1) * sizeof(rv->_x));
157 if (rv == NULL) return NULL;
161 rv->_sign = rv->_wds = 0;
166 _DEFUN (Bfree, (ptr, v), struct _reent *ptr _AND _Bigint * v)
168 _REENT_CHECK_MP(ptr);
171 v->_next = _REENT_MP_FREELIST(ptr)[v->_k];
172 _REENT_MP_FREELIST(ptr)[v->_k] = v;
177 _DEFUN (multadd, (ptr, b, m, a),
178 struct _reent *ptr _AND
197 y = (xi & 0xffff) * m + a;
198 z = (xi >> 16) * m + (y >> 16);
200 *x++ = (z << 16) + (y & 0xffff);
210 if (wds >= b->_maxwds)
212 b1 = Balloc (ptr, b->_k + 1);
224 _DEFUN (s2b, (ptr, s, nd0, nd, y9),
225 struct _reent * ptr _AND
236 for (k = 0, y = 1; x > y; y <<= 1, k++);
242 b = Balloc (ptr, k + 1);
243 b->_x[0] = y9 & 0xffff;
244 b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1;
252 b = multadd (ptr, b, 10, *s++ - '0');
259 b = multadd (ptr, b, 10, *s++ - '0');
265 (x), register __ULong x)
269 if (!(x & 0xffff0000))
274 if (!(x & 0xff000000))
279 if (!(x & 0xf0000000))
284 if (!(x & 0xc0000000))
289 if (!(x & 0x80000000))
292 if (!(x & 0x40000000))
299 _DEFUN (lo0bits, (y), __ULong *y)
302 register __ULong x = *y;
349 _DEFUN (i2b, (ptr, i), struct _reent * ptr _AND int i)
360 _DEFUN (mult, (ptr, a, b), struct _reent * ptr _AND _Bigint * a _AND _Bigint * b)
365 __ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
370 if (a->_wds < b->_wds)
383 for (x = c->_x, xa = x + wc; x < xa; x++)
391 for (; xb < xbe; xb++, xc0++)
393 if ((y = *xb & 0xffff) != 0)
400 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
402 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
404 Storeinc (xc, z2, z);
409 if ((y = *xb >> 16) != 0)
417 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
419 Storeinc (xc, z, z2);
420 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
428 for (; xb < xbe; xc0++)
437 z = *x++ * y + *xc + carry;
446 for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc);
453 (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k)
455 _Bigint *b1, *p5, *p51;
457 static _CONST int p05[3] = {5, 25, 125};
459 if ((i = k & 3) != 0)
460 b = multadd (ptr, b, p05[i - 1], 0);
464 _REENT_CHECK_MP(ptr);
465 if (!(p5 = _REENT_MP_P5S(ptr)))
468 p5 = _REENT_MP_P5S(ptr) = i2b (ptr, 625);
475 b1 = mult (ptr, b, p5);
481 if (!(p51 = p5->_next))
483 p51 = p5->_next = mult (ptr, p5, p5);
492 _DEFUN (lshift, (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k)
496 __ULong *x, *x1, *xe, z;
504 n1 = n + b->_wds + 1;
505 for (i = b->_maxwds; n1 > i; i <<= 1)
507 b1 = Balloc (ptr, k1);
509 for (i = 0; i < n; i++)
534 *x1++ = *x << k & 0xffff | z;
552 _DEFUN (cmp, (a, b), _Bigint * a _AND _Bigint * b)
554 __ULong *xa, *xa0, *xb, *xb0;
560 if (i > 1 && !a->_x[i - 1])
561 Bug ("cmp called with a->_x[a->_wds-1] == 0");
562 if (j > 1 && !b->_x[j - 1])
563 Bug ("cmp called with b->_x[b->_wds-1] == 0");
574 return *xa < *xb ? -1 : 1;
582 _DEFUN (diff, (ptr, a, b), struct _reent * ptr _AND
583 _Bigint * a _AND _Bigint * b)
587 __Long borrow, y; /* We need signed shifts here. */
588 __ULong *xa, *xae, *xb, *xbe, *xc;
610 c = Balloc (ptr, a->_k);
623 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
625 Sign_Extend (borrow, y);
626 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
628 Sign_Extend (borrow, z);
634 y = (*xa & 0xffff) + borrow;
636 Sign_Extend (borrow, y);
637 z = (*xa++ >> 16) + borrow;
639 Sign_Extend (borrow, z);
645 y = *xa++ - *xb++ + borrow;
647 Sign_Extend (borrow, y);
655 Sign_Extend (borrow, y);
666 _DEFUN (ulp, (_x), double _x)
668 union double_union x, a;
673 L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1;
674 #ifndef Sudden_Underflow
682 #ifndef _DOUBLE_IS_32BITS
686 #ifndef Sudden_Underflow
693 word0 (a) = 0x80000 >> L;
694 #ifndef _DOUBLE_IS_32BITS
702 #ifndef _DOUBLE_IS_32BITS
703 word1 (a) = L >= 31 ? 1 : 1 << (31 - L);
713 _Bigint * a _AND int *e)
715 __ULong *xa, *xa0, w, y, z;
717 union double_union d;
730 Bug ("zero y in b2d");
737 d0 = Exp_1 | y >> (Ebits - k);
738 w = xa > xa0 ? *--xa : 0;
739 #ifndef _DOUBLE_IS_32BITS
740 d1 = y << ((32 - Ebits) + k) | w >> (Ebits - k);
744 z = xa > xa0 ? *--xa : 0;
747 d0 = Exp_1 | y << k | z >> (32 - k);
748 y = xa > xa0 ? *--xa : 0;
749 #ifndef _DOUBLE_IS_32BITS
750 d1 = z << k | y >> (32 - k);
756 #ifndef _DOUBLE_IS_32BITS
763 z = xa > xa0 ? *--xa : 0;
764 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
765 w = xa > xa0 ? *--xa : 0;
766 y = xa > xa0 ? *--xa : 0;
767 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
770 z = xa > xa0 ? *--xa : 0;
771 w = xa > xa0 ? *--xa : 0;
773 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
774 y = xa > xa0 ? *--xa : 0;
775 d1 = w << k + 16 | y << k;
779 word0 (d) = d0 >> 16 | d0 << 16;
780 word1 (d) = d1 >> 16 | d1 << 16;
791 struct _reent * ptr _AND
797 union double_union d;
806 d0 = word0 (d) >> 16 | word0 (d) << 16;
807 d1 = word1 (d) >> 16 | word1 (d) << 16;
822 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
823 #ifdef Sudden_Underflow
824 de = (int) (d0 >> Exp_shift);
829 if ((de = (int) (d0 >> Exp_shift)) != 0)
833 #ifndef _DOUBLE_IS_32BITS
840 x[0] = y | z << (32 - k);
845 i = b->_wds = (x[1] = z) ? 2 : 1;
852 Bug ("Zero passed to d2b");
857 #ifndef _DOUBLE_IS_32BITS
869 x[0] = y | z << 32 - k & 0xffff;
870 x[1] = z >> k - 16 & 0xffff;
877 x[1] = y >> 16 | z << 16 - k & 0xffff;
878 x[2] = z >> k & 0xffff;
895 Bug ("Zero passed to d2b");
915 #ifndef Sudden_Underflow
920 *e = (de - Bias - (P - 1) << 2) + k;
921 *bits = 4 * P + 8 - k - hi0bits (word0 (d) & Frac_mask);
923 *e = de - Bias - (P - 1) + k;
926 #ifndef Sudden_Underflow
930 *e = de - Bias - (P - 1) + 1 + k;
932 *bits = 32 * i - hi0bits (x[i - 1]);
934 *bits = (i + 2) * 16 - hi0bits (x[i]);
944 _DEFUN (ratio, (a, b), _Bigint * a _AND _Bigint * b)
947 union double_union da, db;
953 k = ka - kb + 32 * (a->_wds - b->_wds);
955 k = ka - kb + 16 * (a->_wds - b->_wds);
960 word0 (da) += (k >> 2) * Exp_msk1;
967 word0 (db) += (k >> 2) * Exp_msk1;
973 word0 (da) += k * Exp_msk1;
977 word0 (db) += k * Exp_msk1;
987 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
988 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
989 1e20, 1e21, 1e22, 1e23, 1e24
993 #if !defined(_DOUBLE_IS_32BITS) && !defined(__v800)
994 _CONST double bigtens[] =
995 {1e16, 1e32, 1e64, 1e128, 1e256};
997 _CONST double tinytens[] =
998 {1e-16, 1e-32, 1e-64, 1e-128, 1e-256};
1000 _CONST double bigtens[] =
1003 _CONST double tinytens[] =
1009 _DEFUN (_mprec_log10, (dig),