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 */