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(__x86_64__)
151 #elif defined(__ia64)
164 #elif defined(__hppa)
178 #define ULong unsigned Long
182 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
204 #define IEEE_ARITHMETIC
207 #define IEEE_ARITHMETIC
209 #ifdef IEEE_ARITHMETIC
211 #define DBL_MAX_10_EXP 308
212 #define DBL_MAX_EXP 1024
215 #define DBL_MAX 1.7976931348623157e+308
220 #define DBL_MAX_10_EXP 75
221 #define DBL_MAX_EXP 63
224 #define DBL_MAX 7.2370055773322621e+75
229 #define DBL_MAX_10_EXP 38
230 #define DBL_MAX_EXP 127
233 #define DBL_MAX 1.7014118346046923e+38
237 #define LONG_MAX 2147483647
252 #define CONST /* blank */
258 #ifdef Unsigned_Shifts
259 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
261 #define Sign_Extend(a,b) /*no-op*/
264 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
265 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
269 #define word0(x) ((ULong *)&x)[1]
270 #define word1(x) ((ULong *)&x)[0]
272 #define word0(x) ((ULong *)&x)[0]
273 #define word1(x) ((ULong *)&x)[1]
276 /* The following definition of Storeinc is appropriate for MIPS processors.
277 * An alternative that might be better on some machines is
278 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
280 #if defined(IEEE_8087) + defined(VAX)
281 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
282 ((unsigned short *)a)[0] = (unsigned short)c, a++)
284 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
285 ((unsigned short *)a)[1] = (unsigned short)c, a++)
288 /* #define P DBL_MANT_DIG */
289 /* Ten_pmax = floor(P*log(2)/log(5)) */
290 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
291 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
292 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
294 #if defined(IEEE_8087) + defined(IEEE_MC68k)
296 #define Exp_shift1 20
297 #define Exp_msk1 0x100000
298 #define Exp_msk11 0x100000
299 #define Exp_mask 0x7ff00000
304 #define Exp_1 0x3ff00000
305 #define Exp_11 0x3ff00000
307 #define Frac_mask 0xfffff
308 #define Frac_mask1 0xfffff
311 #define Bndry_mask 0xfffff
312 #define Bndry_mask1 0xfffff
314 #define Sign_bit 0x80000000
320 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
322 #undef Sudden_Underflow
323 #define Sudden_Underflow
326 #define Exp_shift1 24
327 #define Exp_msk1 0x1000000
328 #define Exp_msk11 0x1000000
329 #define Exp_mask 0x7f000000
332 #define Exp_1 0x41000000
333 #define Exp_11 0x41000000
334 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
335 #define Frac_mask 0xffffff
336 #define Frac_mask1 0xffffff
339 #define Bndry_mask 0xefffff
340 #define Bndry_mask1 0xffffff
342 #define Sign_bit 0x80000000
344 #define Tiny0 0x100000
351 #define Exp_msk1 0x80
352 #define Exp_msk11 0x800000
353 #define Exp_mask 0x7f80
356 #define Exp_1 0x40800000
357 #define Exp_11 0x4080
359 #define Frac_mask 0x7fffff
360 #define Frac_mask1 0xffff007f
363 #define Bndry_mask 0xffff007f
364 #define Bndry_mask1 0xffff007f
366 #define Sign_bit 0x8000
380 #define rounded_product(a,b) a = rnd_prod(a, b)
381 #define rounded_quotient(a,b) a = rnd_quot(a, b)
383 extern double rnd_prod(), rnd_quot();
385 extern double rnd_prod(double, double), rnd_quot(double, double);
388 #define rounded_product(a,b) a *= b
389 #define rounded_quotient(a,b) a /= b
392 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
393 #define Big1 0xffffffff
396 /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
397 * This makes some inner loops simpler and sometimes saves work
398 * during multiplications, but it often seems to make things slightly
399 * slower. Hence the default is now to store 32 bits per long.
409 extern "C" double bsd_strtod(const char *s00, char **se);
410 extern "C" char *__dtoa(double d, int mode, int ndigits,
411 int *decpt, int *sign, char **rve, char **resultp);
417 int k, maxwds, sign, wds;
421 typedef struct Bigint Bigint;
435 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(Long));
438 rv->sign = rv->wds = 0;
453 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
454 y->wds*sizeof(Long) + 2*sizeof(int))
459 (b, m, a) Bigint *b; int m, a;
461 (Bigint *b, int m, int a) /* multiply by m and add a */
477 y = (xi & 0xffff) * m + a;
478 z = (xi >> 16) * m + (y >> 16);
480 *x++ = (z << 16) + (y & 0xffff);
488 if (wds >= b->maxwds) {
503 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
505 (CONST char *s, int nd0, int nd, ULong y9)
513 for (k = 0, y = 1; x > y; y <<= 1, k++) ;
520 b->x[0] = y9 & 0xffff;
521 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
528 b = multadd(b, 10, *s++ - '0');
534 b = multadd(b, 10, *s++ - '0');
541 (x) register ULong x;
548 if (!(x & 0xffff0000)) {
552 if (!(x & 0xff000000)) {
556 if (!(x & 0xf0000000)) {
560 if (!(x & 0xc0000000)) {
564 if (!(x & 0x80000000)) {
566 if (!(x & 0x40000000))
581 register ULong x = *y;
639 (a, b) Bigint *a, *b;
641 (Bigint *a, Bigint *b)
647 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
652 if (a->wds < b->wds) {
664 for (x = c->x, xa = x + wc; x < xa; x++)
672 for (; xb < xbe; xb++, xc0++) {
673 if ( (y = *xb & 0xffff) ) {
678 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
680 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
686 if ( (y = *xb >> 16) ) {
692 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
695 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
702 for (; xb < xbe; xc0++) {
708 z = *x++ * y + *xc + carry;
716 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
726 (b, k) Bigint *b; int k;
731 Bigint *b1, *p5, *p51;
733 static int p05[3] = { 5, 25, 125 };
736 b = multadd(b, p05[i-1], 0);
753 if (!(p51 = p5->next)) {
754 p51 = p5->next = mult(p5,p5);
765 (b, k) Bigint *b; int k;
772 ULong *x, *x1, *xe, z;
781 for (i = b->maxwds; n1 > i; i <<= 1)
785 for (i = 0; i < n; i++)
805 *x1++ = *x << k & 0xffff | z;
824 (a, b) Bigint *a, *b;
826 (Bigint *a, Bigint *b)
829 ULong *xa, *xa0, *xb, *xb0;
835 if (i > 1 && !a->x[i-1])
836 Bug("cmp called with a->x[a->wds-1] == 0");
837 if (j > 1 && !b->x[j-1])
838 Bug("cmp called with b->x[b->wds-1] == 0");
848 return *xa < *xb ? -1 : 1;
858 (a, b) Bigint *a, *b;
860 (Bigint *a, Bigint *b)
865 Long borrow, y; /* We need signed shifts here. */
866 ULong *xa, *xae, *xb, *xbe, *xc;
897 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
899 Sign_Extend(borrow, y);
900 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
902 Sign_Extend(borrow, z);
906 y = (*xa & 0xffff) + borrow;
908 Sign_Extend(borrow, y);
909 z = (*xa++ >> 16) + borrow;
911 Sign_Extend(borrow, z);
916 y = *xa++ - *xb++ + borrow;
918 Sign_Extend(borrow, y);
924 Sign_Extend(borrow, y);
945 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
946 #ifndef Sudden_Underflow
954 #ifndef Sudden_Underflow
958 word0(a) = 0x80000 >> L;
963 word1(a) = L >= 31 ? 1 : 1 << (31 - L);
973 (a, e) Bigint *a; int *e;
978 ULong *xa, *xa0, w, y, z;
992 if (!y) Bug("zero y in b2d");
998 d0 = Exp_1 | (y >> (Ebits - k));
999 w = xa > xa0 ? *--xa : 0;
1000 d1 = (y << ((32-Ebits) + k)) | (w >> (Ebits - k));
1003 z = xa > xa0 ? *--xa : 0;
1005 d0 = Exp_1 | (y << k) | (z >> (32 - k));
1006 y = xa > xa0 ? *--xa : 0;
1007 d1 = (z << k) | (y >> (32 - k));
1013 if (k < Ebits + 16) {
1014 z = xa > xa0 ? *--xa : 0;
1015 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1016 w = xa > xa0 ? *--xa : 0;
1017 y = xa > xa0 ? *--xa : 0;
1018 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1021 z = xa > xa0 ? *--xa : 0;
1022 w = xa > xa0 ? *--xa : 0;
1024 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1025 y = xa > xa0 ? *--xa : 0;
1026 d1 = w << k + 16 | y << k;
1030 word0(d) = d0 >> 16 | d0 << 16;
1031 word1(d) = d1 >> 16 | d1 << 16;
1042 (d, e, bits) double d; int *e, *bits;
1044 (double d, int *e, int *bits)
1052 d0 = word0(d) >> 16 | word0(d) << 16;
1053 d1 = word1(d) >> 16 | word1(d) << 16;
1067 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1068 #ifdef Sudden_Underflow
1069 de = (int)(d0 >> Exp_shift);
1074 if ( (de = (int)(d0 >> Exp_shift)) )
1079 if ( (k = lo0bits(&y)) ) {
1080 x[0] = y | (z << (32 - k));
1085 i = b->wds = (x[1] = z) ? 2 : 1;
1089 Bug("Zero passed to d2b");
1098 if (k = lo0bits(&y))
1100 x[0] = y | z << 32 - k & 0xffff;
1101 x[1] = z >> k - 16 & 0xffff;
1106 x[1] = y >> 16 | z << 16 - k & 0xffff;
1107 x[2] = z >> k & 0xffff;
1121 Bug("Zero passed to d2b");
1138 #ifndef Sudden_Underflow
1142 *e = (de - Bias - (P-1) << 2) + k;
1143 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1145 *e = de - Bias - (P-1) + k;
1148 #ifndef Sudden_Underflow
1150 *e = de - Bias - (P-1) + 1 + k;
1152 *bits = 32*i - hi0bits(x[i-1]);
1154 *bits = (i+2)*16 - hi0bits(x[i]);
1166 (a, b) Bigint *a, *b;
1168 (Bigint *a, Bigint *b)
1177 k = ka - kb + 32*(a->wds - b->wds);
1179 k = ka - kb + 16*(a->wds - b->wds);
1183 word0(da) += (k >> 2)*Exp_msk1;
1188 word0(db) += (k >> 2)*Exp_msk1;
1194 word0(da) += k*Exp_msk1;
1197 word0(db) += k*Exp_msk1;
1205 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1206 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1215 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1216 static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1220 bigtens[] = { 1e16, 1e32, 1e64 };
1221 static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1224 bigtens[] = { 1e16, 1e32 };
1225 static double tinytens[] = { 1e-16, 1e-32 };
1233 (s00, se) CONST char *s00; char **se;
1235 (CONST char *s00, char **se)
1238 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1239 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1240 CONST char *s, *s0, *s1;
1241 double aadj, aadj1, adj, rv, rv0;
1244 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1245 char decimal_point = '.';
1247 sign = nz0 = nz = 0;
1249 for (s = s00;;s++) switch(*s) {
1261 if (isspace((unsigned char)*s))
1268 while (*++s == '0') ;
1274 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1280 if ((char)c == decimal_point) {
1283 for (; c == '0'; c = *++s)
1285 if (c > '0' && c <= '9') {
1293 for (; c >= '0' && c <= '9'; c = *++s) {
1298 for (i = 1; i < nz; i++)
1301 else if (nd <= DBL_DIG + 1)
1305 else if (nd <= DBL_DIG + 1)
1313 if (c == 'e' || c == 'E') {
1314 if (!nd && !nz && !nz0) {
1326 if (c >= '0' && c <= '9') {
1329 if (c > '0' && c <= '9') {
1332 while ((c = *++s) >= '0' && c <= '9')
1334 if (s - s1 > 8 || L > 19999)
1335 /* Avoid confusion from exponents
1336 * so large that e might overflow.
1338 e = 19999; /* safe for 16 bit ints */
1355 /* Now we have nd0 digits, starting at s0, followed by a
1356 * decimal point, followed by nd-nd0 digits. The number we're
1357 * after is the integer represented by those digits times
1362 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1365 rv = tens[k - 9] * rv + z;
1367 #ifndef RND_PRODQUOT
1374 if (e <= Ten_pmax) {
1376 goto vax_ovfl_check;
1378 /* rv = */ rounded_product(rv, tens[e]);
1383 if (e <= Ten_pmax + i) {
1384 /* A fancier test would sometimes let us do
1385 * this for larger i values.
1390 /* VAX exponent range is so narrow we must
1391 * worry about overflow here...
1394 word0(rv) -= P*Exp_msk1;
1395 /* rv = */ rounded_product(rv, tens[e]);
1396 if ((word0(rv) & Exp_mask)
1397 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1399 word0(rv) += P*Exp_msk1;
1401 /* rv = */ rounded_product(rv, tens[e]);
1406 #ifndef Inaccurate_Divide
1407 else if (e >= -Ten_pmax) {
1408 /* rv = */ rounded_quotient(rv, tens[-e]);
1415 /* Get starting approximation = rv * 10**e1 */
1418 if ( (i = e1 & 15) )
1420 if ( (e1 &= ~15) ) {
1421 if (e1 > DBL_MAX_10_EXP) {
1427 /* Can't trust HUGE_VAL */
1429 word0(rv) = Exp_mask;
1439 for (j = 0; e1 > 1; j++, e1 >>= 1)
1442 /* The last multiplication could overflow. */
1443 word0(rv) -= P*Exp_msk1;
1445 if ((z = word0(rv) & Exp_mask)
1446 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1448 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1449 /* set to largest number */
1450 /* (Can't trust DBL_MAX) */
1455 word0(rv) += P*Exp_msk1;
1458 } else if (e1 < 0) {
1460 if ( (i = e1 & 15) )
1462 if ( (e1 &= ~15) ) {
1464 for (j = 0; e1 > 1; j++, e1 >>= 1)
1467 /* The last multiplication could underflow. */
1481 /* The refinement below will clean
1482 * this approximation up.
1488 /* Now the hard part -- adjusting rv to the correct value.*/
1490 /* Put digits into bd: true value = bd * 10^e */
1492 bd0 = s2b(s0, nd0, nd, y);
1495 bd = Balloc(bd0->k);
1497 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1512 #ifdef Sudden_Underflow
1514 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1519 i = bbe + bbbits - 1; /* logb(rv) */
1520 if (i < Emin) /* denormal */
1527 i = bb2 < bd2 ? bb2 : bd2;
1536 bs = pow5mult(bs, bb5);
1542 bb = lshift(bb, bb2);
1544 bd = pow5mult(bd, bd5);
1546 bd = lshift(bd, bd2);
1548 bs = lshift(bs, bs2);
1549 delta = diff(bb, bd);
1550 dsign = delta->sign;
1554 /* Error is less than half an ulp -- check for
1555 * special case of mantissa a power of two.
1557 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1559 delta = lshift(delta,Log2P);
1560 if (cmp(delta, bs) > 0)
1565 /* exactly half-way between */
1567 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1568 && word1(rv) == 0xffffffff) {
1569 /*boundary case -- increment exponent*/
1570 word0(rv) = (word0(rv) & Exp_mask)
1579 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1581 /* boundary case -- decrement exponent */
1582 #ifdef Sudden_Underflow
1583 L = word0(rv) & Exp_mask;
1592 L = (word0(rv) & Exp_mask) - Exp_msk1;
1594 word0(rv) = L | Bndry_mask1;
1595 word1(rv) = 0xffffffff;
1602 #ifndef ROUND_BIASED
1603 if (!(word1(rv) & LSB))
1608 #ifndef ROUND_BIASED
1611 #ifndef Sudden_Underflow
1619 if ((aadj = ratio(delta, bs)) <= 2.) {
1622 else if (word1(rv) || word0(rv) & Bndry_mask) {
1623 #ifndef Sudden_Underflow
1624 if (word1(rv) == Tiny1 && !word0(rv))
1630 /* special case -- power of FLT_RADIX to be */
1631 /* rounded down... */
1633 if (aadj < 2./FLT_RADIX)
1634 aadj = 1./FLT_RADIX;
1641 aadj1 = dsign ? aadj : -aadj;
1642 #ifdef Check_FLT_ROUNDS
1643 switch(FLT_ROUNDS) {
1644 case 2: /* towards +infinity */
1647 case 0: /* towards 0 */
1648 case 3: /* towards -infinity */
1652 if (FLT_ROUNDS == 0)
1656 y = word0(rv) & Exp_mask;
1658 /* Check for overflow */
1660 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1662 word0(rv) -= P*Exp_msk1;
1663 adj = aadj1 * ulp(rv);
1665 if ((word0(rv) & Exp_mask) >=
1666 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1667 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1673 word0(rv) += P*Exp_msk1;
1675 #ifdef Sudden_Underflow
1676 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1678 word0(rv) += P*Exp_msk1;
1679 adj = aadj1 * ulp(rv);
1682 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1684 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1687 if (word0(rv0) == Tiny0
1688 && word1(rv0) == Tiny1)
1694 word0(rv) -= P*Exp_msk1;
1696 adj = aadj1 * ulp(rv);
1700 /* Compute adj so that the IEEE rounding rules will
1701 * correctly round rv + adj in some half-way cases.
1702 * If rv * ulp(rv) is denormalized (i.e.,
1703 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1704 * trouble from bits lost to denormalization;
1705 * example: 1.2e-307 .
1707 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1708 aadj1 = (double)(int)(aadj + 0.5);
1712 adj = aadj1 * ulp(rv);
1716 z = word0(rv) & Exp_mask;
1718 /* Can we stop now? */
1721 /* The tolerances below are conservative. */
1722 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1723 if (aadj < .4999999 || aadj > .5000001)
1725 } else if (aadj < .4999999/FLT_RADIX)
1742 return sign ? -rv : rv;
1748 (b, S) Bigint *b, *S;
1750 (Bigint *b, Bigint *S)
1756 ULong *bx, *bxe, *sx, *sxe;
1764 /*debug*/ if (b->wds > n)
1765 /*debug*/ Bug("oversize b in quorem");
1773 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1775 /*debug*/ if (q > 9)
1776 /*debug*/ Bug("oversized quotient in quorem");
1784 ys = (si & 0xffff) * q + carry;
1785 zs = (si >> 16) * q + (ys >> 16);
1787 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1789 Sign_Extend(borrow, y);
1790 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1792 Sign_Extend(borrow, z);
1795 ys = *sx++ * q + carry;
1797 y = *bx - (ys & 0xffff) + borrow;
1799 Sign_Extend(borrow, y);
1802 } while (sx <= sxe);
1805 while (--bxe > bx && !*bxe)
1810 if (cmp(b, S) >= 0) {
1819 ys = (si & 0xffff) + carry;
1820 zs = (si >> 16) + (ys >> 16);
1822 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1824 Sign_Extend(borrow, y);
1825 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1827 Sign_Extend(borrow, z);
1832 y = *bx - (ys & 0xffff) + borrow;
1834 Sign_Extend(borrow, y);
1837 } while (sx <= sxe);
1841 while (--bxe > bx && !*bxe)
1849 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1851 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1852 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1855 * 1. Rather than iterating, we use a simple numeric overestimate
1856 * to determine k = floor(log10(d)). We scale relevant
1857 * quantities using O(log2(k)) rather than O(k) multiplications.
1858 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1859 * try to generate digits strictly left to right. Instead, we
1860 * compute with fewer bits and propagate the carry if necessary
1861 * when rounding the final digit up. This is often faster.
1862 * 3. Under the assumption that input will be rounded nearest,
1863 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1864 * That is, we allow equality in stopping tests when the
1865 * round-nearest rule will give the same floating-point value
1866 * as would satisfaction of the stopping test with strict
1868 * 4. We remove common factors of powers of 2 from relevant
1870 * 5. When converting floating-point integers less than 1e16,
1871 * we use floating-point arithmetic rather than resorting
1872 * to multiple-precision integers.
1873 * 6. When asked to produce fewer than 15 digits, we first try
1874 * to get by with floating-point arithmetic; we resort to
1875 * multiple-precision integer arithmetic only if we cannot
1876 * guarantee that the floating-point calculation has given
1877 * the correctly rounded result. For k requested digits and
1878 * "uniformly" distributed input, the probability is
1879 * something like 10^(k-15) that we must resort to the long
1886 (d, mode, ndigits, decpt, sign, rve, resultp)
1887 double d; int mode, ndigits, *decpt, *sign; char **rve, **resultp;
1889 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve,
1893 /* Arguments ndigits, decpt, sign are similar to those
1894 of ecvt and fcvt; trailing zeros are suppressed from
1895 the returned string. If not null, *rve is set to point
1896 to the end of the return value. If d is +-Infinity or NaN,
1897 then *decpt is set to 9999.
1900 0 ==> shortest string that yields d when read in
1901 and rounded to nearest.
1902 1 ==> like 0, but with Steele & White stopping rule;
1903 e.g. with IEEE P754 arithmetic , mode 0 gives
1904 1e23 whereas mode 1 gives 9.999999999999999e22.
1905 2 ==> max(1,ndigits) significant digits. This gives a
1906 return value similar to that of ecvt, except
1907 that trailing zeros are suppressed.
1908 3 ==> through ndigits past the decimal point. This
1909 gives a return value similar to that from fcvt,
1910 except that trailing zeros are suppressed, and
1911 ndigits can be negative.
1912 4-9 should give the same return values as 2-3, i.e.,
1913 4 <= mode <= 9 ==> same return as mode
1914 2 + (mode & 1). These modes are mainly for
1915 debugging; often they run slower but sometimes
1916 faster than modes 2-3.
1917 4,5,8,9 ==> left-to-right digit generation.
1918 6-9 ==> don't try fast floating-point estimate
1921 Values of mode other than 0-9 are treated as mode 0.
1923 Sufficient space is allocated to the return value
1924 to hold the suppressed trailing zeros.
1927 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1928 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1929 spec_case, try_quick;
1931 #ifndef Sudden_Underflow
1935 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1939 if (word0(d) & Sign_bit) {
1940 /* set sign for everything, including 0's and NaNs */
1942 word0(d) &= ~Sign_bit; /* clear sign bit */
1947 #if defined(IEEE_Arith) + defined(VAX)
1949 if ((word0(d) & Exp_mask) == Exp_mask)
1951 if (word0(d) == 0x8000)
1954 /* Infinity or NaN */
1959 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
1962 *resultp = s = malloc (strlen (ss) + 1);
1974 d += 0; /* normalize */
1978 *resultp = s = malloc (2);
1986 b = d2b(d, &be, &bbits);
1987 #ifdef Sudden_Underflow
1988 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1990 if ( (i = (int)((word0(d) >> Exp_shift1) & (Exp_mask>>Exp_shift1))) ) {
1993 word0(d2) &= Frac_mask1;
1994 word0(d2) |= Exp_11;
1996 if ( (j = 11 - hi0bits(word0(d2) & Frac_mask)) )
2000 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2001 * log10(x) = log(x) / log(10)
2002 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2003 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2005 * This suggests computing an approximation k to log10(d) by
2007 * k = (i - Bias)*0.301029995663981
2008 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2010 * We want k to be too large rather than too small.
2011 * The error in the first-order Taylor series approximation
2012 * is in our favor, so we just round up the constant enough
2013 * to compensate for any error in the multiplication of
2014 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2015 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2016 * adding 1e-13 to the constant term more than suffices.
2017 * Hence we adjust the constant term to 0.1760912590558.
2018 * (We could get a more accurate k by invoking log10,
2019 * but this is probably not worthwhile.)
2027 #ifndef Sudden_Underflow
2030 /* d is denormalized */
2032 i = bbits + be + (Bias + (P-1) - 1);
2033 x = i > 32 ? ((word0(d) << (64 - i)) | (word1(d) >> (i - 32)))
2034 : (word1(d) << (32 - i));
2036 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2037 i -= (Bias + (P-1) - 1) + 1;
2041 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2043 if (ds < 0. && ds != k)
2044 k--; /* want k = floor(ds) */
2046 if (k >= 0 && k <= Ten_pmax) {
2068 if (mode < 0 || mode > 9)
2089 ilim = ilim1 = i = ndigits;
2095 i = ndigits + k + 1;
2101 *resultp = (char *) malloc(i + 1);
2104 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2106 /* Try to get by with floating-point arithmetic. */
2112 ieps = 2; /* conservative */
2117 /* prevent overflows */
2119 d /= bigtens[n_bigtens-1];
2122 for (; j; j >>= 1, i++)
2128 } else if ( (j1 = -k) ) {
2129 d *= tens[j1 & 0xf];
2130 for (j = j1 >> 4; j; j >>= 1, i++)
2136 if (k_check && d < 1. && ilim > 0) {
2145 word0(eps) -= (P-1)*Exp_msk1;
2155 #ifndef No_leftright
2157 /* Use Steele & White method of only
2158 * generating digits needed.
2160 eps = 0.5/tens[ilim-1] - eps;
2164 *s++ = '0' + (int)L;
2176 /* Generate ilim digits, then fix them up. */
2177 eps *= tens[ilim-1];
2178 for (i = 1;; i++, d *= 10.) {
2181 *s++ = '0' + (int)L;
2185 else if (d < 0.5 - eps) {
2186 while (*--s == '0');
2193 #ifndef No_leftright
2203 /* Do we have a "small" integer? */
2205 if (be >= 0 && k <= Int_max) {
2208 if (ndigits < 0 && ilim <= 0) {
2210 if (ilim < 0 || d <= 5*ds)
2217 #ifdef Check_FLT_ROUNDS
2218 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2224 *s++ = '0' + (int)L;
2227 if (d > ds || (d == ds && L & 1)) {
2251 #ifndef Sudden_Underflow
2252 denorm ? be + (Bias + (P-1) - 1 + 1) :
2255 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2268 if ((i = ilim) < 0) {
2277 if (m2 > 0 && s2 > 0) {
2278 i = m2 < s2 ? m2 : s2;
2286 mhi = pow5mult(mhi, m5);
2291 if ( (j = b5 - m5) )
2294 b = pow5mult(b, b5);
2298 S = pow5mult(S, s5);
2300 /* Check for special case that d is a normalized power of 2. */
2303 if (!word1(d) && !(word0(d) & Bndry_mask)
2304 #ifndef Sudden_Underflow
2305 && word0(d) & Exp_mask
2308 /* The special case */
2316 /* Arrange for convenient computation of quotients:
2317 * shift left if necessary so divisor has 4 leading 0 bits.
2319 * Perhaps we should just compute leading 28 bits of S once
2320 * and for all and pass them and a shift to quorem, so it
2321 * can do shifts and ors to compute the numerator for q.
2324 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) )
2327 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) )
2348 b = multadd(b, 10, 0); /* we botched the k estimate */
2350 mhi = multadd(mhi, 10, 0);
2354 if (ilim <= 0 && mode > 2) {
2355 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2356 /* no digits, fcvt style */
2368 mhi = lshift(mhi, m2);
2370 /* Compute mlo -- check for special case
2371 * that d is a normalized power of 2.
2376 mhi = Balloc(mhi->k);
2378 mhi = lshift(mhi, Log2P);
2382 dig = quorem(b,S) + '0';
2383 /* Do we yet have the shortest decimal string
2384 * that will round to d?
2387 delta = diff(S, mhi);
2388 j1 = delta->sign ? 1 : cmp(b, delta);
2390 #ifndef ROUND_BIASED
2391 if (j1 == 0 && !mode && !(word1(d) & 1)) {
2400 if (j < 0 || (j == 0 && !mode
2401 #ifndef ROUND_BIASED
2408 if ((j1 > 0 || (j1 == 0 && dig & 1))
2416 if (dig == '9') { /* possible if i == 1 */
2427 b = multadd(b, 10, 0);
2429 mlo = mhi = multadd(mhi, 10, 0);
2431 mlo = multadd(mlo, 10, 0);
2432 mhi = multadd(mhi, 10, 0);
2437 *s++ = dig = quorem(b,S) + '0';
2440 b = multadd(b, 10, 0);
2443 /* Round off last digit */
2447 if (j > 0 || (j == 0 && dig & 1)) {
2457 while (*--s == '0');
2463 if (mlo && mlo != mhi)
2469 if (s == s0) { /* don't return empty string */