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 * $NetBSD: strtod.c,v 1.26 1998/02/03 18:44:21 perry Exp $
21 * $FreeBSD: src/lib/libc/stdlib/netbsd_strtod.c,v 1.2.2.2 2001/03/02 17:14:15 tegge Exp $
22 * $DragonFly: src/lib/libc/stdlib/Attic/netbsd_strtod.c,v 1.3 2003/09/06 08:19:16 asmodai Exp $
26 /* Please send bug reports to
28 AT&T Bell Laboratories, Room 2C-463
30 Murray Hill, NJ 07974-2070
32 dmg@research.att.com or research!dmg
35 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
37 * This strtod returns a nearest machine number to the input decimal
38 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
39 * broken by the IEEE round-even rule. Otherwise ties are broken by
40 * biased rounding (add half and chop).
42 * Inspired loosely by William D. Clinger's paper "How to Read Floating
43 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
47 * 1. We only require IEEE, IBM, or VAX double-precision
48 * arithmetic (not IEEE double-extended).
49 * 2. We get by with floating-point arithmetic in a case that
50 * Clinger missed -- when we're computing d * 10^n
51 * for a small integer d and the integer n is not too
52 * much larger than 22 (the maximum integer k for which
53 * we can represent 10^k exactly), we may be able to
54 * compute (d*10^k) * 10^(e-k) with just one roundoff.
55 * 3. Rather than a bit-at-a-time adjustment of the binary
56 * result in the hard case, we use floating-point
57 * arithmetic to determine the adjustment to within
58 * one bit; only in really hard cases do we need to
59 * compute a second residual.
60 * 4. Because of 3., we don't need a large table of powers of 10
61 * for ten-to-e (just some small tables, e.g. of 10^k
66 * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least
67 * significant byte has the lowest address.
68 * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most
69 * significant byte has the lowest address.
70 * #define Long int on machines with 32-bit ints and 64-bit longs.
71 * #define Sudden_Underflow for IEEE-format machines without gradual
72 * underflow (i.e., that flush to zero on underflow).
73 * #define IBM for IBM mainframe-style floating-point arithmetic.
74 * #define VAX for VAX-style floating-point arithmetic.
75 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
76 * #define No_leftright to omit left-right logic in fast floating-point
77 * computation of dtoa.
78 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
79 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
80 * that use extended-precision instructions to compute rounded
81 * products and quotients) with IBM.
82 * #define ROUND_BIASED for IEEE-format with biased rounding.
83 * #define Inaccurate_Divide for IEEE-format with correctly rounded
84 * products but inaccurate quotients, e.g., for Intel i860.
85 * #define Just_16 to store 16 bits per 32-bit Long when doing high-precision
86 * integer arithmetic. Whether this speeds things up or slows things
87 * down depends on the machine and the number being converted.
88 * #define KR_headers for old-style C function headers.
89 * #define Bad_float_h if your system lacks a float.h or if it does not
90 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
91 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
92 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
93 * if memory is available and otherwise does something you deem
94 * appropriate. If MALLOC is undefined, malloc will be invoked
95 * directly -- and assumed always to succeed.
98 #include <sys/cdefs.h>
100 #if defined(__m68k__) || defined(__sparc__) || defined(__i386__) || \
101 defined(__mips__) || defined(__ns32k__) || defined(__alpha__) || \
103 #include <sys/types.h>
104 #if BYTE_ORDER == BIG_ENDIAN
105 #define IEEE_BIG_ENDIAN
107 #define IEEE_LITTLE_ENDIAN
113 * Although the CPU is little endian the FP has different
114 * byte and word endianness. The byte order is still little endian
115 * but the word order is big endian.
117 #define IEEE_BIG_ENDIAN
125 #define ULong u_int32_t
129 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
145 char *__dtoa (double, int, int, int *, int *, char **, char **);
149 extern char *MALLOC();
151 extern void *MALLOC(size_t);
154 #define MALLOC malloc
162 #ifdef IEEE_BIG_ENDIAN
163 #define IEEE_ARITHMETIC
165 #ifdef IEEE_LITTLE_ENDIAN
166 #define IEEE_ARITHMETIC
169 #ifdef IEEE_ARITHMETIC
171 #define DBL_MAX_10_EXP 308
172 #define DBL_MAX_EXP 1024
175 #define DBL_MAX 1.7976931348623157e+308
180 #define DBL_MAX_10_EXP 75
181 #define DBL_MAX_EXP 63
184 #define DBL_MAX 7.2370055773322621e+75
189 #define DBL_MAX_10_EXP 38
190 #define DBL_MAX_EXP 127
193 #define DBL_MAX 1.7014118346046923e+38
197 #define LONG_MAX 2147483647
212 #define CONST /* blank */
218 #ifdef Unsigned_Shifts
219 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
221 #define Sign_Extend(a,b) /*no-op*/
224 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + \
226 Exactly one of IEEE_LITTLE_ENDIAN IEEE_BIG_ENDIAN, VAX, or
227 IBM should be defined.
230 #ifdef IEEE_LITTLE_ENDIAN
231 #define word0(x) ((ULong *)&x)[1]
232 #define word1(x) ((ULong *)&x)[0]
234 #define word0(x) ((ULong *)&x)[0]
235 #define word1(x) ((ULong *)&x)[1]
238 /* The following definition of Storeinc is appropriate for MIPS processors.
239 * An alternative that might be better on some machines is
240 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
242 #if defined(IEEE_LITTLE_ENDIAN) + defined(VAX) + defined(__arm32__)
243 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
244 ((unsigned short *)a)[0] = (unsigned short)c, a++)
246 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
247 ((unsigned short *)a)[1] = (unsigned short)c, a++)
250 /* #define P DBL_MANT_DIG */
251 /* Ten_pmax = floor(P*log(2)/log(5)) */
252 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
253 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
254 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
256 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN)
258 #define Exp_shift1 20
259 #define Exp_msk1 0x100000
260 #define Exp_msk11 0x100000
261 #define Exp_mask 0x7ff00000
266 #define Exp_1 0x3ff00000
267 #define Exp_11 0x3ff00000
269 #define Frac_mask 0xfffff
270 #define Frac_mask1 0xfffff
273 #define Bndry_mask 0xfffff
274 #define Bndry_mask1 0xfffff
276 #define Sign_bit 0x80000000
282 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
284 #undef Sudden_Underflow
285 #define Sudden_Underflow
288 #define Exp_shift1 24
289 #define Exp_msk1 0x1000000
290 #define Exp_msk11 0x1000000
291 #define Exp_mask 0x7f000000
294 #define Exp_1 0x41000000
295 #define Exp_11 0x41000000
296 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
297 #define Frac_mask 0xffffff
298 #define Frac_mask1 0xffffff
301 #define Bndry_mask 0xefffff
302 #define Bndry_mask1 0xffffff
304 #define Sign_bit 0x80000000
306 #define Tiny0 0x100000
313 #define Exp_msk1 0x80
314 #define Exp_msk11 0x800000
315 #define Exp_mask 0x7f80
318 #define Exp_1 0x40800000
319 #define Exp_11 0x4080
321 #define Frac_mask 0x7fffff
322 #define Frac_mask1 0xffff007f
325 #define Bndry_mask 0xffff007f
326 #define Bndry_mask1 0xffff007f
328 #define Sign_bit 0x8000
342 #define rounded_product(a,b) a = rnd_prod(a, b)
343 #define rounded_quotient(a,b) a = rnd_quot(a, b)
345 extern double rnd_prod(), rnd_quot();
347 extern double rnd_prod(double, double), rnd_quot(double, double);
350 #define rounded_product(a,b) a *= b
351 #define rounded_quotient(a,b) a /= b
354 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
355 #define Big1 0xffffffff
358 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
359 * This makes some inner loops simpler and sometimes saves work
360 * during multiplications, but it often seems to make things slightly
361 * slower. Hence the default is now to store 32 bits per Long.
371 extern "C" double strtod(const char *s00, char **se);
372 extern "C" char *__dtoa(double d, int mode, int ndigits,
373 int *decpt, int *sign, char **rve, char **resultp);
379 int k, maxwds, sign, wds;
383 typedef struct Bigint Bigint;
397 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(Long));
400 rv->sign = rv->wds = 0;
415 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
416 y->wds*sizeof(Long) + 2*sizeof(int))
421 (b, m, a) Bigint *b; int m, a;
423 (Bigint *b, int m, int a) /* multiply by m and add a */
439 y = (xi & 0xffff) * m + a;
440 z = (xi >> 16) * m + (y >> 16);
442 *x++ = (z << 16) + (y & 0xffff);
451 if (wds >= b->maxwds) {
466 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
468 (CONST char *s, int nd0, int nd, ULong y9)
476 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
483 b->x[0] = y9 & 0xffff;
484 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
490 do b = multadd(b, 10, *s++ - '0');
497 b = multadd(b, 10, *s++ - '0');
511 if (!(x & 0xffff0000)) {
515 if (!(x & 0xff000000)) {
519 if (!(x & 0xf0000000)) {
523 if (!(x & 0xc0000000)) {
527 if (!(x & 0x80000000)) {
529 if (!(x & 0x40000000))
602 (a, b) Bigint *a, *b;
604 (Bigint *a, Bigint *b)
610 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
615 if (a->wds < b->wds) {
627 for(x = c->x, xa = x + wc; x < xa; x++)
635 for(; xb < xbe; xb++, xc0++) {
636 if ((y = *xb & 0xffff) != 0) {
641 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
643 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
650 if ((y = *xb >> 16) != 0) {
656 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
659 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
667 for(; xb < xbe; xc0++) {
673 z = *x++ * y + *xc + carry;
682 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
692 (b, k) Bigint *b; int k;
697 Bigint *b1, *p5, *p51;
699 static int p05[3] = { 5, 25, 125 };
701 if ((i = k & 3) != 0)
702 b = multadd(b, p05[i-1], 0);
719 if (!(p51 = p5->next)) {
720 p51 = p5->next = mult(p5,p5);
731 (b, k) Bigint *b; int k;
738 ULong *x, *x1, *xe, z;
747 for(i = b->maxwds; n1 > i; i <<= 1)
751 for(i = 0; i < n; i++)
772 *x1++ = *x << k & 0xffff | z;
791 (a, b) Bigint *a, *b;
793 (Bigint *a, Bigint *b)
796 ULong *xa, *xa0, *xb, *xb0;
802 if (i > 1 && !a->x[i-1])
803 Bug("cmp called with a->x[a->wds-1] == 0");
804 if (j > 1 && !b->x[j-1])
805 Bug("cmp called with b->x[b->wds-1] == 0");
815 return *xa < *xb ? -1 : 1;
825 (a, b) Bigint *a, *b;
827 (Bigint *a, Bigint *b)
832 Long borrow, y; /* We need signed shifts here. */
833 ULong *xa, *xae, *xb, *xbe, *xc;
865 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
867 Sign_Extend(borrow, y);
868 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
870 Sign_Extend(borrow, z);
875 y = (*xa & 0xffff) + borrow;
877 Sign_Extend(borrow, y);
878 z = (*xa++ >> 16) + borrow;
880 Sign_Extend(borrow, z);
885 y = *xa++ - *xb++ + borrow;
887 Sign_Extend(borrow, y);
894 Sign_Extend(borrow, y);
915 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
916 #ifndef Sudden_Underflow
924 #ifndef Sudden_Underflow
929 word0(a) = 0x80000 >> L;
935 word1(a) = L >= 31 ? 1 : 1 << (31 - L);
945 (a, e) Bigint *a; int *e;
950 ULong *xa, *xa0, w, y, z;
964 if (!y) Bug("zero y in b2d");
970 d0 = Exp_1 | y >> (Ebits - k);
971 w = xa > xa0 ? *--xa : 0;
972 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
975 z = xa > xa0 ? *--xa : 0;
977 d0 = Exp_1 | y << k | z >> (32 - k);
978 y = xa > xa0 ? *--xa : 0;
979 d1 = z << k | y >> (32 - k);
986 if (k < Ebits + 16) {
987 z = xa > xa0 ? *--xa : 0;
988 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
989 w = xa > xa0 ? *--xa : 0;
990 y = xa > xa0 ? *--xa : 0;
991 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
994 z = xa > xa0 ? *--xa : 0;
995 w = xa > xa0 ? *--xa : 0;
997 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
998 y = xa > xa0 ? *--xa : 0;
999 d1 = w << k + 16 | y << k;
1003 word0(d) = d0 >> 16 | d0 << 16;
1004 word1(d) = d1 >> 16 | d1 << 16;
1015 (d, e, bits) double d; int *e, *bits;
1017 (double d, int *e, int *bits)
1025 d0 = word0(d) >> 16 | word0(d) << 16;
1026 d1 = word1(d) >> 16 | word1(d) << 16;
1040 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1041 #ifdef Sudden_Underflow
1042 de = (int)(d0 >> Exp_shift);
1047 if ((de = (int)(d0 >> Exp_shift)) != 0)
1051 if ((y = d1) != 0) {
1052 if ((k = lo0bits(&y)) != 0) {
1053 x[0] = y | z << (32 - k);
1058 i = b->wds = (x[1] = z) ? 2 : 1;
1063 Bug("Zero passed to d2b");
1072 if (k = lo0bits(&y))
1074 x[0] = y | z << 32 - k & 0xffff;
1075 x[1] = z >> k - 16 & 0xffff;
1081 x[1] = y >> 16 | z << 16 - k & 0xffff;
1082 x[2] = z >> k & 0xffff;
1097 Bug("Zero passed to d2b");
1115 #ifndef Sudden_Underflow
1119 *e = (de - Bias - (P-1) << 2) + k;
1120 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1122 *e = de - Bias - (P-1) + k;
1125 #ifndef Sudden_Underflow
1128 *e = de - Bias - (P-1) + 1 + k;
1130 *bits = 32*i - hi0bits(x[i-1]);
1132 *bits = (i+2)*16 - hi0bits(x[i]);
1144 (a, b) Bigint *a, *b;
1146 (Bigint *a, Bigint *b)
1155 k = ka - kb + 32*(a->wds - b->wds);
1157 k = ka - kb + 16*(a->wds - b->wds);
1161 word0(da) += (k >> 2)*Exp_msk1;
1167 word0(db) += (k >> 2)*Exp_msk1;
1173 word0(da) += k*Exp_msk1;
1176 word0(db) += k*Exp_msk1;
1184 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1185 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1193 static CONST double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1194 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1198 static CONST double bigtens[] = { 1e16, 1e32, 1e64 };
1199 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1202 static CONST double bigtens[] = { 1e16, 1e32 };
1203 static CONST double tinytens[] = { 1e-16, 1e-32 };
1211 (s00, se) CONST char *s00; char **se;
1213 (CONST char *s00, char **se)
1216 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1217 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1218 CONST char *s, *s0, *s1;
1219 double aadj, aadj1, adj, rv, rv0;
1223 Bigint *bb = NULL, *bd = NULL, *bs = NULL, *delta = NULL;/* pacify gcc */
1226 CONST char decimal_point = localeconv()->decimal_point[0];
1228 CONST char decimal_point = '.';
1231 sign = nz0 = nz = 0;
1235 for(s = s00; isspace((unsigned char) *s); s++)
1241 } else if (*s == '+') {
1252 while(*++s == '0') ;
1258 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1264 if (c == decimal_point) {
1267 for(; c == '0'; c = *++s)
1269 if (c > '0' && c <= '9') {
1277 for(; c >= '0' && c <= '9'; c = *++s) {
1282 for(i = 1; i < nz; i++)
1285 else if (nd <= DBL_DIG + 1)
1289 else if (nd <= DBL_DIG + 1)
1297 if (c == 'e' || c == 'E') {
1298 if (!nd && !nz && !nz0) {
1310 if (c >= '0' && c <= '9') {
1313 if (c > '0' && c <= '9') {
1316 while((c = *++s) >= '0' && c <= '9')
1318 if (s - s1 > 8 || L > 19999)
1319 /* Avoid confusion from exponents
1320 * so large that e might overflow.
1322 e = 19999; /* safe for 16 bit ints */
1341 /* Now we have nd0 digits, starting at s0, followed by a
1342 * decimal point, followed by nd-nd0 digits. The number we're
1343 * after is the integer represented by those digits times
1348 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1351 rv = tens[k - 9] * rv + z;
1354 #ifndef RND_PRODQUOT
1361 if (e <= Ten_pmax) {
1363 goto vax_ovfl_check;
1365 /* rv = */ rounded_product(rv, tens[e]);
1370 if (e <= Ten_pmax + i) {
1371 /* A fancier test would sometimes let us do
1372 * this for larger i values.
1377 /* VAX exponent range is so narrow we must
1378 * worry about overflow here...
1381 word0(rv) -= P*Exp_msk1;
1382 /* rv = */ rounded_product(rv, tens[e]);
1383 if ((word0(rv) & Exp_mask)
1384 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1386 word0(rv) += P*Exp_msk1;
1388 /* rv = */ rounded_product(rv, tens[e]);
1393 #ifndef Inaccurate_Divide
1394 else if (e >= -Ten_pmax) {
1395 /* rv = */ rounded_quotient(rv, tens[-e]);
1402 /* Get starting approximation = rv * 10**e1 */
1405 if ((i = e1 & 15) != 0)
1408 if (e1 > DBL_MAX_10_EXP) {
1414 /* Can't trust HUGE_VAL */
1416 word0(rv) = Exp_mask;
1428 for(j = 0; e1 > 1; j++, e1 >>= 1)
1431 /* The last multiplication could overflow. */
1432 word0(rv) -= P*Exp_msk1;
1434 if ((z = word0(rv) & Exp_mask)
1435 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1437 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1438 /* set to largest number */
1439 /* (Can't trust DBL_MAX) */
1444 word0(rv) += P*Exp_msk1;
1451 if ((i = e1 & 15) != 0)
1455 if (e1 >= 1 << n_bigtens)
1457 for(j = 0; e1 > 1; j++, e1 >>= 1)
1460 /* The last multiplication could underflow. */
1476 /* The refinement below will clean
1477 * this approximation up.
1483 /* Now the hard part -- adjusting rv to the correct value.*/
1485 /* Put digits into bd: true value = bd * 10^e */
1487 bd0 = s2b(s0, nd0, nd, y);
1490 bd = Balloc(bd0->k);
1492 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1508 #ifdef Sudden_Underflow
1510 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1515 i = bbe + bbbits - 1; /* logb(rv) */
1516 if (i < Emin) /* denormal */
1523 i = bb2 < bd2 ? bb2 : bd2;
1532 bs = pow5mult(bs, bb5);
1538 bb = lshift(bb, bb2);
1540 bd = pow5mult(bd, bd5);
1542 bd = lshift(bd, bd2);
1544 bs = lshift(bs, bs2);
1545 delta = diff(bb, bd);
1546 dsign = delta->sign;
1550 /* Error is less than half an ulp -- check for
1551 * special case of mantissa a power of two.
1553 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1555 delta = lshift(delta,Log2P);
1556 if (cmp(delta, bs) > 0)
1561 /* exactly half-way between */
1563 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1564 && word1(rv) == 0xffffffff) {
1565 /*boundary case -- increment exponent*/
1566 word0(rv) = (word0(rv) & Exp_mask)
1576 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1578 /* boundary case -- decrement exponent */
1579 #ifdef Sudden_Underflow
1580 L = word0(rv) & Exp_mask;
1589 L = (word0(rv) & Exp_mask) - Exp_msk1;
1591 word0(rv) = L | Bndry_mask1;
1592 word1(rv) = 0xffffffff;
1599 #ifndef ROUND_BIASED
1600 if (!(word1(rv) & LSB))
1605 #ifndef ROUND_BIASED
1608 #ifndef Sudden_Underflow
1616 if ((aadj = ratio(delta, bs)) <= 2.) {
1619 else if (word1(rv) || word0(rv) & Bndry_mask) {
1620 #ifndef Sudden_Underflow
1621 if (word1(rv) == Tiny1 && !word0(rv))
1628 /* special case -- power of FLT_RADIX to be */
1629 /* rounded down... */
1631 if (aadj < 2./FLT_RADIX)
1632 aadj = 1./FLT_RADIX;
1640 aadj1 = dsign ? aadj : -aadj;
1641 #ifdef Check_FLT_ROUNDS
1642 switch(FLT_ROUNDS) {
1643 case 2: /* towards +infinity */
1646 case 0: /* towards 0 */
1647 case 3: /* towards -infinity */
1651 if (FLT_ROUNDS == 0)
1655 y = word0(rv) & Exp_mask;
1657 /* Check for overflow */
1659 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1661 word0(rv) -= P*Exp_msk1;
1662 adj = aadj1 * ulp(rv);
1664 if ((word0(rv) & Exp_mask) >=
1665 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1666 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1673 word0(rv) += P*Exp_msk1;
1676 #ifdef Sudden_Underflow
1677 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1679 word0(rv) += P*Exp_msk1;
1680 adj = aadj1 * ulp(rv);
1683 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1685 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1688 if (word0(rv0) == Tiny0
1689 && word1(rv0) == Tiny1)
1696 word0(rv) -= P*Exp_msk1;
1699 adj = aadj1 * ulp(rv);
1703 /* Compute adj so that the IEEE rounding rules will
1704 * correctly round rv + adj in some half-way cases.
1705 * If rv * ulp(rv) is denormalized (i.e.,
1706 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1707 * trouble from bits lost to denormalization;
1708 * example: 1.2e-307 .
1710 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1711 aadj1 = (double)(int)(aadj + 0.5);
1715 adj = aadj1 * ulp(rv);
1719 z = word0(rv) & Exp_mask;
1721 /* Can we stop now? */
1724 /* The tolerances below are conservative. */
1725 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1726 if (aadj < .4999999 || aadj > .5000001)
1729 else if (aadj < .4999999/FLT_RADIX)
1747 return sign ? -rv : rv;
1753 (b, S) Bigint *b, *S;
1755 (Bigint *b, Bigint *S)
1761 ULong *bx, *bxe, *sx, *sxe;
1769 /*debug*/ if (b->wds > n)
1770 /*debug*/ Bug("oversize b in quorem");
1778 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1780 /*debug*/ if (q > 9)
1781 /*debug*/ Bug("oversized quotient in quorem");
1789 ys = (si & 0xffff) * q + carry;
1790 zs = (si >> 16) * q + (ys >> 16);
1792 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1794 Sign_Extend(borrow, y);
1795 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1797 Sign_Extend(borrow, z);
1800 ys = *sx++ * q + carry;
1802 y = *bx - (ys & 0xffff) + borrow;
1804 Sign_Extend(borrow, y);
1811 while(--bxe > bx && !*bxe)
1816 if (cmp(b, S) >= 0) {
1825 ys = (si & 0xffff) + carry;
1826 zs = (si >> 16) + (ys >> 16);
1828 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1830 Sign_Extend(borrow, y);
1831 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1833 Sign_Extend(borrow, z);
1838 y = *bx - (ys & 0xffff) + borrow;
1840 Sign_Extend(borrow, y);
1848 while(--bxe > bx && !*bxe)
1856 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1858 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1859 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1862 * 1. Rather than iterating, we use a simple numeric overestimate
1863 * to determine k = floor(log10(d)). We scale relevant
1864 * quantities using O(log2(k)) rather than O(k) multiplications.
1865 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1866 * try to generate digits strictly left to right. Instead, we
1867 * compute with fewer bits and propagate the carry if necessary
1868 * when rounding the final digit up. This is often faster.
1869 * 3. Under the assumption that input will be rounded nearest,
1870 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1871 * That is, we allow equality in stopping tests when the
1872 * round-nearest rule will give the same floating-point value
1873 * as would satisfaction of the stopping test with strict
1875 * 4. We remove common factors of powers of 2 from relevant
1877 * 5. When converting floating-point integers less than 1e16,
1878 * we use floating-point arithmetic rather than resorting
1879 * to multiple-precision integers.
1880 * 6. When asked to produce fewer than 15 digits, we first try
1881 * to get by with floating-point arithmetic; we resort to
1882 * multiple-precision integer arithmetic only if we cannot
1883 * guarantee that the floating-point calculation has given
1884 * the correctly rounded result. For k requested digits and
1885 * "uniformly" distributed input, the probability is
1886 * something like 10^(k-15) that we must resort to the Long
1893 (d, mode, ndigits, decpt, sign, rve, resultp)
1894 double d; int mode, ndigits, *decpt, *sign; char **rve, **resultp;
1896 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve,
1900 /* Arguments ndigits, decpt, sign are similar to those
1901 of ecvt and fcvt; trailing zeros are suppressed from
1902 the returned string. If not null, *rve is set to point
1903 to the end of the return value. If d is +-Infinity or NaN,
1904 then *decpt is set to 9999.
1907 0 ==> shortest string that yields d when read in
1908 and rounded to nearest.
1909 1 ==> like 0, but with Steele & White stopping rule;
1910 e.g. with IEEE P754 arithmetic , mode 0 gives
1911 1e23 whereas mode 1 gives 9.999999999999999e22.
1912 2 ==> max(1,ndigits) significant digits. This gives a
1913 return value similar to that of ecvt, except
1914 that trailing zeros are suppressed.
1915 3 ==> through ndigits past the decimal point. This
1916 gives a return value similar to that from fcvt,
1917 except that trailing zeros are suppressed, and
1918 ndigits can be negative.
1919 4-9 should give the same return values as 2-3, i.e.,
1920 4 <= mode <= 9 ==> same return as mode
1921 2 + (mode & 1). These modes are mainly for
1922 debugging; often they run slower but sometimes
1923 faster than modes 2-3.
1924 4,5,8,9 ==> left-to-right digit generation.
1925 6-9 ==> don't try fast floating-point estimate
1928 Values of mode other than 0-9 are treated as mode 0.
1930 Sufficient space is allocated to the return value
1931 to hold the suppressed trailing zeros.
1934 int bbits, b2, b5, be, dig, i, ieps, ilim0,
1935 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1937 int ilim = 0, ilim1 = 0, spec_case = 0; /* pacify gcc */
1939 #ifndef Sudden_Underflow
1943 Bigint *b, *b1, *delta, *mhi, *S;
1944 Bigint *mlo = NULL; /* pacify gcc */
1948 if (word0(d) & Sign_bit) {
1949 /* set sign for everything, including 0's and NaNs */
1951 word0(d) &= ~Sign_bit; /* clear sign bit */
1956 #if defined(IEEE_Arith) + defined(VAX)
1958 if ((word0(d) & Exp_mask) == Exp_mask)
1960 if (word0(d) == 0x8000)
1963 /* Infinity or NaN */
1967 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
1980 d += 0; /* normalize */
1990 b = d2b(d, &be, &bbits);
1991 #ifdef Sudden_Underflow
1992 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1994 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) {
1997 word0(d2) &= Frac_mask1;
1998 word0(d2) |= Exp_11;
2000 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
2004 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2005 * log10(x) = log(x) / log(10)
2006 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2007 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2009 * This suggests computing an approximation k to log10(d) by
2011 * k = (i - Bias)*0.301029995663981
2012 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2014 * We want k to be too large rather than too small.
2015 * The error in the first-order Taylor series approximation
2016 * is in our favor, so we just round up the constant enough
2017 * to compensate for any error in the multiplication of
2018 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2019 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2020 * adding 1e-13 to the constant term more than suffices.
2021 * Hence we adjust the constant term to 0.1760912590558.
2022 * (We could get a more accurate k by invoking log10,
2023 * but this is probably not worthwhile.)
2031 #ifndef Sudden_Underflow
2035 /* d is denormalized */
2037 i = bbits + be + (Bias + (P-1) - 1);
2038 x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
2039 : word1(d) << (32 - i);
2041 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2042 i -= (Bias + (P-1) - 1) + 1;
2046 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2048 if (ds < 0. && ds != k)
2049 k--; /* want k = floor(ds) */
2051 if (k >= 0 && k <= Ten_pmax) {
2075 if (mode < 0 || mode > 9)
2096 ilim = ilim1 = i = ndigits;
2102 i = ndigits + k + 1;
2108 *resultp = (char *) malloc(i + 1);
2111 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2113 /* Try to get by with floating-point arithmetic. */
2119 ieps = 2; /* conservative */
2124 /* prevent overflows */
2126 d /= bigtens[n_bigtens-1];
2129 for(; j; j >>= 1, i++)
2136 else if ((j1 = -k) != 0) {
2137 d *= tens[j1 & 0xf];
2138 for(j = j1 >> 4; j; j >>= 1, i++)
2144 if (k_check && d < 1. && ilim > 0) {
2153 word0(eps) -= (P-1)*Exp_msk1;
2163 #ifndef No_leftright
2165 /* Use Steele & White method of only
2166 * generating digits needed.
2168 eps = 0.5/tens[ilim-1] - eps;
2172 *s++ = '0' + (int)L;
2185 /* Generate ilim digits, then fix them up. */
2186 eps *= tens[ilim-1];
2187 for(i = 1;; i++, d *= 10.) {
2190 *s++ = '0' + (int)L;
2194 else if (d < 0.5 - eps) {
2202 #ifndef No_leftright
2212 /* Do we have a "small" integer? */
2214 if (be >= 0 && k <= Int_max) {
2217 if (ndigits < 0 && ilim <= 0) {
2219 if (ilim < 0 || d <= 5*ds)
2226 #ifdef Check_FLT_ROUNDS
2227 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2233 *s++ = '0' + (int)L;
2236 if (d > ds || (d == ds && L & 1)) {
2260 #ifndef Sudden_Underflow
2261 denorm ? be + (Bias + (P-1) - 1 + 1) :
2264 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2278 if ((i = ilim) < 0) {
2287 if (m2 > 0 && s2 > 0) {
2288 i = m2 < s2 ? m2 : s2;
2296 mhi = pow5mult(mhi, m5);
2301 if ((j = b5 - m5) != 0)
2305 b = pow5mult(b, b5);
2309 S = pow5mult(S, s5);
2311 /* Check for special case that d is a normalized power of 2. */
2314 if (!word1(d) && !(word0(d) & Bndry_mask)
2315 #ifndef Sudden_Underflow
2316 && word0(d) & Exp_mask
2319 /* The special case */
2328 /* Arrange for convenient computation of quotients:
2329 * shift left if necessary so divisor has 4 leading 0 bits.
2331 * Perhaps we should just compute leading 28 bits of S once
2332 * and for all and pass them and a shift to quorem, so it
2333 * can do shifts and ors to compute the numerator for q.
2336 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0)
2339 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
2361 b = multadd(b, 10, 0); /* we botched the k estimate */
2363 mhi = multadd(mhi, 10, 0);
2367 if (ilim <= 0 && mode > 2) {
2368 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2369 /* no digits, fcvt style */
2381 mhi = lshift(mhi, m2);
2383 /* Compute mlo -- check for special case
2384 * that d is a normalized power of 2.
2389 mhi = Balloc(mhi->k);
2391 mhi = lshift(mhi, Log2P);
2395 dig = quorem(b,S) + '0';
2396 /* Do we yet have the shortest decimal string
2397 * that will round to d?
2400 delta = diff(S, mhi);
2401 j1 = delta->sign ? 1 : cmp(b, delta);
2403 #ifndef ROUND_BIASED
2404 if (j1 == 0 && !mode && !(word1(d) & 1)) {
2413 if (j < 0 || (j == 0 && !mode
2414 #ifndef ROUND_BIASED
2421 if ((j1 > 0 || (j1 == 0 && dig & 1))
2429 if (dig == '9') { /* possible if i == 1 */
2440 b = multadd(b, 10, 0);
2442 mlo = mhi = multadd(mhi, 10, 0);
2444 mlo = multadd(mlo, 10, 0);
2445 mhi = multadd(mhi, 10, 0);
2451 *s++ = dig = quorem(b,S) + '0';
2454 b = multadd(b, 10, 0);
2457 /* Round off last digit */
2461 if (j > 0 || (j == 0 && dig & 1)) {
2478 if (mlo && mlo != mhi)
2484 if (s == s0) { /* don't return empty string */