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)
149 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
171 #define IEEE_ARITHMETIC
174 #define IEEE_ARITHMETIC
176 #ifdef IEEE_ARITHMETIC
178 #define DBL_MAX_10_EXP 308
179 #define DBL_MAX_EXP 1024
182 #define DBL_MAX 1.7976931348623157e+308
187 #define DBL_MAX_10_EXP 75
188 #define DBL_MAX_EXP 63
191 #define DBL_MAX 7.2370055773322621e+75
196 #define DBL_MAX_10_EXP 38
197 #define DBL_MAX_EXP 127
200 #define DBL_MAX 1.7014118346046923e+38
204 #define LONG_MAX 2147483647
219 #define CONST /* blank */
225 #ifdef Unsigned_Shifts
226 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
228 #define Sign_Extend(a,b) /*no-op*/
231 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
232 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
236 #define word0(x) ((unsigned long *)&x)[1]
237 #define word1(x) ((unsigned long *)&x)[0]
239 #define word0(x) ((unsigned long *)&x)[0]
240 #define word1(x) ((unsigned long *)&x)[1]
243 /* The following definition of Storeinc is appropriate for MIPS processors.
244 * An alternative that might be better on some machines is
245 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
247 #if defined(IEEE_8087) + defined(VAX)
248 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
249 ((unsigned short *)a)[0] = (unsigned short)c, a++)
251 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
252 ((unsigned short *)a)[1] = (unsigned short)c, a++)
255 /* #define P DBL_MANT_DIG */
256 /* Ten_pmax = floor(P*log(2)/log(5)) */
257 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
258 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
259 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
261 #if defined(IEEE_8087) + defined(IEEE_MC68k)
263 #define Exp_shift1 20
264 #define Exp_msk1 0x100000
265 #define Exp_msk11 0x100000
266 #define Exp_mask 0x7ff00000
271 #define Exp_1 0x3ff00000
272 #define Exp_11 0x3ff00000
274 #define Frac_mask 0xfffff
275 #define Frac_mask1 0xfffff
278 #define Bndry_mask 0xfffff
279 #define Bndry_mask1 0xfffff
281 #define Sign_bit 0x80000000
287 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
289 #undef Sudden_Underflow
290 #define Sudden_Underflow
293 #define Exp_shift1 24
294 #define Exp_msk1 0x1000000
295 #define Exp_msk11 0x1000000
296 #define Exp_mask 0x7f000000
299 #define Exp_1 0x41000000
300 #define Exp_11 0x41000000
301 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
302 #define Frac_mask 0xffffff
303 #define Frac_mask1 0xffffff
306 #define Bndry_mask 0xefffff
307 #define Bndry_mask1 0xffffff
309 #define Sign_bit 0x80000000
311 #define Tiny0 0x100000
318 #define Exp_msk1 0x80
319 #define Exp_msk11 0x800000
320 #define Exp_mask 0x7f80
323 #define Exp_1 0x40800000
324 #define Exp_11 0x4080
326 #define Frac_mask 0x7fffff
327 #define Frac_mask1 0xffff007f
330 #define Bndry_mask 0xffff007f
331 #define Bndry_mask1 0xffff007f
333 #define Sign_bit 0x8000
347 #define rounded_product(a,b) a = rnd_prod(a, b)
348 #define rounded_quotient(a,b) a = rnd_quot(a, b)
350 extern double rnd_prod(), rnd_quot();
352 extern double rnd_prod(double, double), rnd_quot(double, double);
355 #define rounded_product(a,b) a *= b
356 #define rounded_quotient(a,b) a /= b
359 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
360 #define Big1 0xffffffff
363 /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
364 * This makes some inner loops simpler and sometimes saves work
365 * during multiplications, but it often seems to make things slightly
366 * slower. Hence the default is now to store 32 bits per long.
376 extern "C" double bsd_strtod(const char *s00, char **se);
377 extern "C" char *__dtoa(double d, int mode, int ndigits,
378 int *decpt, int *sign, char **rve, char **resultp);
384 int k, maxwds, sign, wds;
388 typedef struct Bigint Bigint;
402 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(long));
405 rv->sign = rv->wds = 0;
420 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
421 y->wds*sizeof(long) + 2*sizeof(int))
426 (b, m, a) Bigint *b; int m, a;
428 (Bigint *b, int m, int a) /* multiply by m and add a */
444 y = (xi & 0xffff) * m + a;
445 z = (xi >> 16) * m + (y >> 16);
447 *x++ = (z << 16) + (y & 0xffff);
455 if (wds >= b->maxwds) {
470 (s, nd0, nd, y9) CONST char *s; int nd0, nd; unsigned long y9;
472 (CONST char *s, int nd0, int nd, unsigned long y9)
480 for (k = 0, y = 1; x > y; y <<= 1, k++) ;
487 b->x[0] = y9 & 0xffff;
488 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
495 b = multadd(b, 10, *s++ - '0');
501 b = multadd(b, 10, *s++ - '0');
508 (x) register unsigned long x;
510 (register unsigned long x)
515 if (!(x & 0xffff0000)) {
519 if (!(x & 0xff000000)) {
523 if (!(x & 0xf0000000)) {
527 if (!(x & 0xc0000000)) {
531 if (!(x & 0x80000000)) {
533 if (!(x & 0x40000000))
542 (y) unsigned long *y;
548 register unsigned long x = *y;
606 (a, b) Bigint *a, *b;
608 (Bigint *a, Bigint *b)
613 unsigned long carry, y, z;
614 unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
619 if (a->wds < b->wds) {
631 for (x = c->x, xa = x + wc; x < xa; x++)
639 for (; xb < xbe; xb++, xc0++) {
640 if ( (y = *xb & 0xffff) ) {
645 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
647 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
653 if ( (y = *xb >> 16) ) {
659 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
662 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
669 for (; xb < xbe; xc0++) {
675 z = *x++ * y + *xc + carry;
683 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
693 (b, k) Bigint *b; int k;
698 Bigint *b1, *p5, *p51;
700 static int p05[3] = { 5, 25, 125 };
703 b = multadd(b, p05[i-1], 0);
720 if (!(p51 = p5->next)) {
721 p51 = p5->next = mult(p5,p5);
732 (b, k) Bigint *b; int k;
739 unsigned long *x, *x1, *xe, z;
748 for (i = b->maxwds; n1 > i; i <<= 1)
752 for (i = 0; i < n; i++)
772 *x1++ = *x << k & 0xffff | z;
791 (a, b) Bigint *a, *b;
793 (Bigint *a, Bigint *b)
796 unsigned long *xa, *xa0, *xb, *xb0;
802 if (i > 1 && !a->x[i-1])
803 Bug("cmp called with a->x[a->wds-1] == 0");
804 if (j > 1 && !b->x[j-1])
805 Bug("cmp called with b->x[b->wds-1] == 0");
815 return *xa < *xb ? -1 : 1;
825 (a, b) Bigint *a, *b;
827 (Bigint *a, Bigint *b)
832 long borrow, y; /* We need signed shifts here. */
833 unsigned long *xa, *xae, *xb, *xbe, *xc;
864 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
866 Sign_Extend(borrow, y);
867 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
869 Sign_Extend(borrow, z);
873 y = (*xa & 0xffff) + borrow;
875 Sign_Extend(borrow, y);
876 z = (*xa++ >> 16) + borrow;
878 Sign_Extend(borrow, z);
883 y = *xa++ - *xb++ + borrow;
885 Sign_Extend(borrow, y);
891 Sign_Extend(borrow, y);
912 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
913 #ifndef Sudden_Underflow
921 #ifndef Sudden_Underflow
925 word0(a) = 0x80000 >> L;
930 word1(a) = L >= 31 ? 1 : 1 << (31 - L);
940 (a, e) Bigint *a; int *e;
945 unsigned long *xa, *xa0, w, y, z;
949 unsigned long d0, d1;
959 if (!y) Bug("zero y in b2d");
965 d0 = Exp_1 | (y >> (Ebits - k));
966 w = xa > xa0 ? *--xa : 0;
967 d1 = (y << ((32-Ebits) + k)) | (w >> (Ebits - k));
970 z = xa > xa0 ? *--xa : 0;
972 d0 = Exp_1 | (y << k) | (z >> (32 - k));
973 y = xa > xa0 ? *--xa : 0;
974 d1 = (z << k) | (y >> (32 - k));
980 if (k < Ebits + 16) {
981 z = xa > xa0 ? *--xa : 0;
982 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
983 w = xa > xa0 ? *--xa : 0;
984 y = xa > xa0 ? *--xa : 0;
985 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
988 z = xa > xa0 ? *--xa : 0;
989 w = xa > xa0 ? *--xa : 0;
991 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
992 y = xa > xa0 ? *--xa : 0;
993 d1 = w << k + 16 | y << k;
997 word0(d) = d0 >> 16 | d0 << 16;
998 word1(d) = d1 >> 16 | d1 << 16;
1009 (d, e, bits) double d; int *e, *bits;
1011 (double d, int *e, int *bits)
1016 unsigned long *x, y, z;
1018 unsigned long d0, d1;
1019 d0 = word0(d) >> 16 | word0(d) << 16;
1020 d1 = word1(d) >> 16 | word1(d) << 16;
1034 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1035 #ifdef Sudden_Underflow
1036 de = (int)(d0 >> Exp_shift);
1041 if ( (de = (int)(d0 >> Exp_shift)) )
1046 if ( (k = lo0bits(&y)) ) {
1047 x[0] = y | (z << (32 - k));
1052 i = b->wds = (x[1] = z) ? 2 : 1;
1056 Bug("Zero passed to d2b");
1065 if (k = lo0bits(&y))
1067 x[0] = y | z << 32 - k & 0xffff;
1068 x[1] = z >> k - 16 & 0xffff;
1073 x[1] = y >> 16 | z << 16 - k & 0xffff;
1074 x[2] = z >> k & 0xffff;
1088 Bug("Zero passed to d2b");
1105 #ifndef Sudden_Underflow
1109 *e = (de - Bias - (P-1) << 2) + k;
1110 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1112 *e = de - Bias - (P-1) + k;
1115 #ifndef Sudden_Underflow
1117 *e = de - Bias - (P-1) + 1 + k;
1119 *bits = 32*i - hi0bits(x[i-1]);
1121 *bits = (i+2)*16 - hi0bits(x[i]);
1133 (a, b) Bigint *a, *b;
1135 (Bigint *a, Bigint *b)
1144 k = ka - kb + 32*(a->wds - b->wds);
1146 k = ka - kb + 16*(a->wds - b->wds);
1150 word0(da) += (k >> 2)*Exp_msk1;
1155 word0(db) += (k >> 2)*Exp_msk1;
1161 word0(da) += k*Exp_msk1;
1164 word0(db) += k*Exp_msk1;
1172 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1173 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1182 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1183 static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1187 bigtens[] = { 1e16, 1e32, 1e64 };
1188 static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1191 bigtens[] = { 1e16, 1e32 };
1192 static double tinytens[] = { 1e-16, 1e-32 };
1200 (s00, se) CONST char *s00; char **se;
1202 (CONST char *s00, char **se)
1205 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1206 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1207 CONST char *s, *s0, *s1;
1208 double aadj, aadj1, adj, rv, rv0;
1211 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1212 char decimal_point = '.';
1214 sign = nz0 = nz = 0;
1216 for (s = s00;;s++) switch(*s) {
1228 if (isspace((unsigned char)*s))
1235 while (*++s == '0') ;
1241 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1247 if ((char)c == decimal_point) {
1250 for (; c == '0'; c = *++s)
1252 if (c > '0' && c <= '9') {
1260 for (; c >= '0' && c <= '9'; c = *++s) {
1265 for (i = 1; i < nz; i++)
1268 else if (nd <= DBL_DIG + 1)
1272 else if (nd <= DBL_DIG + 1)
1280 if (c == 'e' || c == 'E') {
1281 if (!nd && !nz && !nz0) {
1293 if (c >= '0' && c <= '9') {
1296 if (c > '0' && c <= '9') {
1299 while ((c = *++s) >= '0' && c <= '9')
1301 if (s - s1 > 8 || L > 19999)
1302 /* Avoid confusion from exponents
1303 * so large that e might overflow.
1305 e = 19999; /* safe for 16 bit ints */
1322 /* Now we have nd0 digits, starting at s0, followed by a
1323 * decimal point, followed by nd-nd0 digits. The number we're
1324 * after is the integer represented by those digits times
1329 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1332 rv = tens[k - 9] * rv + z;
1334 #ifndef RND_PRODQUOT
1341 if (e <= Ten_pmax) {
1343 goto vax_ovfl_check;
1345 /* rv = */ rounded_product(rv, tens[e]);
1350 if (e <= Ten_pmax + i) {
1351 /* A fancier test would sometimes let us do
1352 * this for larger i values.
1357 /* VAX exponent range is so narrow we must
1358 * worry about overflow here...
1361 word0(rv) -= P*Exp_msk1;
1362 /* rv = */ rounded_product(rv, tens[e]);
1363 if ((word0(rv) & Exp_mask)
1364 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1366 word0(rv) += P*Exp_msk1;
1368 /* rv = */ rounded_product(rv, tens[e]);
1373 #ifndef Inaccurate_Divide
1374 else if (e >= -Ten_pmax) {
1375 /* rv = */ rounded_quotient(rv, tens[-e]);
1382 /* Get starting approximation = rv * 10**e1 */
1385 if ( (i = e1 & 15) )
1387 if ( (e1 &= ~15) ) {
1388 if (e1 > DBL_MAX_10_EXP) {
1394 /* Can't trust HUGE_VAL */
1396 word0(rv) = Exp_mask;
1406 for (j = 0; e1 > 1; j++, e1 >>= 1)
1409 /* The last multiplication could overflow. */
1410 word0(rv) -= P*Exp_msk1;
1412 if ((z = word0(rv) & Exp_mask)
1413 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1415 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1416 /* set to largest number */
1417 /* (Can't trust DBL_MAX) */
1422 word0(rv) += P*Exp_msk1;
1425 } else if (e1 < 0) {
1427 if ( (i = e1 & 15) )
1429 if ( (e1 &= ~15) ) {
1431 for (j = 0; e1 > 1; j++, e1 >>= 1)
1434 /* The last multiplication could underflow. */
1448 /* The refinement below will clean
1449 * this approximation up.
1455 /* Now the hard part -- adjusting rv to the correct value.*/
1457 /* Put digits into bd: true value = bd * 10^e */
1459 bd0 = s2b(s0, nd0, nd, y);
1462 bd = Balloc(bd0->k);
1464 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1479 #ifdef Sudden_Underflow
1481 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1486 i = bbe + bbbits - 1; /* logb(rv) */
1487 if (i < Emin) /* denormal */
1494 i = bb2 < bd2 ? bb2 : bd2;
1503 bs = pow5mult(bs, bb5);
1509 bb = lshift(bb, bb2);
1511 bd = pow5mult(bd, bd5);
1513 bd = lshift(bd, bd2);
1515 bs = lshift(bs, bs2);
1516 delta = diff(bb, bd);
1517 dsign = delta->sign;
1521 /* Error is less than half an ulp -- check for
1522 * special case of mantissa a power of two.
1524 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1526 delta = lshift(delta,Log2P);
1527 if (cmp(delta, bs) > 0)
1532 /* exactly half-way between */
1534 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1535 && word1(rv) == 0xffffffff) {
1536 /*boundary case -- increment exponent*/
1537 word0(rv) = (word0(rv) & Exp_mask)
1546 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1548 /* boundary case -- decrement exponent */
1549 #ifdef Sudden_Underflow
1550 L = word0(rv) & Exp_mask;
1559 L = (word0(rv) & Exp_mask) - Exp_msk1;
1561 word0(rv) = L | Bndry_mask1;
1562 word1(rv) = 0xffffffff;
1569 #ifndef ROUND_BIASED
1570 if (!(word1(rv) & LSB))
1575 #ifndef ROUND_BIASED
1578 #ifndef Sudden_Underflow
1586 if ((aadj = ratio(delta, bs)) <= 2.) {
1589 else if (word1(rv) || word0(rv) & Bndry_mask) {
1590 #ifndef Sudden_Underflow
1591 if (word1(rv) == Tiny1 && !word0(rv))
1597 /* special case -- power of FLT_RADIX to be */
1598 /* rounded down... */
1600 if (aadj < 2./FLT_RADIX)
1601 aadj = 1./FLT_RADIX;
1608 aadj1 = dsign ? aadj : -aadj;
1609 #ifdef Check_FLT_ROUNDS
1610 switch(FLT_ROUNDS) {
1611 case 2: /* towards +infinity */
1614 case 0: /* towards 0 */
1615 case 3: /* towards -infinity */
1619 if (FLT_ROUNDS == 0)
1623 y = word0(rv) & Exp_mask;
1625 /* Check for overflow */
1627 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1629 word0(rv) -= P*Exp_msk1;
1630 adj = aadj1 * ulp(rv);
1632 if ((word0(rv) & Exp_mask) >=
1633 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1634 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1640 word0(rv) += P*Exp_msk1;
1642 #ifdef Sudden_Underflow
1643 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1645 word0(rv) += P*Exp_msk1;
1646 adj = aadj1 * ulp(rv);
1649 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1651 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1654 if (word0(rv0) == Tiny0
1655 && word1(rv0) == Tiny1)
1661 word0(rv) -= P*Exp_msk1;
1663 adj = aadj1 * ulp(rv);
1667 /* Compute adj so that the IEEE rounding rules will
1668 * correctly round rv + adj in some half-way cases.
1669 * If rv * ulp(rv) is denormalized (i.e.,
1670 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1671 * trouble from bits lost to denormalization;
1672 * example: 1.2e-307 .
1674 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1675 aadj1 = (double)(int)(aadj + 0.5);
1679 adj = aadj1 * ulp(rv);
1683 z = word0(rv) & Exp_mask;
1685 /* Can we stop now? */
1688 /* The tolerances below are conservative. */
1689 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1690 if (aadj < .4999999 || aadj > .5000001)
1692 } else if (aadj < .4999999/FLT_RADIX)
1709 return sign ? -rv : rv;
1715 (b, S) Bigint *b, *S;
1717 (Bigint *b, Bigint *S)
1722 unsigned long carry, q, ys;
1723 unsigned long *bx, *bxe, *sx, *sxe;
1726 unsigned long si, zs;
1731 /*debug*/ if (b->wds > n)
1732 /*debug*/ Bug("oversize b in quorem");
1740 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1742 /*debug*/ if (q > 9)
1743 /*debug*/ Bug("oversized quotient in quorem");
1751 ys = (si & 0xffff) * q + carry;
1752 zs = (si >> 16) * q + (ys >> 16);
1754 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1756 Sign_Extend(borrow, y);
1757 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1759 Sign_Extend(borrow, z);
1762 ys = *sx++ * q + carry;
1764 y = *bx - (ys & 0xffff) + borrow;
1766 Sign_Extend(borrow, y);
1769 } while (sx <= sxe);
1772 while (--bxe > bx && !*bxe)
1777 if (cmp(b, S) >= 0) {
1786 ys = (si & 0xffff) + carry;
1787 zs = (si >> 16) + (ys >> 16);
1789 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1791 Sign_Extend(borrow, y);
1792 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1794 Sign_Extend(borrow, z);
1799 y = *bx - (ys & 0xffff) + borrow;
1801 Sign_Extend(borrow, y);
1804 } while (sx <= sxe);
1808 while (--bxe > bx && !*bxe)
1816 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1818 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1819 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1822 * 1. Rather than iterating, we use a simple numeric overestimate
1823 * to determine k = floor(log10(d)). We scale relevant
1824 * quantities using O(log2(k)) rather than O(k) multiplications.
1825 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1826 * try to generate digits strictly left to right. Instead, we
1827 * compute with fewer bits and propagate the carry if necessary
1828 * when rounding the final digit up. This is often faster.
1829 * 3. Under the assumption that input will be rounded nearest,
1830 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1831 * That is, we allow equality in stopping tests when the
1832 * round-nearest rule will give the same floating-point value
1833 * as would satisfaction of the stopping test with strict
1835 * 4. We remove common factors of powers of 2 from relevant
1837 * 5. When converting floating-point integers less than 1e16,
1838 * we use floating-point arithmetic rather than resorting
1839 * to multiple-precision integers.
1840 * 6. When asked to produce fewer than 15 digits, we first try
1841 * to get by with floating-point arithmetic; we resort to
1842 * multiple-precision integer arithmetic only if we cannot
1843 * guarantee that the floating-point calculation has given
1844 * the correctly rounded result. For k requested digits and
1845 * "uniformly" distributed input, the probability is
1846 * something like 10^(k-15) that we must resort to the long
1853 (d, mode, ndigits, decpt, sign, rve, resultp)
1854 double d; int mode, ndigits, *decpt, *sign; char **rve, **resultp;
1856 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve,
1860 /* Arguments ndigits, decpt, sign are similar to those
1861 of ecvt and fcvt; trailing zeros are suppressed from
1862 the returned string. If not null, *rve is set to point
1863 to the end of the return value. If d is +-Infinity or NaN,
1864 then *decpt is set to 9999.
1867 0 ==> shortest string that yields d when read in
1868 and rounded to nearest.
1869 1 ==> like 0, but with Steele & White stopping rule;
1870 e.g. with IEEE P754 arithmetic , mode 0 gives
1871 1e23 whereas mode 1 gives 9.999999999999999e22.
1872 2 ==> max(1,ndigits) significant digits. This gives a
1873 return value similar to that of ecvt, except
1874 that trailing zeros are suppressed.
1875 3 ==> through ndigits past the decimal point. This
1876 gives a return value similar to that from fcvt,
1877 except that trailing zeros are suppressed, and
1878 ndigits can be negative.
1879 4-9 should give the same return values as 2-3, i.e.,
1880 4 <= mode <= 9 ==> same return as mode
1881 2 + (mode & 1). These modes are mainly for
1882 debugging; often they run slower but sometimes
1883 faster than modes 2-3.
1884 4,5,8,9 ==> left-to-right digit generation.
1885 6-9 ==> don't try fast floating-point estimate
1888 Values of mode other than 0-9 are treated as mode 0.
1890 Sufficient space is allocated to the return value
1891 to hold the suppressed trailing zeros.
1894 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1895 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1896 spec_case, try_quick;
1898 #ifndef Sudden_Underflow
1902 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1906 if (word0(d) & Sign_bit) {
1907 /* set sign for everything, including 0's and NaNs */
1909 word0(d) &= ~Sign_bit; /* clear sign bit */
1914 #if defined(IEEE_Arith) + defined(VAX)
1916 if ((word0(d) & Exp_mask) == Exp_mask)
1918 if (word0(d) == 0x8000)
1921 /* Infinity or NaN */
1926 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
1929 *resultp = s = malloc (strlen (ss) + 1);
1941 d += 0; /* normalize */
1945 *resultp = s = malloc (2);
1953 b = d2b(d, &be, &bbits);
1954 #ifdef Sudden_Underflow
1955 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1957 if ( (i = (int)((word0(d) >> Exp_shift1) & (Exp_mask>>Exp_shift1))) ) {
1960 word0(d2) &= Frac_mask1;
1961 word0(d2) |= Exp_11;
1963 if ( (j = 11 - hi0bits(word0(d2) & Frac_mask)) )
1967 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1968 * log10(x) = log(x) / log(10)
1969 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1970 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1972 * This suggests computing an approximation k to log10(d) by
1974 * k = (i - Bias)*0.301029995663981
1975 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1977 * We want k to be too large rather than too small.
1978 * The error in the first-order Taylor series approximation
1979 * is in our favor, so we just round up the constant enough
1980 * to compensate for any error in the multiplication of
1981 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1982 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1983 * adding 1e-13 to the constant term more than suffices.
1984 * Hence we adjust the constant term to 0.1760912590558.
1985 * (We could get a more accurate k by invoking log10,
1986 * but this is probably not worthwhile.)
1994 #ifndef Sudden_Underflow
1997 /* d is denormalized */
1999 i = bbits + be + (Bias + (P-1) - 1);
2000 x = i > 32 ? ((word0(d) << (64 - i)) | (word1(d) >> (i - 32)))
2001 : (word1(d) << (32 - i));
2003 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2004 i -= (Bias + (P-1) - 1) + 1;
2008 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2010 if (ds < 0. && ds != k)
2011 k--; /* want k = floor(ds) */
2013 if (k >= 0 && k <= Ten_pmax) {
2035 if (mode < 0 || mode > 9)
2056 ilim = ilim1 = i = ndigits;
2062 i = ndigits + k + 1;
2068 *resultp = (char *) malloc(i + 1);
2071 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2073 /* Try to get by with floating-point arithmetic. */
2079 ieps = 2; /* conservative */
2084 /* prevent overflows */
2086 d /= bigtens[n_bigtens-1];
2089 for (; j; j >>= 1, i++)
2095 } else if ( (j1 = -k) ) {
2096 d *= tens[j1 & 0xf];
2097 for (j = j1 >> 4; j; j >>= 1, i++)
2103 if (k_check && d < 1. && ilim > 0) {
2112 word0(eps) -= (P-1)*Exp_msk1;
2122 #ifndef No_leftright
2124 /* Use Steele & White method of only
2125 * generating digits needed.
2127 eps = 0.5/tens[ilim-1] - eps;
2131 *s++ = '0' + (int)L;
2143 /* Generate ilim digits, then fix them up. */
2144 eps *= tens[ilim-1];
2145 for (i = 1;; i++, d *= 10.) {
2148 *s++ = '0' + (int)L;
2152 else if (d < 0.5 - eps) {
2153 while (*--s == '0');
2160 #ifndef No_leftright
2170 /* Do we have a "small" integer? */
2172 if (be >= 0 && k <= Int_max) {
2175 if (ndigits < 0 && ilim <= 0) {
2177 if (ilim < 0 || d <= 5*ds)
2184 #ifdef Check_FLT_ROUNDS
2185 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2191 *s++ = '0' + (int)L;
2194 if (d > ds || (d == ds && L & 1)) {
2218 #ifndef Sudden_Underflow
2219 denorm ? be + (Bias + (P-1) - 1 + 1) :
2222 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2235 if ((i = ilim) < 0) {
2244 if (m2 > 0 && s2 > 0) {
2245 i = m2 < s2 ? m2 : s2;
2253 mhi = pow5mult(mhi, m5);
2258 if ( (j = b5 - m5) )
2261 b = pow5mult(b, b5);
2265 S = pow5mult(S, s5);
2267 /* Check for special case that d is a normalized power of 2. */
2270 if (!word1(d) && !(word0(d) & Bndry_mask)
2271 #ifndef Sudden_Underflow
2272 && word0(d) & Exp_mask
2275 /* The special case */
2283 /* Arrange for convenient computation of quotients:
2284 * shift left if necessary so divisor has 4 leading 0 bits.
2286 * Perhaps we should just compute leading 28 bits of S once
2287 * and for all and pass them and a shift to quorem, so it
2288 * can do shifts and ors to compute the numerator for q.
2291 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) )
2294 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) )
2315 b = multadd(b, 10, 0); /* we botched the k estimate */
2317 mhi = multadd(mhi, 10, 0);
2321 if (ilim <= 0 && mode > 2) {
2322 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2323 /* no digits, fcvt style */
2335 mhi = lshift(mhi, m2);
2337 /* Compute mlo -- check for special case
2338 * that d is a normalized power of 2.
2343 mhi = Balloc(mhi->k);
2345 mhi = lshift(mhi, Log2P);
2349 dig = quorem(b,S) + '0';
2350 /* Do we yet have the shortest decimal string
2351 * that will round to d?
2354 delta = diff(S, mhi);
2355 j1 = delta->sign ? 1 : cmp(b, delta);
2357 #ifndef ROUND_BIASED
2358 if (j1 == 0 && !mode && !(word1(d) & 1)) {
2367 if (j < 0 || (j == 0 && !mode
2368 #ifndef ROUND_BIASED
2375 if ((j1 > 0 || (j1 == 0 && dig & 1))
2383 if (dig == '9') { /* possible if i == 1 */
2394 b = multadd(b, 10, 0);
2396 mlo = mhi = multadd(mhi, 10, 0);
2398 mlo = multadd(mlo, 10, 0);
2399 mhi = multadd(mhi, 10, 0);
2404 *s++ = dig = quorem(b,S) + '0';
2407 b = multadd(b, 10, 0);
2410 /* Round off last digit */
2414 if (j > 0 || (j == 0 && dig & 1)) {
2424 while (*--s == '0');
2430 if (mlo && mlo != mhi)
2436 if (s == s0) { /* don't return empty string */