3 * The Regents of the University of California. All rights reserved.
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 * 3. All advertising materials mentioning features or use of this software
14 * must display the following acknowledgement:
15 * This product includes software developed by the University of
16 * California, Berkeley and its contributors.
17 * 4. Neither the name of the University nor the names of its contributors
18 * may be used to endorse or promote products derived from this software
19 * without specific prior written permission.
21 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
22 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
25 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
30 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
33 * $FreeBSD: src/lib/libc/stdlib/strtod.c,v 1.3.8.4 2002/08/31 22:26:35 dwmalone Exp $
34 * $DragonFly: src/lib/libc/stdlib/strtod.c,v 1.3 2003/09/06 08:10:46 asmodai Exp $
36 * @(#)strtod.c 8.1 (Berkeley) 6/4/93
39 /****************************************************************
41 * The author of this software is David M. Gay.
43 * Copyright (c) 1991 by AT&T.
45 * Permission to use, copy, modify, and distribute this software for any
46 * purpose without fee is hereby granted, provided that this entire notice
47 * is included in all copies of any software which is or includes a copy
48 * or modification of this software and in all copies of the supporting
49 * documentation for such software.
51 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
52 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
53 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
54 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
56 ***************************************************************/
58 /* Please send bug reports to
60 AT&T Bell Laboratories, Room 2C-463
62 Murray Hill, NJ 07974-2070
64 dmg@research.att.com or research!dmg
67 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
69 * This strtod returns a nearest machine number to the input decimal
70 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
71 * broken by the IEEE round-even rule. Otherwise ties are broken by
72 * biased rounding (add half and chop).
74 * Inspired loosely by William D. Clinger's paper "How to Read Floating
75 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
79 * 1. We only require IEEE, IBM, or VAX double-precision
80 * arithmetic (not IEEE double-extended).
81 * 2. We get by with floating-point arithmetic in a case that
82 * Clinger missed -- when we're computing d * 10^n
83 * for a small integer d and the integer n is not too
84 * much larger than 22 (the maximum integer k for which
85 * we can represent 10^k exactly), we may be able to
86 * compute (d*10^k) * 10^(e-k) with just one roundoff.
87 * 3. Rather than a bit-at-a-time adjustment of the binary
88 * result in the hard case, we use floating-point
89 * arithmetic to determine the adjustment to within
90 * one bit; only in really hard cases do we need to
91 * compute a second residual.
92 * 4. Because of 3., we don't need a large table of powers of 10
93 * for ten-to-e (just some small tables, e.g. of 10^k
98 * #define IEEE_8087 for IEEE-arithmetic machines where the least
99 * significant byte has the lowest address.
100 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
101 * significant byte has the lowest address.
102 * #define Sudden_Underflow for IEEE-format machines without gradual
103 * underflow (i.e., that flush to zero on underflow).
104 * #define IBM for IBM mainframe-style floating-point arithmetic.
105 * #define VAX for VAX-style floating-point arithmetic.
106 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
107 * #define No_leftright to omit left-right logic in fast floating-point
108 * computation of dtoa.
109 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
110 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
111 * that use extended-precision instructions to compute rounded
112 * products and quotients) with IBM.
113 * #define ROUND_BIASED for IEEE-format with biased rounding.
114 * #define Inaccurate_Divide for IEEE-format with correctly rounded
115 * products but inaccurate quotients, e.g., for Intel i860.
116 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
117 * integer arithmetic. Whether this speeds things up or slows things
118 * down depends on the machine and the number being converted.
119 * #define KR_headers for old-style C function headers.
120 * #define Bad_float_h if your system lacks a float.h or if it does not
121 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
122 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
125 #if defined(i386) || defined(mips) && defined(MIPSEL)
133 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
155 #define IEEE_ARITHMETIC
158 #define IEEE_ARITHMETIC
160 #ifdef IEEE_ARITHMETIC
162 #define DBL_MAX_10_EXP 308
163 #define DBL_MAX_EXP 1024
166 #define DBL_MAX 1.7976931348623157e+308
171 #define DBL_MAX_10_EXP 75
172 #define DBL_MAX_EXP 63
175 #define DBL_MAX 7.2370055773322621e+75
180 #define DBL_MAX_10_EXP 38
181 #define DBL_MAX_EXP 127
184 #define DBL_MAX 1.7014118346046923e+38
188 #define LONG_MAX 2147483647
203 #define CONST /* blank */
209 #ifdef Unsigned_Shifts
210 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
212 #define Sign_Extend(a,b) /*no-op*/
215 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
216 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
219 union doubleasulongs {
224 #define word0(x) (((union doubleasulongs *)&x)->w)[1]
225 #define word1(x) (((union doubleasulongs *)&x)->w)[0]
227 #define word0(x) (((union doubleasulongs *)&x)->w)[0]
228 #define word1(x) (((union doubleasulongs *)&x)->w)[1]
231 /* The following definition of Storeinc is appropriate for MIPS processors.
232 * An alternative that might be better on some machines is
233 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
235 #if defined(IEEE_8087) + defined(VAX)
236 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
237 ((unsigned short *)a)[0] = (unsigned short)c, a++)
239 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
240 ((unsigned short *)a)[1] = (unsigned short)c, a++)
243 /* #define P DBL_MANT_DIG */
244 /* Ten_pmax = floor(P*log(2)/log(5)) */
245 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
246 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
247 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
249 #if defined(IEEE_8087) + defined(IEEE_MC68k)
251 #define Exp_shift1 20
252 #define Exp_msk1 0x100000
253 #define Exp_msk11 0x100000
254 #define Exp_mask 0x7ff00000
259 #define Exp_1 0x3ff00000
260 #define Exp_11 0x3ff00000
262 #define Frac_mask 0xfffff
263 #define Frac_mask1 0xfffff
266 #define Bndry_mask 0xfffff
267 #define Bndry_mask1 0xfffff
269 #define Sign_bit 0x80000000
275 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
277 #undef Sudden_Underflow
278 #define Sudden_Underflow
281 #define Exp_shift1 24
282 #define Exp_msk1 0x1000000
283 #define Exp_msk11 0x1000000
284 #define Exp_mask 0x7f000000
287 #define Exp_1 0x41000000
288 #define Exp_11 0x41000000
289 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
290 #define Frac_mask 0xffffff
291 #define Frac_mask1 0xffffff
294 #define Bndry_mask 0xefffff
295 #define Bndry_mask1 0xffffff
297 #define Sign_bit 0x80000000
299 #define Tiny0 0x100000
306 #define Exp_msk1 0x80
307 #define Exp_msk11 0x800000
308 #define Exp_mask 0x7f80
311 #define Exp_1 0x40800000
312 #define Exp_11 0x4080
314 #define Frac_mask 0x7fffff
315 #define Frac_mask1 0xffff007f
318 #define Bndry_mask 0xffff007f
319 #define Bndry_mask1 0xffff007f
321 #define Sign_bit 0x8000
335 #define rounded_product(a,b) a = rnd_prod(a, b)
336 #define rounded_quotient(a,b) a = rnd_quot(a, b)
338 extern double rnd_prod(), rnd_quot();
340 extern double rnd_prod(double, double), rnd_quot(double, double);
343 #define rounded_product(a,b) a *= b
344 #define rounded_quotient(a,b) a /= b
347 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
348 #define Big1 0xffffffff
351 /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
352 * This makes some inner loops simpler and sometimes saves work
353 * during multiplications, but it often seems to make things slightly
354 * slower. Hence the default is now to store 32 bits per long.
364 extern "C" double strtod(const char *s00, char **se);
365 extern "C" char *__dtoa(double d, int mode, int ndigits,
366 int *decpt, int *sign, char **rve, char **resultp);
372 int k, maxwds, sign, wds;
376 typedef struct Bigint Bigint;
390 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(long));
393 rv->sign = rv->wds = 0;
408 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
409 y->wds*sizeof(long) + 2*sizeof(int))
414 (b, m, a) Bigint *b; int m, a;
416 (Bigint *b, int m, int a) /* multiply by m and add a */
432 y = (xi & 0xffff) * m + a;
433 z = (xi >> 16) * m + (y >> 16);
435 *x++ = (z << 16) + (y & 0xffff);
443 if (wds >= b->maxwds) {
458 (s, nd0, nd, y9) CONST char *s; int nd0, nd; unsigned long y9;
460 (CONST char *s, int nd0, int nd, unsigned long y9)
468 for (k = 0, y = 1; x > y; y <<= 1, k++) ;
475 b->x[0] = y9 & 0xffff;
476 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
483 b = multadd(b, 10, *s++ - '0');
489 b = multadd(b, 10, *s++ - '0');
503 if (!(x & 0xffff0000)) {
507 if (!(x & 0xff000000)) {
511 if (!(x & 0xf0000000)) {
515 if (!(x & 0xc0000000)) {
519 if (!(x & 0x80000000)) {
521 if (!(x & 0x40000000))
530 (y) unsigned long *y;
536 unsigned long x = *y;
594 (a, b) Bigint *a, *b;
596 (Bigint *a, Bigint *b)
601 unsigned long carry, y, z;
602 unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
607 if (a->wds < b->wds) {
619 for (x = c->x, xa = x + wc; x < xa; x++)
627 for (; xb < xbe; xb++, xc0++) {
628 if ( (y = *xb & 0xffff) ) {
633 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
635 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
641 if ( (y = *xb >> 16) ) {
647 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
650 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
657 for (; xb < xbe; xc0++) {
663 z = *x++ * y + *xc + carry;
671 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
681 (b, k) Bigint *b; int k;
686 Bigint *b1, *p5, *p51;
688 static int p05[3] = { 5, 25, 125 };
691 b = multadd(b, p05[i-1], 0);
708 if (!(p51 = p5->next)) {
709 p51 = p5->next = mult(p5,p5);
720 (b, k) Bigint *b; int k;
727 unsigned long *x, *x1, *xe, z;
736 for (i = b->maxwds; n1 > i; i <<= 1)
740 for (i = 0; i < n; i++)
760 *x1++ = *x << k & 0xffff | z;
779 (a, b) Bigint *a, *b;
781 (Bigint *a, Bigint *b)
784 unsigned long *xa, *xa0, *xb, *xb0;
790 if (i > 1 && !a->x[i-1])
791 Bug("cmp called with a->x[a->wds-1] == 0");
792 if (j > 1 && !b->x[j-1])
793 Bug("cmp called with b->x[b->wds-1] == 0");
803 return *xa < *xb ? -1 : 1;
813 (a, b) Bigint *a, *b;
815 (Bigint *a, Bigint *b)
820 long borrow, y; /* We need signed shifts here. */
821 unsigned long *xa, *xae, *xb, *xbe, *xc;
852 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
854 Sign_Extend(borrow, y);
855 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
857 Sign_Extend(borrow, z);
861 y = (*xa & 0xffff) + borrow;
863 Sign_Extend(borrow, y);
864 z = (*xa++ >> 16) + borrow;
866 Sign_Extend(borrow, z);
871 y = *xa++ - *xb++ + borrow;
873 Sign_Extend(borrow, y);
879 Sign_Extend(borrow, y);
900 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
901 #ifndef Sudden_Underflow
909 #ifndef Sudden_Underflow
913 word0(a) = 0x80000 >> L;
918 word1(a) = L >= 31 ? 1 : 1 << (31 - L);
928 (a, e) Bigint *a; int *e;
933 unsigned long *xa, *xa0, w, y, z;
937 unsigned long d0, d1;
947 if (!y) Bug("zero y in b2d");
953 d0 = Exp_1 | (y >> (Ebits - k));
954 w = xa > xa0 ? *--xa : 0;
955 d1 = (y << ((32-Ebits) + k)) | (w >> (Ebits - k));
958 z = xa > xa0 ? *--xa : 0;
960 d0 = Exp_1 | (y << k) | (z >> (32 - k));
961 y = xa > xa0 ? *--xa : 0;
962 d1 = (z << k) | (y >> (32 - k));
968 if (k < Ebits + 16) {
969 z = xa > xa0 ? *--xa : 0;
970 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
971 w = xa > xa0 ? *--xa : 0;
972 y = xa > xa0 ? *--xa : 0;
973 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
976 z = xa > xa0 ? *--xa : 0;
977 w = xa > xa0 ? *--xa : 0;
979 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
980 y = xa > xa0 ? *--xa : 0;
981 d1 = w << k + 16 | y << k;
985 word0(d) = d0 >> 16 | d0 << 16;
986 word1(d) = d1 >> 16 | d1 << 16;
997 (d, e, bits) double d; int *e, *bits;
999 (double d, int *e, int *bits)
1004 unsigned long *x, y, z;
1006 unsigned long d0, d1;
1007 d0 = word0(d) >> 16 | word0(d) << 16;
1008 d1 = word1(d) >> 16 | word1(d) << 16;
1022 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1023 #ifdef Sudden_Underflow
1024 de = (int)(d0 >> Exp_shift);
1029 if ( (de = (int)(d0 >> Exp_shift)) )
1034 if ( (k = lo0bits(&y)) ) {
1035 x[0] = y | (z << (32 - k));
1040 i = b->wds = (x[1] = z) ? 2 : 1;
1044 Bug("Zero passed to d2b");
1053 if (k = lo0bits(&y))
1055 x[0] = y | z << 32 - k & 0xffff;
1056 x[1] = z >> k - 16 & 0xffff;
1061 x[1] = y >> 16 | z << 16 - k & 0xffff;
1062 x[2] = z >> k & 0xffff;
1076 Bug("Zero passed to d2b");
1093 #ifndef Sudden_Underflow
1097 *e = (de - Bias - (P-1) << 2) + k;
1098 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1100 *e = de - Bias - (P-1) + k;
1103 #ifndef Sudden_Underflow
1105 *e = de - Bias - (P-1) + 1 + k;
1107 *bits = 32*i - hi0bits(x[i-1]);
1109 *bits = (i+2)*16 - hi0bits(x[i]);
1121 (a, b) Bigint *a, *b;
1123 (Bigint *a, Bigint *b)
1132 k = ka - kb + 32*(a->wds - b->wds);
1134 k = ka - kb + 16*(a->wds - b->wds);
1138 word0(da) += (k >> 2)*Exp_msk1;
1143 word0(db) += (k >> 2)*Exp_msk1;
1149 word0(da) += k*Exp_msk1;
1152 word0(db) += k*Exp_msk1;
1160 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1161 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1170 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1171 static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1175 bigtens[] = { 1e16, 1e32, 1e64 };
1176 static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1179 bigtens[] = { 1e16, 1e32 };
1180 static double tinytens[] = { 1e-16, 1e-32 };
1188 (s00, se) CONST char *s00; char **se;
1190 (CONST char *s00, char **se)
1193 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1194 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1195 CONST char *s, *s0, *s1;
1196 double aadj, aadj1, adj, rv, rv0;
1199 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1200 char decimal_point = localeconv()->decimal_point[0];
1202 sign = nz0 = nz = 0;
1204 for (s = s00;;s++) switch(*s) {
1216 if (isspace((unsigned char)*s))
1223 while (*++s == '0') ;
1229 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1235 if ((char)c == decimal_point) {
1238 for (; c == '0'; c = *++s)
1240 if (c > '0' && c <= '9') {
1248 for (; c >= '0' && c <= '9'; c = *++s) {
1253 for (i = 1; i < nz; i++)
1256 else if (nd <= DBL_DIG + 1)
1260 else if (nd <= DBL_DIG + 1)
1268 if (c == 'e' || c == 'E') {
1269 if (!nd && !nz && !nz0) {
1281 if (c >= '0' && c <= '9') {
1284 if (c > '0' && c <= '9') {
1287 while ((c = *++s) >= '0' && c <= '9')
1289 if (s - s1 > 8 || L > 19999)
1290 /* Avoid confusion from exponents
1291 * so large that e might overflow.
1293 e = 19999; /* safe for 16 bit ints */
1310 /* Now we have nd0 digits, starting at s0, followed by a
1311 * decimal point, followed by nd-nd0 digits. The number we're
1312 * after is the integer represented by those digits times
1317 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1320 rv = tens[k - 9] * rv + z;
1322 #ifndef RND_PRODQUOT
1329 if (e <= Ten_pmax) {
1331 goto vax_ovfl_check;
1333 /* rv = */ rounded_product(rv, tens[e]);
1338 if (e <= Ten_pmax + i) {
1339 /* A fancier test would sometimes let us do
1340 * this for larger i values.
1345 /* VAX exponent range is so narrow we must
1346 * worry about overflow here...
1349 word0(rv) -= P*Exp_msk1;
1350 /* rv = */ rounded_product(rv, tens[e]);
1351 if ((word0(rv) & Exp_mask)
1352 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1354 word0(rv) += P*Exp_msk1;
1356 /* rv = */ rounded_product(rv, tens[e]);
1361 #ifndef Inaccurate_Divide
1362 else if (e >= -Ten_pmax) {
1363 /* rv = */ rounded_quotient(rv, tens[-e]);
1370 /* Get starting approximation = rv * 10**e1 */
1373 if ( (i = e1 & 15) )
1375 if ( (e1 &= ~15) ) {
1376 if (e1 > DBL_MAX_10_EXP) {
1382 /* Can't trust HUGE_VAL */
1384 word0(rv) = Exp_mask;
1394 for (j = 0; e1 > 1; j++, e1 >>= 1)
1397 /* The last multiplication could overflow. */
1398 word0(rv) -= P*Exp_msk1;
1400 if ((z = word0(rv) & Exp_mask)
1401 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1403 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1404 /* set to largest number */
1405 /* (Can't trust DBL_MAX) */
1410 word0(rv) += P*Exp_msk1;
1413 } else if (e1 < 0) {
1415 if ( (i = e1 & 15) )
1417 if ( (e1 &= ~15) ) {
1419 for (j = 0; e1 > 1; j++, e1 >>= 1)
1422 /* The last multiplication could underflow. */
1436 /* The refinement below will clean
1437 * this approximation up.
1443 /* Now the hard part -- adjusting rv to the correct value.*/
1445 /* Put digits into bd: true value = bd * 10^e */
1447 bd0 = s2b(s0, nd0, nd, y);
1450 bd = Balloc(bd0->k);
1452 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1467 #ifdef Sudden_Underflow
1469 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1474 i = bbe + bbbits - 1; /* logb(rv) */
1475 if (i < Emin) /* denormal */
1482 i = bb2 < bd2 ? bb2 : bd2;
1491 bs = pow5mult(bs, bb5);
1497 bb = lshift(bb, bb2);
1499 bd = pow5mult(bd, bd5);
1501 bd = lshift(bd, bd2);
1503 bs = lshift(bs, bs2);
1504 delta = diff(bb, bd);
1505 dsign = delta->sign;
1509 /* Error is less than half an ulp -- check for
1510 * special case of mantissa a power of two.
1512 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1514 delta = lshift(delta,Log2P);
1515 if (cmp(delta, bs) > 0)
1520 /* exactly half-way between */
1522 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1523 && word1(rv) == 0xffffffff) {
1524 /*boundary case -- increment exponent*/
1525 word0(rv) = (word0(rv) & Exp_mask)
1534 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1536 /* boundary case -- decrement exponent */
1537 #ifdef Sudden_Underflow
1538 L = word0(rv) & Exp_mask;
1547 L = (word0(rv) & Exp_mask) - Exp_msk1;
1549 word0(rv) = L | Bndry_mask1;
1550 word1(rv) = 0xffffffff;
1557 #ifndef ROUND_BIASED
1558 if (!(word1(rv) & LSB))
1563 #ifndef ROUND_BIASED
1566 #ifndef Sudden_Underflow
1574 if ((aadj = ratio(delta, bs)) <= 2.) {
1577 else if (word1(rv) || word0(rv) & Bndry_mask) {
1578 #ifndef Sudden_Underflow
1579 if (word1(rv) == Tiny1 && !word0(rv))
1585 /* special case -- power of FLT_RADIX to be */
1586 /* rounded down... */
1588 if (aadj < 2./FLT_RADIX)
1589 aadj = 1./FLT_RADIX;
1596 aadj1 = dsign ? aadj : -aadj;
1597 #ifdef Check_FLT_ROUNDS
1598 switch(FLT_ROUNDS) {
1599 case 2: /* towards +infinity */
1602 case 0: /* towards 0 */
1603 case 3: /* towards -infinity */
1607 if (FLT_ROUNDS == 0)
1611 y = word0(rv) & Exp_mask;
1613 /* Check for overflow */
1615 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1617 word0(rv) -= P*Exp_msk1;
1618 adj = aadj1 * ulp(rv);
1620 if ((word0(rv) & Exp_mask) >=
1621 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1622 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1628 word0(rv) += P*Exp_msk1;
1630 #ifdef Sudden_Underflow
1631 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1633 word0(rv) += P*Exp_msk1;
1634 adj = aadj1 * ulp(rv);
1637 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1639 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1642 if (word0(rv0) == Tiny0
1643 && word1(rv0) == Tiny1)
1649 word0(rv) -= P*Exp_msk1;
1651 adj = aadj1 * ulp(rv);
1655 /* Compute adj so that the IEEE rounding rules will
1656 * correctly round rv + adj in some half-way cases.
1657 * If rv * ulp(rv) is denormalized (i.e.,
1658 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1659 * trouble from bits lost to denormalization;
1660 * example: 1.2e-307 .
1662 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1663 aadj1 = (double)(int)(aadj + 0.5);
1667 adj = aadj1 * ulp(rv);
1671 z = word0(rv) & Exp_mask;
1673 /* Can we stop now? */
1676 /* The tolerances below are conservative. */
1677 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1678 if (aadj < .4999999 || aadj > .5000001)
1680 } else if (aadj < .4999999/FLT_RADIX)
1697 return sign ? -rv : rv;
1703 (b, S) Bigint *b, *S;
1705 (Bigint *b, Bigint *S)
1710 unsigned long carry, q, ys;
1711 unsigned long *bx, *bxe, *sx, *sxe;
1714 unsigned long si, zs;
1719 /*debug*/ if (b->wds > n)
1720 /*debug*/ Bug("oversize b in quorem");
1728 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1730 /*debug*/ if (q > 9)
1731 /*debug*/ Bug("oversized quotient in quorem");
1739 ys = (si & 0xffff) * q + carry;
1740 zs = (si >> 16) * q + (ys >> 16);
1742 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1744 Sign_Extend(borrow, y);
1745 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1747 Sign_Extend(borrow, z);
1750 ys = *sx++ * q + carry;
1752 y = *bx - (ys & 0xffff) + borrow;
1754 Sign_Extend(borrow, y);
1757 } while (sx <= sxe);
1760 while (--bxe > bx && !*bxe)
1765 if (cmp(b, S) >= 0) {
1774 ys = (si & 0xffff) + carry;
1775 zs = (si >> 16) + (ys >> 16);
1777 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1779 Sign_Extend(borrow, y);
1780 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1782 Sign_Extend(borrow, z);
1787 y = *bx - (ys & 0xffff) + borrow;
1789 Sign_Extend(borrow, y);
1792 } while (sx <= sxe);
1796 while (--bxe > bx && !*bxe)
1804 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1806 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1807 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1810 * 1. Rather than iterating, we use a simple numeric overestimate
1811 * to determine k = floor(log10(d)). We scale relevant
1812 * quantities using O(log2(k)) rather than O(k) multiplications.
1813 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1814 * try to generate digits strictly left to right. Instead, we
1815 * compute with fewer bits and propagate the carry if necessary
1816 * when rounding the final digit up. This is often faster.
1817 * 3. Under the assumption that input will be rounded nearest,
1818 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1819 * That is, we allow equality in stopping tests when the
1820 * round-nearest rule will give the same floating-point value
1821 * as would satisfaction of the stopping test with strict
1823 * 4. We remove common factors of powers of 2 from relevant
1825 * 5. When converting floating-point integers less than 1e16,
1826 * we use floating-point arithmetic rather than resorting
1827 * to multiple-precision integers.
1828 * 6. When asked to produce fewer than 15 digits, we first try
1829 * to get by with floating-point arithmetic; we resort to
1830 * multiple-precision integer arithmetic only if we cannot
1831 * guarantee that the floating-point calculation has given
1832 * the correctly rounded result. For k requested digits and
1833 * "uniformly" distributed input, the probability is
1834 * something like 10^(k-15) that we must resort to the long
1841 (d, mode, ndigits, decpt, sign, rve, resultp)
1842 double d; int mode, ndigits, *decpt, *sign; char **rve, **resultp;
1844 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve,
1848 /* Arguments ndigits, decpt, sign are similar to those
1849 of ecvt and fcvt; trailing zeros are suppressed from
1850 the returned string. If not null, *rve is set to point
1851 to the end of the return value. If d is +-Infinity or NaN,
1852 then *decpt is set to 9999.
1855 0 ==> shortest string that yields d when read in
1856 and rounded to nearest.
1857 1 ==> like 0, but with Steele & White stopping rule;
1858 e.g. with IEEE P754 arithmetic , mode 0 gives
1859 1e23 whereas mode 1 gives 9.999999999999999e22.
1860 2 ==> max(1,ndigits) significant digits. This gives a
1861 return value similar to that of ecvt, except
1862 that trailing zeros are suppressed.
1863 3 ==> through ndigits past the decimal point. This
1864 gives a return value similar to that from fcvt,
1865 except that trailing zeros are suppressed, and
1866 ndigits can be negative.
1867 4-9 should give the same return values as 2-3, i.e.,
1868 4 <= mode <= 9 ==> same return as mode
1869 2 + (mode & 1). These modes are mainly for
1870 debugging; often they run slower but sometimes
1871 faster than modes 2-3.
1872 4,5,8,9 ==> left-to-right digit generation.
1873 6-9 ==> don't try fast floating-point estimate
1876 Values of mode other than 0-9 are treated as mode 0.
1878 Sufficient space is allocated to the return value
1879 to hold the suppressed trailing zeros.
1882 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1883 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1884 spec_case, try_quick;
1886 #ifndef Sudden_Underflow
1890 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1894 if (word0(d) & Sign_bit) {
1895 /* set sign for everything, including 0's and NaNs */
1897 word0(d) &= ~Sign_bit; /* clear sign bit */
1902 #if defined(IEEE_Arith) + defined(VAX)
1904 if ((word0(d) & Exp_mask) == Exp_mask)
1906 if (word0(d) == 0x8000)
1909 /* Infinity or NaN */
1913 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
1926 d += 0; /* normalize */
1936 b = d2b(d, &be, &bbits);
1937 #ifdef Sudden_Underflow
1938 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1940 if ( (i = (int)((word0(d) >> Exp_shift1) & (Exp_mask>>Exp_shift1))) ) {
1943 word0(d2) &= Frac_mask1;
1944 word0(d2) |= Exp_11;
1946 if ( (j = 11 - hi0bits(word0(d2) & Frac_mask)) )
1950 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1951 * log10(x) = log(x) / log(10)
1952 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1953 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1955 * This suggests computing an approximation k to log10(d) by
1957 * k = (i - Bias)*0.301029995663981
1958 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1960 * We want k to be too large rather than too small.
1961 * The error in the first-order Taylor series approximation
1962 * is in our favor, so we just round up the constant enough
1963 * to compensate for any error in the multiplication of
1964 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1965 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1966 * adding 1e-13 to the constant term more than suffices.
1967 * Hence we adjust the constant term to 0.1760912590558.
1968 * (We could get a more accurate k by invoking log10,
1969 * but this is probably not worthwhile.)
1977 #ifndef Sudden_Underflow
1980 /* d is denormalized */
1982 i = bbits + be + (Bias + (P-1) - 1);
1983 x = i > 32 ? ((word0(d) << (64 - i)) | (word1(d) >> (i - 32)))
1984 : (word1(d) << (32 - i));
1986 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
1987 i -= (Bias + (P-1) - 1) + 1;
1991 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1993 if (ds < 0. && ds != k)
1994 k--; /* want k = floor(ds) */
1996 if (k >= 0 && k <= Ten_pmax) {
2018 if (mode < 0 || mode > 9)
2039 ilim = ilim1 = i = ndigits;
2045 i = ndigits + k + 1;
2051 *resultp = (char *) malloc(i + 1);
2054 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2056 /* Try to get by with floating-point arithmetic. */
2062 ieps = 2; /* conservative */
2067 /* prevent overflows */
2069 d /= bigtens[n_bigtens-1];
2072 for (; j; j >>= 1, i++)
2078 } else if ( (j1 = -k) ) {
2079 d *= tens[j1 & 0xf];
2080 for (j = j1 >> 4; j; j >>= 1, i++)
2086 if (k_check && d < 1. && ilim > 0) {
2095 word0(eps) -= (P-1)*Exp_msk1;
2105 #ifndef No_leftright
2107 /* Use Steele & White method of only
2108 * generating digits needed.
2110 eps = 0.5/tens[ilim-1] - eps;
2114 *s++ = '0' + (int)L;
2126 /* Generate ilim digits, then fix them up. */
2127 eps *= tens[ilim-1];
2128 for (i = 1;; i++, d *= 10.) {
2131 *s++ = '0' + (int)L;
2135 else if (d < 0.5 - eps) {
2136 while (*--s == '0');
2143 #ifndef No_leftright
2153 /* Do we have a "small" integer? */
2155 if (be >= 0 && k <= Int_max) {
2158 if (ndigits < 0 && ilim <= 0) {
2160 if (ilim < 0 || d <= 5*ds)
2167 #ifdef Check_FLT_ROUNDS
2168 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2174 *s++ = '0' + (int)L;
2177 if (d > ds || (d == ds && L & 1)) {
2201 #ifndef Sudden_Underflow
2202 denorm ? be + (Bias + (P-1) - 1 + 1) :
2205 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2218 if ((i = ilim) < 0) {
2227 if (m2 > 0 && s2 > 0) {
2228 i = m2 < s2 ? m2 : s2;
2236 mhi = pow5mult(mhi, m5);
2241 if ( (j = b5 - m5) )
2244 b = pow5mult(b, b5);
2248 S = pow5mult(S, s5);
2250 /* Check for special case that d is a normalized power of 2. */
2253 if (!word1(d) && !(word0(d) & Bndry_mask)
2254 #ifndef Sudden_Underflow
2255 && word0(d) & Exp_mask
2258 /* The special case */
2266 /* Arrange for convenient computation of quotients:
2267 * shift left if necessary so divisor has 4 leading 0 bits.
2269 * Perhaps we should just compute leading 28 bits of S once
2270 * and for all and pass them and a shift to quorem, so it
2271 * can do shifts and ors to compute the numerator for q.
2274 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) )
2277 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) )
2298 b = multadd(b, 10, 0); /* we botched the k estimate */
2300 mhi = multadd(mhi, 10, 0);
2304 if (ilim <= 0 && mode > 2) {
2305 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2306 /* no digits, fcvt style */
2318 mhi = lshift(mhi, m2);
2320 /* Compute mlo -- check for special case
2321 * that d is a normalized power of 2.
2326 mhi = Balloc(mhi->k);
2328 mhi = lshift(mhi, Log2P);
2332 dig = quorem(b,S) + '0';
2333 /* Do we yet have the shortest decimal string
2334 * that will round to d?
2337 delta = diff(S, mhi);
2338 j1 = delta->sign ? 1 : cmp(b, delta);
2340 #ifndef ROUND_BIASED
2341 if (j1 == 0 && !mode && !(word1(d) & 1)) {
2350 if (j < 0 || (j == 0 && !mode
2351 #ifndef ROUND_BIASED
2358 if ((j1 > 0 || (j1 == 0 && dig & 1))
2366 if (dig == '9') { /* possible if i == 1 */
2377 b = multadd(b, 10, 0);
2379 mlo = mhi = multadd(mhi, 10, 0);
2381 mlo = multadd(mlo, 10, 0);
2382 mhi = multadd(mhi, 10, 0);
2387 *s++ = dig = quorem(b,S) + '0';
2390 b = multadd(b, 10, 0);
2393 /* Round off last digit */
2397 if (j > 0 || (j == 0 && dig & 1)) {
2407 while (*--s == '0');
2413 if (mlo && mlo != mhi)
2419 if (s == s0) { /* don't return empty string */