3 * The Regents of the University of California. All rights reserved.
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 * 3. All advertising materials mentioning features or use of this software
14 * must display the following acknowledgement:
15 * This product includes software developed by the University of
16 * California, Berkeley and its contributors.
17 * 4. Neither the name of the University nor the names of its contributors
18 * may be used to endorse or promote products derived from this software
19 * without specific prior written permission.
21 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
22 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
25 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
30 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
34 #if defined(LIBC_SCCS) && !defined(lint)
35 static char sccsid[] = "@(#)strtod.c 8.1 (Berkeley) 6/4/93";
36 #endif /* LIBC_SCCS and not lint */
37 #include <sys/cdefs.h>
38 __FBSDID("$FreeBSD$");
40 /****************************************************************
42 * The author of this software is David M. Gay.
44 * Copyright (c) 1991 by AT&T.
46 * Permission to use, copy, modify, and distribute this software for any
47 * purpose without fee is hereby granted, provided that this entire notice
48 * is included in all copies of any software which is or includes a copy
49 * or modification of this software and in all copies of the supporting
50 * documentation for such software.
52 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
53 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
54 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
55 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
57 ***************************************************************/
59 /* Please send bug reports to
61 AT&T Bell Laboratories, Room 2C-463
63 Murray Hill, NJ 07974-2070
65 dmg@research.att.com or research!dmg
68 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
70 * This strtod returns a nearest machine number to the input decimal
71 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
72 * broken by the IEEE round-even rule. Otherwise ties are broken by
73 * biased rounding (add half and chop).
75 * Inspired loosely by William D. Clinger's paper "How to Read Floating
76 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
80 * 1. We only require IEEE, IBM, or VAX double-precision
81 * arithmetic (not IEEE double-extended).
82 * 2. We get by with floating-point arithmetic in a case that
83 * Clinger missed -- when we're computing d * 10^n
84 * for a small integer d and the integer n is not too
85 * much larger than 22 (the maximum integer k for which
86 * we can represent 10^k exactly), we may be able to
87 * compute (d*10^k) * 10^(e-k) with just one roundoff.
88 * 3. Rather than a bit-at-a-time adjustment of the binary
89 * result in the hard case, we use floating-point
90 * arithmetic to determine the adjustment to within
91 * one bit; only in really hard cases do we need to
92 * compute a second residual.
93 * 4. Because of 3., we don't need a large table of powers of 10
94 * for ten-to-e (just some small tables, e.g. of 10^k
99 * #define Sudden_Underflow for IEEE-format machines without gradual
100 * underflow (i.e., that flush to zero on underflow).
101 * #define IBM for IBM mainframe-style floating-point arithmetic.
102 * #define VAX for VAX-style floating-point arithmetic.
103 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
104 * #define No_leftright to omit left-right logic in fast floating-point
105 * computation of dtoa.
106 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
107 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
108 * that use extended-precision instructions to compute rounded
109 * products and quotients) with IBM.
110 * #define ROUND_BIASED for IEEE-format with biased rounding.
111 * #define Inaccurate_Divide for IEEE-format with correctly rounded
112 * products but inaccurate quotients, e.g., for Intel i860.
113 * #define Just_16 to store 16 bits per 32-bit Long when doing high-precision
114 * integer arithmetic. Whether this speeds things up or slows things
115 * down depends on the machine and the number being converted.
116 * #define KR_headers for old-style C function headers.
117 * #define Bad_float_h if your system lacks a float.h or if it does not
118 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
119 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
122 #if defined(__i386__) || defined(__ia64__) || defined(__alpha__) || \
123 defined(__sparc64__) || defined(__powerpc__)
124 #include <sys/types.h>
125 #if BYTE_ORDER == BIG_ENDIAN
126 #define IEEE_BIG_ENDIAN
128 #define IEEE_LITTLE_ENDIAN
130 #endif /* defined(__i386__) ... */
132 typedef int32_t Long;
133 typedef u_int32_t ULong;
137 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
159 #ifdef IEEE_BIG_ENDIAN
160 #define IEEE_ARITHMETIC
162 #ifdef IEEE_LITTLE_ENDIAN
163 #define IEEE_ARITHMETIC
165 #ifdef IEEE_ARITHMETIC
167 #define DBL_MAX_10_EXP 308
168 #define DBL_MAX_EXP 1024
171 #define DBL_MAX 1.7976931348623157e+308
176 #define DBL_MAX_10_EXP 75
177 #define DBL_MAX_EXP 63
180 #define DBL_MAX 7.2370055773322621e+75
185 #define DBL_MAX_10_EXP 38
186 #define DBL_MAX_EXP 127
189 #define DBL_MAX 1.7014118346046923e+38
193 #define LONG_MAX 2147483647
208 #define CONST /* blank */
214 #ifdef Unsigned_Shifts
215 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
217 #define Sign_Extend(a,b) /*no-op*/
220 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + \
222 Only one of IEEE_LITTLE_ENDIAN, IEEE_BIG_ENDIAN, VAX, or IBM should be defined.
225 #ifdef IEEE_LITTLE_ENDIAN
226 #define word0(x) ((ULong *)&x)[1]
227 #define word1(x) ((ULong *)&x)[0]
229 #define word0(x) ((ULong *)&x)[0]
230 #define word1(x) ((ULong *)&x)[1]
233 /* The following definition of Storeinc is appropriate for MIPS processors.
234 * An alternative that might be better on some machines is
235 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
237 #if defined(IEEE_LITTLE_ENDIAN) + defined(VAX)
238 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
239 ((unsigned short *)a)[0] = (unsigned short)c, a++)
241 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
242 ((unsigned short *)a)[1] = (unsigned short)c, a++)
245 /* #define P DBL_MANT_DIG */
246 /* Ten_pmax = floor(P*log(2)/log(5)) */
247 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
248 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
249 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
251 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN)
253 #define Exp_shift1 20
254 #define Exp_msk1 0x100000
255 #define Exp_msk11 0x100000
256 #define Exp_mask 0x7ff00000
261 #define Exp_1 0x3ff00000
262 #define Exp_11 0x3ff00000
264 #define Frac_mask 0xfffff
265 #define Frac_mask1 0xfffff
268 #define Bndry_mask 0xfffff
269 #define Bndry_mask1 0xfffff
271 #define Sign_bit 0x80000000
277 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
279 #undef Sudden_Underflow
280 #define Sudden_Underflow
283 #define Exp_shift1 24
284 #define Exp_msk1 0x1000000
285 #define Exp_msk11 0x1000000
286 #define Exp_mask 0x7f000000
289 #define Exp_1 0x41000000
290 #define Exp_11 0x41000000
291 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
292 #define Frac_mask 0xffffff
293 #define Frac_mask1 0xffffff
296 #define Bndry_mask 0xefffff
297 #define Bndry_mask1 0xffffff
299 #define Sign_bit 0x80000000
301 #define Tiny0 0x100000
308 #define Exp_msk1 0x80
309 #define Exp_msk11 0x800000
310 #define Exp_mask 0x7f80
313 #define Exp_1 0x40800000
314 #define Exp_11 0x4080
316 #define Frac_mask 0x7fffff
317 #define Frac_mask1 0xffff007f
320 #define Bndry_mask 0xffff007f
321 #define Bndry_mask1 0xffff007f
323 #define Sign_bit 0x8000
337 #define rounded_product(a,b) a = rnd_prod(a, b)
338 #define rounded_quotient(a,b) a = rnd_quot(a, b)
340 extern double rnd_prod(), rnd_quot();
342 extern double rnd_prod(double, double), rnd_quot(double, double);
345 #define rounded_product(a,b) a *= b
346 #define rounded_quotient(a,b) a /= b
349 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
350 #define Big1 0xffffffff
353 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
354 * This makes some inner loops simpler and sometimes saves work
355 * during multiplications, but it often seems to make things slightly
356 * slower. Hence the default is now to store 32 bits per Long.
366 extern "C" double strtod(const char *s00, char **se);
367 extern "C" char *__dtoa(double d, int mode, int ndigits,
368 int *decpt, int *sign, char **rve, char **resultp);
374 int k, maxwds, sign, wds;
378 typedef struct Bigint Bigint;
392 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(Long));
395 rv->sign = rv->wds = 0;
410 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
411 y->wds*sizeof(Long) + 2*sizeof(int))
416 (b, m, a) Bigint *b; int m, a;
418 (Bigint *b, int m, int a) /* multiply by m and add a */
434 y = (xi & 0xffff) * m + a;
435 z = (xi >> 16) * m + (y >> 16);
437 *x++ = (z << 16) + (y & 0xffff);
445 if (wds >= b->maxwds) {
460 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
462 (CONST char *s, int nd0, int nd, ULong y9)
470 for (k = 0, y = 1; x > y; y <<= 1, k++) ;
477 b->x[0] = y9 & 0xffff;
478 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
485 b = multadd(b, 10, *s++ - '0');
491 b = multadd(b, 10, *s++ - '0');
505 if (!(x & 0xffff0000)) {
509 if (!(x & 0xff000000)) {
513 if (!(x & 0xf0000000)) {
517 if (!(x & 0xc0000000)) {
521 if (!(x & 0x80000000)) {
523 if (!(x & 0x40000000))
596 (a, b) Bigint *a, *b;
598 (Bigint *a, Bigint *b)
604 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
609 if (a->wds < b->wds) {
621 for (x = c->x, xa = x + wc; x < xa; x++)
629 for (; xb < xbe; xb++, xc0++) {
630 if ( (y = *xb & 0xffff) ) {
635 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
637 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
643 if ( (y = *xb >> 16) ) {
649 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
652 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
659 for (; xb < xbe; xc0++) {
665 z = *x++ * y + *xc + carry;
673 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
683 (b, k) Bigint *b; int k;
688 Bigint *b1, *p5, *p51;
690 static int p05[3] = { 5, 25, 125 };
693 b = multadd(b, p05[i-1], 0);
710 if (!(p51 = p5->next)) {
711 p51 = p5->next = mult(p5,p5);
722 (b, k) Bigint *b; int k;
729 ULong *x, *x1, *xe, z;
738 for (i = b->maxwds; n1 > i; i <<= 1)
742 for (i = 0; i < n; i++)
762 *x1++ = *x << k & 0xffff | z;
781 (a, b) Bigint *a, *b;
783 (Bigint *a, Bigint *b)
786 ULong *xa, *xa0, *xb, *xb0;
792 if (i > 1 && !a->x[i-1])
793 Bug("cmp called with a->x[a->wds-1] == 0");
794 if (j > 1 && !b->x[j-1])
795 Bug("cmp called with b->x[b->wds-1] == 0");
805 return *xa < *xb ? -1 : 1;
815 (a, b) Bigint *a, *b;
817 (Bigint *a, Bigint *b)
822 Long borrow, y; /* We need signed shifts here. */
823 ULong *xa, *xae, *xb, *xbe, *xc;
854 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
856 Sign_Extend(borrow, y);
857 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
859 Sign_Extend(borrow, z);
863 y = (*xa & 0xffff) + borrow;
865 Sign_Extend(borrow, y);
866 z = (*xa++ >> 16) + borrow;
868 Sign_Extend(borrow, z);
873 y = *xa++ - *xb++ + borrow;
875 Sign_Extend(borrow, y);
881 Sign_Extend(borrow, y);
902 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
903 #ifndef Sudden_Underflow
911 #ifndef Sudden_Underflow
915 word0(a) = 0x80000 >> L;
920 word1(a) = L >= 31 ? 1 : 1 << (31 - L);
930 (a, e) Bigint *a; int *e;
935 ULong *xa, *xa0, w, y, z;
949 if (!y) Bug("zero y in b2d");
955 d0 = Exp_1 | (y >> (Ebits - k));
956 w = xa > xa0 ? *--xa : 0;
957 d1 = (y << ((32-Ebits) + k)) | (w >> (Ebits - k));
960 z = xa > xa0 ? *--xa : 0;
962 d0 = Exp_1 | (y << k) | (z >> (32 - k));
963 y = xa > xa0 ? *--xa : 0;
964 d1 = (z << k) | (y >> (32 - k));
970 if (k < Ebits + 16) {
971 z = xa > xa0 ? *--xa : 0;
972 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
973 w = xa > xa0 ? *--xa : 0;
974 y = xa > xa0 ? *--xa : 0;
975 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
978 z = xa > xa0 ? *--xa : 0;
979 w = xa > xa0 ? *--xa : 0;
981 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
982 y = xa > xa0 ? *--xa : 0;
983 d1 = w << k + 16 | y << k;
987 word0(d) = d0 >> 16 | d0 << 16;
988 word1(d) = d1 >> 16 | d1 << 16;
999 (d, e, bits) double d; int *e, *bits;
1001 (double d, int *e, int *bits)
1009 d0 = word0(d) >> 16 | word0(d) << 16;
1010 d1 = word1(d) >> 16 | word1(d) << 16;
1024 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1025 #ifdef Sudden_Underflow
1026 de = (int)(d0 >> Exp_shift);
1031 if ( (de = (int)(d0 >> Exp_shift)) )
1036 if ( (k = lo0bits(&y)) ) {
1037 x[0] = y | (z << (32 - k));
1042 i = b->wds = (x[1] = z) ? 2 : 1;
1046 Bug("Zero passed to d2b");
1055 if (k = lo0bits(&y))
1057 x[0] = y | z << 32 - k & 0xffff;
1058 x[1] = z >> k - 16 & 0xffff;
1063 x[1] = y >> 16 | z << 16 - k & 0xffff;
1064 x[2] = z >> k & 0xffff;
1078 Bug("Zero passed to d2b");
1095 #ifndef Sudden_Underflow
1099 *e = (de - Bias - (P-1) << 2) + k;
1100 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1102 *e = de - Bias - (P-1) + k;
1105 #ifndef Sudden_Underflow
1107 *e = de - Bias - (P-1) + 1 + k;
1109 *bits = 32*i - hi0bits(x[i-1]);
1111 *bits = (i+2)*16 - hi0bits(x[i]);
1123 (a, b) Bigint *a, *b;
1125 (Bigint *a, Bigint *b)
1134 k = ka - kb + 32*(a->wds - b->wds);
1136 k = ka - kb + 16*(a->wds - b->wds);
1140 word0(da) += (k >> 2)*Exp_msk1;
1145 word0(db) += (k >> 2)*Exp_msk1;
1151 word0(da) += k*Exp_msk1;
1154 word0(db) += k*Exp_msk1;
1162 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1163 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1172 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1173 static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1177 bigtens[] = { 1e16, 1e32, 1e64 };
1178 static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1181 bigtens[] = { 1e16, 1e32 };
1182 static double tinytens[] = { 1e-16, 1e-32 };
1190 (s00, se) CONST char *s00; char **se;
1192 (CONST char *s00, char **se)
1195 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1196 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1197 CONST char *s, *s0, *s1;
1198 double aadj, aadj1, adj, rv, rv0;
1201 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1202 char decimal_point = localeconv()->decimal_point[0];
1204 sign = nz0 = nz = 0;
1206 for (s = s00;;s++) switch(*s) {
1218 if (isspace((unsigned char)*s))
1225 while (*++s == '0') ;
1231 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1237 if ((char)c == decimal_point) {
1240 for (; c == '0'; c = *++s)
1242 if (c > '0' && c <= '9') {
1250 for (; c >= '0' && c <= '9'; c = *++s) {
1255 for (i = 1; i < nz; i++)
1258 else if (nd <= DBL_DIG + 1)
1262 else if (nd <= DBL_DIG + 1)
1270 if (c == 'e' || c == 'E') {
1271 if (!nd && !nz && !nz0) {
1283 if (c >= '0' && c <= '9') {
1286 if (c > '0' && c <= '9') {
1289 while ((c = *++s) >= '0' && c <= '9')
1291 if (s - s1 > 8 || L > 19999)
1292 /* Avoid confusion from exponents
1293 * so large that e might overflow.
1295 e = 19999; /* safe for 16 bit ints */
1312 /* Now we have nd0 digits, starting at s0, followed by a
1313 * decimal point, followed by nd-nd0 digits. The number we're
1314 * after is the integer represented by those digits times
1319 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1322 rv = tens[k - 9] * rv + z;
1324 #ifndef RND_PRODQUOT
1331 if (e <= Ten_pmax) {
1333 goto vax_ovfl_check;
1335 /* rv = */ rounded_product(rv, tens[e]);
1340 if (e <= Ten_pmax + i) {
1341 /* A fancier test would sometimes let us do
1342 * this for larger i values.
1347 /* VAX exponent range is so narrow we must
1348 * worry about overflow here...
1351 word0(rv) -= P*Exp_msk1;
1352 /* rv = */ rounded_product(rv, tens[e]);
1353 if ((word0(rv) & Exp_mask)
1354 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1356 word0(rv) += P*Exp_msk1;
1358 /* rv = */ rounded_product(rv, tens[e]);
1363 #ifndef Inaccurate_Divide
1364 else if (e >= -Ten_pmax) {
1365 /* rv = */ rounded_quotient(rv, tens[-e]);
1372 /* Get starting approximation = rv * 10**e1 */
1375 if ( (i = e1 & 15) )
1377 if ( (e1 &= ~15) ) {
1378 if (e1 > DBL_MAX_10_EXP) {
1385 for (j = 0; e1 > 1; j++, e1 >>= 1)
1388 /* The last multiplication could overflow. */
1389 word0(rv) -= P*Exp_msk1;
1391 if ((z = word0(rv) & Exp_mask)
1392 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1394 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1395 /* set to largest number */
1396 /* (Can't trust DBL_MAX) */
1401 word0(rv) += P*Exp_msk1;
1404 } else if (e1 < 0) {
1406 if ( (i = e1 & 15) )
1408 if ( (e1 &= ~15) ) {
1410 for (j = 0; e1 > 1; j++, e1 >>= 1)
1413 /* The last multiplication could underflow. */
1427 /* The refinement below will clean
1428 * this approximation up.
1434 /* Now the hard part -- adjusting rv to the correct value.*/
1436 /* Put digits into bd: true value = bd * 10^e */
1438 bd0 = s2b(s0, nd0, nd, y);
1441 bd = Balloc(bd0->k);
1443 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1458 #ifdef Sudden_Underflow
1460 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1465 i = bbe + bbbits - 1; /* logb(rv) */
1466 if (i < Emin) /* denormal */
1473 i = bb2 < bd2 ? bb2 : bd2;
1482 bs = pow5mult(bs, bb5);
1488 bb = lshift(bb, bb2);
1490 bd = pow5mult(bd, bd5);
1492 bd = lshift(bd, bd2);
1494 bs = lshift(bs, bs2);
1495 delta = diff(bb, bd);
1496 dsign = delta->sign;
1500 /* Error is less than half an ulp -- check for
1501 * special case of mantissa a power of two.
1503 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1505 delta = lshift(delta,Log2P);
1506 if (cmp(delta, bs) > 0)
1511 /* exactly half-way between */
1513 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1514 && word1(rv) == 0xffffffff) {
1515 /*boundary case -- increment exponent*/
1516 word0(rv) = (word0(rv) & Exp_mask)
1525 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1527 /* boundary case -- decrement exponent */
1528 #ifdef Sudden_Underflow
1529 L = word0(rv) & Exp_mask;
1538 L = (word0(rv) & Exp_mask) - Exp_msk1;
1540 word0(rv) = L | Bndry_mask1;
1541 word1(rv) = 0xffffffff;
1548 #ifndef ROUND_BIASED
1549 if (!(word1(rv) & LSB))
1554 #ifndef ROUND_BIASED
1557 #ifndef Sudden_Underflow
1565 if ((aadj = ratio(delta, bs)) <= 2.) {
1568 else if (word1(rv) || word0(rv) & Bndry_mask) {
1569 #ifndef Sudden_Underflow
1570 if (word1(rv) == Tiny1 && !word0(rv))
1576 /* special case -- power of FLT_RADIX to be */
1577 /* rounded down... */
1579 if (aadj < 2./FLT_RADIX)
1580 aadj = 1./FLT_RADIX;
1587 aadj1 = dsign ? aadj : -aadj;
1588 #ifdef Check_FLT_ROUNDS
1589 switch(FLT_ROUNDS) {
1590 case 2: /* towards +infinity */
1593 case 0: /* towards 0 */
1594 case 3: /* towards -infinity */
1598 if (FLT_ROUNDS == 0)
1602 y = word0(rv) & Exp_mask;
1604 /* Check for overflow */
1606 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1608 word0(rv) -= P*Exp_msk1;
1609 adj = aadj1 * ulp(rv);
1611 if ((word0(rv) & Exp_mask) >=
1612 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1613 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1619 word0(rv) += P*Exp_msk1;
1621 #ifdef Sudden_Underflow
1622 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1624 word0(rv) += P*Exp_msk1;
1625 adj = aadj1 * ulp(rv);
1628 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1630 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1633 if (word0(rv0) == Tiny0
1634 && word1(rv0) == Tiny1)
1640 word0(rv) -= P*Exp_msk1;
1642 adj = aadj1 * ulp(rv);
1646 /* Compute adj so that the IEEE rounding rules will
1647 * correctly round rv + adj in some half-way cases.
1648 * If rv * ulp(rv) is denormalized (i.e.,
1649 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1650 * trouble from bits lost to denormalization;
1651 * example: 1.2e-307 .
1653 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1654 aadj1 = (double)(int)(aadj + 0.5);
1658 adj = aadj1 * ulp(rv);
1662 z = word0(rv) & Exp_mask;
1664 /* Can we stop now? */
1667 /* The tolerances below are conservative. */
1668 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1669 if (aadj < .4999999 || aadj > .5000001)
1671 } else if (aadj < .4999999/FLT_RADIX)
1688 return sign ? -rv : rv;
1694 (b, S) Bigint *b, *S;
1696 (Bigint *b, Bigint *S)
1702 ULong *bx, *bxe, *sx, *sxe;
1710 /*debug*/ if (b->wds > n)
1711 /*debug*/ Bug("oversize b in quorem");
1719 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1721 /*debug*/ if (q > 9)
1722 /*debug*/ Bug("oversized quotient in quorem");
1730 ys = (si & 0xffff) * q + carry;
1731 zs = (si >> 16) * q + (ys >> 16);
1733 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1735 Sign_Extend(borrow, y);
1736 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1738 Sign_Extend(borrow, z);
1741 ys = *sx++ * q + carry;
1743 y = *bx - (ys & 0xffff) + borrow;
1745 Sign_Extend(borrow, y);
1748 } while (sx <= sxe);
1751 while (--bxe > bx && !*bxe)
1756 if (cmp(b, S) >= 0) {
1765 ys = (si & 0xffff) + carry;
1766 zs = (si >> 16) + (ys >> 16);
1768 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1770 Sign_Extend(borrow, y);
1771 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1773 Sign_Extend(borrow, z);
1778 y = *bx - (ys & 0xffff) + borrow;
1780 Sign_Extend(borrow, y);
1783 } while (sx <= sxe);
1787 while (--bxe > bx && !*bxe)
1795 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1797 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1798 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1801 * 1. Rather than iterating, we use a simple numeric overestimate
1802 * to determine k = floor(log10(d)). We scale relevant
1803 * quantities using O(log2(k)) rather than O(k) multiplications.
1804 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1805 * try to generate digits strictly left to right. Instead, we
1806 * compute with fewer bits and propagate the carry if necessary
1807 * when rounding the final digit up. This is often faster.
1808 * 3. Under the assumption that input will be rounded nearest,
1809 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1810 * That is, we allow equality in stopping tests when the
1811 * round-nearest rule will give the same floating-point value
1812 * as would satisfaction of the stopping test with strict
1814 * 4. We remove common factors of powers of 2 from relevant
1816 * 5. When converting floating-point integers less than 1e16,
1817 * we use floating-point arithmetic rather than resorting
1818 * to multiple-precision integers.
1819 * 6. When asked to produce fewer than 15 digits, we first try
1820 * to get by with floating-point arithmetic; we resort to
1821 * multiple-precision integer arithmetic only if we cannot
1822 * guarantee that the floating-point calculation has given
1823 * the correctly rounded result. For k requested digits and
1824 * "uniformly" distributed input, the probability is
1825 * something like 10^(k-15) that we must resort to the Long
1832 (d, mode, ndigits, decpt, sign, rve, resultp)
1833 double d; int mode, ndigits, *decpt, *sign; char **rve, **resultp;
1835 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve,
1839 /* Arguments ndigits, decpt, sign are similar to those
1840 of ecvt and fcvt; trailing zeros are suppressed from
1841 the returned string. If not null, *rve is set to point
1842 to the end of the return value. If d is +-Infinity or NaN,
1843 then *decpt is set to 9999.
1846 0 ==> shortest string that yields d when read in
1847 and rounded to nearest.
1848 1 ==> like 0, but with Steele & White stopping rule;
1849 e.g. with IEEE P754 arithmetic , mode 0 gives
1850 1e23 whereas mode 1 gives 9.999999999999999e22.
1851 2 ==> max(1,ndigits) significant digits. This gives a
1852 return value similar to that of ecvt, except
1853 that trailing zeros are suppressed.
1854 3 ==> through ndigits past the decimal point. This
1855 gives a return value similar to that from fcvt,
1856 except that trailing zeros are suppressed, and
1857 ndigits can be negative.
1858 4-9 should give the same return values as 2-3, i.e.,
1859 4 <= mode <= 9 ==> same return as mode
1860 2 + (mode & 1). These modes are mainly for
1861 debugging; often they run slower but sometimes
1862 faster than modes 2-3.
1863 4,5,8,9 ==> left-to-right digit generation.
1864 6-9 ==> don't try fast floating-point estimate
1867 Values of mode other than 0-9 are treated as mode 0.
1869 Sufficient space is allocated to the return value
1870 to hold the suppressed trailing zeros.
1873 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1874 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1875 spec_case, try_quick;
1877 #ifndef Sudden_Underflow
1881 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1885 if (word0(d) & Sign_bit) {
1886 /* set sign for everything, including 0's and NaNs */
1888 word0(d) &= ~Sign_bit; /* clear sign bit */
1893 #if defined(IEEE_Arith) + defined(VAX)
1895 if ((word0(d) & Exp_mask) == Exp_mask)
1897 if (word0(d) == 0x8000)
1900 /* Infinity or NaN */
1904 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
1917 d += 0; /* normalize */
1927 b = d2b(d, &be, &bbits);
1928 #ifdef Sudden_Underflow
1929 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1931 if ( (i = (int)((word0(d) >> Exp_shift1) & (Exp_mask>>Exp_shift1))) ) {
1934 word0(d2) &= Frac_mask1;
1935 word0(d2) |= Exp_11;
1937 if ( (j = 11 - hi0bits(word0(d2) & Frac_mask)) )
1941 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1942 * log10(x) = log(x) / log(10)
1943 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1944 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1946 * This suggests computing an approximation k to log10(d) by
1948 * k = (i - Bias)*0.301029995663981
1949 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1951 * We want k to be too large rather than too small.
1952 * The error in the first-order Taylor series approximation
1953 * is in our favor, so we just round up the constant enough
1954 * to compensate for any error in the multiplication of
1955 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1956 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1957 * adding 1e-13 to the constant term more than suffices.
1958 * Hence we adjust the constant term to 0.1760912590558.
1959 * (We could get a more accurate k by invoking log10,
1960 * but this is probably not worthwhile.)
1968 #ifndef Sudden_Underflow
1971 /* d is denormalized */
1973 i = bbits + be + (Bias + (P-1) - 1);
1974 x = i > 32 ? ((word0(d) << (64 - i)) | (word1(d) >> (i - 32)))
1975 : (word1(d) << (32 - i));
1977 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
1978 i -= (Bias + (P-1) - 1) + 1;
1982 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1984 if (ds < 0. && ds != k)
1985 k--; /* want k = floor(ds) */
1987 if (k >= 0 && k <= Ten_pmax) {
2009 if (mode < 0 || mode > 9)
2030 ilim = ilim1 = i = ndigits;
2036 i = ndigits + k + 1;
2042 *resultp = (char *) malloc(i + 1);
2045 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2047 /* Try to get by with floating-point arithmetic. */
2053 ieps = 2; /* conservative */
2058 /* prevent overflows */
2060 d /= bigtens[n_bigtens-1];
2063 for (; j; j >>= 1, i++)
2069 } else if ( (j1 = -k) ) {
2070 d *= tens[j1 & 0xf];
2071 for (j = j1 >> 4; j; j >>= 1, i++)
2077 if (k_check && d < 1. && ilim > 0) {
2086 word0(eps) -= (P-1)*Exp_msk1;
2096 #ifndef No_leftright
2098 /* Use Steele & White method of only
2099 * generating digits needed.
2101 eps = 0.5/tens[ilim-1] - eps;
2105 *s++ = '0' + (int)L;
2117 /* Generate ilim digits, then fix them up. */
2118 eps *= tens[ilim-1];
2119 for (i = 1;; i++, d *= 10.) {
2122 *s++ = '0' + (int)L;
2126 else if (d < 0.5 - eps) {
2127 while (*--s == '0');
2134 #ifndef No_leftright
2144 /* Do we have a "small" integer? */
2146 if (be >= 0 && k <= Int_max) {
2149 if (ndigits < 0 && ilim <= 0) {
2151 if (ilim < 0 || d <= 5*ds)
2158 #ifdef Check_FLT_ROUNDS
2159 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2165 *s++ = '0' + (int)L;
2168 if (d > ds || (d == ds && L & 1)) {
2192 #ifndef Sudden_Underflow
2193 denorm ? be + (Bias + (P-1) - 1 + 1) :
2196 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2209 if ((i = ilim) < 0) {
2218 if (m2 > 0 && s2 > 0) {
2219 i = m2 < s2 ? m2 : s2;
2227 mhi = pow5mult(mhi, m5);
2232 if ( (j = b5 - m5) )
2235 b = pow5mult(b, b5);
2239 S = pow5mult(S, s5);
2241 /* Check for special case that d is a normalized power of 2. */
2244 if (!word1(d) && !(word0(d) & Bndry_mask)
2245 #ifndef Sudden_Underflow
2246 && word0(d) & Exp_mask
2249 /* The special case */
2257 /* Arrange for convenient computation of quotients:
2258 * shift left if necessary so divisor has 4 leading 0 bits.
2260 * Perhaps we should just compute leading 28 bits of S once
2261 * and for all and pass them and a shift to quorem, so it
2262 * can do shifts and ors to compute the numerator for q.
2265 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) )
2268 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) )
2289 b = multadd(b, 10, 0); /* we botched the k estimate */
2291 mhi = multadd(mhi, 10, 0);
2295 if (ilim <= 0 && mode > 2) {
2296 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2297 /* no digits, fcvt style */
2309 mhi = lshift(mhi, m2);
2311 /* Compute mlo -- check for special case
2312 * that d is a normalized power of 2.
2317 mhi = Balloc(mhi->k);
2319 mhi = lshift(mhi, Log2P);
2323 dig = quorem(b,S) + '0';
2324 /* Do we yet have the shortest decimal string
2325 * that will round to d?
2328 delta = diff(S, mhi);
2329 j1 = delta->sign ? 1 : cmp(b, delta);
2331 #ifndef ROUND_BIASED
2332 if (j1 == 0 && !mode && !(word1(d) & 1)) {
2341 if (j < 0 || (j == 0 && !mode
2342 #ifndef ROUND_BIASED
2349 if ((j1 > 0 || (j1 == 0 && dig & 1))
2357 if (dig == '9') { /* possible if i == 1 */
2368 b = multadd(b, 10, 0);
2370 mlo = mhi = multadd(mhi, 10, 0);
2372 mlo = multadd(mlo, 10, 0);
2373 mhi = multadd(mhi, 10, 0);
2378 *s++ = dig = quorem(b,S) + '0';
2381 b = multadd(b, 10, 0);
2384 /* Round off last digit */
2388 if (j > 0 || (j == 0 && dig & 1)) {
2398 while (*--s == '0');
2404 if (mlo && mlo != mhi)
2410 if (s == s0) { /* don't return empty string */