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 $
34 * $DragonFly: src/lib/libc/stdlib/strtod.c,v 1.6 2006/09/10 21:22:32 swildner Exp $
36 * @(#)strtod.c 8.1 (Berkeley) 6/4/93
39 /****************************************************************
41 * The author of this software is David M. Gay.
43 * Copyright (c) 1991 by AT&T.
45 * Permission to use, copy, modify, and distribute this software for any
46 * purpose without fee is hereby granted, provided that this entire notice
47 * is included in all copies of any software which is or includes a copy
48 * or modification of this software and in all copies of the supporting
49 * documentation for such software.
51 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
52 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
53 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
54 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
56 ***************************************************************/
58 /* Please send bug reports to
60 AT&T Bell Laboratories, Room 2C-463
62 Murray Hill, NJ 07974-2070
64 dmg@research.att.com or research!dmg
67 /* strtod for IEEE--arithmetic machines.
69 * This strtod returns a nearest machine number to the input decimal
70 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
71 * broken by the IEEE round-even rule. Otherwise ties are broken by
72 * biased rounding (add half and chop).
74 * Inspired loosely by William D. Clinger's paper "How to Read Floating
75 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
79 * 1. We only require IEEE double-precision arithmetic
80 * (not IEEE double-extended).
81 * 2. We get by with floating-point arithmetic in a case that
82 * Clinger missed -- when we're computing d * 10^n
83 * for a small integer d and the integer n is not too
84 * much larger than 22 (the maximum integer k for which
85 * we can represent 10^k exactly), we may be able to
86 * compute (d*10^k) * 10^(e-k) with just one roundoff.
87 * 3. Rather than a bit-at-a-time adjustment of the binary
88 * result in the hard case, we use floating-point
89 * arithmetic to determine the adjustment to within
90 * one bit; only in really hard cases do we need to
91 * compute a second residual.
92 * 4. Because of 3., we don't need a large table of powers of 10
93 * for ten-to-e (just some small tables, e.g. of 10^k
98 * #define IEEE_8087 for IEEE-arithmetic machines where the least
99 * significant byte has the lowest address.
100 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
101 * significant byte has the lowest address.
102 * #define Sudden_Underflow for IEEE-format machines without gradual
103 * underflow (i.e., that flush to zero on underflow).
104 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
105 * #define No_leftright to omit left-right logic in fast floating-point
106 * computation of dtoa.
107 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
108 * #define ROUND_BIASED for IEEE-format with biased rounding.
109 * #define Inaccurate_Divide for IEEE-format with correctly rounded
110 * products but inaccurate quotients, e.g., for Intel i860.
111 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
112 * integer arithmetic. Whether this speeds things up or slows things
113 * down depends on the machine and the number being converted.
116 #if defined(i386) || defined(mips) && defined(MIPSEL)
124 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
152 #ifdef Unsigned_Shifts
153 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
155 #define Sign_Extend(a,b) /*no-op*/
158 #if defined(IEEE_8087) + defined(IEEE_MC68k) != 1
159 Exactly one of IEEE_8087 or IEEE_MC68k should be defined.
162 union doubleasulongs {
167 #define word0(x) (((union doubleasulongs *)&x)->w)[1]
168 #define word1(x) (((union doubleasulongs *)&x)->w)[0]
170 #define word0(x) (((union doubleasulongs *)&x)->w)[0]
171 #define word1(x) (((union doubleasulongs *)&x)->w)[1]
174 /* The following definition of Storeinc is appropriate for MIPS processors.
175 * An alternative that might be better on some machines is
176 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
178 #if defined(IEEE_8087)
179 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
180 ((unsigned short *)a)[0] = (unsigned short)c, a++)
182 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
183 ((unsigned short *)a)[1] = (unsigned short)c, a++)
186 /* #define P DBL_MANT_DIG */
187 /* Ten_pmax = floor(P*log(2)/log(5)) */
188 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
189 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
190 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
192 #if defined(IEEE_8087) + defined(IEEE_MC68k)
194 #define Exp_shift1 20
195 #define Exp_msk1 0x100000
196 #define Exp_msk11 0x100000
197 #define Exp_mask 0x7ff00000
202 #define Exp_1 0x3ff00000
203 #define Exp_11 0x3ff00000
205 #define Frac_mask 0xfffff
206 #define Frac_mask1 0xfffff
209 #define Bndry_mask 0xfffff
210 #define Bndry_mask1 0xfffff
212 #define Sign_bit 0x80000000
218 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
225 #define rounded_product(a,b) a *= b
226 #define rounded_quotient(a,b) a /= b
228 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
229 #define Big1 0xffffffff
232 /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
233 * This makes some inner loops simpler and sometimes saves work
234 * during multiplications, but it often seems to make things slightly
235 * slower. Hence the default is now to store 32 bits per long.
245 extern "C" double strtod(const char *s00, char **se);
246 extern "C" char *__dtoa(double d, int mode, int ndigits,
247 int *decpt, int *sign, char **rve, char **resultp);
253 int k, maxwds, sign, wds;
257 typedef struct Bigint Bigint;
266 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(long));
269 rv->sign = rv->wds = 0;
279 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
280 y->wds*sizeof(long) + 2*sizeof(int))
283 multadd(Bigint *b, int m, int a) /* multiply by m and add a */
298 y = (xi & 0xffff) * m + a;
299 z = (xi >> 16) * m + (y >> 16);
301 *x++ = (z << 16) + (y & 0xffff);
309 if (wds >= b->maxwds) {
322 s2b(CONST char *s, int nd0, int nd, unsigned long y9)
329 for (k = 0, y = 1; x > y; y <<= 1, k++) ;
336 b->x[0] = y9 & 0xffff;
337 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
344 b = multadd(b, 10, *s++ - '0');
350 b = multadd(b, 10, *s++ - '0');
355 hi0bits(unsigned long x)
359 if (!(x & 0xffff0000)) {
363 if (!(x & 0xff000000)) {
367 if (!(x & 0xf0000000)) {
371 if (!(x & 0xc0000000)) {
375 if (!(x & 0x80000000)) {
377 if (!(x & 0x40000000))
384 lo0bits(unsigned long *y)
387 unsigned long x = *y;
438 mult(Bigint *a, Bigint *b)
442 unsigned long carry, y, z;
443 unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
448 if (a->wds < b->wds) {
460 for (x = c->x, xa = x + wc; x < xa; x++)
468 for (; xb < xbe; xb++, xc0++) {
469 if ( (y = *xb & 0xffff) ) {
474 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
476 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
482 if ( (y = *xb >> 16) ) {
488 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
491 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
498 for (; xb < xbe; xc0++) {
504 z = *x++ * y + *xc + carry;
512 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
520 pow5mult(Bigint *b, int k)
522 Bigint *b1, *p5, *p51;
524 static int p05[3] = { 5, 25, 125 };
527 b = multadd(b, p05[i-1], 0);
544 if (!(p51 = p5->next)) {
545 p51 = p5->next = mult(p5,p5);
554 lshift(Bigint *b, int k)
558 unsigned long *x, *x1, *xe, z;
567 for (i = b->maxwds; n1 > i; i <<= 1)
571 for (i = 0; i < n; i++)
591 *x1++ = *x << k & 0xffff | z;
608 cmp(Bigint *a, Bigint *b)
610 unsigned long *xa, *xa0, *xb, *xb0;
616 if (i > 1 && !a->x[i-1])
617 Bug("cmp called with a->x[a->wds-1] == 0");
618 if (j > 1 && !b->x[j-1])
619 Bug("cmp called with b->x[b->wds-1] == 0");
629 return *xa < *xb ? -1 : 1;
637 diff(Bigint *a, Bigint *b)
641 long borrow, y; /* We need signed shifts here. */
642 unsigned long *xa, *xae, *xb, *xbe, *xc;
673 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
675 Sign_Extend(borrow, y);
676 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
678 Sign_Extend(borrow, z);
682 y = (*xa & 0xffff) + borrow;
684 Sign_Extend(borrow, y);
685 z = (*xa++ >> 16) + borrow;
687 Sign_Extend(borrow, z);
692 y = *xa++ - *xb++ + borrow;
694 Sign_Extend(borrow, y);
700 Sign_Extend(borrow, y);
716 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
717 #ifndef Sudden_Underflow
722 #ifndef Sudden_Underflow
726 word0(a) = 0x80000 >> L;
731 word1(a) = L >= 31 ? 1 : 1 << (31 - L);
739 b2d(Bigint *a, int *e)
741 unsigned long *xa, *xa0, w, y, z;
751 if (!y) Bug("zero y in b2d");
757 d0 = Exp_1 | (y >> (Ebits - k));
758 w = xa > xa0 ? *--xa : 0;
759 d1 = (y << ((32-Ebits) + k)) | (w >> (Ebits - k));
762 z = xa > xa0 ? *--xa : 0;
764 d0 = Exp_1 | (y << k) | (z >> (32 - k));
765 y = xa > xa0 ? *--xa : 0;
766 d1 = (z << k) | (y >> (32 - k));
772 if (k < Ebits + 16) {
773 z = xa > xa0 ? *--xa : 0;
774 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
775 w = xa > xa0 ? *--xa : 0;
776 y = xa > xa0 ? *--xa : 0;
777 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
780 z = xa > xa0 ? *--xa : 0;
781 w = xa > xa0 ? *--xa : 0;
783 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
784 y = xa > xa0 ? *--xa : 0;
785 d1 = w << k + 16 | y << k;
794 d2b(double d, int *e, int *bits)
798 unsigned long *x, y, z;
810 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
811 #ifdef Sudden_Underflow
812 de = (int)(d0 >> Exp_shift);
815 if ( (de = (int)(d0 >> Exp_shift)) )
820 if ( (k = lo0bits(&y)) ) {
821 x[0] = y | (z << (32 - k));
826 i = b->wds = (x[1] = z) ? 2 : 1;
830 Bug("Zero passed to d2b");
841 x[0] = y | z << 32 - k & 0xffff;
842 x[1] = z >> k - 16 & 0xffff;
847 x[1] = y >> 16 | z << 16 - k & 0xffff;
848 x[2] = z >> k & 0xffff;
862 Bug("Zero passed to d2b");
879 #ifndef Sudden_Underflow
882 *e = de - Bias - (P-1) + k;
884 #ifndef Sudden_Underflow
886 *e = de - Bias - (P-1) + 1 + k;
888 *bits = 32*i - hi0bits(x[i-1]);
890 *bits = (i+2)*16 - hi0bits(x[i]);
900 ratio(Bigint *a, Bigint *b)
908 k = ka - kb + 32*(a->wds - b->wds);
910 k = ka - kb + 16*(a->wds - b->wds);
913 word0(da) += k*Exp_msk1;
916 word0(db) += k*Exp_msk1;
923 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
924 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
930 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
931 static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
934 bigtens[] = { 1e16, 1e32 };
935 static double tinytens[] = { 1e-16, 1e-32 };
940 strtod(CONST char *s00, char **se)
942 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
943 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
944 CONST char *s, *s0, *s1;
945 double aadj, aadj1, adj, rv, rv0;
948 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
949 char decimal_point = localeconv()->decimal_point[0];
953 for (s = s00;;s++) switch(*s) {
965 if (isspace((unsigned char)*s))
972 while (*++s == '0') ;
978 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
984 if ((char)c == decimal_point) {
987 for (; c == '0'; c = *++s)
989 if (c > '0' && c <= '9') {
997 for (; c >= '0' && c <= '9'; c = *++s) {
1002 for (i = 1; i < nz; i++)
1005 else if (nd <= DBL_DIG + 1)
1009 else if (nd <= DBL_DIG + 1)
1017 if (c == 'e' || c == 'E') {
1018 if (!nd && !nz && !nz0) {
1030 if (c >= '0' && c <= '9') {
1033 if (c > '0' && c <= '9') {
1036 while ((c = *++s) >= '0' && c <= '9')
1038 if (s - s1 > 8 || L > 19999)
1039 /* Avoid confusion from exponents
1040 * so large that e might overflow.
1042 e = 19999; /* safe for 16 bit ints */
1059 /* Now we have nd0 digits, starting at s0, followed by a
1060 * decimal point, followed by nd-nd0 digits. The number we're
1061 * after is the integer represented by those digits times
1066 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1069 rv = tens[k - 9] * rv + z;
1070 if (nd <= DBL_DIG) {
1074 if (e <= Ten_pmax) {
1075 /* rv = */ rounded_product(rv, tens[e]);
1079 if (e <= Ten_pmax + i) {
1080 /* A fancier test would sometimes let us do
1081 * this for larger i values.
1085 /* rv = */ rounded_product(rv, tens[e]);
1089 #ifndef Inaccurate_Divide
1090 else if (e >= -Ten_pmax) {
1091 /* rv = */ rounded_quotient(rv, tens[-e]);
1098 /* Get starting approximation = rv * 10**e1 */
1101 if ( (i = e1 & 15) )
1103 if ( (e1 &= ~15) ) {
1104 if (e1 > DBL_MAX_10_EXP) {
1111 for (j = 0; e1 > 1; j++, e1 >>= 1)
1114 /* The last multiplication could overflow. */
1115 word0(rv) -= P*Exp_msk1;
1117 if ((z = word0(rv) & Exp_mask)
1118 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1120 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1121 /* set to largest number */
1122 /* (Can't trust DBL_MAX) */
1127 word0(rv) += P*Exp_msk1;
1130 } else if (e1 < 0) {
1132 if ( (i = e1 & 15) )
1134 if ( (e1 &= ~15) ) {
1136 for (j = 0; e1 > 1; j++, e1 >>= 1)
1139 /* The last multiplication could underflow. */
1153 /* The refinement below will clean
1154 * this approximation up.
1160 /* Now the hard part -- adjusting rv to the correct value.*/
1162 /* Put digits into bd: true value = bd * 10^e */
1164 bd0 = s2b(s0, nd0, nd, y);
1167 bd = Balloc(bd0->k);
1169 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1184 #ifdef Sudden_Underflow
1187 i = bbe + bbbits - 1; /* logb(rv) */
1188 if (i < Emin) /* denormal */
1195 i = bb2 < bd2 ? bb2 : bd2;
1204 bs = pow5mult(bs, bb5);
1210 bb = lshift(bb, bb2);
1212 bd = pow5mult(bd, bd5);
1214 bd = lshift(bd, bd2);
1216 bs = lshift(bs, bs2);
1217 delta = diff(bb, bd);
1218 dsign = delta->sign;
1222 /* Error is less than half an ulp -- check for
1223 * special case of mantissa a power of two.
1225 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1227 delta = lshift(delta,Log2P);
1228 if (cmp(delta, bs) > 0)
1233 /* exactly half-way between */
1235 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1236 && word1(rv) == 0xffffffff) {
1237 /*boundary case -- increment exponent*/
1238 word0(rv) = (word0(rv) & Exp_mask)
1243 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1245 /* boundary case -- decrement exponent */
1246 #ifdef Sudden_Underflow
1247 L = word0(rv) & Exp_mask;
1252 L = (word0(rv) & Exp_mask) - Exp_msk1;
1254 word0(rv) = L | Bndry_mask1;
1255 word1(rv) = 0xffffffff;
1258 #ifndef ROUND_BIASED
1259 if (!(word1(rv) & LSB))
1264 #ifndef ROUND_BIASED
1267 #ifndef Sudden_Underflow
1275 if ((aadj = ratio(delta, bs)) <= 2.) {
1278 else if (word1(rv) || word0(rv) & Bndry_mask) {
1279 #ifndef Sudden_Underflow
1280 if (word1(rv) == Tiny1 && !word0(rv))
1286 /* special case -- power of FLT_RADIX to be */
1287 /* rounded down... */
1289 if (aadj < 2./FLT_RADIX)
1290 aadj = 1./FLT_RADIX;
1297 aadj1 = dsign ? aadj : -aadj;
1298 #ifdef Check_FLT_ROUNDS
1299 switch(FLT_ROUNDS) {
1300 case 2: /* towards +infinity */
1303 case 0: /* towards 0 */
1304 case 3: /* towards -infinity */
1308 if (FLT_ROUNDS == 0)
1312 y = word0(rv) & Exp_mask;
1314 /* Check for overflow */
1316 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1318 word0(rv) -= P*Exp_msk1;
1319 adj = aadj1 * ulp(rv);
1321 if ((word0(rv) & Exp_mask) >=
1322 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1323 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1329 word0(rv) += P*Exp_msk1;
1331 #ifdef Sudden_Underflow
1332 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1334 word0(rv) += P*Exp_msk1;
1335 adj = aadj1 * ulp(rv);
1337 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1339 if (word0(rv0) == Tiny0
1340 && word1(rv0) == Tiny1)
1346 word0(rv) -= P*Exp_msk1;
1348 adj = aadj1 * ulp(rv);
1352 /* Compute adj so that the IEEE rounding rules will
1353 * correctly round rv + adj in some half-way cases.
1354 * If rv * ulp(rv) is denormalized (i.e.,
1355 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1356 * trouble from bits lost to denormalization;
1357 * example: 1.2e-307 .
1359 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1360 aadj1 = (double)(int)(aadj + 0.5);
1364 adj = aadj1 * ulp(rv);
1368 z = word0(rv) & Exp_mask;
1370 /* Can we stop now? */
1373 /* The tolerances below are conservative. */
1374 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1375 if (aadj < .4999999 || aadj > .5000001)
1377 } else if (aadj < .4999999/FLT_RADIX)
1394 return sign ? -rv : rv;
1398 quorem(Bigint *b, Bigint *S)
1402 unsigned long carry, q, ys;
1403 unsigned long *bx, *bxe, *sx, *sxe;
1406 unsigned long si, zs;
1411 /*debug*/ if (b->wds > n)
1412 /*debug*/ Bug("oversize b in quorem");
1420 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1422 /*debug*/ if (q > 9)
1423 /*debug*/ Bug("oversized quotient in quorem");
1431 ys = (si & 0xffff) * q + carry;
1432 zs = (si >> 16) * q + (ys >> 16);
1434 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1436 Sign_Extend(borrow, y);
1437 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1439 Sign_Extend(borrow, z);
1442 ys = *sx++ * q + carry;
1444 y = *bx - (ys & 0xffff) + borrow;
1446 Sign_Extend(borrow, y);
1449 } while (sx <= sxe);
1452 while (--bxe > bx && !*bxe)
1457 if (cmp(b, S) >= 0) {
1466 ys = (si & 0xffff) + carry;
1467 zs = (si >> 16) + (ys >> 16);
1469 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1471 Sign_Extend(borrow, y);
1472 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1474 Sign_Extend(borrow, z);
1479 y = *bx - (ys & 0xffff) + borrow;
1481 Sign_Extend(borrow, y);
1484 } while (sx <= sxe);
1488 while (--bxe > bx && !*bxe)
1496 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1498 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1499 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1502 * 1. Rather than iterating, we use a simple numeric overestimate
1503 * to determine k = floor(log10(d)). We scale relevant
1504 * quantities using O(log2(k)) rather than O(k) multiplications.
1505 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1506 * try to generate digits strictly left to right. Instead, we
1507 * compute with fewer bits and propagate the carry if necessary
1508 * when rounding the final digit up. This is often faster.
1509 * 3. Under the assumption that input will be rounded nearest,
1510 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1511 * That is, we allow equality in stopping tests when the
1512 * round-nearest rule will give the same floating-point value
1513 * as would satisfaction of the stopping test with strict
1515 * 4. We remove common factors of powers of 2 from relevant
1517 * 5. When converting floating-point integers less than 1e16,
1518 * we use floating-point arithmetic rather than resorting
1519 * to multiple-precision integers.
1520 * 6. When asked to produce fewer than 15 digits, we first try
1521 * to get by with floating-point arithmetic; we resort to
1522 * multiple-precision integer arithmetic only if we cannot
1523 * guarantee that the floating-point calculation has given
1524 * the correctly rounded result. For k requested digits and
1525 * "uniformly" distributed input, the probability is
1526 * something like 10^(k-15) that we must resort to the long
1531 __dtoa(double d, int mode, int ndigits, int *decpt, int *sign, char **rve,
1534 /* Arguments ndigits, decpt, sign are similar to those
1535 of ecvt and fcvt; trailing zeros are suppressed from
1536 the returned string. If not null, *rve is set to point
1537 to the end of the return value. If d is +-Infinity or NaN,
1538 then *decpt is set to 9999.
1541 0 ==> shortest string that yields d when read in
1542 and rounded to nearest.
1543 1 ==> like 0, but with Steele & White stopping rule;
1544 e.g. with IEEE P754 arithmetic , mode 0 gives
1545 1e23 whereas mode 1 gives 9.999999999999999e22.
1546 2 ==> max(1,ndigits) significant digits. This gives a
1547 return value similar to that of ecvt, except
1548 that trailing zeros are suppressed.
1549 3 ==> through ndigits past the decimal point. This
1550 gives a return value similar to that from fcvt,
1551 except that trailing zeros are suppressed, and
1552 ndigits can be negative.
1553 4-9 should give the same return values as 2-3, i.e.,
1554 4 <= mode <= 9 ==> same return as mode
1555 2 + (mode & 1). These modes are mainly for
1556 debugging; often they run slower but sometimes
1557 faster than modes 2-3.
1558 4,5,8,9 ==> left-to-right digit generation.
1559 6-9 ==> don't try fast floating-point estimate
1562 Values of mode other than 0-9 are treated as mode 0.
1564 Sufficient space is allocated to the return value
1565 to hold the suppressed trailing zeros.
1568 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1569 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1570 spec_case, try_quick;
1572 #ifndef Sudden_Underflow
1576 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1581 * XXX initialize to silence GCC warnings
1588 if (word0(d) & Sign_bit) {
1589 /* set sign for everything, including 0's and NaNs */
1591 word0(d) &= ~Sign_bit; /* clear sign bit */
1597 if ((word0(d) & Exp_mask) == Exp_mask)
1599 /* Infinity or NaN */
1601 s = !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" : "NaN";
1603 *rve = s[3] ? s + 8 : s + 3;
1615 b = d2b(d, &be, &bbits);
1616 #ifdef Sudden_Underflow
1617 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1619 if ( (i = (int)((word0(d) >> Exp_shift1) & (Exp_mask>>Exp_shift1))) ) {
1622 word0(d2) &= Frac_mask1;
1623 word0(d2) |= Exp_11;
1625 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1626 * log10(x) = log(x) / log(10)
1627 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1628 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1630 * This suggests computing an approximation k to log10(d) by
1632 * k = (i - Bias)*0.301029995663981
1633 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1635 * We want k to be too large rather than too small.
1636 * The error in the first-order Taylor series approximation
1637 * is in our favor, so we just round up the constant enough
1638 * to compensate for any error in the multiplication of
1639 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1640 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1641 * adding 1e-13 to the constant term more than suffices.
1642 * Hence we adjust the constant term to 0.1760912590558.
1643 * (We could get a more accurate k by invoking log10,
1644 * but this is probably not worthwhile.)
1648 #ifndef Sudden_Underflow
1651 /* d is denormalized */
1653 i = bbits + be + (Bias + (P-1) - 1);
1654 x = i > 32 ? ((word0(d) << (64 - i)) | (word1(d) >> (i - 32)))
1655 : (word1(d) << (32 - i));
1657 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
1658 i -= (Bias + (P-1) - 1) + 1;
1662 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1664 if (ds < 0. && ds != k)
1665 k--; /* want k = floor(ds) */
1667 if (k >= 0 && k <= Ten_pmax) {
1689 if (mode < 0 || mode > 9)
1710 ilim = ilim1 = i = ndigits;
1716 i = ndigits + k + 1;
1722 *resultp = (char *) malloc(i + 1);
1725 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
1727 /* Try to get by with floating-point arithmetic. */
1733 ieps = 2; /* conservative */
1738 /* prevent overflows */
1740 d /= bigtens[n_bigtens-1];
1743 for (; j; j >>= 1, i++)
1749 } else if ( (j1 = -k) ) {
1750 d *= tens[j1 & 0xf];
1751 for (j = j1 >> 4; j; j >>= 1, i++)
1757 if (k_check && d < 1. && ilim > 0) {
1766 word0(eps) -= (P-1)*Exp_msk1;
1776 #ifndef No_leftright
1778 /* Use Steele & White method of only
1779 * generating digits needed.
1781 eps = 0.5/tens[ilim-1] - eps;
1785 *s++ = '0' + (int)L;
1797 /* Generate ilim digits, then fix them up. */
1798 eps *= tens[ilim-1];
1799 for (i = 1;; i++, d *= 10.) {
1802 *s++ = '0' + (int)L;
1806 else if (d < 0.5 - eps) {
1807 while (*--s == '0');
1814 #ifndef No_leftright
1824 /* Do we have a "small" integer? */
1826 if (be >= 0 && k <= Int_max) {
1829 if (ndigits < 0 && ilim <= 0) {
1831 if (ilim < 0 || d <= 5*ds)
1838 #ifdef Check_FLT_ROUNDS
1839 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
1845 *s++ = '0' + (int)L;
1848 if (d > ds || (d == ds && L & 1)) {
1872 #ifndef Sudden_Underflow
1873 denorm ? be + (Bias + (P-1) - 1 + 1) :
1885 if ((i = ilim) < 0) {
1894 if (m2 > 0 && s2 > 0) {
1895 i = m2 < s2 ? m2 : s2;
1903 mhi = pow5mult(mhi, m5);
1908 if ( (j = b5 - m5) )
1911 b = pow5mult(b, b5);
1915 S = pow5mult(S, s5);
1917 /* Check for special case that d is a normalized power of 2. */
1920 if (!word1(d) && !(word0(d) & Bndry_mask)
1921 #ifndef Sudden_Underflow
1922 && word0(d) & Exp_mask
1925 /* The special case */
1933 /* Arrange for convenient computation of quotients:
1934 * shift left if necessary so divisor has 4 leading 0 bits.
1936 * Perhaps we should just compute leading 28 bits of S once
1937 * and for all and pass them and a shift to quorem, so it
1938 * can do shifts and ors to compute the numerator for q.
1941 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) )
1944 if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) )
1965 b = multadd(b, 10, 0); /* we botched the k estimate */
1967 mhi = multadd(mhi, 10, 0);
1971 if (ilim <= 0 && mode > 2) {
1972 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
1973 /* no digits, fcvt style */
1985 mhi = lshift(mhi, m2);
1987 /* Compute mlo -- check for special case
1988 * that d is a normalized power of 2.
1993 mhi = Balloc(mhi->k);
1995 mhi = lshift(mhi, Log2P);
1999 dig = quorem(b,S) + '0';
2000 /* Do we yet have the shortest decimal string
2001 * that will round to d?
2004 delta = diff(S, mhi);
2005 j1 = delta->sign ? 1 : cmp(b, delta);
2007 #ifndef ROUND_BIASED
2008 if (j1 == 0 && !mode && !(word1(d) & 1)) {
2017 if (j < 0 || (j == 0 && !mode
2018 #ifndef ROUND_BIASED
2025 if ((j1 > 0 || (j1 == 0 && dig & 1))
2033 if (dig == '9') { /* possible if i == 1 */
2044 b = multadd(b, 10, 0);
2046 mlo = mhi = multadd(mhi, 10, 0);
2048 mlo = multadd(mlo, 10, 0);
2049 mhi = multadd(mhi, 10, 0);
2054 *s++ = dig = quorem(b,S) + '0';
2057 b = multadd(b, 10, 0);
2060 /* Round off last digit */
2064 if (j > 0 || (j == 0 && dig & 1)) {
2074 while (*--s == '0');
2080 if (mlo && mlo != mhi)
2086 if (s == s0) { /* don't return empty string */