2 * This strtod has been modified to not use values from the locale,
3 * but to hardcode the `.' as the separator. Our class libraries will
4 * make sure that only the dot is passed.
6 * This is so we do not call `setlocale' from our runtime before doing
7 * a strtod, because this could have unwanted effects in code that is
8 * co-hosted with the Mono runtime
10 * The entry point has been renamed `bsd_strtod'.
12 * Taken from the FreeBSD distribution.
20 * The Regents of the University of California. All rights reserved.
22 * Redistribution and use in source and binary forms, with or without
23 * modification, are permitted provided that the following conditions
25 * 1. Redistributions of source code must retain the above copyright
26 * notice, this list of conditions and the following disclaimer.
27 * 2. Redistributions in binary form must reproduce the above copyright
28 * notice, this list of conditions and the following disclaimer in the
29 * documentation and/or other materials provided with the distribution.
30 * 3. All advertising materials mentioning features or use of this software
31 * must display the following acknowledgement:
32 * This product includes software developed by the University of
33 * California, Berkeley and its contributors.
34 * 4. Neither the name of the University nor the names of its contributors
35 * may be used to endorse or promote products derived from this software
36 * without specific prior written permission.
38 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
39 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
40 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
41 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
42 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
43 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
44 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
45 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
46 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
47 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
50 * $FreeBSD: src/lib/libc/stdlib/strtod.c,v 1.3.8.3 2002/04/17 12:01:21 ache Exp $
53 #if defined(LIBC_SCCS) && !defined(lint)
54 static char sccsid[] = "@(#)strtod.c 8.1 (Berkeley) 6/4/93";
55 #endif /* LIBC_SCCS and not lint */
57 /****************************************************************
59 * The author of this software is David M. Gay.
61 * Copyright (c) 1991 by AT&T.
63 * Permission to use, copy, modify, and distribute this software for any
64 * purpose without fee is hereby granted, provided that this entire notice
65 * is included in all copies of any software which is or includes a copy
66 * or modification of this software and in all copies of the supporting
67 * documentation for such software.
69 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
70 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
71 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
72 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
74 ***************************************************************/
76 /* Please send bug reports to
78 AT&T Bell Laboratories, Room 2C-463
80 Murray Hill, NJ 07974-2070
82 dmg@research.att.com or research!dmg
85 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
87 * This strtod returns a nearest machine number to the input decimal
88 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
89 * broken by the IEEE round-even rule. Otherwise ties are broken by
90 * biased rounding (add half and chop).
92 * Inspired loosely by William D. Clinger's paper "How to Read Floating
93 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
97 * 1. We only require IEEE, IBM, or VAX double-precision
98 * arithmetic (not IEEE double-extended).
99 * 2. We get by with floating-point arithmetic in a case that
100 * Clinger missed -- when we're computing d * 10^n
101 * for a small integer d and the integer n is not too
102 * much larger than 22 (the maximum integer k for which
103 * we can represent 10^k exactly), we may be able to
104 * compute (d*10^k) * 10^(e-k) with just one roundoff.
105 * 3. Rather than a bit-at-a-time adjustment of the binary
106 * result in the hard case, we use floating-point
107 * arithmetic to determine the adjustment to within
108 * one bit; only in really hard cases do we need to
109 * compute a second residual.
110 * 4. Because of 3., we don't need a large table of powers of 10
111 * for ten-to-e (just some small tables, e.g. of 10^k
116 * #define IEEE_8087 for IEEE-arithmetic machines where the least
117 * significant byte has the lowest address.
118 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
119 * significant byte has the lowest address.
120 * #define Sudden_Underflow for IEEE-format machines without gradual
121 * underflow (i.e., that flush to zero on underflow).
122 * #define IBM for IBM mainframe-style floating-point arithmetic.
123 * #define VAX for VAX-style floating-point arithmetic.
124 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
125 * #define No_leftright to omit left-right logic in fast floating-point
126 * computation of dtoa.
127 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
128 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
129 * that use extended-precision instructions to compute rounded
130 * products and quotients) with IBM.
131 * #define ROUND_BIASED for IEEE-format with biased rounding.
132 * #define Inaccurate_Divide for IEEE-format with correctly rounded
133 * products but inaccurate quotients, e.g., for Intel i860.
134 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
135 * integer arithmetic. Whether this speeds things up or slows things
136 * down depends on the machine and the number being converted.
137 * #define KR_headers for old-style C function headers.
138 * #define Bad_float_h if your system lacks a float.h or if it does not
139 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
140 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
143 #if defined(i386) || defined(mips) && defined(MIPSEL) || defined (__arm__)
147 #elif defined(__x86_64__)
151 #elif defined(__ia64)
159 #elif defined(__hppa)
168 #define ULong guint32
172 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
194 #define IEEE_ARITHMETIC
197 #define IEEE_ARITHMETIC
199 #ifdef IEEE_ARITHMETIC
201 #define DBL_MAX_10_EXP 308
202 #define DBL_MAX_EXP 1024
205 #define DBL_MAX 1.7976931348623157e+308
210 #define DBL_MAX_10_EXP 75
211 #define DBL_MAX_EXP 63
214 #define DBL_MAX 7.2370055773322621e+75
219 #define DBL_MAX_10_EXP 38
220 #define DBL_MAX_EXP 127
223 #define DBL_MAX 1.7014118346046923e+38
227 #define LONG_MAX 2147483647
242 #define CONST /* blank */
248 #ifdef Unsigned_Shifts
249 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
251 #define Sign_Extend(a,b) /*no-op*/
254 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
255 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
259 #define word0(x) ((ULong *)&x)[1]
260 #define word1(x) ((ULong *)&x)[0]
262 #define word0(x) ((ULong *)&x)[0]
263 #define word1(x) ((ULong *)&x)[1]
266 /* The following definition of Storeinc is appropriate for MIPS processors.
267 * An alternative that might be better on some machines is
268 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
270 #if defined(IEEE_8087) + defined(VAX)
271 #define Storeinc(a,b,c) do { (((unsigned short *)a)[1] = (unsigned short)b, \
272 ((unsigned short *)a)[0] = (unsigned short)c); a ++; } while (0)
274 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
275 ((unsigned short *)a)[1] = (unsigned short)c, a++)
278 /* #define P DBL_MANT_DIG */
279 /* Ten_pmax = floor(P*log(2)/log(5)) */
280 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
281 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
282 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
284 #if defined(IEEE_8087) + defined(IEEE_MC68k)
286 #define Exp_shift1 20
287 #define Exp_msk1 0x100000
288 #define Exp_msk11 0x100000
289 #define Exp_mask 0x7ff00000
294 #define Exp_1 0x3ff00000
295 #define Exp_11 0x3ff00000
297 #define Frac_mask 0xfffff
298 #define Frac_mask1 0xfffff
301 #define Bndry_mask 0xfffff
302 #define Bndry_mask1 0xfffff
304 #define Sign_bit 0x80000000
310 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
312 #undef Sudden_Underflow
313 #define Sudden_Underflow
316 #define Exp_shift1 24
317 #define Exp_msk1 0x1000000
318 #define Exp_msk11 0x1000000
319 #define Exp_mask 0x7f000000
322 #define Exp_1 0x41000000
323 #define Exp_11 0x41000000
324 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
325 #define Frac_mask 0xffffff
326 #define Frac_mask1 0xffffff
329 #define Bndry_mask 0xefffff
330 #define Bndry_mask1 0xffffff
332 #define Sign_bit 0x80000000
334 #define Tiny0 0x100000
341 #define Exp_msk1 0x80
342 #define Exp_msk11 0x800000
343 #define Exp_mask 0x7f80
346 #define Exp_1 0x40800000
347 #define Exp_11 0x4080
349 #define Frac_mask 0x7fffff
350 #define Frac_mask1 0xffff007f
353 #define Bndry_mask 0xffff007f
354 #define Bndry_mask1 0xffff007f
356 #define Sign_bit 0x8000
370 #define rounded_product(a,b) a = rnd_prod(a, b)
371 #define rounded_quotient(a,b) a = rnd_quot(a, b)
373 extern double rnd_prod(), rnd_quot();
375 extern double rnd_prod(double, double), rnd_quot(double, double);
378 #define rounded_product(a,b) a *= b
379 #define rounded_quotient(a,b) a /= b
382 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
383 #define Big1 0xffffffff
386 /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
387 * This makes some inner loops simpler and sometimes saves work
388 * during multiplications, but it often seems to make things slightly
389 * slower. Hence the default is now to store 32 bits per long.
399 extern "C" double bsd_strtod(const char *s00, char **se);
400 extern "C" char *__dtoa(double d, int mode, int ndigits,
401 int *decpt, int *sign, char **rve, char **resultp);
407 int k, maxwds, sign, wds;
411 typedef struct Bigint Bigint;
425 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(Long));
428 rv->sign = rv->wds = 0;
443 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
444 y->wds*sizeof(Long) + 2*sizeof(int))
449 (b, m, a) Bigint *b; int m, a;
451 (Bigint *b, int m, int a) /* multiply by m and add a */
467 y = (xi & 0xffff) * m + a;
468 z = (xi >> 16) * m + (y >> 16);
470 *x++ = (z << 16) + (y & 0xffff);
478 if (wds >= b->maxwds) {
493 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
495 (CONST char *s, int nd0, int nd, ULong y9)
503 for (k = 0, y = 1; x > y; y <<= 1, k++) ;
510 b->x[0] = y9 & 0xffff;
511 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
518 b = multadd(b, 10, *s++ - '0');
524 b = multadd(b, 10, *s++ - '0');
531 (x) register ULong x;
538 if (!(x & 0xffff0000)) {
542 if (!(x & 0xff000000)) {
546 if (!(x & 0xf0000000)) {
550 if (!(x & 0xc0000000)) {
554 if (!(x & 0x80000000)) {
556 if (!(x & 0x40000000))
571 register ULong x = *y;
629 (a, b) Bigint *a, *b;
631 (Bigint *a, Bigint *b)
637 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
642 if (a->wds < b->wds) {
654 for (x = c->x, xa = x + wc; x < xa; x++)
663 for (; xb < xbe; xb++, xc0++) {
664 if ( (y = *xb & 0xffff) ) {
669 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
671 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
677 if ( (y = *xb >> 16) ) {
683 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
686 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
693 for (; xb < xbe; xc0++) {
699 z = *x++ * y + *xc + carry;
707 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
717 (b, k) Bigint *b; int k;
722 Bigint *b1, *p5, *p51;
724 static int p05[3] = { 5, 25, 125 };
727 b = multadd(b, p05[i-1], 0);
744 if (!(p51 = p5->next)) {
745 p51 = p5->next = mult(p5,p5);
756 (b, k) Bigint *b; int k;
763 ULong *x, *x1, *xe, z;
772 for (i = b->maxwds; n1 > i; i <<= 1)
776 for (i = 0; i < n; i++)
796 *x1++ = *x << k & 0xffff | z;
815 (a, b) Bigint *a, *b;
817 (Bigint *a, Bigint *b)
820 ULong *xa, *xa0, *xb, *xb0;
826 if (i > 1 && !a->x[i-1])
827 Bug("cmp called with a->x[a->wds-1] == 0");
828 if (j > 1 && !b->x[j-1])
829 Bug("cmp called with b->x[b->wds-1] == 0");
839 return *xa < *xb ? -1 : 1;
849 (a, b) Bigint *a, *b;
851 (Bigint *a, Bigint *b)
856 Long borrow, y; /* We need signed shifts here. */
857 ULong *xa, *xae, *xb, *xbe, *xc;
888 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
890 Sign_Extend(borrow, y);
891 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
893 Sign_Extend(borrow, z);
897 y = (*xa & 0xffff) + borrow;
899 Sign_Extend(borrow, y);
900 z = (*xa++ >> 16) + borrow;
902 Sign_Extend(borrow, z);
907 y = *xa++ - *xb++ + borrow;
909 Sign_Extend(borrow, y);
915 Sign_Extend(borrow, y);
936 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
937 #ifndef Sudden_Underflow
945 #ifndef Sudden_Underflow
949 word0(a) = 0x80000 >> L;
954 word1(a) = L >= 31 ? 1 : 1 << (31 - L);
964 (a, e) Bigint *a; int *e;
969 ULong *xa, *xa0, w, y, z;
983 if (!y) Bug("zero y in b2d");
989 d0 = Exp_1 | (y >> (Ebits - k));
990 w = xa > xa0 ? *--xa : 0;
991 d1 = (y << ((32-Ebits) + k)) | (w >> (Ebits - k));
994 z = xa > xa0 ? *--xa : 0;
996 d0 = Exp_1 | (y << k) | (z >> (32 - k));
997 y = xa > xa0 ? *--xa : 0;
998 d1 = (z << k) | (y >> (32 - k));
1004 if (k < Ebits + 16) {
1005 z = xa > xa0 ? *--xa : 0;
1006 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1007 w = xa > xa0 ? *--xa : 0;
1008 y = xa > xa0 ? *--xa : 0;
1009 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1012 z = xa > xa0 ? *--xa : 0;
1013 w = xa > xa0 ? *--xa : 0;
1015 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1016 y = xa > xa0 ? *--xa : 0;
1017 d1 = w << k + 16 | y << k;
1021 word0(d) = d0 >> 16 | d0 << 16;
1022 word1(d) = d1 >> 16 | d1 << 16;
1033 (d, e, bits) double d; int *e, *bits;
1035 (double d, int *e, int *bits)
1043 d0 = word0(d) >> 16 | word0(d) << 16;
1044 d1 = word1(d) >> 16 | word1(d) << 16;
1058 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1059 #ifdef Sudden_Underflow
1060 de = (int)(d0 >> Exp_shift);
1065 if ( (de = (int)(d0 >> Exp_shift)) )
1070 if ( (k = lo0bits(&y)) ) {
1071 x[0] = y | (z << (32 - k));
1076 i = b->wds = (x[1] = z) ? 2 : 1;
1080 Bug("Zero passed to d2b");
1089 if (k = lo0bits(&y))
1091 x[0] = y | z << 32 - k & 0xffff;
1092 x[1] = z >> k - 16 & 0xffff;
1097 x[1] = y >> 16 | z << 16 - k & 0xffff;
1098 x[2] = z >> k & 0xffff;
1112 Bug("Zero passed to d2b");
1129 #ifndef Sudden_Underflow
1133 *e = (de - Bias - (P-1) << 2) + k;
1134 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1136 *e = de - Bias - (P-1) + k;
1139 #ifndef Sudden_Underflow
1141 *e = de - Bias - (P-1) + 1 + k;
1143 *bits = 32*i - hi0bits(x[i-1]);
1145 *bits = (i+2)*16 - hi0bits(x[i]);
1157 (a, b) Bigint *a, *b;
1159 (Bigint *a, Bigint *b)
1168 k = ka - kb + 32*(a->wds - b->wds);
1170 k = ka - kb + 16*(a->wds - b->wds);
1174 word0(da) += (k >> 2)*Exp_msk1;
1179 word0(db) += (k >> 2)*Exp_msk1;
1185 word0(da) += k*Exp_msk1;
1188 word0(db) += k*Exp_msk1;
1196 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1197 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1206 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1207 static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1211 bigtens[] = { 1e16, 1e32, 1e64 };
1212 static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1215 bigtens[] = { 1e16, 1e32 };
1216 static double tinytens[] = { 1e-16, 1e-32 };
1224 (s00, se) CONST char *s00; char **se;
1226 (CONST char *s00, char **se)
1229 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1230 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1231 CONST char *s, *s0, *s1;
1232 double aadj, aadj1, adj, rv, rv0;
1235 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1236 char decimal_point = '.';
1238 sign = nz0 = nz = 0;
1240 for (s = s00;;s++) switch(*s) {
1252 if (isspace((unsigned char)*s))
1259 while (*++s == '0') ;
1265 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1271 if ((char)c == decimal_point) {
1274 for (; c == '0'; c = *++s)
1276 if (c > '0' && c <= '9') {
1284 for (; c >= '0' && c <= '9'; c = *++s) {
1289 for (i = 1; i < nz; i++)
1292 else if (nd <= DBL_DIG + 1)
1296 else if (nd <= DBL_DIG + 1)
1304 if (c == 'e' || c == 'E') {
1305 if (!nd && !nz && !nz0) {
1317 if (c >= '0' && c <= '9') {
1320 if (c > '0' && c <= '9') {
1323 while ((c = *++s) >= '0' && c <= '9')
1325 if (s - s1 > 8 || L > 19999)
1326 /* Avoid confusion from exponents
1327 * so large that e might overflow.
1329 e = 19999; /* safe for 16 bit ints */
1346 /* Now we have nd0 digits, starting at s0, followed by a
1347 * decimal point, followed by nd-nd0 digits. The number we're
1348 * after is the integer represented by those digits times
1353 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1356 rv = tens[k - 9] * rv + z;
1358 #ifndef RND_PRODQUOT
1365 if (e <= Ten_pmax) {
1367 goto vax_ovfl_check;
1369 /* rv = */ rounded_product(rv, tens[e]);
1374 if (e <= Ten_pmax + i) {
1375 /* A fancier test would sometimes let us do
1376 * this for larger i values.
1381 /* VAX exponent range is so narrow we must
1382 * worry about overflow here...
1385 word0(rv) -= P*Exp_msk1;
1386 /* rv = */ rounded_product(rv, tens[e]);
1387 if ((word0(rv) & Exp_mask)
1388 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1390 word0(rv) += P*Exp_msk1;
1392 /* rv = */ rounded_product(rv, tens[e]);
1397 #ifndef Inaccurate_Divide
1398 else if (e >= -Ten_pmax) {
1399 /* rv = */ rounded_quotient(rv, tens[-e]);
1406 /* Get starting approximation = rv * 10**e1 */
1409 if ( (i = e1 & 15) )
1411 if ( (e1 &= ~15) ) {
1412 if (e1 > DBL_MAX_10_EXP) {
1418 /* Can't trust HUGE_VAL */
1420 word0(rv) = Exp_mask;
1430 for (j = 0; e1 > 1; j++, e1 >>= 1)
1433 /* The last multiplication could overflow. */
1434 word0(rv) -= P*Exp_msk1;
1436 if ((z = word0(rv) & Exp_mask)
1437 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1439 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1440 /* set to largest number */
1441 /* (Can't trust DBL_MAX) */
1446 word0(rv) += P*Exp_msk1;
1449 } else if (e1 < 0) {
1451 if ( (i = e1 & 15) )
1453 if ( (e1 &= ~15) ) {
1455 for (j = 0; e1 > 1; j++, e1 >>= 1)
1458 /* The last multiplication could underflow. */
1472 /* The refinement below will clean
1473 * this approximation up.
1479 /* Now the hard part -- adjusting rv to the correct value.*/
1481 /* Put digits into bd: true value = bd * 10^e */
1483 bd0 = s2b(s0, nd0, nd, y);
1486 bd = Balloc(bd0->k);
1488 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1503 #ifdef Sudden_Underflow
1505 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1510 i = bbe + bbbits - 1; /* logb(rv) */
1511 if (i < Emin) /* denormal */
1518 i = bb2 < bd2 ? bb2 : bd2;
1527 bs = pow5mult(bs, bb5);
1533 bb = lshift(bb, bb2);
1535 bd = pow5mult(bd, bd5);
1537 bd = lshift(bd, bd2);
1539 bs = lshift(bs, bs2);
1540 delta = diff(bb, bd);
1541 dsign = delta->sign;
1545 /* Error is less than half an ulp -- check for
1546 * special case of mantissa a power of two.
1548 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1550 delta = lshift(delta,Log2P);
1551 if (cmp(delta, bs) > 0)
1556 /* exactly half-way between */
1558 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1559 && word1(rv) == 0xffffffff) {
1560 /*boundary case -- increment exponent*/
1561 word0(rv) = (word0(rv) & Exp_mask)
1570 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1572 /* boundary case -- decrement exponent */
1573 #ifdef Sudden_Underflow
1574 L = word0(rv) & Exp_mask;
1583 L = (word0(rv) & Exp_mask) - Exp_msk1;
1585 word0(rv) = L | Bndry_mask1;
1586 word1(rv) = 0xffffffff;
1593 #ifndef ROUND_BIASED
1594 if (!(word1(rv) & LSB))
1599 #ifndef ROUND_BIASED
1602 #ifndef Sudden_Underflow
1610 if ((aadj = ratio(delta, bs)) <= 2.) {
1613 else if (word1(rv) || word0(rv) & Bndry_mask) {
1614 #ifndef Sudden_Underflow
1615 if (word1(rv) == Tiny1 && !word0(rv))
1621 /* special case -- power of FLT_RADIX to be */
1622 /* rounded down... */
1624 if (aadj < 2./FLT_RADIX)
1625 aadj = 1./FLT_RADIX;
1632 aadj1 = dsign ? aadj : -aadj;
1633 #ifdef Check_FLT_ROUNDS
1634 switch(FLT_ROUNDS) {
1635 case 2: /* towards +infinity */
1638 case 0: /* towards 0 */
1639 case 3: /* towards -infinity */
1643 if (FLT_ROUNDS == 0)
1647 y = word0(rv) & Exp_mask;
1649 /* Check for overflow */
1651 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1653 word0(rv) -= P*Exp_msk1;
1654 adj = aadj1 * ulp(rv);
1656 if ((word0(rv) & Exp_mask) >=
1657 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1658 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1664 word0(rv) += P*Exp_msk1;
1666 #ifdef Sudden_Underflow
1667 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1669 word0(rv) += P*Exp_msk1;
1670 adj = aadj1 * ulp(rv);
1673 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1675 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1678 if (word0(rv0) == Tiny0
1679 && word1(rv0) == Tiny1)
1685 word0(rv) -= P*Exp_msk1;
1687 adj = aadj1 * ulp(rv);
1691 /* Compute adj so that the IEEE rounding rules will
1692 * correctly round rv + adj in some half-way cases.
1693 * If rv * ulp(rv) is denormalized (i.e.,
1694 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1695 * trouble from bits lost to denormalization;
1696 * example: 1.2e-307 .
1698 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1699 aadj1 = (double)(int)(aadj + 0.5);
1703 adj = aadj1 * ulp(rv);
1707 z = word0(rv) & Exp_mask;
1709 /* Can we stop now? */
1712 /* The tolerances below are conservative. */
1713 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1714 if (aadj < .4999999 || aadj > .5000001)
1716 } else if (aadj < .4999999/FLT_RADIX)
1733 return sign ? -rv : rv;
1739 (b, S) Bigint *b, *S;
1741 (Bigint *b, Bigint *S)
1747 ULong *bx, *bxe, *sx, *sxe;
1755 /*debug*/ if (b->wds > n)
1756 /*debug*/ Bug("oversize b in quorem");
1764 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1766 /*debug*/ if (q > 9)
1767 /*debug*/ Bug("oversized quotient in quorem");
1775 ys = (si & 0xffff) * q + carry;
1776 zs = (si >> 16) * q + (ys >> 16);
1778 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1780 Sign_Extend(borrow, y);
1781 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1783 Sign_Extend(borrow, z);
1786 ys = *sx++ * q + carry;
1788 y = *bx - (ys & 0xffff) + borrow;
1790 Sign_Extend(borrow, y);
1793 } while (sx <= sxe);
1796 while (--bxe > bx && !*bxe)
1801 if (cmp(b, S) >= 0) {
1810 ys = (si & 0xffff) + carry;
1811 zs = (si >> 16) + (ys >> 16);
1813 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1815 Sign_Extend(borrow, y);
1816 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1818 Sign_Extend(borrow, z);
1823 y = *bx - (ys & 0xffff) + borrow;
1825 Sign_Extend(borrow, y);
1828 } while (sx <= sxe);
1832 while (--bxe > bx && !*bxe)
1840 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1842 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1843 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1846 * 1. Rather than iterating, we use a simple numeric overestimate
1847 * to determine k = floor(log10(d)). We scale relevant
1848 * quantities using O(log2(k)) rather than O(k) multiplications.
1849 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1850 * try to generate digits strictly left to right. Instead, we
1851 * compute with fewer bits and propagate the carry if necessary
1852 * when rounding the final digit up. This is often faster.
1853 * 3. Under the assumption that input will be rounded nearest,
1854 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1855 * That is, we allow equality in stopping tests when the
1856 * round-nearest rule will give the same floating-point value
1857 * as would satisfaction of the stopping test with strict
1859 * 4. We remove common factors of powers of 2 from relevant
1861 * 5. When converting floating-point integers less than 1e16,
1862 * we use floating-point arithmetic rather than resorting
1863 * to multiple-precision integers.
1864 * 6. When asked to produce fewer than 15 digits, we first try
1865 * to get by with floating-point arithmetic; we resort to
1866 * multiple-precision integer arithmetic only if we cannot
1867 * guarantee that the floating-point calculation has given
1868 * the correctly rounded result. For k requested digits and
1869 * "uniformly" distributed input, the probability is
1870 * something like 10^(k-15) that we must resort to the long
1877 (d, mode, ndigits, decpt, sign, rve, resultp)
1878 double d; int mode, ndigits, *decpt, *sign; char **rve, **resultp;
1880 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve,
1884 /* Arguments ndigits, decpt, sign are similar to those
1885 of ecvt and fcvt; trailing zeros are suppressed from
1886 the returned string. If not null, *rve is set to point
1887 to the end of the return value. If d is +-Infinity or NaN,
1888 then *decpt is set to 9999.
1891 0 ==> shortest string that yields d when read in
1892 and rounded to nearest.
1893 1 ==> like 0, but with Steele & White stopping rule;
1894 e.g. with IEEE P754 arithmetic , mode 0 gives
1895 1e23 whereas mode 1 gives 9.999999999999999e22.
1896 2 ==> max(1,ndigits) significant digits. This gives a
1897 return value similar to that of ecvt, except
1898 that trailing zeros are suppressed.
1899 3 ==> through ndigits past the decimal point. This
1900 gives a return value similar to that from fcvt,
1901 except that trailing zeros are suppressed, and
1902 ndigits can be negative.
1903 4-9 should give the same return values as 2-3, i.e.,
1904 4 <= mode <= 9 ==> same return as mode
1905 2 + (mode & 1). These modes are mainly for
1906 debugging; often they run slower but sometimes
1907 faster than modes 2-3.
1908 4,5,8,9 ==> left-to-right digit generation.
1909 6-9 ==> don't try fast floating-point estimate
1912 Values of mode other than 0-9 are treated as mode 0.
1914 Sufficient space is allocated to the return value
1915 to hold the suppressed trailing zeros.
1918 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1919 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1920 spec_case, try_quick;
1922 #ifndef Sudden_Underflow
1926 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1930 if (word0(d) & Sign_bit) {
1931 /* set sign for everything, including 0's and NaNs */
1933 word0(d) &= ~Sign_bit; /* clear sign bit */
1938 #if defined(IEEE_Arith) + defined(VAX)
1940 if ((word0(d) & Exp_mask) == Exp_mask)
1942 if (word0(d) == 0x8000)
1945 /* Infinity or NaN */
1950 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
1953 *resultp = s = malloc (strlen (ss) + 1);
1965 d += 0; /* normalize */
1969 *resultp = s = malloc (2);
1977 b = d2b(d, &be, &bbits);
1978 #ifdef Sudden_Underflow
1979 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1981 if ( (i = (int)((word0(d) >> Exp_shift1) & (Exp_mask>>Exp_shift1))) ) {
1984 word0(d2) &= Frac_mask1;
1985 word0(d2) |= Exp_11;
1987 if ( (j = 11 - hi0bits(word0(d2) & Frac_mask)) )
1991 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1992 * log10(x) = log(x) / log(10)
1993 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1994 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1996 * This suggests computing an approximation k to log10(d) by
1998 * k = (i - Bias)*0.301029995663981
1999 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2001 * We want k to be too large rather than too small.
2002 * The error in the first-order Taylor series approximation
2003 * is in our favor, so we just round up the constant enough
2004 * to compensate for any error in the multiplication of
2005 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2006 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2007 * adding 1e-13 to the constant term more than suffices.
2008 * Hence we adjust the constant term to 0.1760912590558.
2009 * (We could get a more accurate k by invoking log10,
2010 * but this is probably not worthwhile.)
2018 #ifndef Sudden_Underflow
2021 /* d is denormalized */
2023 i = bbits + be + (Bias + (P-1) - 1);
2024 x = i > 32 ? ((word0(d) << (64 - i)) | (word1(d) >> (i - 32)))
2025 : (word1(d) << (32 - i));
2027 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2028 i -= (Bias + (P-1) - 1) + 1;
2032 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2034 if (ds < 0. && ds != k)
2035 k--; /* want k = floor(ds) */
2037 if (k >= 0 && k <= Ten_pmax) {
2059 if (mode < 0 || mode > 9)
2080 ilim = ilim1 = i = ndigits;
2086 i = ndigits + k + 1;
2092 *resultp = (char *) malloc(i + 1);
2095 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2097 /* Try to get by with floating-point arithmetic. */
2103 ieps = 2; /* conservative */
2108 /* prevent overflows */
2110 d /= bigtens[n_bigtens-1];
2113 for (; j; j >>= 1, i++)
2119 } else if ( (j1 = -k) ) {
2120 d *= tens[j1 & 0xf];
2121 for (j = j1 >> 4; j; j >>= 1, i++)
2127 if (k_check && d < 1. && ilim > 0) {
2136 word0(eps) -= (P-1)*Exp_msk1;
2146 #ifndef No_leftright
2148 /* Use Steele & White method of only
2149 * generating digits needed.
2151 eps = 0.5/tens[ilim-1] - eps;
2155 *s++ = '0' + (int)L;
2167 /* Generate ilim digits, then fix them up. */
2168 eps *= tens[ilim-1];
2169 for (i = 1;; i++, d *= 10.) {
2172 *s++ = '0' + (int)L;
2176 else if (d < 0.5 - eps) {
2177 while (*--s == '0');
2184 #ifndef No_leftright
2194 /* Do we have a "small" integer? */
2196 if (be >= 0 && k <= Int_max) {
2199 if (ndigits < 0 && ilim <= 0) {
2201 if (ilim < 0 || d <= 5*ds)
2208 #ifdef Check_FLT_ROUNDS
2209 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2215 *s++ = '0' + (int)L;
2218 if (d > ds || (d == ds && L & 1)) {
2242 #ifndef Sudden_Underflow
2243 denorm ? be + (Bias + (P-1) - 1 + 1) :
2246 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2259 if ((i = ilim) < 0) {
2268 if (m2 > 0 && s2 > 0) {
2269 i = m2 < s2 ? m2 : s2;
2277 mhi = pow5mult(mhi, m5);
2282 if ( (j = b5 - m5) )
2285 b = pow5mult(b, b5);
2289 S = pow5mult(S, s5);
2291 /* Check for special case that d is a normalized power of 2. */
2294 if (!word1(d) && !(word0(d) & Bndry_mask)
2295 #ifndef Sudden_Underflow
2296 && word0(d) & Exp_mask
2299 /* The special case */
2307 /* Arrange for convenient computation of quotients:
2308 * shift left if necessary so divisor has 4 leading 0 bits.
2310 * Perhaps we should just compute leading 28 bits of S once
2311 * and for all and pass them and a shift to quorem, so it
2312 * can do shifts and ors to compute the numerator for q.
2315 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) )
2318 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) )
2339 b = multadd(b, 10, 0); /* we botched the k estimate */
2341 mhi = multadd(mhi, 10, 0);
2345 if (ilim <= 0 && mode > 2) {
2346 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2347 /* no digits, fcvt style */
2359 mhi = lshift(mhi, m2);
2361 /* Compute mlo -- check for special case
2362 * that d is a normalized power of 2.
2367 mhi = Balloc(mhi->k);
2369 mhi = lshift(mhi, Log2P);
2373 dig = quorem(b,S) + '0';
2374 /* Do we yet have the shortest decimal string
2375 * that will round to d?
2378 delta = diff(S, mhi);
2379 j1 = delta->sign ? 1 : cmp(b, delta);
2381 #ifndef ROUND_BIASED
2382 if (j1 == 0 && !mode && !(word1(d) & 1)) {
2391 if (j < 0 || (j == 0 && !mode
2392 #ifndef ROUND_BIASED
2399 if ((j1 > 0 || (j1 == 0 && dig & 1))
2407 if (dig == '9') { /* possible if i == 1 */
2418 b = multadd(b, 10, 0);
2420 mlo = mhi = multadd(mhi, 10, 0);
2422 mlo = multadd(mlo, 10, 0);
2423 mhi = multadd(mhi, 10, 0);
2428 *s++ = dig = quorem(b,S) + '0';
2431 b = multadd(b, 10, 0);
2434 /* Round off last digit */
2438 if (j > 0 || (j == 0 && dig & 1)) {
2448 while (*--s == '0');
2454 if (mlo && mlo != mhi)
2460 if (s == s0) { /* don't return empty string */