2 * This strtod has been modified to not use values from the locale,
3 * but to hardcode the `.' as the separator. Our class libraries will
4 * make sure that only the dot is passed.
6 * This is so we do not call `setlocale' from our runtime before doing
7 * a strtod, because this could have unwanted effects in code that is
8 * co-hosted with the Mono runtime
10 * The entry point has been renamed `bsd_strtod'.
12 * Taken from the FreeBSD distribution.
20 * The Regents of the University of California. All rights reserved.
22 * Redistribution and use in source and binary forms, with or without
23 * modification, are permitted provided that the following conditions
25 * 1. Redistributions of source code must retain the above copyright
26 * notice, this list of conditions and the following disclaimer.
27 * 2. Redistributions in binary form must reproduce the above copyright
28 * notice, this list of conditions and the following disclaimer in the
29 * documentation and/or other materials provided with the distribution.
30 * 3. All advertising materials mentioning features or use of this software
31 * must display the following acknowledgement:
32 * This product includes software developed by the University of
33 * California, Berkeley and its contributors.
34 * 4. Neither the name of the University nor the names of its contributors
35 * may be used to endorse or promote products derived from this software
36 * without specific prior written permission.
38 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
39 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
40 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
41 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
42 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
43 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
44 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
45 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
46 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
47 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
50 * $FreeBSD: src/lib/libc/stdlib/strtod.c,v 1.3.8.3 2002/04/17 12:01:21 ache Exp $
53 #if defined(LIBC_SCCS) && !defined(lint)
54 static char sccsid[] = "@(#)strtod.c 8.1 (Berkeley) 6/4/93";
55 #endif /* LIBC_SCCS and not lint */
57 /****************************************************************
59 * The author of this software is David M. Gay.
61 * Copyright (c) 1991 by AT&T.
63 * Permission to use, copy, modify, and distribute this software for any
64 * purpose without fee is hereby granted, provided that this entire notice
65 * is included in all copies of any software which is or includes a copy
66 * or modification of this software and in all copies of the supporting
67 * documentation for such software.
69 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
70 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
71 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
72 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
74 ***************************************************************/
76 /* Please send bug reports to
78 AT&T Bell Laboratories, Room 2C-463
80 Murray Hill, NJ 07974-2070
82 dmg@research.att.com or research!dmg
85 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
87 * This strtod returns a nearest machine number to the input decimal
88 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
89 * broken by the IEEE round-even rule. Otherwise ties are broken by
90 * biased rounding (add half and chop).
92 * Inspired loosely by William D. Clinger's paper "How to Read Floating
93 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
97 * 1. We only require IEEE, IBM, or VAX double-precision
98 * arithmetic (not IEEE double-extended).
99 * 2. We get by with floating-point arithmetic in a case that
100 * Clinger missed -- when we're computing d * 10^n
101 * for a small integer d and the integer n is not too
102 * much larger than 22 (the maximum integer k for which
103 * we can represent 10^k exactly), we may be able to
104 * compute (d*10^k) * 10^(e-k) with just one roundoff.
105 * 3. Rather than a bit-at-a-time adjustment of the binary
106 * result in the hard case, we use floating-point
107 * arithmetic to determine the adjustment to within
108 * one bit; only in really hard cases do we need to
109 * compute a second residual.
110 * 4. Because of 3., we don't need a large table of powers of 10
111 * for ten-to-e (just some small tables, e.g. of 10^k
116 * #define IEEE_8087 for IEEE-arithmetic machines where the least
117 * significant byte has the lowest address.
118 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
119 * significant byte has the lowest address.
120 * #define Sudden_Underflow for IEEE-format machines without gradual
121 * underflow (i.e., that flush to zero on underflow).
122 * #define IBM for IBM mainframe-style floating-point arithmetic.
123 * #define VAX for VAX-style floating-point arithmetic.
124 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
125 * #define No_leftright to omit left-right logic in fast floating-point
126 * computation of dtoa.
127 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
128 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
129 * that use extended-precision instructions to compute rounded
130 * products and quotients) with IBM.
131 * #define ROUND_BIASED for IEEE-format with biased rounding.
132 * #define Inaccurate_Divide for IEEE-format with correctly rounded
133 * products but inaccurate quotients, e.g., for Intel i860.
134 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
135 * integer arithmetic. Whether this speeds things up or slows things
136 * down depends on the machine and the number being converted.
137 * #define KR_headers for old-style C function headers.
138 * #define Bad_float_h if your system lacks a float.h or if it does not
139 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
140 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
143 #if defined(i386) || defined(mips) && defined(MIPSEL) || defined (__arm__)
147 #elif defined(__x86_64__)
151 #elif defined(__ia64)
159 #elif defined(__hppa)
168 #define ULong guint32
172 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
194 #define IEEE_ARITHMETIC
197 #define IEEE_ARITHMETIC
199 #ifdef IEEE_ARITHMETIC
201 #define DBL_MAX_10_EXP 308
202 #define DBL_MAX_EXP 1024
205 #define DBL_MAX 1.7976931348623157e+308
210 #define DBL_MAX_10_EXP 75
211 #define DBL_MAX_EXP 63
214 #define DBL_MAX 7.2370055773322621e+75
219 #define DBL_MAX_10_EXP 38
220 #define DBL_MAX_EXP 127
223 #define DBL_MAX 1.7014118346046923e+38
227 #define LONG_MAX 2147483647
242 #define CONST /* blank */
248 #ifdef Unsigned_Shifts
249 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
251 #define Sign_Extend(a,b) /*no-op*/
254 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
255 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
259 #define word0(x) ((ULong *)&x)[1]
260 #define word1(x) ((ULong *)&x)[0]
262 #define word0(x) ((ULong *)&x)[0]
263 #define word1(x) ((ULong *)&x)[1]
266 /* The following definition of Storeinc is appropriate for MIPS processors.
267 * An alternative that might be better on some machines is
268 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
270 #if defined(IEEE_8087) + defined(VAX)
271 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
272 ((unsigned short *)a)[0] = (unsigned short)c, a++)
274 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
275 ((unsigned short *)a)[1] = (unsigned short)c, a++)
278 /* #define P DBL_MANT_DIG */
279 /* Ten_pmax = floor(P*log(2)/log(5)) */
280 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
281 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
282 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
284 #if defined(IEEE_8087) + defined(IEEE_MC68k)
286 #define Exp_shift1 20
287 #define Exp_msk1 0x100000
288 #define Exp_msk11 0x100000
289 #define Exp_mask 0x7ff00000
294 #define Exp_1 0x3ff00000
295 #define Exp_11 0x3ff00000
297 #define Frac_mask 0xfffff
298 #define Frac_mask1 0xfffff
301 #define Bndry_mask 0xfffff
302 #define Bndry_mask1 0xfffff
304 #define Sign_bit 0x80000000
310 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
312 #undef Sudden_Underflow
313 #define Sudden_Underflow
316 #define Exp_shift1 24
317 #define Exp_msk1 0x1000000
318 #define Exp_msk11 0x1000000
319 #define Exp_mask 0x7f000000
322 #define Exp_1 0x41000000
323 #define Exp_11 0x41000000
324 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
325 #define Frac_mask 0xffffff
326 #define Frac_mask1 0xffffff
329 #define Bndry_mask 0xefffff
330 #define Bndry_mask1 0xffffff
332 #define Sign_bit 0x80000000
334 #define Tiny0 0x100000
341 #define Exp_msk1 0x80
342 #define Exp_msk11 0x800000
343 #define Exp_mask 0x7f80
346 #define Exp_1 0x40800000
347 #define Exp_11 0x4080
349 #define Frac_mask 0x7fffff
350 #define Frac_mask1 0xffff007f
353 #define Bndry_mask 0xffff007f
354 #define Bndry_mask1 0xffff007f
356 #define Sign_bit 0x8000
370 #define rounded_product(a,b) a = rnd_prod(a, b)
371 #define rounded_quotient(a,b) a = rnd_quot(a, b)
373 extern double rnd_prod(), rnd_quot();
375 extern double rnd_prod(double, double), rnd_quot(double, double);
378 #define rounded_product(a,b) a *= b
379 #define rounded_quotient(a,b) a /= b
382 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
383 #define Big1 0xffffffff
386 /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
387 * This makes some inner loops simpler and sometimes saves work
388 * during multiplications, but it often seems to make things slightly
389 * slower. Hence the default is now to store 32 bits per long.
399 extern "C" double bsd_strtod(const char *s00, char **se);
400 extern "C" char *__dtoa(double d, int mode, int ndigits,
401 int *decpt, int *sign, char **rve, char **resultp);
407 int k, maxwds, sign, wds;
411 typedef struct Bigint Bigint;
425 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(Long));
428 rv->sign = rv->wds = 0;
443 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
444 y->wds*sizeof(Long) + 2*sizeof(int))
449 (b, m, a) Bigint *b; int m, a;
451 (Bigint *b, int m, int a) /* multiply by m and add a */
467 y = (xi & 0xffff) * m + a;
468 z = (xi >> 16) * m + (y >> 16);
470 *x++ = (z << 16) + (y & 0xffff);
478 if (wds >= b->maxwds) {
493 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
495 (CONST char *s, int nd0, int nd, ULong y9)
503 for (k = 0, y = 1; x > y; y <<= 1, k++) ;
510 b->x[0] = y9 & 0xffff;
511 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
518 b = multadd(b, 10, *s++ - '0');
524 b = multadd(b, 10, *s++ - '0');
531 (x) register ULong x;
538 if (!(x & 0xffff0000)) {
542 if (!(x & 0xff000000)) {
546 if (!(x & 0xf0000000)) {
550 if (!(x & 0xc0000000)) {
554 if (!(x & 0x80000000)) {
556 if (!(x & 0x40000000))
571 register ULong x = *y;
629 (a, b) Bigint *a, *b;
631 (Bigint *a, Bigint *b)
637 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
642 if (a->wds < b->wds) {
654 for (x = c->x, xa = x + wc; x < xa; x++)
662 for (; xb < xbe; xb++, xc0++) {
663 if ( (y = *xb & 0xffff) ) {
668 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
670 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
676 if ( (y = *xb >> 16) ) {
682 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
685 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
692 for (; xb < xbe; xc0++) {
698 z = *x++ * y + *xc + carry;
706 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
716 (b, k) Bigint *b; int k;
721 Bigint *b1, *p5, *p51;
723 static int p05[3] = { 5, 25, 125 };
726 b = multadd(b, p05[i-1], 0);
743 if (!(p51 = p5->next)) {
744 p51 = p5->next = mult(p5,p5);
755 (b, k) Bigint *b; int k;
762 ULong *x, *x1, *xe, z;
771 for (i = b->maxwds; n1 > i; i <<= 1)
775 for (i = 0; i < n; i++)
795 *x1++ = *x << k & 0xffff | z;
814 (a, b) Bigint *a, *b;
816 (Bigint *a, Bigint *b)
819 ULong *xa, *xa0, *xb, *xb0;
825 if (i > 1 && !a->x[i-1])
826 Bug("cmp called with a->x[a->wds-1] == 0");
827 if (j > 1 && !b->x[j-1])
828 Bug("cmp called with b->x[b->wds-1] == 0");
838 return *xa < *xb ? -1 : 1;
848 (a, b) Bigint *a, *b;
850 (Bigint *a, Bigint *b)
855 Long borrow, y; /* We need signed shifts here. */
856 ULong *xa, *xae, *xb, *xbe, *xc;
887 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
889 Sign_Extend(borrow, y);
890 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
892 Sign_Extend(borrow, z);
896 y = (*xa & 0xffff) + borrow;
898 Sign_Extend(borrow, y);
899 z = (*xa++ >> 16) + borrow;
901 Sign_Extend(borrow, z);
906 y = *xa++ - *xb++ + borrow;
908 Sign_Extend(borrow, y);
914 Sign_Extend(borrow, y);
935 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
936 #ifndef Sudden_Underflow
944 #ifndef Sudden_Underflow
948 word0(a) = 0x80000 >> L;
953 word1(a) = L >= 31 ? 1 : 1 << (31 - L);
963 (a, e) Bigint *a; int *e;
968 ULong *xa, *xa0, w, y, z;
982 if (!y) Bug("zero y in b2d");
988 d0 = Exp_1 | (y >> (Ebits - k));
989 w = xa > xa0 ? *--xa : 0;
990 d1 = (y << ((32-Ebits) + k)) | (w >> (Ebits - k));
993 z = xa > xa0 ? *--xa : 0;
995 d0 = Exp_1 | (y << k) | (z >> (32 - k));
996 y = xa > xa0 ? *--xa : 0;
997 d1 = (z << k) | (y >> (32 - k));
1003 if (k < Ebits + 16) {
1004 z = xa > xa0 ? *--xa : 0;
1005 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1006 w = xa > xa0 ? *--xa : 0;
1007 y = xa > xa0 ? *--xa : 0;
1008 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1011 z = xa > xa0 ? *--xa : 0;
1012 w = xa > xa0 ? *--xa : 0;
1014 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1015 y = xa > xa0 ? *--xa : 0;
1016 d1 = w << k + 16 | y << k;
1020 word0(d) = d0 >> 16 | d0 << 16;
1021 word1(d) = d1 >> 16 | d1 << 16;
1032 (d, e, bits) double d; int *e, *bits;
1034 (double d, int *e, int *bits)
1042 d0 = word0(d) >> 16 | word0(d) << 16;
1043 d1 = word1(d) >> 16 | word1(d) << 16;
1057 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1058 #ifdef Sudden_Underflow
1059 de = (int)(d0 >> Exp_shift);
1064 if ( (de = (int)(d0 >> Exp_shift)) )
1069 if ( (k = lo0bits(&y)) ) {
1070 x[0] = y | (z << (32 - k));
1075 i = b->wds = (x[1] = z) ? 2 : 1;
1079 Bug("Zero passed to d2b");
1088 if (k = lo0bits(&y))
1090 x[0] = y | z << 32 - k & 0xffff;
1091 x[1] = z >> k - 16 & 0xffff;
1096 x[1] = y >> 16 | z << 16 - k & 0xffff;
1097 x[2] = z >> k & 0xffff;
1111 Bug("Zero passed to d2b");
1128 #ifndef Sudden_Underflow
1132 *e = (de - Bias - (P-1) << 2) + k;
1133 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1135 *e = de - Bias - (P-1) + k;
1138 #ifndef Sudden_Underflow
1140 *e = de - Bias - (P-1) + 1 + k;
1142 *bits = 32*i - hi0bits(x[i-1]);
1144 *bits = (i+2)*16 - hi0bits(x[i]);
1156 (a, b) Bigint *a, *b;
1158 (Bigint *a, Bigint *b)
1167 k = ka - kb + 32*(a->wds - b->wds);
1169 k = ka - kb + 16*(a->wds - b->wds);
1173 word0(da) += (k >> 2)*Exp_msk1;
1178 word0(db) += (k >> 2)*Exp_msk1;
1184 word0(da) += k*Exp_msk1;
1187 word0(db) += k*Exp_msk1;
1195 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1196 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1205 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1206 static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1210 bigtens[] = { 1e16, 1e32, 1e64 };
1211 static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1214 bigtens[] = { 1e16, 1e32 };
1215 static double tinytens[] = { 1e-16, 1e-32 };
1223 (s00, se) CONST char *s00; char **se;
1225 (CONST char *s00, char **se)
1228 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1229 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1230 CONST char *s, *s0, *s1;
1231 double aadj, aadj1, adj, rv, rv0;
1234 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1235 char decimal_point = '.';
1237 sign = nz0 = nz = 0;
1239 for (s = s00;;s++) switch(*s) {
1251 if (isspace((unsigned char)*s))
1258 while (*++s == '0') ;
1264 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1270 if ((char)c == decimal_point) {
1273 for (; c == '0'; c = *++s)
1275 if (c > '0' && c <= '9') {
1283 for (; c >= '0' && c <= '9'; c = *++s) {
1288 for (i = 1; i < nz; i++)
1291 else if (nd <= DBL_DIG + 1)
1295 else if (nd <= DBL_DIG + 1)
1303 if (c == 'e' || c == 'E') {
1304 if (!nd && !nz && !nz0) {
1316 if (c >= '0' && c <= '9') {
1319 if (c > '0' && c <= '9') {
1322 while ((c = *++s) >= '0' && c <= '9')
1324 if (s - s1 > 8 || L > 19999)
1325 /* Avoid confusion from exponents
1326 * so large that e might overflow.
1328 e = 19999; /* safe for 16 bit ints */
1345 /* Now we have nd0 digits, starting at s0, followed by a
1346 * decimal point, followed by nd-nd0 digits. The number we're
1347 * after is the integer represented by those digits times
1352 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1355 rv = tens[k - 9] * rv + z;
1357 #ifndef RND_PRODQUOT
1364 if (e <= Ten_pmax) {
1366 goto vax_ovfl_check;
1368 /* rv = */ rounded_product(rv, tens[e]);
1373 if (e <= Ten_pmax + i) {
1374 /* A fancier test would sometimes let us do
1375 * this for larger i values.
1380 /* VAX exponent range is so narrow we must
1381 * worry about overflow here...
1384 word0(rv) -= P*Exp_msk1;
1385 /* rv = */ rounded_product(rv, tens[e]);
1386 if ((word0(rv) & Exp_mask)
1387 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1389 word0(rv) += P*Exp_msk1;
1391 /* rv = */ rounded_product(rv, tens[e]);
1396 #ifndef Inaccurate_Divide
1397 else if (e >= -Ten_pmax) {
1398 /* rv = */ rounded_quotient(rv, tens[-e]);
1405 /* Get starting approximation = rv * 10**e1 */
1408 if ( (i = e1 & 15) )
1410 if ( (e1 &= ~15) ) {
1411 if (e1 > DBL_MAX_10_EXP) {
1417 /* Can't trust HUGE_VAL */
1419 word0(rv) = Exp_mask;
1429 for (j = 0; e1 > 1; j++, e1 >>= 1)
1432 /* The last multiplication could overflow. */
1433 word0(rv) -= P*Exp_msk1;
1435 if ((z = word0(rv) & Exp_mask)
1436 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1438 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1439 /* set to largest number */
1440 /* (Can't trust DBL_MAX) */
1445 word0(rv) += P*Exp_msk1;
1448 } else if (e1 < 0) {
1450 if ( (i = e1 & 15) )
1452 if ( (e1 &= ~15) ) {
1454 for (j = 0; e1 > 1; j++, e1 >>= 1)
1457 /* The last multiplication could underflow. */
1471 /* The refinement below will clean
1472 * this approximation up.
1478 /* Now the hard part -- adjusting rv to the correct value.*/
1480 /* Put digits into bd: true value = bd * 10^e */
1482 bd0 = s2b(s0, nd0, nd, y);
1485 bd = Balloc(bd0->k);
1487 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1502 #ifdef Sudden_Underflow
1504 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1509 i = bbe + bbbits - 1; /* logb(rv) */
1510 if (i < Emin) /* denormal */
1517 i = bb2 < bd2 ? bb2 : bd2;
1526 bs = pow5mult(bs, bb5);
1532 bb = lshift(bb, bb2);
1534 bd = pow5mult(bd, bd5);
1536 bd = lshift(bd, bd2);
1538 bs = lshift(bs, bs2);
1539 delta = diff(bb, bd);
1540 dsign = delta->sign;
1544 /* Error is less than half an ulp -- check for
1545 * special case of mantissa a power of two.
1547 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1549 delta = lshift(delta,Log2P);
1550 if (cmp(delta, bs) > 0)
1555 /* exactly half-way between */
1557 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1558 && word1(rv) == 0xffffffff) {
1559 /*boundary case -- increment exponent*/
1560 word0(rv) = (word0(rv) & Exp_mask)
1569 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1571 /* boundary case -- decrement exponent */
1572 #ifdef Sudden_Underflow
1573 L = word0(rv) & Exp_mask;
1582 L = (word0(rv) & Exp_mask) - Exp_msk1;
1584 word0(rv) = L | Bndry_mask1;
1585 word1(rv) = 0xffffffff;
1592 #ifndef ROUND_BIASED
1593 if (!(word1(rv) & LSB))
1598 #ifndef ROUND_BIASED
1601 #ifndef Sudden_Underflow
1609 if ((aadj = ratio(delta, bs)) <= 2.) {
1612 else if (word1(rv) || word0(rv) & Bndry_mask) {
1613 #ifndef Sudden_Underflow
1614 if (word1(rv) == Tiny1 && !word0(rv))
1620 /* special case -- power of FLT_RADIX to be */
1621 /* rounded down... */
1623 if (aadj < 2./FLT_RADIX)
1624 aadj = 1./FLT_RADIX;
1631 aadj1 = dsign ? aadj : -aadj;
1632 #ifdef Check_FLT_ROUNDS
1633 switch(FLT_ROUNDS) {
1634 case 2: /* towards +infinity */
1637 case 0: /* towards 0 */
1638 case 3: /* towards -infinity */
1642 if (FLT_ROUNDS == 0)
1646 y = word0(rv) & Exp_mask;
1648 /* Check for overflow */
1650 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1652 word0(rv) -= P*Exp_msk1;
1653 adj = aadj1 * ulp(rv);
1655 if ((word0(rv) & Exp_mask) >=
1656 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1657 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1663 word0(rv) += P*Exp_msk1;
1665 #ifdef Sudden_Underflow
1666 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1668 word0(rv) += P*Exp_msk1;
1669 adj = aadj1 * ulp(rv);
1672 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1674 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1677 if (word0(rv0) == Tiny0
1678 && word1(rv0) == Tiny1)
1684 word0(rv) -= P*Exp_msk1;
1686 adj = aadj1 * ulp(rv);
1690 /* Compute adj so that the IEEE rounding rules will
1691 * correctly round rv + adj in some half-way cases.
1692 * If rv * ulp(rv) is denormalized (i.e.,
1693 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1694 * trouble from bits lost to denormalization;
1695 * example: 1.2e-307 .
1697 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1698 aadj1 = (double)(int)(aadj + 0.5);
1702 adj = aadj1 * ulp(rv);
1706 z = word0(rv) & Exp_mask;
1708 /* Can we stop now? */
1711 /* The tolerances below are conservative. */
1712 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1713 if (aadj < .4999999 || aadj > .5000001)
1715 } else if (aadj < .4999999/FLT_RADIX)
1732 return sign ? -rv : rv;
1738 (b, S) Bigint *b, *S;
1740 (Bigint *b, Bigint *S)
1746 ULong *bx, *bxe, *sx, *sxe;
1754 /*debug*/ if (b->wds > n)
1755 /*debug*/ Bug("oversize b in quorem");
1763 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1765 /*debug*/ if (q > 9)
1766 /*debug*/ Bug("oversized quotient in quorem");
1774 ys = (si & 0xffff) * q + carry;
1775 zs = (si >> 16) * q + (ys >> 16);
1777 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1779 Sign_Extend(borrow, y);
1780 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1782 Sign_Extend(borrow, z);
1785 ys = *sx++ * q + carry;
1787 y = *bx - (ys & 0xffff) + borrow;
1789 Sign_Extend(borrow, y);
1792 } while (sx <= sxe);
1795 while (--bxe > bx && !*bxe)
1800 if (cmp(b, S) >= 0) {
1809 ys = (si & 0xffff) + carry;
1810 zs = (si >> 16) + (ys >> 16);
1812 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1814 Sign_Extend(borrow, y);
1815 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1817 Sign_Extend(borrow, z);
1822 y = *bx - (ys & 0xffff) + borrow;
1824 Sign_Extend(borrow, y);
1827 } while (sx <= sxe);
1831 while (--bxe > bx && !*bxe)
1839 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1841 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1842 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1845 * 1. Rather than iterating, we use a simple numeric overestimate
1846 * to determine k = floor(log10(d)). We scale relevant
1847 * quantities using O(log2(k)) rather than O(k) multiplications.
1848 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1849 * try to generate digits strictly left to right. Instead, we
1850 * compute with fewer bits and propagate the carry if necessary
1851 * when rounding the final digit up. This is often faster.
1852 * 3. Under the assumption that input will be rounded nearest,
1853 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1854 * That is, we allow equality in stopping tests when the
1855 * round-nearest rule will give the same floating-point value
1856 * as would satisfaction of the stopping test with strict
1858 * 4. We remove common factors of powers of 2 from relevant
1860 * 5. When converting floating-point integers less than 1e16,
1861 * we use floating-point arithmetic rather than resorting
1862 * to multiple-precision integers.
1863 * 6. When asked to produce fewer than 15 digits, we first try
1864 * to get by with floating-point arithmetic; we resort to
1865 * multiple-precision integer arithmetic only if we cannot
1866 * guarantee that the floating-point calculation has given
1867 * the correctly rounded result. For k requested digits and
1868 * "uniformly" distributed input, the probability is
1869 * something like 10^(k-15) that we must resort to the long
1876 (d, mode, ndigits, decpt, sign, rve, resultp)
1877 double d; int mode, ndigits, *decpt, *sign; char **rve, **resultp;
1879 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve,
1883 /* Arguments ndigits, decpt, sign are similar to those
1884 of ecvt and fcvt; trailing zeros are suppressed from
1885 the returned string. If not null, *rve is set to point
1886 to the end of the return value. If d is +-Infinity or NaN,
1887 then *decpt is set to 9999.
1890 0 ==> shortest string that yields d when read in
1891 and rounded to nearest.
1892 1 ==> like 0, but with Steele & White stopping rule;
1893 e.g. with IEEE P754 arithmetic , mode 0 gives
1894 1e23 whereas mode 1 gives 9.999999999999999e22.
1895 2 ==> max(1,ndigits) significant digits. This gives a
1896 return value similar to that of ecvt, except
1897 that trailing zeros are suppressed.
1898 3 ==> through ndigits past the decimal point. This
1899 gives a return value similar to that from fcvt,
1900 except that trailing zeros are suppressed, and
1901 ndigits can be negative.
1902 4-9 should give the same return values as 2-3, i.e.,
1903 4 <= mode <= 9 ==> same return as mode
1904 2 + (mode & 1). These modes are mainly for
1905 debugging; often they run slower but sometimes
1906 faster than modes 2-3.
1907 4,5,8,9 ==> left-to-right digit generation.
1908 6-9 ==> don't try fast floating-point estimate
1911 Values of mode other than 0-9 are treated as mode 0.
1913 Sufficient space is allocated to the return value
1914 to hold the suppressed trailing zeros.
1917 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1918 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1919 spec_case, try_quick;
1921 #ifndef Sudden_Underflow
1925 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1929 if (word0(d) & Sign_bit) {
1930 /* set sign for everything, including 0's and NaNs */
1932 word0(d) &= ~Sign_bit; /* clear sign bit */
1937 #if defined(IEEE_Arith) + defined(VAX)
1939 if ((word0(d) & Exp_mask) == Exp_mask)
1941 if (word0(d) == 0x8000)
1944 /* Infinity or NaN */
1949 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
1952 *resultp = s = malloc (strlen (ss) + 1);
1964 d += 0; /* normalize */
1968 *resultp = s = malloc (2);
1976 b = d2b(d, &be, &bbits);
1977 #ifdef Sudden_Underflow
1978 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1980 if ( (i = (int)((word0(d) >> Exp_shift1) & (Exp_mask>>Exp_shift1))) ) {
1983 word0(d2) &= Frac_mask1;
1984 word0(d2) |= Exp_11;
1986 if ( (j = 11 - hi0bits(word0(d2) & Frac_mask)) )
1990 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1991 * log10(x) = log(x) / log(10)
1992 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1993 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1995 * This suggests computing an approximation k to log10(d) by
1997 * k = (i - Bias)*0.301029995663981
1998 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2000 * We want k to be too large rather than too small.
2001 * The error in the first-order Taylor series approximation
2002 * is in our favor, so we just round up the constant enough
2003 * to compensate for any error in the multiplication of
2004 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2005 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2006 * adding 1e-13 to the constant term more than suffices.
2007 * Hence we adjust the constant term to 0.1760912590558.
2008 * (We could get a more accurate k by invoking log10,
2009 * but this is probably not worthwhile.)
2017 #ifndef Sudden_Underflow
2020 /* d is denormalized */
2022 i = bbits + be + (Bias + (P-1) - 1);
2023 x = i > 32 ? ((word0(d) << (64 - i)) | (word1(d) >> (i - 32)))
2024 : (word1(d) << (32 - i));
2026 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2027 i -= (Bias + (P-1) - 1) + 1;
2031 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2033 if (ds < 0. && ds != k)
2034 k--; /* want k = floor(ds) */
2036 if (k >= 0 && k <= Ten_pmax) {
2058 if (mode < 0 || mode > 9)
2079 ilim = ilim1 = i = ndigits;
2085 i = ndigits + k + 1;
2091 *resultp = (char *) malloc(i + 1);
2094 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2096 /* Try to get by with floating-point arithmetic. */
2102 ieps = 2; /* conservative */
2107 /* prevent overflows */
2109 d /= bigtens[n_bigtens-1];
2112 for (; j; j >>= 1, i++)
2118 } else if ( (j1 = -k) ) {
2119 d *= tens[j1 & 0xf];
2120 for (j = j1 >> 4; j; j >>= 1, i++)
2126 if (k_check && d < 1. && ilim > 0) {
2135 word0(eps) -= (P-1)*Exp_msk1;
2145 #ifndef No_leftright
2147 /* Use Steele & White method of only
2148 * generating digits needed.
2150 eps = 0.5/tens[ilim-1] - eps;
2154 *s++ = '0' + (int)L;
2166 /* Generate ilim digits, then fix them up. */
2167 eps *= tens[ilim-1];
2168 for (i = 1;; i++, d *= 10.) {
2171 *s++ = '0' + (int)L;
2175 else if (d < 0.5 - eps) {
2176 while (*--s == '0');
2183 #ifndef No_leftright
2193 /* Do we have a "small" integer? */
2195 if (be >= 0 && k <= Int_max) {
2198 if (ndigits < 0 && ilim <= 0) {
2200 if (ilim < 0 || d <= 5*ds)
2207 #ifdef Check_FLT_ROUNDS
2208 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2214 *s++ = '0' + (int)L;
2217 if (d > ds || (d == ds && L & 1)) {
2241 #ifndef Sudden_Underflow
2242 denorm ? be + (Bias + (P-1) - 1 + 1) :
2245 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2258 if ((i = ilim) < 0) {
2267 if (m2 > 0 && s2 > 0) {
2268 i = m2 < s2 ? m2 : s2;
2276 mhi = pow5mult(mhi, m5);
2281 if ( (j = b5 - m5) )
2284 b = pow5mult(b, b5);
2288 S = pow5mult(S, s5);
2290 /* Check for special case that d is a normalized power of 2. */
2293 if (!word1(d) && !(word0(d) & Bndry_mask)
2294 #ifndef Sudden_Underflow
2295 && word0(d) & Exp_mask
2298 /* The special case */
2306 /* Arrange for convenient computation of quotients:
2307 * shift left if necessary so divisor has 4 leading 0 bits.
2309 * Perhaps we should just compute leading 28 bits of S once
2310 * and for all and pass them and a shift to quorem, so it
2311 * can do shifts and ors to compute the numerator for q.
2314 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) )
2317 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) )
2338 b = multadd(b, 10, 0); /* we botched the k estimate */
2340 mhi = multadd(mhi, 10, 0);
2344 if (ilim <= 0 && mode > 2) {
2345 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2346 /* no digits, fcvt style */
2358 mhi = lshift(mhi, m2);
2360 /* Compute mlo -- check for special case
2361 * that d is a normalized power of 2.
2366 mhi = Balloc(mhi->k);
2368 mhi = lshift(mhi, Log2P);
2372 dig = quorem(b,S) + '0';
2373 /* Do we yet have the shortest decimal string
2374 * that will round to d?
2377 delta = diff(S, mhi);
2378 j1 = delta->sign ? 1 : cmp(b, delta);
2380 #ifndef ROUND_BIASED
2381 if (j1 == 0 && !mode && !(word1(d) & 1)) {
2390 if (j < 0 || (j == 0 && !mode
2391 #ifndef ROUND_BIASED
2398 if ((j1 > 0 || (j1 == 0 && dig & 1))
2406 if (dig == '9') { /* possible if i == 1 */
2417 b = multadd(b, 10, 0);
2419 mlo = mhi = multadd(mhi, 10, 0);
2421 mlo = multadd(mlo, 10, 0);
2422 mhi = multadd(mhi, 10, 0);
2427 *s++ = dig = quorem(b,S) + '0';
2430 b = multadd(b, 10, 0);
2433 /* Round off last digit */
2437 if (j > 0 || (j == 0 && dig & 1)) {
2447 while (*--s == '0');
2453 if (mlo && mlo != mhi)
2459 if (s == s0) { /* don't return empty string */