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
33 * $FreeBSD: src/lib/libc/stdlib/strtod.c,v 1.3.8.4 2002/08/31 22:26:35 dwmalone Exp $
36 #if defined(LIBC_SCCS) && !defined(lint)
37 static char sccsid[] = "@(#)strtod.c 8.1 (Berkeley) 6/4/93";
38 #endif /* LIBC_SCCS and not lint */
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 IEEE_8087 for IEEE-arithmetic machines where the least
100 * significant byte has the lowest address.
101 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
102 * significant byte has the lowest address.
103 * #define Sudden_Underflow for IEEE-format machines without gradual
104 * underflow (i.e., that flush to zero on underflow).
105 * #define IBM for IBM mainframe-style floating-point arithmetic.
106 * #define VAX for VAX-style floating-point arithmetic.
107 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
108 * #define No_leftright to omit left-right logic in fast floating-point
109 * computation of dtoa.
110 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
111 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
112 * that use extended-precision instructions to compute rounded
113 * products and quotients) with IBM.
114 * #define ROUND_BIASED for IEEE-format with biased rounding.
115 * #define Inaccurate_Divide for IEEE-format with correctly rounded
116 * products but inaccurate quotients, e.g., for Intel i860.
117 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
118 * integer arithmetic. Whether this speeds things up or slows things
119 * down depends on the machine and the number being converted.
120 * #define KR_headers for old-style C function headers.
121 * #define Bad_float_h if your system lacks a float.h or if it does not
122 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
123 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
126 #if defined(i386) || defined(mips) && defined(MIPSEL)
134 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
156 #define IEEE_ARITHMETIC
159 #define IEEE_ARITHMETIC
161 #ifdef IEEE_ARITHMETIC
163 #define DBL_MAX_10_EXP 308
164 #define DBL_MAX_EXP 1024
167 #define DBL_MAX 1.7976931348623157e+308
172 #define DBL_MAX_10_EXP 75
173 #define DBL_MAX_EXP 63
176 #define DBL_MAX 7.2370055773322621e+75
181 #define DBL_MAX_10_EXP 38
182 #define DBL_MAX_EXP 127
185 #define DBL_MAX 1.7014118346046923e+38
189 #define LONG_MAX 2147483647
204 #define CONST /* blank */
210 #ifdef Unsigned_Shifts
211 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
213 #define Sign_Extend(a,b) /*no-op*/
216 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
217 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
220 union doubleasulongs {
225 #define word0(x) (((union doubleasulongs *)&x)->w)[1]
226 #define word1(x) (((union doubleasulongs *)&x)->w)[0]
228 #define word0(x) (((union doubleasulongs *)&x)->w)[0]
229 #define word1(x) (((union doubleasulongs *)&x)->w)[1]
232 /* The following definition of Storeinc is appropriate for MIPS processors.
233 * An alternative that might be better on some machines is
234 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
236 #if defined(IEEE_8087) + defined(VAX)
237 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
238 ((unsigned short *)a)[0] = (unsigned short)c, a++)
240 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
241 ((unsigned short *)a)[1] = (unsigned short)c, a++)
244 /* #define P DBL_MANT_DIG */
245 /* Ten_pmax = floor(P*log(2)/log(5)) */
246 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
247 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
248 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
250 #if defined(IEEE_8087) + defined(IEEE_MC68k)
252 #define Exp_shift1 20
253 #define Exp_msk1 0x100000
254 #define Exp_msk11 0x100000
255 #define Exp_mask 0x7ff00000
260 #define Exp_1 0x3ff00000
261 #define Exp_11 0x3ff00000
263 #define Frac_mask 0xfffff
264 #define Frac_mask1 0xfffff
267 #define Bndry_mask 0xfffff
268 #define Bndry_mask1 0xfffff
270 #define Sign_bit 0x80000000
276 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
278 #undef Sudden_Underflow
279 #define Sudden_Underflow
282 #define Exp_shift1 24
283 #define Exp_msk1 0x1000000
284 #define Exp_msk11 0x1000000
285 #define Exp_mask 0x7f000000
288 #define Exp_1 0x41000000
289 #define Exp_11 0x41000000
290 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
291 #define Frac_mask 0xffffff
292 #define Frac_mask1 0xffffff
295 #define Bndry_mask 0xefffff
296 #define Bndry_mask1 0xffffff
298 #define Sign_bit 0x80000000
300 #define Tiny0 0x100000
307 #define Exp_msk1 0x80
308 #define Exp_msk11 0x800000
309 #define Exp_mask 0x7f80
312 #define Exp_1 0x40800000
313 #define Exp_11 0x4080
315 #define Frac_mask 0x7fffff
316 #define Frac_mask1 0xffff007f
319 #define Bndry_mask 0xffff007f
320 #define Bndry_mask1 0xffff007f
322 #define Sign_bit 0x8000
336 #define rounded_product(a,b) a = rnd_prod(a, b)
337 #define rounded_quotient(a,b) a = rnd_quot(a, b)
339 extern double rnd_prod(), rnd_quot();
341 extern double rnd_prod(double, double), rnd_quot(double, double);
344 #define rounded_product(a,b) a *= b
345 #define rounded_quotient(a,b) a /= b
348 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
349 #define Big1 0xffffffff
352 /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
353 * This makes some inner loops simpler and sometimes saves work
354 * during multiplications, but it often seems to make things slightly
355 * slower. Hence the default is now to store 32 bits per long.
365 extern "C" double strtod(const char *s00, char **se);
366 extern "C" char *__dtoa(double d, int mode, int ndigits,
367 int *decpt, int *sign, char **rve, char **resultp);
373 int k, maxwds, sign, wds;
377 typedef struct Bigint Bigint;
391 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(long));
394 rv->sign = rv->wds = 0;
409 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
410 y->wds*sizeof(long) + 2*sizeof(int))
415 (b, m, a) Bigint *b; int m, a;
417 (Bigint *b, int m, int a) /* multiply by m and add a */
433 y = (xi & 0xffff) * m + a;
434 z = (xi >> 16) * m + (y >> 16);
436 *x++ = (z << 16) + (y & 0xffff);
444 if (wds >= b->maxwds) {
459 (s, nd0, nd, y9) CONST char *s; int nd0, nd; unsigned long y9;
461 (CONST char *s, int nd0, int nd, unsigned long y9)
469 for (k = 0, y = 1; x > y; y <<= 1, k++) ;
476 b->x[0] = y9 & 0xffff;
477 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
484 b = multadd(b, 10, *s++ - '0');
490 b = multadd(b, 10, *s++ - '0');
497 (x) register unsigned long x;
499 (register unsigned long x)
504 if (!(x & 0xffff0000)) {
508 if (!(x & 0xff000000)) {
512 if (!(x & 0xf0000000)) {
516 if (!(x & 0xc0000000)) {
520 if (!(x & 0x80000000)) {
522 if (!(x & 0x40000000))
531 (y) unsigned long *y;
537 register unsigned long x = *y;
595 (a, b) Bigint *a, *b;
597 (Bigint *a, Bigint *b)
602 unsigned long carry, y, z;
603 unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
608 if (a->wds < b->wds) {
620 for (x = c->x, xa = x + wc; x < xa; x++)
628 for (; xb < xbe; xb++, xc0++) {
629 if ( (y = *xb & 0xffff) ) {
634 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
636 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
642 if ( (y = *xb >> 16) ) {
648 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
651 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
658 for (; xb < xbe; xc0++) {
664 z = *x++ * y + *xc + carry;
672 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
682 (b, k) Bigint *b; int k;
687 Bigint *b1, *p5, *p51;
689 static int p05[3] = { 5, 25, 125 };
692 b = multadd(b, p05[i-1], 0);
709 if (!(p51 = p5->next)) {
710 p51 = p5->next = mult(p5,p5);
721 (b, k) Bigint *b; int k;
728 unsigned long *x, *x1, *xe, z;
737 for (i = b->maxwds; n1 > i; i <<= 1)
741 for (i = 0; i < n; i++)
761 *x1++ = *x << k & 0xffff | z;
780 (a, b) Bigint *a, *b;
782 (Bigint *a, Bigint *b)
785 unsigned long *xa, *xa0, *xb, *xb0;
791 if (i > 1 && !a->x[i-1])
792 Bug("cmp called with a->x[a->wds-1] == 0");
793 if (j > 1 && !b->x[j-1])
794 Bug("cmp called with b->x[b->wds-1] == 0");
804 return *xa < *xb ? -1 : 1;
814 (a, b) Bigint *a, *b;
816 (Bigint *a, Bigint *b)
821 long borrow, y; /* We need signed shifts here. */
822 unsigned long *xa, *xae, *xb, *xbe, *xc;
853 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
855 Sign_Extend(borrow, y);
856 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
858 Sign_Extend(borrow, z);
862 y = (*xa & 0xffff) + borrow;
864 Sign_Extend(borrow, y);
865 z = (*xa++ >> 16) + borrow;
867 Sign_Extend(borrow, z);
872 y = *xa++ - *xb++ + borrow;
874 Sign_Extend(borrow, y);
880 Sign_Extend(borrow, y);
901 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
902 #ifndef Sudden_Underflow
910 #ifndef Sudden_Underflow
914 word0(a) = 0x80000 >> L;
919 word1(a) = L >= 31 ? 1 : 1 << (31 - L);
929 (a, e) Bigint *a; int *e;
934 unsigned long *xa, *xa0, w, y, z;
938 unsigned long d0, d1;
948 if (!y) Bug("zero y in b2d");
954 d0 = Exp_1 | (y >> (Ebits - k));
955 w = xa > xa0 ? *--xa : 0;
956 d1 = (y << ((32-Ebits) + k)) | (w >> (Ebits - k));
959 z = xa > xa0 ? *--xa : 0;
961 d0 = Exp_1 | (y << k) | (z >> (32 - k));
962 y = xa > xa0 ? *--xa : 0;
963 d1 = (z << k) | (y >> (32 - k));
969 if (k < Ebits + 16) {
970 z = xa > xa0 ? *--xa : 0;
971 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
972 w = xa > xa0 ? *--xa : 0;
973 y = xa > xa0 ? *--xa : 0;
974 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
977 z = xa > xa0 ? *--xa : 0;
978 w = xa > xa0 ? *--xa : 0;
980 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
981 y = xa > xa0 ? *--xa : 0;
982 d1 = w << k + 16 | y << k;
986 word0(d) = d0 >> 16 | d0 << 16;
987 word1(d) = d1 >> 16 | d1 << 16;
998 (d, e, bits) double d; int *e, *bits;
1000 (double d, int *e, int *bits)
1005 unsigned long *x, y, z;
1007 unsigned long d0, d1;
1008 d0 = word0(d) >> 16 | word0(d) << 16;
1009 d1 = word1(d) >> 16 | word1(d) << 16;
1023 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1024 #ifdef Sudden_Underflow
1025 de = (int)(d0 >> Exp_shift);
1030 if ( (de = (int)(d0 >> Exp_shift)) )
1035 if ( (k = lo0bits(&y)) ) {
1036 x[0] = y | (z << (32 - k));
1041 i = b->wds = (x[1] = z) ? 2 : 1;
1045 Bug("Zero passed to d2b");
1054 if (k = lo0bits(&y))
1056 x[0] = y | z << 32 - k & 0xffff;
1057 x[1] = z >> k - 16 & 0xffff;
1062 x[1] = y >> 16 | z << 16 - k & 0xffff;
1063 x[2] = z >> k & 0xffff;
1077 Bug("Zero passed to d2b");
1094 #ifndef Sudden_Underflow
1098 *e = (de - Bias - (P-1) << 2) + k;
1099 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1101 *e = de - Bias - (P-1) + k;
1104 #ifndef Sudden_Underflow
1106 *e = de - Bias - (P-1) + 1 + k;
1108 *bits = 32*i - hi0bits(x[i-1]);
1110 *bits = (i+2)*16 - hi0bits(x[i]);
1122 (a, b) Bigint *a, *b;
1124 (Bigint *a, Bigint *b)
1133 k = ka - kb + 32*(a->wds - b->wds);
1135 k = ka - kb + 16*(a->wds - b->wds);
1139 word0(da) += (k >> 2)*Exp_msk1;
1144 word0(db) += (k >> 2)*Exp_msk1;
1150 word0(da) += k*Exp_msk1;
1153 word0(db) += k*Exp_msk1;
1161 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1162 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1171 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1172 static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1176 bigtens[] = { 1e16, 1e32, 1e64 };
1177 static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1180 bigtens[] = { 1e16, 1e32 };
1181 static double tinytens[] = { 1e-16, 1e-32 };
1189 (s00, se) CONST char *s00; char **se;
1191 (CONST char *s00, char **se)
1194 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1195 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1196 CONST char *s, *s0, *s1;
1197 double aadj, aadj1, adj, rv, rv0;
1200 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1201 char decimal_point = localeconv()->decimal_point[0];
1203 sign = nz0 = nz = 0;
1205 for (s = s00;;s++) switch(*s) {
1217 if (isspace((unsigned char)*s))
1224 while (*++s == '0') ;
1230 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1236 if ((char)c == decimal_point) {
1239 for (; c == '0'; c = *++s)
1241 if (c > '0' && c <= '9') {
1249 for (; c >= '0' && c <= '9'; c = *++s) {
1254 for (i = 1; i < nz; i++)
1257 else if (nd <= DBL_DIG + 1)
1261 else if (nd <= DBL_DIG + 1)
1269 if (c == 'e' || c == 'E') {
1270 if (!nd && !nz && !nz0) {
1282 if (c >= '0' && c <= '9') {
1285 if (c > '0' && c <= '9') {
1288 while ((c = *++s) >= '0' && c <= '9')
1290 if (s - s1 > 8 || L > 19999)
1291 /* Avoid confusion from exponents
1292 * so large that e might overflow.
1294 e = 19999; /* safe for 16 bit ints */
1311 /* Now we have nd0 digits, starting at s0, followed by a
1312 * decimal point, followed by nd-nd0 digits. The number we're
1313 * after is the integer represented by those digits times
1318 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1321 rv = tens[k - 9] * rv + z;
1323 #ifndef RND_PRODQUOT
1330 if (e <= Ten_pmax) {
1332 goto vax_ovfl_check;
1334 /* rv = */ rounded_product(rv, tens[e]);
1339 if (e <= Ten_pmax + i) {
1340 /* A fancier test would sometimes let us do
1341 * this for larger i values.
1346 /* VAX exponent range is so narrow we must
1347 * worry about overflow here...
1350 word0(rv) -= P*Exp_msk1;
1351 /* rv = */ rounded_product(rv, tens[e]);
1352 if ((word0(rv) & Exp_mask)
1353 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1355 word0(rv) += P*Exp_msk1;
1357 /* rv = */ rounded_product(rv, tens[e]);
1362 #ifndef Inaccurate_Divide
1363 else if (e >= -Ten_pmax) {
1364 /* rv = */ rounded_quotient(rv, tens[-e]);
1371 /* Get starting approximation = rv * 10**e1 */
1374 if ( (i = e1 & 15) )
1376 if ( (e1 &= ~15) ) {
1377 if (e1 > DBL_MAX_10_EXP) {
1383 /* Can't trust HUGE_VAL */
1385 word0(rv) = Exp_mask;
1395 for (j = 0; e1 > 1; j++, e1 >>= 1)
1398 /* The last multiplication could overflow. */
1399 word0(rv) -= P*Exp_msk1;
1401 if ((z = word0(rv) & Exp_mask)
1402 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1404 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1405 /* set to largest number */
1406 /* (Can't trust DBL_MAX) */
1411 word0(rv) += P*Exp_msk1;
1414 } else if (e1 < 0) {
1416 if ( (i = e1 & 15) )
1418 if ( (e1 &= ~15) ) {
1420 for (j = 0; e1 > 1; j++, e1 >>= 1)
1423 /* The last multiplication could underflow. */
1437 /* The refinement below will clean
1438 * this approximation up.
1444 /* Now the hard part -- adjusting rv to the correct value.*/
1446 /* Put digits into bd: true value = bd * 10^e */
1448 bd0 = s2b(s0, nd0, nd, y);
1451 bd = Balloc(bd0->k);
1453 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1468 #ifdef Sudden_Underflow
1470 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1475 i = bbe + bbbits - 1; /* logb(rv) */
1476 if (i < Emin) /* denormal */
1483 i = bb2 < bd2 ? bb2 : bd2;
1492 bs = pow5mult(bs, bb5);
1498 bb = lshift(bb, bb2);
1500 bd = pow5mult(bd, bd5);
1502 bd = lshift(bd, bd2);
1504 bs = lshift(bs, bs2);
1505 delta = diff(bb, bd);
1506 dsign = delta->sign;
1510 /* Error is less than half an ulp -- check for
1511 * special case of mantissa a power of two.
1513 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1515 delta = lshift(delta,Log2P);
1516 if (cmp(delta, bs) > 0)
1521 /* exactly half-way between */
1523 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1524 && word1(rv) == 0xffffffff) {
1525 /*boundary case -- increment exponent*/
1526 word0(rv) = (word0(rv) & Exp_mask)
1535 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1537 /* boundary case -- decrement exponent */
1538 #ifdef Sudden_Underflow
1539 L = word0(rv) & Exp_mask;
1548 L = (word0(rv) & Exp_mask) - Exp_msk1;
1550 word0(rv) = L | Bndry_mask1;
1551 word1(rv) = 0xffffffff;
1558 #ifndef ROUND_BIASED
1559 if (!(word1(rv) & LSB))
1564 #ifndef ROUND_BIASED
1567 #ifndef Sudden_Underflow
1575 if ((aadj = ratio(delta, bs)) <= 2.) {
1578 else if (word1(rv) || word0(rv) & Bndry_mask) {
1579 #ifndef Sudden_Underflow
1580 if (word1(rv) == Tiny1 && !word0(rv))
1586 /* special case -- power of FLT_RADIX to be */
1587 /* rounded down... */
1589 if (aadj < 2./FLT_RADIX)
1590 aadj = 1./FLT_RADIX;
1597 aadj1 = dsign ? aadj : -aadj;
1598 #ifdef Check_FLT_ROUNDS
1599 switch(FLT_ROUNDS) {
1600 case 2: /* towards +infinity */
1603 case 0: /* towards 0 */
1604 case 3: /* towards -infinity */
1608 if (FLT_ROUNDS == 0)
1612 y = word0(rv) & Exp_mask;
1614 /* Check for overflow */
1616 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1618 word0(rv) -= P*Exp_msk1;
1619 adj = aadj1 * ulp(rv);
1621 if ((word0(rv) & Exp_mask) >=
1622 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1623 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1629 word0(rv) += P*Exp_msk1;
1631 #ifdef Sudden_Underflow
1632 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1634 word0(rv) += P*Exp_msk1;
1635 adj = aadj1 * ulp(rv);
1638 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1640 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1643 if (word0(rv0) == Tiny0
1644 && word1(rv0) == Tiny1)
1650 word0(rv) -= P*Exp_msk1;
1652 adj = aadj1 * ulp(rv);
1656 /* Compute adj so that the IEEE rounding rules will
1657 * correctly round rv + adj in some half-way cases.
1658 * If rv * ulp(rv) is denormalized (i.e.,
1659 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1660 * trouble from bits lost to denormalization;
1661 * example: 1.2e-307 .
1663 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1664 aadj1 = (double)(int)(aadj + 0.5);
1668 adj = aadj1 * ulp(rv);
1672 z = word0(rv) & Exp_mask;
1674 /* Can we stop now? */
1677 /* The tolerances below are conservative. */
1678 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1679 if (aadj < .4999999 || aadj > .5000001)
1681 } else if (aadj < .4999999/FLT_RADIX)
1698 return sign ? -rv : rv;
1704 (b, S) Bigint *b, *S;
1706 (Bigint *b, Bigint *S)
1711 unsigned long carry, q, ys;
1712 unsigned long *bx, *bxe, *sx, *sxe;
1715 unsigned long si, zs;
1720 /*debug*/ if (b->wds > n)
1721 /*debug*/ Bug("oversize b in quorem");
1729 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1731 /*debug*/ if (q > 9)
1732 /*debug*/ Bug("oversized quotient in quorem");
1740 ys = (si & 0xffff) * q + carry;
1741 zs = (si >> 16) * q + (ys >> 16);
1743 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1745 Sign_Extend(borrow, y);
1746 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1748 Sign_Extend(borrow, z);
1751 ys = *sx++ * q + carry;
1753 y = *bx - (ys & 0xffff) + borrow;
1755 Sign_Extend(borrow, y);
1758 } while (sx <= sxe);
1761 while (--bxe > bx && !*bxe)
1766 if (cmp(b, S) >= 0) {
1775 ys = (si & 0xffff) + carry;
1776 zs = (si >> 16) + (ys >> 16);
1778 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1780 Sign_Extend(borrow, y);
1781 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1783 Sign_Extend(borrow, z);
1788 y = *bx - (ys & 0xffff) + borrow;
1790 Sign_Extend(borrow, y);
1793 } while (sx <= sxe);
1797 while (--bxe > bx && !*bxe)
1805 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1807 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1808 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1811 * 1. Rather than iterating, we use a simple numeric overestimate
1812 * to determine k = floor(log10(d)). We scale relevant
1813 * quantities using O(log2(k)) rather than O(k) multiplications.
1814 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1815 * try to generate digits strictly left to right. Instead, we
1816 * compute with fewer bits and propagate the carry if necessary
1817 * when rounding the final digit up. This is often faster.
1818 * 3. Under the assumption that input will be rounded nearest,
1819 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1820 * That is, we allow equality in stopping tests when the
1821 * round-nearest rule will give the same floating-point value
1822 * as would satisfaction of the stopping test with strict
1824 * 4. We remove common factors of powers of 2 from relevant
1826 * 5. When converting floating-point integers less than 1e16,
1827 * we use floating-point arithmetic rather than resorting
1828 * to multiple-precision integers.
1829 * 6. When asked to produce fewer than 15 digits, we first try
1830 * to get by with floating-point arithmetic; we resort to
1831 * multiple-precision integer arithmetic only if we cannot
1832 * guarantee that the floating-point calculation has given
1833 * the correctly rounded result. For k requested digits and
1834 * "uniformly" distributed input, the probability is
1835 * something like 10^(k-15) that we must resort to the long
1842 (d, mode, ndigits, decpt, sign, rve, resultp)
1843 double d; int mode, ndigits, *decpt, *sign; char **rve, **resultp;
1845 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve,
1849 /* Arguments ndigits, decpt, sign are similar to those
1850 of ecvt and fcvt; trailing zeros are suppressed from
1851 the returned string. If not null, *rve is set to point
1852 to the end of the return value. If d is +-Infinity or NaN,
1853 then *decpt is set to 9999.
1856 0 ==> shortest string that yields d when read in
1857 and rounded to nearest.
1858 1 ==> like 0, but with Steele & White stopping rule;
1859 e.g. with IEEE P754 arithmetic , mode 0 gives
1860 1e23 whereas mode 1 gives 9.999999999999999e22.
1861 2 ==> max(1,ndigits) significant digits. This gives a
1862 return value similar to that of ecvt, except
1863 that trailing zeros are suppressed.
1864 3 ==> through ndigits past the decimal point. This
1865 gives a return value similar to that from fcvt,
1866 except that trailing zeros are suppressed, and
1867 ndigits can be negative.
1868 4-9 should give the same return values as 2-3, i.e.,
1869 4 <= mode <= 9 ==> same return as mode
1870 2 + (mode & 1). These modes are mainly for
1871 debugging; often they run slower but sometimes
1872 faster than modes 2-3.
1873 4,5,8,9 ==> left-to-right digit generation.
1874 6-9 ==> don't try fast floating-point estimate
1877 Values of mode other than 0-9 are treated as mode 0.
1879 Sufficient space is allocated to the return value
1880 to hold the suppressed trailing zeros.
1883 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1884 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1885 spec_case, try_quick;
1887 #ifndef Sudden_Underflow
1891 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1895 if (word0(d) & Sign_bit) {
1896 /* set sign for everything, including 0's and NaNs */
1898 word0(d) &= ~Sign_bit; /* clear sign bit */
1903 #if defined(IEEE_Arith) + defined(VAX)
1905 if ((word0(d) & Exp_mask) == Exp_mask)
1907 if (word0(d) == 0x8000)
1910 /* Infinity or NaN */
1914 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
1927 d += 0; /* normalize */
1937 b = d2b(d, &be, &bbits);
1938 #ifdef Sudden_Underflow
1939 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1941 if ( (i = (int)((word0(d) >> Exp_shift1) & (Exp_mask>>Exp_shift1))) ) {
1944 word0(d2) &= Frac_mask1;
1945 word0(d2) |= Exp_11;
1947 if ( (j = 11 - hi0bits(word0(d2) & Frac_mask)) )
1951 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1952 * log10(x) = log(x) / log(10)
1953 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1954 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1956 * This suggests computing an approximation k to log10(d) by
1958 * k = (i - Bias)*0.301029995663981
1959 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1961 * We want k to be too large rather than too small.
1962 * The error in the first-order Taylor series approximation
1963 * is in our favor, so we just round up the constant enough
1964 * to compensate for any error in the multiplication of
1965 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1966 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1967 * adding 1e-13 to the constant term more than suffices.
1968 * Hence we adjust the constant term to 0.1760912590558.
1969 * (We could get a more accurate k by invoking log10,
1970 * but this is probably not worthwhile.)
1978 #ifndef Sudden_Underflow
1981 /* d is denormalized */
1983 i = bbits + be + (Bias + (P-1) - 1);
1984 x = i > 32 ? ((word0(d) << (64 - i)) | (word1(d) >> (i - 32)))
1985 : (word1(d) << (32 - i));
1987 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
1988 i -= (Bias + (P-1) - 1) + 1;
1992 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1994 if (ds < 0. && ds != k)
1995 k--; /* want k = floor(ds) */
1997 if (k >= 0 && k <= Ten_pmax) {
2019 if (mode < 0 || mode > 9)
2040 ilim = ilim1 = i = ndigits;
2046 i = ndigits + k + 1;
2052 *resultp = (char *) malloc(i + 1);
2055 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2057 /* Try to get by with floating-point arithmetic. */
2063 ieps = 2; /* conservative */
2068 /* prevent overflows */
2070 d /= bigtens[n_bigtens-1];
2073 for (; j; j >>= 1, i++)
2079 } else if ( (j1 = -k) ) {
2080 d *= tens[j1 & 0xf];
2081 for (j = j1 >> 4; j; j >>= 1, i++)
2087 if (k_check && d < 1. && ilim > 0) {
2096 word0(eps) -= (P-1)*Exp_msk1;
2106 #ifndef No_leftright
2108 /* Use Steele & White method of only
2109 * generating digits needed.
2111 eps = 0.5/tens[ilim-1] - eps;
2115 *s++ = '0' + (int)L;
2127 /* Generate ilim digits, then fix them up. */
2128 eps *= tens[ilim-1];
2129 for (i = 1;; i++, d *= 10.) {
2132 *s++ = '0' + (int)L;
2136 else if (d < 0.5 - eps) {
2137 while (*--s == '0');
2144 #ifndef No_leftright
2154 /* Do we have a "small" integer? */
2156 if (be >= 0 && k <= Int_max) {
2159 if (ndigits < 0 && ilim <= 0) {
2161 if (ilim < 0 || d <= 5*ds)
2168 #ifdef Check_FLT_ROUNDS
2169 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2175 *s++ = '0' + (int)L;
2178 if (d > ds || (d == ds && L & 1)) {
2202 #ifndef Sudden_Underflow
2203 denorm ? be + (Bias + (P-1) - 1 + 1) :
2206 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2219 if ((i = ilim) < 0) {
2228 if (m2 > 0 && s2 > 0) {
2229 i = m2 < s2 ? m2 : s2;
2237 mhi = pow5mult(mhi, m5);
2242 if ( (j = b5 - m5) )
2245 b = pow5mult(b, b5);
2249 S = pow5mult(S, s5);
2251 /* Check for special case that d is a normalized power of 2. */
2254 if (!word1(d) && !(word0(d) & Bndry_mask)
2255 #ifndef Sudden_Underflow
2256 && word0(d) & Exp_mask
2259 /* The special case */
2267 /* Arrange for convenient computation of quotients:
2268 * shift left if necessary so divisor has 4 leading 0 bits.
2270 * Perhaps we should just compute leading 28 bits of S once
2271 * and for all and pass them and a shift to quorem, so it
2272 * can do shifts and ors to compute the numerator for q.
2275 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) )
2278 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) )
2299 b = multadd(b, 10, 0); /* we botched the k estimate */
2301 mhi = multadd(mhi, 10, 0);
2305 if (ilim <= 0 && mode > 2) {
2306 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2307 /* no digits, fcvt style */
2319 mhi = lshift(mhi, m2);
2321 /* Compute mlo -- check for special case
2322 * that d is a normalized power of 2.
2327 mhi = Balloc(mhi->k);
2329 mhi = lshift(mhi, Log2P);
2333 dig = quorem(b,S) + '0';
2334 /* Do we yet have the shortest decimal string
2335 * that will round to d?
2338 delta = diff(S, mhi);
2339 j1 = delta->sign ? 1 : cmp(b, delta);
2341 #ifndef ROUND_BIASED
2342 if (j1 == 0 && !mode && !(word1(d) & 1)) {
2351 if (j < 0 || (j == 0 && !mode
2352 #ifndef ROUND_BIASED
2359 if ((j1 > 0 || (j1 == 0 && dig & 1))
2367 if (dig == '9') { /* possible if i == 1 */
2378 b = multadd(b, 10, 0);
2380 mlo = mhi = multadd(mhi, 10, 0);
2382 mlo = multadd(mlo, 10, 0);
2383 mhi = multadd(mhi, 10, 0);
2388 *s++ = dig = quorem(b,S) + '0';
2391 b = multadd(b, 10, 0);
2394 /* Round off last digit */
2398 if (j > 0 || (j == 0 && dig & 1)) {
2408 while (*--s == '0');
2414 if (mlo && mlo != mhi)
2420 if (s == s0) { /* don't return empty string */