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.
18 * The Regents of the University of California. All rights reserved.
20 * Redistribution and use in source and binary forms, with or without
21 * modification, are permitted provided that the following conditions
23 * 1. Redistributions of source code must retain the above copyright
24 * notice, this list of conditions and the following disclaimer.
25 * 2. Redistributions in binary form must reproduce the above copyright
26 * notice, this list of conditions and the following disclaimer in the
27 * documentation and/or other materials provided with the distribution.
28 * 3. All advertising materials mentioning features or use of this software
29 * must display the following acknowledgement:
30 * This product includes software developed by the University of
31 * California, Berkeley and its contributors.
32 * 4. Neither the name of the University nor the names of its contributors
33 * may be used to endorse or promote products derived from this software
34 * without specific prior written permission.
36 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
37 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
38 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
39 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
40 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
41 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
42 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
43 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
44 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
45 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
48 * $FreeBSD: src/lib/libc/stdlib/strtod.c,v 1.3.8.3 2002/04/17 12:01:21 ache Exp $
51 #if defined(LIBC_SCCS) && !defined(lint)
52 static char sccsid[] = "@(#)strtod.c 8.1 (Berkeley) 6/4/93";
53 #endif /* LIBC_SCCS and not lint */
55 /****************************************************************
57 * The author of this software is David M. Gay.
59 * Copyright (c) 1991 by AT&T.
61 * Permission to use, copy, modify, and distribute this software for any
62 * purpose without fee is hereby granted, provided that this entire notice
63 * is included in all copies of any software which is or includes a copy
64 * or modification of this software and in all copies of the supporting
65 * documentation for such software.
67 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
68 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
69 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
70 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
72 ***************************************************************/
74 /* Please send bug reports to
76 AT&T Bell Laboratories, Room 2C-463
78 Murray Hill, NJ 07974-2070
80 dmg@research.att.com or research!dmg
83 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
85 * This strtod returns a nearest machine number to the input decimal
86 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
87 * broken by the IEEE round-even rule. Otherwise ties are broken by
88 * biased rounding (add half and chop).
90 * Inspired loosely by William D. Clinger's paper "How to Read Floating
91 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
95 * 1. We only require IEEE, IBM, or VAX double-precision
96 * arithmetic (not IEEE double-extended).
97 * 2. We get by with floating-point arithmetic in a case that
98 * Clinger missed -- when we're computing d * 10^n
99 * for a small integer d and the integer n is not too
100 * much larger than 22 (the maximum integer k for which
101 * we can represent 10^k exactly), we may be able to
102 * compute (d*10^k) * 10^(e-k) with just one roundoff.
103 * 3. Rather than a bit-at-a-time adjustment of the binary
104 * result in the hard case, we use floating-point
105 * arithmetic to determine the adjustment to within
106 * one bit; only in really hard cases do we need to
107 * compute a second residual.
108 * 4. Because of 3., we don't need a large table of powers of 10
109 * for ten-to-e (just some small tables, e.g. of 10^k
114 * #define IEEE_8087 for IEEE-arithmetic machines where the least
115 * significant byte has the lowest address.
116 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
117 * significant byte has the lowest address.
118 * #define Sudden_Underflow for IEEE-format machines without gradual
119 * underflow (i.e., that flush to zero on underflow).
120 * #define IBM for IBM mainframe-style floating-point arithmetic.
121 * #define VAX for VAX-style floating-point arithmetic.
122 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
123 * #define No_leftright to omit left-right logic in fast floating-point
124 * computation of dtoa.
125 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
126 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
127 * that use extended-precision instructions to compute rounded
128 * products and quotients) with IBM.
129 * #define ROUND_BIASED for IEEE-format with biased rounding.
130 * #define Inaccurate_Divide for IEEE-format with correctly rounded
131 * products but inaccurate quotients, e.g., for Intel i860.
132 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
133 * integer arithmetic. Whether this speeds things up or slows things
134 * down depends on the machine and the number being converted.
135 * #define KR_headers for old-style C function headers.
136 * #define Bad_float_h if your system lacks a float.h or if it does not
137 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
138 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
141 #if defined(i386) || defined(mips) && defined(MIPSEL) || defined (__arm__)
146 #elif defined(__ia64)
159 #elif defined(__hppa)
173 #define ULong unsigned Long
177 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
199 #define IEEE_ARITHMETIC
202 #define IEEE_ARITHMETIC
204 #ifdef IEEE_ARITHMETIC
206 #define DBL_MAX_10_EXP 308
207 #define DBL_MAX_EXP 1024
210 #define DBL_MAX 1.7976931348623157e+308
215 #define DBL_MAX_10_EXP 75
216 #define DBL_MAX_EXP 63
219 #define DBL_MAX 7.2370055773322621e+75
224 #define DBL_MAX_10_EXP 38
225 #define DBL_MAX_EXP 127
228 #define DBL_MAX 1.7014118346046923e+38
232 #define LONG_MAX 2147483647
247 #define CONST /* blank */
253 #ifdef Unsigned_Shifts
254 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
256 #define Sign_Extend(a,b) /*no-op*/
259 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
260 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
264 #define word0(x) ((ULong *)&x)[1]
265 #define word1(x) ((ULong *)&x)[0]
267 #define word0(x) ((ULong *)&x)[0]
268 #define word1(x) ((ULong *)&x)[1]
271 /* The following definition of Storeinc is appropriate for MIPS processors.
272 * An alternative that might be better on some machines is
273 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
275 #if defined(IEEE_8087) + defined(VAX)
276 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
277 ((unsigned short *)a)[0] = (unsigned short)c, a++)
279 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
280 ((unsigned short *)a)[1] = (unsigned short)c, a++)
283 /* #define P DBL_MANT_DIG */
284 /* Ten_pmax = floor(P*log(2)/log(5)) */
285 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
286 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
287 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
289 #if defined(IEEE_8087) + defined(IEEE_MC68k)
291 #define Exp_shift1 20
292 #define Exp_msk1 0x100000
293 #define Exp_msk11 0x100000
294 #define Exp_mask 0x7ff00000
299 #define Exp_1 0x3ff00000
300 #define Exp_11 0x3ff00000
302 #define Frac_mask 0xfffff
303 #define Frac_mask1 0xfffff
306 #define Bndry_mask 0xfffff
307 #define Bndry_mask1 0xfffff
309 #define Sign_bit 0x80000000
315 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
317 #undef Sudden_Underflow
318 #define Sudden_Underflow
321 #define Exp_shift1 24
322 #define Exp_msk1 0x1000000
323 #define Exp_msk11 0x1000000
324 #define Exp_mask 0x7f000000
327 #define Exp_1 0x41000000
328 #define Exp_11 0x41000000
329 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
330 #define Frac_mask 0xffffff
331 #define Frac_mask1 0xffffff
334 #define Bndry_mask 0xefffff
335 #define Bndry_mask1 0xffffff
337 #define Sign_bit 0x80000000
339 #define Tiny0 0x100000
346 #define Exp_msk1 0x80
347 #define Exp_msk11 0x800000
348 #define Exp_mask 0x7f80
351 #define Exp_1 0x40800000
352 #define Exp_11 0x4080
354 #define Frac_mask 0x7fffff
355 #define Frac_mask1 0xffff007f
358 #define Bndry_mask 0xffff007f
359 #define Bndry_mask1 0xffff007f
361 #define Sign_bit 0x8000
375 #define rounded_product(a,b) a = rnd_prod(a, b)
376 #define rounded_quotient(a,b) a = rnd_quot(a, b)
378 extern double rnd_prod(), rnd_quot();
380 extern double rnd_prod(double, double), rnd_quot(double, double);
383 #define rounded_product(a,b) a *= b
384 #define rounded_quotient(a,b) a /= b
387 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
388 #define Big1 0xffffffff
391 /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
392 * This makes some inner loops simpler and sometimes saves work
393 * during multiplications, but it often seems to make things slightly
394 * slower. Hence the default is now to store 32 bits per long.
404 extern "C" double bsd_strtod(const char *s00, char **se);
405 extern "C" char *__dtoa(double d, int mode, int ndigits,
406 int *decpt, int *sign, char **rve, char **resultp);
412 int k, maxwds, sign, wds;
416 typedef struct Bigint Bigint;
430 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(Long));
433 rv->sign = rv->wds = 0;
448 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
449 y->wds*sizeof(Long) + 2*sizeof(int))
454 (b, m, a) Bigint *b; int m, a;
456 (Bigint *b, int m, int a) /* multiply by m and add a */
472 y = (xi & 0xffff) * m + a;
473 z = (xi >> 16) * m + (y >> 16);
475 *x++ = (z << 16) + (y & 0xffff);
483 if (wds >= b->maxwds) {
498 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
500 (CONST char *s, int nd0, int nd, ULong y9)
508 for (k = 0, y = 1; x > y; y <<= 1, k++) ;
515 b->x[0] = y9 & 0xffff;
516 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
523 b = multadd(b, 10, *s++ - '0');
529 b = multadd(b, 10, *s++ - '0');
536 (x) register ULong x;
543 if (!(x & 0xffff0000)) {
547 if (!(x & 0xff000000)) {
551 if (!(x & 0xf0000000)) {
555 if (!(x & 0xc0000000)) {
559 if (!(x & 0x80000000)) {
561 if (!(x & 0x40000000))
576 register ULong x = *y;
634 (a, b) Bigint *a, *b;
636 (Bigint *a, Bigint *b)
642 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
647 if (a->wds < b->wds) {
659 for (x = c->x, xa = x + wc; x < xa; x++)
667 for (; xb < xbe; xb++, xc0++) {
668 if ( (y = *xb & 0xffff) ) {
673 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
675 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
681 if ( (y = *xb >> 16) ) {
687 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
690 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
697 for (; xb < xbe; xc0++) {
703 z = *x++ * y + *xc + carry;
711 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
721 (b, k) Bigint *b; int k;
726 Bigint *b1, *p5, *p51;
728 static int p05[3] = { 5, 25, 125 };
731 b = multadd(b, p05[i-1], 0);
748 if (!(p51 = p5->next)) {
749 p51 = p5->next = mult(p5,p5);
760 (b, k) Bigint *b; int k;
767 ULong *x, *x1, *xe, z;
776 for (i = b->maxwds; n1 > i; i <<= 1)
780 for (i = 0; i < n; i++)
800 *x1++ = *x << k & 0xffff | z;
819 (a, b) Bigint *a, *b;
821 (Bigint *a, Bigint *b)
824 ULong *xa, *xa0, *xb, *xb0;
830 if (i > 1 && !a->x[i-1])
831 Bug("cmp called with a->x[a->wds-1] == 0");
832 if (j > 1 && !b->x[j-1])
833 Bug("cmp called with b->x[b->wds-1] == 0");
843 return *xa < *xb ? -1 : 1;
853 (a, b) Bigint *a, *b;
855 (Bigint *a, Bigint *b)
860 Long borrow, y; /* We need signed shifts here. */
861 ULong *xa, *xae, *xb, *xbe, *xc;
892 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
894 Sign_Extend(borrow, y);
895 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
897 Sign_Extend(borrow, z);
901 y = (*xa & 0xffff) + borrow;
903 Sign_Extend(borrow, y);
904 z = (*xa++ >> 16) + borrow;
906 Sign_Extend(borrow, z);
911 y = *xa++ - *xb++ + borrow;
913 Sign_Extend(borrow, y);
919 Sign_Extend(borrow, y);
940 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
941 #ifndef Sudden_Underflow
949 #ifndef Sudden_Underflow
953 word0(a) = 0x80000 >> L;
958 word1(a) = L >= 31 ? 1 : 1 << (31 - L);
968 (a, e) Bigint *a; int *e;
973 ULong *xa, *xa0, w, y, z;
987 if (!y) Bug("zero y in b2d");
993 d0 = Exp_1 | (y >> (Ebits - k));
994 w = xa > xa0 ? *--xa : 0;
995 d1 = (y << ((32-Ebits) + k)) | (w >> (Ebits - k));
998 z = xa > xa0 ? *--xa : 0;
1000 d0 = Exp_1 | (y << k) | (z >> (32 - k));
1001 y = xa > xa0 ? *--xa : 0;
1002 d1 = (z << k) | (y >> (32 - k));
1008 if (k < Ebits + 16) {
1009 z = xa > xa0 ? *--xa : 0;
1010 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1011 w = xa > xa0 ? *--xa : 0;
1012 y = xa > xa0 ? *--xa : 0;
1013 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1016 z = xa > xa0 ? *--xa : 0;
1017 w = xa > xa0 ? *--xa : 0;
1019 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1020 y = xa > xa0 ? *--xa : 0;
1021 d1 = w << k + 16 | y << k;
1025 word0(d) = d0 >> 16 | d0 << 16;
1026 word1(d) = d1 >> 16 | d1 << 16;
1037 (d, e, bits) double d; int *e, *bits;
1039 (double d, int *e, int *bits)
1047 d0 = word0(d) >> 16 | word0(d) << 16;
1048 d1 = word1(d) >> 16 | word1(d) << 16;
1062 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1063 #ifdef Sudden_Underflow
1064 de = (int)(d0 >> Exp_shift);
1069 if ( (de = (int)(d0 >> Exp_shift)) )
1074 if ( (k = lo0bits(&y)) ) {
1075 x[0] = y | (z << (32 - k));
1080 i = b->wds = (x[1] = z) ? 2 : 1;
1084 Bug("Zero passed to d2b");
1093 if (k = lo0bits(&y))
1095 x[0] = y | z << 32 - k & 0xffff;
1096 x[1] = z >> k - 16 & 0xffff;
1101 x[1] = y >> 16 | z << 16 - k & 0xffff;
1102 x[2] = z >> k & 0xffff;
1116 Bug("Zero passed to d2b");
1133 #ifndef Sudden_Underflow
1137 *e = (de - Bias - (P-1) << 2) + k;
1138 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1140 *e = de - Bias - (P-1) + k;
1143 #ifndef Sudden_Underflow
1145 *e = de - Bias - (P-1) + 1 + k;
1147 *bits = 32*i - hi0bits(x[i-1]);
1149 *bits = (i+2)*16 - hi0bits(x[i]);
1161 (a, b) Bigint *a, *b;
1163 (Bigint *a, Bigint *b)
1172 k = ka - kb + 32*(a->wds - b->wds);
1174 k = ka - kb + 16*(a->wds - b->wds);
1178 word0(da) += (k >> 2)*Exp_msk1;
1183 word0(db) += (k >> 2)*Exp_msk1;
1189 word0(da) += k*Exp_msk1;
1192 word0(db) += k*Exp_msk1;
1200 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1201 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1210 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1211 static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1215 bigtens[] = { 1e16, 1e32, 1e64 };
1216 static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1219 bigtens[] = { 1e16, 1e32 };
1220 static double tinytens[] = { 1e-16, 1e-32 };
1228 (s00, se) CONST char *s00; char **se;
1230 (CONST char *s00, char **se)
1233 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1234 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1235 CONST char *s, *s0, *s1;
1236 double aadj, aadj1, adj, rv, rv0;
1239 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1240 char decimal_point = '.';
1242 sign = nz0 = nz = 0;
1244 for (s = s00;;s++) switch(*s) {
1256 if (isspace((unsigned char)*s))
1263 while (*++s == '0') ;
1269 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1275 if ((char)c == decimal_point) {
1278 for (; c == '0'; c = *++s)
1280 if (c > '0' && c <= '9') {
1288 for (; c >= '0' && c <= '9'; c = *++s) {
1293 for (i = 1; i < nz; i++)
1296 else if (nd <= DBL_DIG + 1)
1300 else if (nd <= DBL_DIG + 1)
1308 if (c == 'e' || c == 'E') {
1309 if (!nd && !nz && !nz0) {
1321 if (c >= '0' && c <= '9') {
1324 if (c > '0' && c <= '9') {
1327 while ((c = *++s) >= '0' && c <= '9')
1329 if (s - s1 > 8 || L > 19999)
1330 /* Avoid confusion from exponents
1331 * so large that e might overflow.
1333 e = 19999; /* safe for 16 bit ints */
1350 /* Now we have nd0 digits, starting at s0, followed by a
1351 * decimal point, followed by nd-nd0 digits. The number we're
1352 * after is the integer represented by those digits times
1357 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1360 rv = tens[k - 9] * rv + z;
1362 #ifndef RND_PRODQUOT
1369 if (e <= Ten_pmax) {
1371 goto vax_ovfl_check;
1373 /* rv = */ rounded_product(rv, tens[e]);
1378 if (e <= Ten_pmax + i) {
1379 /* A fancier test would sometimes let us do
1380 * this for larger i values.
1385 /* VAX exponent range is so narrow we must
1386 * worry about overflow here...
1389 word0(rv) -= P*Exp_msk1;
1390 /* rv = */ rounded_product(rv, tens[e]);
1391 if ((word0(rv) & Exp_mask)
1392 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1394 word0(rv) += P*Exp_msk1;
1396 /* rv = */ rounded_product(rv, tens[e]);
1401 #ifndef Inaccurate_Divide
1402 else if (e >= -Ten_pmax) {
1403 /* rv = */ rounded_quotient(rv, tens[-e]);
1410 /* Get starting approximation = rv * 10**e1 */
1413 if ( (i = e1 & 15) )
1415 if ( (e1 &= ~15) ) {
1416 if (e1 > DBL_MAX_10_EXP) {
1422 /* Can't trust HUGE_VAL */
1424 word0(rv) = Exp_mask;
1434 for (j = 0; e1 > 1; j++, e1 >>= 1)
1437 /* The last multiplication could overflow. */
1438 word0(rv) -= P*Exp_msk1;
1440 if ((z = word0(rv) & Exp_mask)
1441 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1443 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1444 /* set to largest number */
1445 /* (Can't trust DBL_MAX) */
1450 word0(rv) += P*Exp_msk1;
1453 } else if (e1 < 0) {
1455 if ( (i = e1 & 15) )
1457 if ( (e1 &= ~15) ) {
1459 for (j = 0; e1 > 1; j++, e1 >>= 1)
1462 /* 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 */
1507 #ifdef Sudden_Underflow
1509 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1514 i = bbe + bbbits - 1; /* logb(rv) */
1515 if (i < Emin) /* denormal */
1522 i = bb2 < bd2 ? bb2 : bd2;
1531 bs = pow5mult(bs, bb5);
1537 bb = lshift(bb, bb2);
1539 bd = pow5mult(bd, bd5);
1541 bd = lshift(bd, bd2);
1543 bs = lshift(bs, bs2);
1544 delta = diff(bb, bd);
1545 dsign = delta->sign;
1549 /* Error is less than half an ulp -- check for
1550 * special case of mantissa a power of two.
1552 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1554 delta = lshift(delta,Log2P);
1555 if (cmp(delta, bs) > 0)
1560 /* exactly half-way between */
1562 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1563 && word1(rv) == 0xffffffff) {
1564 /*boundary case -- increment exponent*/
1565 word0(rv) = (word0(rv) & Exp_mask)
1574 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1576 /* boundary case -- decrement exponent */
1577 #ifdef Sudden_Underflow
1578 L = word0(rv) & Exp_mask;
1587 L = (word0(rv) & Exp_mask) - Exp_msk1;
1589 word0(rv) = L | Bndry_mask1;
1590 word1(rv) = 0xffffffff;
1597 #ifndef ROUND_BIASED
1598 if (!(word1(rv) & LSB))
1603 #ifndef ROUND_BIASED
1606 #ifndef Sudden_Underflow
1614 if ((aadj = ratio(delta, bs)) <= 2.) {
1617 else if (word1(rv) || word0(rv) & Bndry_mask) {
1618 #ifndef Sudden_Underflow
1619 if (word1(rv) == Tiny1 && !word0(rv))
1625 /* special case -- power of FLT_RADIX to be */
1626 /* rounded down... */
1628 if (aadj < 2./FLT_RADIX)
1629 aadj = 1./FLT_RADIX;
1636 aadj1 = dsign ? aadj : -aadj;
1637 #ifdef Check_FLT_ROUNDS
1638 switch(FLT_ROUNDS) {
1639 case 2: /* towards +infinity */
1642 case 0: /* towards 0 */
1643 case 3: /* towards -infinity */
1647 if (FLT_ROUNDS == 0)
1651 y = word0(rv) & Exp_mask;
1653 /* Check for overflow */
1655 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1657 word0(rv) -= P*Exp_msk1;
1658 adj = aadj1 * ulp(rv);
1660 if ((word0(rv) & Exp_mask) >=
1661 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1662 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1668 word0(rv) += P*Exp_msk1;
1670 #ifdef Sudden_Underflow
1671 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1673 word0(rv) += P*Exp_msk1;
1674 adj = aadj1 * ulp(rv);
1677 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1679 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1682 if (word0(rv0) == Tiny0
1683 && word1(rv0) == Tiny1)
1689 word0(rv) -= P*Exp_msk1;
1691 adj = aadj1 * ulp(rv);
1695 /* Compute adj so that the IEEE rounding rules will
1696 * correctly round rv + adj in some half-way cases.
1697 * If rv * ulp(rv) is denormalized (i.e.,
1698 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1699 * trouble from bits lost to denormalization;
1700 * example: 1.2e-307 .
1702 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1703 aadj1 = (double)(int)(aadj + 0.5);
1707 adj = aadj1 * ulp(rv);
1711 z = word0(rv) & Exp_mask;
1713 /* Can we stop now? */
1716 /* The tolerances below are conservative. */
1717 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1718 if (aadj < .4999999 || aadj > .5000001)
1720 } else if (aadj < .4999999/FLT_RADIX)
1737 return sign ? -rv : rv;
1743 (b, S) Bigint *b, *S;
1745 (Bigint *b, Bigint *S)
1751 ULong *bx, *bxe, *sx, *sxe;
1759 /*debug*/ if (b->wds > n)
1760 /*debug*/ Bug("oversize b in quorem");
1768 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1770 /*debug*/ if (q > 9)
1771 /*debug*/ Bug("oversized quotient in quorem");
1779 ys = (si & 0xffff) * q + carry;
1780 zs = (si >> 16) * q + (ys >> 16);
1782 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1784 Sign_Extend(borrow, y);
1785 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1787 Sign_Extend(borrow, z);
1790 ys = *sx++ * q + carry;
1792 y = *bx - (ys & 0xffff) + borrow;
1794 Sign_Extend(borrow, y);
1797 } while (sx <= sxe);
1800 while (--bxe > bx && !*bxe)
1805 if (cmp(b, S) >= 0) {
1814 ys = (si & 0xffff) + carry;
1815 zs = (si >> 16) + (ys >> 16);
1817 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1819 Sign_Extend(borrow, y);
1820 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1822 Sign_Extend(borrow, z);
1827 y = *bx - (ys & 0xffff) + borrow;
1829 Sign_Extend(borrow, y);
1832 } while (sx <= sxe);
1836 while (--bxe > bx && !*bxe)
1844 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1846 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1847 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1850 * 1. Rather than iterating, we use a simple numeric overestimate
1851 * to determine k = floor(log10(d)). We scale relevant
1852 * quantities using O(log2(k)) rather than O(k) multiplications.
1853 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1854 * try to generate digits strictly left to right. Instead, we
1855 * compute with fewer bits and propagate the carry if necessary
1856 * when rounding the final digit up. This is often faster.
1857 * 3. Under the assumption that input will be rounded nearest,
1858 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1859 * That is, we allow equality in stopping tests when the
1860 * round-nearest rule will give the same floating-point value
1861 * as would satisfaction of the stopping test with strict
1863 * 4. We remove common factors of powers of 2 from relevant
1865 * 5. When converting floating-point integers less than 1e16,
1866 * we use floating-point arithmetic rather than resorting
1867 * to multiple-precision integers.
1868 * 6. When asked to produce fewer than 15 digits, we first try
1869 * to get by with floating-point arithmetic; we resort to
1870 * multiple-precision integer arithmetic only if we cannot
1871 * guarantee that the floating-point calculation has given
1872 * the correctly rounded result. For k requested digits and
1873 * "uniformly" distributed input, the probability is
1874 * something like 10^(k-15) that we must resort to the long
1881 (d, mode, ndigits, decpt, sign, rve, resultp)
1882 double d; int mode, ndigits, *decpt, *sign; char **rve, **resultp;
1884 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve,
1888 /* Arguments ndigits, decpt, sign are similar to those
1889 of ecvt and fcvt; trailing zeros are suppressed from
1890 the returned string. If not null, *rve is set to point
1891 to the end of the return value. If d is +-Infinity or NaN,
1892 then *decpt is set to 9999.
1895 0 ==> shortest string that yields d when read in
1896 and rounded to nearest.
1897 1 ==> like 0, but with Steele & White stopping rule;
1898 e.g. with IEEE P754 arithmetic , mode 0 gives
1899 1e23 whereas mode 1 gives 9.999999999999999e22.
1900 2 ==> max(1,ndigits) significant digits. This gives a
1901 return value similar to that of ecvt, except
1902 that trailing zeros are suppressed.
1903 3 ==> through ndigits past the decimal point. This
1904 gives a return value similar to that from fcvt,
1905 except that trailing zeros are suppressed, and
1906 ndigits can be negative.
1907 4-9 should give the same return values as 2-3, i.e.,
1908 4 <= mode <= 9 ==> same return as mode
1909 2 + (mode & 1). These modes are mainly for
1910 debugging; often they run slower but sometimes
1911 faster than modes 2-3.
1912 4,5,8,9 ==> left-to-right digit generation.
1913 6-9 ==> don't try fast floating-point estimate
1916 Values of mode other than 0-9 are treated as mode 0.
1918 Sufficient space is allocated to the return value
1919 to hold the suppressed trailing zeros.
1922 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1923 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1924 spec_case, try_quick;
1926 #ifndef Sudden_Underflow
1930 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1934 if (word0(d) & Sign_bit) {
1935 /* set sign for everything, including 0's and NaNs */
1937 word0(d) &= ~Sign_bit; /* clear sign bit */
1942 #if defined(IEEE_Arith) + defined(VAX)
1944 if ((word0(d) & Exp_mask) == Exp_mask)
1946 if (word0(d) == 0x8000)
1949 /* Infinity or NaN */
1954 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
1957 *resultp = s = malloc (strlen (ss) + 1);
1969 d += 0; /* normalize */
1973 *resultp = s = malloc (2);
1981 b = d2b(d, &be, &bbits);
1982 #ifdef Sudden_Underflow
1983 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1985 if ( (i = (int)((word0(d) >> Exp_shift1) & (Exp_mask>>Exp_shift1))) ) {
1988 word0(d2) &= Frac_mask1;
1989 word0(d2) |= Exp_11;
1991 if ( (j = 11 - hi0bits(word0(d2) & Frac_mask)) )
1995 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1996 * log10(x) = log(x) / log(10)
1997 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1998 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2000 * This suggests computing an approximation k to log10(d) by
2002 * k = (i - Bias)*0.301029995663981
2003 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2005 * We want k to be too large rather than too small.
2006 * The error in the first-order Taylor series approximation
2007 * is in our favor, so we just round up the constant enough
2008 * to compensate for any error in the multiplication of
2009 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2010 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2011 * adding 1e-13 to the constant term more than suffices.
2012 * Hence we adjust the constant term to 0.1760912590558.
2013 * (We could get a more accurate k by invoking log10,
2014 * but this is probably not worthwhile.)
2022 #ifndef Sudden_Underflow
2025 /* d is denormalized */
2027 i = bbits + be + (Bias + (P-1) - 1);
2028 x = i > 32 ? ((word0(d) << (64 - i)) | (word1(d) >> (i - 32)))
2029 : (word1(d) << (32 - i));
2031 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2032 i -= (Bias + (P-1) - 1) + 1;
2036 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2038 if (ds < 0. && ds != k)
2039 k--; /* want k = floor(ds) */
2041 if (k >= 0 && k <= Ten_pmax) {
2063 if (mode < 0 || mode > 9)
2084 ilim = ilim1 = i = ndigits;
2090 i = ndigits + k + 1;
2096 *resultp = (char *) malloc(i + 1);
2099 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2101 /* Try to get by with floating-point arithmetic. */
2107 ieps = 2; /* conservative */
2112 /* prevent overflows */
2114 d /= bigtens[n_bigtens-1];
2117 for (; j; j >>= 1, i++)
2123 } else if ( (j1 = -k) ) {
2124 d *= tens[j1 & 0xf];
2125 for (j = j1 >> 4; j; j >>= 1, i++)
2131 if (k_check && d < 1. && ilim > 0) {
2140 word0(eps) -= (P-1)*Exp_msk1;
2150 #ifndef No_leftright
2152 /* Use Steele & White method of only
2153 * generating digits needed.
2155 eps = 0.5/tens[ilim-1] - eps;
2159 *s++ = '0' + (int)L;
2171 /* Generate ilim digits, then fix them up. */
2172 eps *= tens[ilim-1];
2173 for (i = 1;; i++, d *= 10.) {
2176 *s++ = '0' + (int)L;
2180 else if (d < 0.5 - eps) {
2181 while (*--s == '0');
2188 #ifndef No_leftright
2198 /* Do we have a "small" integer? */
2200 if (be >= 0 && k <= Int_max) {
2203 if (ndigits < 0 && ilim <= 0) {
2205 if (ilim < 0 || d <= 5*ds)
2212 #ifdef Check_FLT_ROUNDS
2213 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2219 *s++ = '0' + (int)L;
2222 if (d > ds || (d == ds && L & 1)) {
2246 #ifndef Sudden_Underflow
2247 denorm ? be + (Bias + (P-1) - 1 + 1) :
2250 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2263 if ((i = ilim) < 0) {
2272 if (m2 > 0 && s2 > 0) {
2273 i = m2 < s2 ? m2 : s2;
2281 mhi = pow5mult(mhi, m5);
2286 if ( (j = b5 - m5) )
2289 b = pow5mult(b, b5);
2293 S = pow5mult(S, s5);
2295 /* Check for special case that d is a normalized power of 2. */
2298 if (!word1(d) && !(word0(d) & Bndry_mask)
2299 #ifndef Sudden_Underflow
2300 && word0(d) & Exp_mask
2303 /* The special case */
2311 /* Arrange for convenient computation of quotients:
2312 * shift left if necessary so divisor has 4 leading 0 bits.
2314 * Perhaps we should just compute leading 28 bits of S once
2315 * and for all and pass them and a shift to quorem, so it
2316 * can do shifts and ors to compute the numerator for q.
2319 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) )
2322 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) )
2343 b = multadd(b, 10, 0); /* we botched the k estimate */
2345 mhi = multadd(mhi, 10, 0);
2349 if (ilim <= 0 && mode > 2) {
2350 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2351 /* no digits, fcvt style */
2363 mhi = lshift(mhi, m2);
2365 /* Compute mlo -- check for special case
2366 * that d is a normalized power of 2.
2371 mhi = Balloc(mhi->k);
2373 mhi = lshift(mhi, Log2P);
2377 dig = quorem(b,S) + '0';
2378 /* Do we yet have the shortest decimal string
2379 * that will round to d?
2382 delta = diff(S, mhi);
2383 j1 = delta->sign ? 1 : cmp(b, delta);
2385 #ifndef ROUND_BIASED
2386 if (j1 == 0 && !mode && !(word1(d) & 1)) {
2395 if (j < 0 || (j == 0 && !mode
2396 #ifndef ROUND_BIASED
2403 if ((j1 > 0 || (j1 == 0 && dig & 1))
2411 if (dig == '9') { /* possible if i == 1 */
2422 b = multadd(b, 10, 0);
2424 mlo = mhi = multadd(mhi, 10, 0);
2426 mlo = multadd(mlo, 10, 0);
2427 mhi = multadd(mhi, 10, 0);
2432 *s++ = dig = quorem(b,S) + '0';
2435 b = multadd(b, 10, 0);
2438 /* Round off last digit */
2442 if (j > 0 || (j == 0 && dig & 1)) {
2452 while (*--s == '0');
2458 if (mlo && mlo != mhi)
2464 if (s == s0) { /* don't return empty string */