1 /****************************************************************
3 * The author of this software is David M. Gay.
5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
7 * Permission to use, copy, modify, and distribute this software for any
8 * purpose without fee is hereby granted, provided that this entire notice
9 * is included in all copies of any software which is or includes a copy
10 * or modification of this software and in all copies of the supporting
11 * documentation for such software.
13 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
15 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
18 ***************************************************************/
20 /* Please send bug reports to David M. Gay (dmg at acm dot org,
21 * with " at " changed at "@" and " dot " changed to "."). */
23 /* On a machine with IEEE extended-precision registers, it is
24 * necessary to specify double-precision (53-bit) rounding precision
25 * before invoking strtod or dtoa. If the machine uses (the equivalent
26 * of) Intel 80x87 arithmetic, the call
27 * _control87(PC_53, MCW_PC);
28 * does this with many compilers. Whether this or another call is
29 * appropriate depends on the compiler; for this to work, it may be
30 * necessary to #include "float.h" or another system-dependent header
34 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
36 * This strtod returns a nearest machine number to the input decimal
37 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
38 * broken by the IEEE round-even rule. Otherwise ties are broken by
39 * biased rounding (add half and chop).
41 * Inspired loosely by William D. Clinger's paper "How to Read Floating
42 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
46 * 1. We only require IEEE, IBM, or VAX double-precision
47 * arithmetic (not IEEE double-extended).
48 * 2. We get by with floating-point arithmetic in a case that
49 * Clinger missed -- when we're computing d * 10^n
50 * for a small integer d and the integer n is not too
51 * much larger than 22 (the maximum integer k for which
52 * we can represent 10^k exactly), we may be able to
53 * compute (d*10^k) * 10^(e-k) with just one roundoff.
54 * 3. Rather than a bit-at-a-time adjustment of the binary
55 * result in the hard case, we use floating-point
56 * arithmetic to determine the adjustment to within
57 * one bit; only in really hard cases do we need to
58 * compute a second residual.
59 * 4. Because of 3., we don't need a large table of powers of 10
60 * for ten-to-e (just some small tables, e.g. of 10^k
65 * #define IEEE_8087 for IEEE-arithmetic machines where the least
66 * significant byte has the lowest address.
67 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
68 * significant byte has the lowest address.
69 * #define Long int on machines with 32-bit ints and 64-bit longs.
70 * #define IBM for IBM mainframe-style floating-point arithmetic.
71 * #define VAX for VAX-style floating-point arithmetic (D_floating).
72 * #define No_leftright to omit left-right logic in fast floating-point
73 * computation of dtoa. This will cause dtoa modes 4 and 5 to be
74 * treated the same as modes 2 and 3 for some inputs.
75 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
76 * and strtod and dtoa should round accordingly. Unless Trust_FLT_ROUNDS
77 * is also #defined, fegetround() will be queried for the rounding mode.
78 * Note that both FLT_ROUNDS and fegetround() are specified by the C99
79 * standard (and are specified to be consistent, with fesetround()
80 * affecting the value of FLT_ROUNDS), but that some (Linux) systems
81 * do not work correctly in this regard, so using fegetround() is more
82 * portable than using FLT_ROUNDS directly.
83 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
84 * and Honor_FLT_ROUNDS is not #defined.
85 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
86 * that use extended-precision instructions to compute rounded
87 * products and quotients) with IBM.
88 * #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic
89 * that rounds toward +Infinity.
90 * #define ROUND_BIASED_without_Round_Up for IEEE-format with biased
91 * rounding when the underlying floating-point arithmetic uses
92 * unbiased rounding. This prevent using ordinary floating-point
93 * arithmetic when the result could be computed with one rounding error.
94 * #define Inaccurate_Divide for IEEE-format with correctly rounded
95 * products but inaccurate quotients, e.g., for Intel i860.
96 * #define NO_LONG_LONG on machines that do not have a "long long"
97 * integer type (of >= 64 bits). On such machines, you can
98 * #define Just_16 to store 16 bits per 32-bit Long when doing
99 * high-precision integer arithmetic. Whether this speeds things
100 * up or slows things down depends on the machine and the number
101 * being converted. If long long is available and the name is
102 * something other than "long long", #define Llong to be the name,
103 * and if "unsigned Llong" does not work as an unsigned version of
104 * Llong, #define #ULLong to be the corresponding unsigned type.
105 * #define KR_headers for old-style C function headers.
106 * #define Bad_float_h if your system lacks a float.h or if it does not
107 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
108 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
109 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
110 * if memory is available and otherwise does something you deem
111 * appropriate. If MALLOC is undefined, malloc will be invoked
112 * directly -- and assumed always to succeed. Similarly, if you
113 * want something other than the system's free() to be called to
114 * recycle memory acquired from MALLOC, #define FREE to be the
115 * name of the alternate routine. (FREE or free is only called in
116 * pathological cases, e.g., in a dtoa call after a dtoa return in
117 * mode 3 with thousands of digits requested.)
118 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
119 * memory allocations from a private pool of memory when possible.
120 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
121 * unless #defined to be a different length. This default length
122 * suffices to get rid of MALLOC calls except for unusual cases,
123 * such as decimal-to-binary conversion of a very long string of
124 * digits. The longest string dtoa can return is about 751 bytes
125 * long. For conversions by strtod of strings of 800 digits and
126 * all dtoa conversions in single-threaded executions with 8-byte
127 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
128 * pointers, PRIVATE_MEM >= 7112 appears adequate.
129 * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
130 * #defined automatically on IEEE systems. On such systems,
131 * when INFNAN_CHECK is #defined, strtod checks
132 * for Infinity and NaN (case insensitively). On some systems
133 * (e.g., some HP systems), it may be necessary to #define NAN_WORD0
134 * appropriately -- to the most significant word of a quiet NaN.
135 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
136 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
137 * strtod also accepts (case insensitively) strings of the form
138 * NaN(x), where x is a string of hexadecimal digits and spaces;
139 * if there is only one string of hexadecimal digits, it is taken
140 * for the 52 fraction bits of the resulting NaN; if there are two
141 * or more strings of hex digits, the first is for the high 20 bits,
142 * the second and subsequent for the low 32 bits, with intervening
143 * white space ignored; but if this results in none of the 52
144 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
145 * and NAN_WORD1 are used instead.
146 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
147 * multiple threads. In this case, you must provide (or suitably
148 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
149 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
150 * in pow5mult, ensures lazy evaluation of only one copy of high
151 * powers of 5; omitting this lock would introduce a small
152 * probability of wasting memory, but would otherwise be harmless.)
153 * You must also invoke freedtoa(s) to free the value s returned by
154 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
155 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
156 * avoids underflows on inputs whose result does not underflow.
157 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
158 * floating-point numbers and flushes underflows to zero rather
159 * than implementing gradual underflow, then you must also #define
161 * #define USE_LOCALE to use the current locale's decimal_point value.
162 * #define SET_INEXACT if IEEE arithmetic is being used and extra
163 * computation should be done to set the inexact flag when the
164 * result is inexact and avoid setting inexact when the result
165 * is exact. In this case, dtoa.c must be compiled in
166 * an environment, perhaps provided by #include "dtoa.c" in a
167 * suitable wrapper, that defines two functions,
168 * int get_inexact(void);
169 * void clear_inexact(void);
170 * such that get_inexact() returns a nonzero value if the
171 * inexact bit is already set, and clear_inexact() sets the
172 * inexact bit to 0. When SET_INEXACT is #defined, strtod
173 * also does extra computations to set the underflow and overflow
174 * flags when appropriate (i.e., when the result is tiny and
175 * inexact or when it is a numeric value rounded to +-infinity).
176 * #define NO_ERRNO if strtod should not assign errno = ERANGE when
177 * the result overflows to +-Infinity or underflows to 0.
178 * #define NO_HEX_FP to omit recognition of hexadecimal floating-point
180 * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now)
181 * to disable logic for "fast" testing of very long input strings
182 * to strtod. This testing proceeds by initially truncating the
183 * input string, then if necessary comparing the whole string with
184 * a decimal expansion to decide close cases. This logic is only
185 * used for input more than STRTOD_DIGLIM digits long (default 40).
192 typedef unsigned Long ULong
;
197 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
207 #ifdef Honor_FLT_ROUNDS
208 #ifndef Trust_FLT_ROUNDS
215 extern char *MALLOC();
217 extern void *MALLOC(size_t);
220 #define MALLOC malloc
223 #ifndef Omit_Private_Memory
225 #define PRIVATE_MEM 2304
227 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
228 static double private_mem
[PRIVATE_mem
], *pmem_next
= private_mem
;
232 #undef Avoid_Underflow
241 #ifndef NO_INFNAN_CHECK
247 #define NO_STRTOD_BIGCOMP
256 #define DBL_MAX_10_EXP 308
257 #define DBL_MAX_EXP 1024
259 #endif /*IEEE_Arith*/
263 #define DBL_MAX_10_EXP 75
264 #define DBL_MAX_EXP 63
266 #define DBL_MAX 7.2370055773322621e+75
271 #define DBL_MAX_10_EXP 38
272 #define DBL_MAX_EXP 127
274 #define DBL_MAX 1.7014118346046923e+38
278 #define LONG_MAX 2147483647
281 #else /* ifndef Bad_float_h */
283 #endif /* Bad_float_h */
295 #define CONST /* blank */
301 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
302 Exactly one of IEEE_8087
, IEEE_MC68k
, VAX
, or IBM should be defined
.
311 #define word0(x) (x)->L[1]
312 #define word1(x) (x)->L[0]
314 #define word0(x) (x)->L[0]
315 #define word1(x) (x)->L[1]
317 #define dval(x) (x)->d
319 #ifndef STRTOD_DIGLIM
320 #define STRTOD_DIGLIM 40
324 extern int strtod_diglim
;
326 #define strtod_diglim STRTOD_DIGLIM
329 /* The following definition of Storeinc is appropriate for MIPS processors.
330 * An alternative that might be better on some machines is
331 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
333 #if defined(IEEE_8087) + defined(VAX)
334 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
335 ((unsigned short *)a)[0] = (unsigned short)c, a++)
337 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
338 ((unsigned short *)a)[1] = (unsigned short)c, a++)
341 /* #define P DBL_MANT_DIG */
342 /* Ten_pmax = floor(P*log(2)/log(5)) */
343 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
344 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
345 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
349 #define Exp_shift1 20
350 #define Exp_msk1 0x100000
351 #define Exp_msk11 0x100000
352 #define Exp_mask 0x7ff00000
358 #define Exp_1 0x3ff00000
359 #define Exp_11 0x3ff00000
361 #define Frac_mask 0xfffff
362 #define Frac_mask1 0xfffff
365 #define Bndry_mask 0xfffff
366 #define Bndry_mask1 0xfffff
368 #define Sign_bit 0x80000000
374 #ifndef NO_IEEE_Scale
375 #define Avoid_Underflow
376 #ifdef Flush_Denorm /* debugging option */
377 #undef Sudden_Underflow
383 #define Flt_Rounds FLT_ROUNDS
387 #endif /*Flt_Rounds*/
389 #ifdef Honor_FLT_ROUNDS
390 #undef Check_FLT_ROUNDS
391 #define Check_FLT_ROUNDS
393 #define Rounding Flt_Rounds
396 #else /* ifndef IEEE_Arith */
397 #undef Check_FLT_ROUNDS
398 #undef Honor_FLT_ROUNDS
400 #undef Sudden_Underflow
401 #define Sudden_Underflow
406 #define Exp_shift1 24
407 #define Exp_msk1 0x1000000
408 #define Exp_msk11 0x1000000
409 #define Exp_mask 0x7f000000
415 #define Exp_1 0x41000000
416 #define Exp_11 0x41000000
417 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
418 #define Frac_mask 0xffffff
419 #define Frac_mask1 0xffffff
422 #define Bndry_mask 0xefffff
423 #define Bndry_mask1 0xffffff
425 #define Sign_bit 0x80000000
427 #define Tiny0 0x100000
436 #define Exp_msk1 0x80
437 #define Exp_msk11 0x800000
438 #define Exp_mask 0x7f80
444 #define Exp_1 0x40800000
445 #define Exp_11 0x4080
447 #define Frac_mask 0x7fffff
448 #define Frac_mask1 0xffff007f
451 #define Bndry_mask 0xffff007f
452 #define Bndry_mask1 0xffff007f
454 #define Sign_bit 0x8000
460 #endif /* IBM, VAX */
461 #endif /* IEEE_Arith */
466 #ifdef ROUND_BIASED_without_Round_Up
473 #define rounded_product(a,b) a = rnd_prod(a, b)
474 #define rounded_quotient(a,b) a = rnd_quot(a, b)
476 extern double rnd_prod(), rnd_quot();
478 extern double rnd_prod(double, double), rnd_quot(double, double);
481 #define rounded_product(a,b) a *= b
482 #define rounded_quotient(a,b) a /= b
485 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
486 #define Big1 0xffffffff
492 typedef struct BCinfo BCinfo
;
495 int dp0
, dp1
, dplen
, dsign
, e0
, inexact
, nd
, nd0
, rounding
, scale
, uflchk
;
499 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
501 #define FFFFFFFF 0xffffffffUL
508 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
509 * This makes some inner loops simpler and sometimes saves work
510 * during multiplications, but it often seems to make things slightly
511 * slower. Hence the default is now to store 32 bits per Long.
514 #else /* long long available */
516 #define Llong long long
519 #define ULLong unsigned Llong
521 #endif /* NO_LONG_LONG */
523 #ifndef MULTIPLE_THREADS
524 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
525 #define FREE_DTOA_LOCK(n) /*nothing*/
531 extern "C" double strtod(const char *s00
, char **se
);
532 extern "C" char *dtoa(double d
, int mode
, int ndigits
,
533 int *decpt
, int *sign
, char **rve
);
539 int k
, maxwds
, sign
, wds
;
543 typedef struct Bigint Bigint
;
545 static Bigint
*freelist
[Kmax
+1];
557 #ifndef Omit_Private_Memory
561 ACQUIRE_DTOA_LOCK(0);
562 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
563 /* but this case seems very unlikely. */
564 if (k
<= Kmax
&& (rv
= freelist
[k
])) {
565 freelist
[k
] = rv
->next
;
569 #ifdef Omit_Private_Memory
570 rv
= (Bigint
*)MALLOC(sizeof(Bigint
) + (x
-1)*sizeof(ULong
));
572 len
= (sizeof(Bigint
) + (x
-1)*sizeof(ULong
) + sizeof(double) - 1)
574 if (k
<= Kmax
&& pmem_next
- private_mem
+ len
<= PRIVATE_mem
) {
575 rv
= (Bigint
*)pmem_next
;
579 rv
= (Bigint
*)MALLOC(len
*sizeof(double));
586 rv
->sign
= rv
->wds
= 0;
606 ACQUIRE_DTOA_LOCK(0);
607 v
->next
= freelist
[v
->k
];
614 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
615 y->wds*sizeof(Long) + 2*sizeof(int))
620 (b
, m
, a
) Bigint
*b
; int m
, a
;
622 (Bigint
*b
, int m
, int a
) /* multiply by m and add a */
643 y
= *x
* (ULLong
)m
+ carry
;
649 y
= (xi
& 0xffff) * m
+ carry
;
650 z
= (xi
>> 16) * m
+ (y
>> 16);
652 *x
++ = (z
<< 16) + (y
& 0xffff);
662 if (wds
>= b
->maxwds
) {
677 (s
, nd0
, nd
, y9
, dplen
) CONST
char *s
; int nd0
, nd
, dplen
; ULong y9
;
679 (const char *s
, int nd0
, int nd
, ULong y9
, int dplen
)
687 for(k
= 0, y
= 1; x
> y
; y
<<= 1, k
++) ;
694 b
->x
[0] = y9
& 0xffff;
695 b
->wds
= (b
->x
[1] = y9
>> 16) ? 2 : 1;
702 b
= multadd(b
, 10, *s
++ - '0');
711 b
= multadd(b
, 10, *s
++ - '0');
726 if (!(x
& 0xffff0000)) {
730 if (!(x
& 0xff000000)) {
734 if (!(x
& 0xf0000000)) {
738 if (!(x
& 0xc0000000)) {
742 if (!(x
& 0x80000000)) {
744 if (!(x
& 0x40000000)) {
820 (a
, b
) Bigint
*a
, *b
;
822 (Bigint
*a
, Bigint
*b
)
827 ULong
*x
, *xa
, *xae
, *xb
, *xbe
, *xc
, *xc0
;
838 if (a
->wds
< b
->wds
) {
847 if (wc
> a
->maxwds
) {
851 for(x
= c
->x
, xa
= x
+ wc
; x
< xa
; x
++) {
860 for(; xb
< xbe
; xc0
++) {
866 z
= *x
++ * (ULLong
)y
+ *xc
+ carry
;
868 *xc
++ = z
& FFFFFFFF
;
876 for(; xb
< xbe
; xb
++, xc0
++) {
877 if (y
= *xb
& 0xffff) {
882 z
= (*x
& 0xffff) * y
+ (*xc
& 0xffff) + carry
;
884 z2
= (*x
++ >> 16) * y
+ (*xc
>> 16) + carry
;
897 z
= (*x
& 0xffff) * y
+ (*xc
>> 16) + carry
;
900 z2
= (*x
++ >> 16) * y
+ (*xc
& 0xffff) + carry
;
908 for(; xb
< xbe
; xc0
++) {
914 z
= *x
++ * y
+ *xc
+ carry
;
924 for(xc0
= c
->x
, xc
= xc0
+ wc
; wc
> 0 && !*--xc
; --wc
) ;
934 (b
, k
) Bigint
*b
; int k
;
939 Bigint
*b1
, *p5
, *p51
;
941 static int p05
[3] = { 5, 25, 125 };
944 b
= multadd(b
, p05
[i
-1], 0);
952 #ifdef MULTIPLE_THREADS
953 ACQUIRE_DTOA_LOCK(1);
973 if (!(p51
= p5
->next
)) {
974 #ifdef MULTIPLE_THREADS
975 ACQUIRE_DTOA_LOCK(1);
976 if (!(p51
= p5
->next
)) {
977 p51
= p5
->next
= mult(p5
,p5
);
982 p51
= p5
->next
= mult(p5
,p5
);
994 (b
, k
) Bigint
*b
; int k
;
1001 ULong
*x
, *x1
, *xe
, z
;
1009 n1
= n
+ b
->wds
+ 1;
1010 for(i
= b
->maxwds
; n1
> i
; i
<<= 1) {
1015 for(i
= 0; i
< n
; i
++) {
1025 *x1
++ = *x
<< k
| z
;
1038 *x1
++ = *x
<< k
& 0xffff | z
;
1059 (a
, b
) Bigint
*a
, *b
;
1061 (Bigint
*a
, Bigint
*b
)
1064 ULong
*xa
, *xa0
, *xb
, *xb0
;
1070 if (i
> 1 && !a
->x
[i
-1]) {
1071 Bug("cmp called with a->x[a->wds-1] == 0");
1073 if (j
> 1 && !b
->x
[j
-1]) {
1074 Bug("cmp called with b->x[b->wds-1] == 0");
1085 if (*--xa
!= *--xb
) {
1086 return *xa
< *xb
? -1 : 1;
1098 (a
, b
) Bigint
*a
, *b
;
1100 (Bigint
*a
, Bigint
*b
)
1105 ULong
*xa
, *xae
, *xb
, *xbe
, *xc
;
1143 y
= (ULLong
)*xa
++ - *xb
++ - borrow
;
1144 borrow
= y
>> 32 & (ULong
)1;
1145 *xc
++ = y
& FFFFFFFF
;
1150 borrow
= y
>> 32 & (ULong
)1;
1151 *xc
++ = y
& FFFFFFFF
;
1156 y
= (*xa
& 0xffff) - (*xb
& 0xffff) - borrow
;
1157 borrow
= (y
& 0x10000) >> 16;
1158 z
= (*xa
++ >> 16) - (*xb
++ >> 16) - borrow
;
1159 borrow
= (z
& 0x10000) >> 16;
1164 y
= (*xa
& 0xffff) - borrow
;
1165 borrow
= (y
& 0x10000) >> 16;
1166 z
= (*xa
++ >> 16) - borrow
;
1167 borrow
= (z
& 0x10000) >> 16;
1172 y
= *xa
++ - *xb
++ - borrow
;
1173 borrow
= (y
& 0x10000) >> 16;
1179 borrow
= (y
& 0x10000) >> 16;
1202 L
= (word0(x
) & Exp_mask
) - (P
-1)*Exp_msk1
;
1203 #ifndef Avoid_Underflow
1204 #ifndef Sudden_Underflow
1213 #ifndef Avoid_Underflow
1214 #ifndef Sudden_Underflow
1217 L
= -L
>> Exp_shift
;
1218 if (L
< Exp_shift
) {
1219 word0(&u
) = 0x80000 >> L
;
1225 word1(&u
) = L
>= 31 ? 1 : 1 << 31 - L
;
1236 (a
, e
) Bigint
*a
; int *e
;
1241 ULong
*xa
, *xa0
, w
, y
, z
;
1247 #define d0 word0(&d)
1248 #define d1 word1(&d)
1256 Bug("zero y in b2d");
1263 d0
= Exp_1
| y
>> (Ebits
- k
);
1264 w
= xa
> xa0
? *--xa
: 0;
1265 d1
= y
<< ((32-Ebits
) + k
) | w
>> (Ebits
- k
);
1268 z
= xa
> xa0
? *--xa
: 0;
1270 d0
= Exp_1
| y
<< k
| z
>> (32 - k
);
1271 y
= xa
> xa0
? *--xa
: 0;
1272 d1
= z
<< k
| y
>> (32 - k
);
1279 if (k
< Ebits
+ 16) {
1280 z
= xa
> xa0
? *--xa
: 0;
1281 d0
= Exp_1
| y
<< k
- Ebits
| z
>> Ebits
+ 16 - k
;
1282 w
= xa
> xa0
? *--xa
: 0;
1283 y
= xa
> xa0
? *--xa
: 0;
1284 d1
= z
<< k
+ 16 - Ebits
| w
<< k
- Ebits
| y
>> 16 + Ebits
- k
;
1287 z
= xa
> xa0
? *--xa
: 0;
1288 w
= xa
> xa0
? *--xa
: 0;
1290 d0
= Exp_1
| y
<< k
+ 16 | z
<< k
| w
>> 16 - k
;
1291 y
= xa
> xa0
? *--xa
: 0;
1292 d1
= w
<< k
+ 16 | y
<< k
;
1296 word0(&d
) = d0
>> 16 | d0
<< 16;
1297 word1(&d
) = d1
>> 16 | d1
<< 16;
1308 (d
, e
, bits
) U
*d
; int *e
, *bits
;
1310 (U
*d
, int *e
, int *bits
)
1316 #ifndef Sudden_Underflow
1321 d0
= word0(d
) >> 16 | word0(d
) << 16;
1322 d1
= word1(d
) >> 16 | word1(d
) << 16;
1336 d0
&= 0x7fffffff; /* clear sign bit, which we ignore */
1337 #ifdef Sudden_Underflow
1338 de
= (int)(d0
>> Exp_shift
);
1343 if ((de
= (int)(d0
>> Exp_shift
))) {
1349 if ((k
= lo0bits(&y
))) {
1350 x
[0] = y
| z
<< (32 - k
);
1356 #ifndef Sudden_Underflow
1359 b
->wds
= (x
[1] = z
) ? 2 : 1;
1364 #ifndef Sudden_Underflow
1372 if (k
= lo0bits(&y
))
1374 x
[0] = y
| z
<< 32 - k
& 0xffff;
1375 x
[1] = z
>> k
- 16 & 0xffff;
1381 x
[1] = y
>> 16 | z
<< 16 - k
& 0xffff;
1382 x
[2] = z
>> k
& 0xffff;
1397 Bug("Zero passed to d2b");
1417 #ifndef Sudden_Underflow
1421 *e
= (de
- Bias
- (P
-1) << 2) + k
;
1422 *bits
= 4*P
+ 8 - k
- hi0bits(word0(d
) & Frac_mask
);
1424 *e
= de
- Bias
- (P
-1) + k
;
1427 #ifndef Sudden_Underflow
1430 *e
= de
- Bias
- (P
-1) + 1 + k
;
1432 *bits
= 32*i
- hi0bits(x
[i
-1]);
1434 *bits
= (i
+2)*16 - hi0bits(x
[i
]);
1446 (a
, b
) Bigint
*a
, *b
;
1448 (Bigint
*a
, Bigint
*b
)
1454 dval(&da
) = b2d(a
, &ka
);
1455 dval(&db
) = b2d(b
, &kb
);
1457 k
= ka
- kb
+ 32*(a
->wds
- b
->wds
);
1459 k
= ka
- kb
+ 16*(a
->wds
- b
->wds
);
1463 word0(&da
) += (k
>> 2)*Exp_msk1
;
1465 dval(&da
) *= 1 << k
;
1470 word0(&db
) += (k
>> 2)*Exp_msk1
;
1472 dval(&db
) *= 1 << k
;
1477 word0(&da
) += k
*Exp_msk1
;
1481 word0(&db
) += k
*Exp_msk1
;
1484 return dval(&da
) / dval(&db
);
1489 1e0
, 1e1
, 1e2
, 1e3
, 1e4
, 1e5
, 1e6
, 1e7
, 1e8
, 1e9
,
1490 1e10
, 1e11
, 1e12
, 1e13
, 1e14
, 1e15
, 1e16
, 1e17
, 1e18
, 1e19
,
1499 bigtens
[] = { 1e16
, 1e32
, 1e64
, 1e128
, 1e256
};
1500 static CONST
double tinytens
[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1501 #ifdef Avoid_Underflow
1502 9007199254740992.*9007199254740992.e
-256
1503 /* = 2^106 * 1e-256 */
1508 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1509 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1510 #define Scale_Bit 0x10
1514 bigtens
[] = { 1e16
, 1e32
, 1e64
};
1515 static CONST
double tinytens
[] = { 1e-16, 1e-32, 1e-64 };
1518 bigtens
[] = { 1e16
, 1e32
};
1519 static CONST
double tinytens
[] = { 1e-16, 1e-32 };
1537 #ifdef Need_Hexdig /*{*/
1538 static unsigned char hexdig
[256];
1542 htinit(h
, s
, inc
) unsigned char *h
; unsigned char *s
; int inc
;
1544 htinit(unsigned char *h
, unsigned char *s
, int inc
)
1548 for(i
= 0; (j
= s
[i
]) !=0; i
++) {
1560 #define USC (unsigned char *)
1561 htinit(hexdig
, USC
"0123456789", 0x10);
1562 htinit(hexdig
, USC
"abcdef", 0x10 + 10);
1563 htinit(hexdig
, USC
"ABCDEF", 0x10 + 10);
1565 #endif /* } Need_Hexdig */
1570 #define NAN_WORD0 0x7ff80000
1580 (sp
, t
) char **sp
, *t
;
1582 (const char **sp
, const char *t
)
1586 CONST
char *s
= *sp
;
1589 if ((c
= *++s
) >= 'A' && c
<= 'Z') {
1604 (rvp
, sp
) U
*rvp
; CONST
char **sp
;
1606 (U
*rvp
, const char **sp
)
1611 int c1
, havedig
, udx0
, xshift
;
1617 havedig
= xshift
= 0;
1620 /* allow optional initial 0x or 0X */
1621 while((c
= *(CONST
unsigned char*)(s
+1)) && c
<= ' ') {
1624 if (s
[1] == '0' && (s
[2] == 'x' || s
[2] == 'X')) {
1627 while((c
= *(CONST
unsigned char*)++s
)) {
1628 if ((c1
= hexdig
[c
])) {
1631 else if (c
<= ' ') {
1632 if (udx0
&& havedig
) {
1638 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
1639 else if (/*(*/ c
== ')' && havedig
) {
1644 return; /* invalid form: don't change *sp */
1649 if (/*(*/ c
== ')') {
1653 } while((c
= *++s
));
1664 x
[0] = (x
[0] << 4) | (x
[1] >> 28);
1666 x
[1] = (x
[1] << 4) | c
;
1668 if ((x
[0] &= 0xfffff) || x
[1]) {
1669 word0(rvp
) = Exp_mask
| x
[0];
1673 #endif /*No_Hex_NaN*/
1674 #endif /* INFNAN_CHECK */
1686 #if !defined(NO_HEX_FP) || defined(Honor_FLT_ROUNDS) /*{*/
1689 increment(b
) Bigint
*b
;
1691 increment(Bigint
*b
)
1700 if (*x
< (ULong
)0xffffffffL
) {
1707 if (b
->wds
>= b
->maxwds
) {
1708 b1
= Balloc(b
->k
+1);
1720 #ifndef NO_HEX_FP /*{*/
1724 rshift(b
, k
) Bigint
*b
; int k
;
1726 rshift(Bigint
*b
, int k
)
1729 ULong
*x
, *x1
, *xe
, y
;
1741 *x1
++ = (y
| (*x
<< n
)) & 0xffffffff;
1744 if ((*x1
= y
) !=0) {
1753 if ((b
->wds
= x1
- b
->x
) == 0) {
1760 any_on(b
, k
) Bigint
*b
; int k
;
1762 any_on(Bigint
*b
, int k
)
1766 ULong
*x
, *x0
, x1
, x2
;
1774 else if (n
< nwds
&& (k
&= kmask
)) {
1791 enum { /* rounding values: same as FLT_ROUNDS */
1800 gethex(sp
, rvp
, rounding
, sign
)
1801 CONST
char **sp
; U
*rvp
; int rounding
, sign
;
1803 gethex( CONST
char **sp
, U
*rvp
, int rounding
, int sign
)
1807 CONST
unsigned char *decpt
, *s0
, *s
, *s1
;
1809 ULong L
, lostbits
, *x
;
1810 int big
, denorm
, esign
, havedig
, k
, n
, nbits
, up
, zret
;
1815 #ifdef IEEE_Arith /*{{*/
1816 emax
= 0x7fe - Bias
- P
+ 1,
1821 emax
= 0x7ff - Bias
- P
+ 1
1824 emax
= 0x7f - Bias
- P
1830 #ifdef NO_LOCALE_CACHE
1831 const unsigned char *decimalpoint
= (unsigned char*)
1832 localeconv()->decimal_point
;
1834 const unsigned char *decimalpoint
;
1835 static unsigned char *decimalpoint_cache
;
1836 if (!(s0
= decimalpoint_cache
)) {
1837 s0
= (unsigned char*)localeconv()->decimal_point
;
1838 if ((decimalpoint_cache
= (unsigned char*)
1839 MALLOC(strlen((CONST
char*)s0
) + 1))) {
1840 strcpy((char*)decimalpoint_cache
, (CONST
char*)s0
);
1841 s0
= decimalpoint_cache
;
1852 s0
= *(CONST
unsigned char **)sp
+ 2;
1853 while(s0
[havedig
] == '0') {
1867 for(i
= 0; decimalpoint
[i
]; ++i
) {
1868 if (s
[i
] != decimalpoint
[i
]) {
1895 if (*s
== *decimalpoint
&& !decpt
) {
1896 for(i
= 1; decimalpoint
[i
]; ++i
) {
1897 if (s
[i
] != decimalpoint
[i
]) {
1903 if (*s
== '.' && !decpt
) {
1911 e
= -(((Long
)(s
-decpt
)) << 2);
1926 if ((n
= hexdig
[*s
]) == 0 || n
> 0x19) {
1931 while((n
= hexdig
[*++s
]) !=0 && n
<= 0x19) {
1932 if (e1
& 0xf8000000) {
1935 e1
= 10*e1
+ n
- 0x10;
1944 *sp
= (char*)s0
- 1;
1974 #endif /* IEEE_Arith */
1996 for(k
= 0; n
> (1 << (kshift
-2)) - 1; n
>>= 1) {
2004 for(i
= 0; decimalpoint
[i
+1]; ++i
);
2008 if (*--s1
== decimalpoint
[i
]) {
2022 L
|= (hexdig
[*s1
] & 0x0f) << n
;
2026 b
->wds
= n
= x
- b
->x
;
2027 n
= ULbits
*n
- hi0bits(L
);
2036 if (x
[k
>>kshift
] & 1 << (k
& kmask
)) {
2038 if (k
> 0 && any_on(b
,k
)) {
2046 else if (n
< nbits
) {
2059 word0(rvp
) = Exp_mask
;
2068 #ifdef IEEE_Arith /*{*/
2071 if (n
== nbits
&& (n
< 2 || any_on(b
,n
-1))) {
2085 #endif /* } IEEE_Arith */
2100 lostbits
= any_on(b
,k
);
2102 if (x
[k
>>kshift
] & 1 << (k
& kmask
)) {
2116 && (lostbits
& 1) | (x
[0] & 1)) {
2132 if (nbits
== Nbits
- 1
2133 && x
[nbits
>> kshift
] & 1 << (nbits
& kmask
)) {
2134 denorm
= 0; /* not currently used */
2139 || ((n
= nbits
& kmask
) !=0
2140 && hi0bits(x
[k
-1]) < 32-n
)) {
2150 word0(rvp
) = b
->wds
> 1 ? b
->x
[1] & ~0x100000 : 0;
2153 word0(rvp
) = (b
->x
[1] & ~0x100000) | ((e
+ 0x3ff + 52) << 20);
2155 word1(rvp
) = b
->x
[0];
2159 k
= b
->x
[0] & ((1 << j
) - 1);
2175 if (k
& j
&& ((k
& (j
-1)) | lostbits
)) {
2182 word0(rvp
) = b
->x
[1] | ((e
+ 65 + 13) << 24);
2183 word1(rvp
) = b
->x
[0];
2186 /* The next two lines ignore swap of low- and high-order 2 bytes. */
2187 /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */
2188 /* word1(rvp) = b->x[0]; */
2189 word0(rvp
) = ((b
->x
[1] & ~0x800000) >> 16) | ((e
+ 129 + 55) << 7) | (b
->x
[1] << 16);
2190 word1(rvp
) = (b
->x
[0] >> 16) | (b
->x
[0] << 16);
2194 #endif /*!NO_HEX_FP}*/
2198 dshift(b
, p2
) Bigint
*b
; int p2
;
2200 dshift(Bigint
*b
, int p2
)
2203 int rv
= hi0bits(b
->x
[b
->wds
-1]) - 4;
2213 (b
, S
) Bigint
*b
, *S
;
2215 (Bigint
*b
, Bigint
*S
)
2219 ULong
*bx
, *bxe
, q
, *sx
, *sxe
;
2221 ULLong borrow
, carry
, y
, ys
;
2223 ULong borrow
, carry
, y
, ys
;
2231 /*debug*/ if (b
->wds
> n
)
2233 Bug("oversize b in quorem");
2243 q
= *bxe
/ (*sxe
+ 1); /* ensure q <= true quotient */
2245 #ifdef NO_STRTOD_BIGCOMP
2246 /*debug*/ if (q
> 9)
2248 /* An oversized q is possible when quorem is called from bigcomp and */
2249 /* the input is near, e.g., twice the smallest denormalized number. */
2250 /*debug*/ if (q
> 15)
2252 /*debug*/ Bug("oversized quotient in quorem");
2259 ys
= *sx
++ * (ULLong
)q
+ carry
;
2261 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
2262 borrow
= y
>> 32 & (ULong
)1;
2263 *bx
++ = y
& FFFFFFFF
;
2267 ys
= (si
& 0xffff) * q
+ carry
;
2268 zs
= (si
>> 16) * q
+ (ys
>> 16);
2270 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
2271 borrow
= (y
& 0x10000) >> 16;
2272 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
2273 borrow
= (z
& 0x10000) >> 16;
2276 ys
= *sx
++ * q
+ carry
;
2278 y
= *bx
- (ys
& 0xffff) - borrow
;
2279 borrow
= (y
& 0x10000) >> 16;
2287 while(--bxe
> bx
&& !*bxe
) {
2293 if (cmp(b
, S
) >= 0) {
2303 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
2304 borrow
= y
>> 32 & (ULong
)1;
2305 *bx
++ = y
& FFFFFFFF
;
2309 ys
= (si
& 0xffff) + carry
;
2310 zs
= (si
>> 16) + (ys
>> 16);
2312 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
2313 borrow
= (y
& 0x10000) >> 16;
2314 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
2315 borrow
= (z
& 0x10000) >> 16;
2320 y
= *bx
- (ys
& 0xffff) - borrow
;
2321 borrow
= (y
& 0x10000) >> 16;
2330 while(--bxe
> bx
&& !*bxe
) {
2339 #if defined(Avoid_Underflow) || !defined(NO_STRTOD_BIGCOMP) /*{*/
2343 (x
, bc
) U
*x
; BCinfo
*bc
;
2353 if (!bc
->scale
|| (i
= 2*P
+ 1 - ((word0(x
) & Exp_mask
) >> Exp_shift
)) <= 0) {
2354 return rv
; /* Is there an example where i <= 0 ? */
2356 word0(&u
) = Exp_1
+ (i
<< Exp_shift
);
2362 #ifndef NO_STRTOD_BIGCOMP
2367 U
*rv
; CONST
char *s0
; BCinfo
*bc
;
2369 (U
*rv
, const char *s0
, BCinfo
*bc
)
2373 int b2
, bbits
, d2
, dd
, dig
, dsign
, i
, j
, nd
, nd0
, p2
, p5
, speccase
;
2378 p5
= nd
+ bc
->e0
- 1;
2380 #ifndef Sudden_Underflow
2381 if (rv
->d
== 0.) { /* special case: value near underflow-to-zero */
2382 /* threshold was rounded to zero */
2386 #ifdef Avoid_Underflow
2387 word0(rv
) = (P
+2) << Exp_shift
;
2392 #ifdef Honor_FLT_ROUNDS
2393 if (bc
->rounding
== 1)
2404 b
= d2b(rv
, &p2
, &bbits
);
2405 #ifdef Avoid_Underflow
2408 /* floor(log2(rv)) == bbits - 1 + p2 */
2409 /* Check for denormal case. */
2411 if (i
> (j
= P
- Emin
- 1 + p2
)) {
2412 #ifdef Sudden_Underflow
2417 #ifdef Avoid_Underflow
2418 word0(rv
) = (1 + bc
->scale
) << Exp_shift
;
2420 word0(rv
) = Exp_msk1
;
2427 #ifdef Honor_FLT_ROUNDS
2428 if (bc
->rounding
!= 1) {
2442 #ifndef Sudden_Underflow
2447 /* Arrange for convenient computation of quotients:
2448 * shift left if necessary so divisor has 4 leading 0 bits.
2451 d
= pow5mult(d
, p5
);
2454 b
= pow5mult(b
, -p5
);
2465 if ((b2
+= i
) > 0) {
2468 if ((d2
+= i
) > 0) {
2472 /* Now b/d = exactly half-way between the two floating-point values */
2473 /* on either side of the input string. Compute first digit of b/d. */
2475 if (!(dig
= quorem(b
,d
))) {
2476 b
= multadd(b
, 10, 0); /* very unlikely */
2480 /* Compare b/d with s0 */
2482 for(i
= 0; i
< nd0
; ) {
2483 if ((dd
= s0
[i
++] - '0' - dig
)) {
2486 if (!b
->x
[0] && b
->wds
== 1) {
2492 b
= multadd(b
, 10, 0);
2495 for(j
= bc
->dp1
; i
++ < nd
;) {
2496 if ((dd
= s0
[j
++] - '0' - dig
)) {
2499 if (!b
->x
[0] && b
->wds
== 1) {
2505 b
= multadd(b
, 10, 0);
2508 if (b
->x
[0] || b
->wds
> 1 || dig
> 0) {
2514 #ifdef Honor_FLT_ROUNDS
2515 if (bc
->rounding
!= 1) {
2517 if (bc
->rounding
== 0) {
2527 if (bc
->rounding
== 0) {
2536 dval(rv
) += 2.*sulp(rv
,bc
);
2553 if (!dsign
) /* does not happen for round-near */
2555 dval(rv
) -= sulp(rv
,bc
);
2560 dval(rv
) += sulp(rv
,bc
);
2564 /* Exact half-way case: apply round-even rule. */
2565 if ((j
= ((word0(rv
) & Exp_mask
) >> Exp_shift
) - bc
->scale
) <= 0) {
2568 if (word1(rv
) & (0x1 << i
)) {
2572 else if (word0(rv
) & (0x1 << (i
-32))) {
2576 else if (word1(rv
) & 1) {
2585 #ifdef Honor_FLT_ROUNDS
2590 #endif /* NO_STRTOD_BIGCOMP */
2595 (s00
, se
) CONST
char *s00
; char **se
;
2597 (const char *s00
, char **se
)
2600 int bb2
, bb5
, bbe
, bd2
, bd5
, bbbits
, bs2
, c
, e
, e1
;
2601 int esign
, i
, j
, k
, nd
, nd0
, nf
, nz
, nz0
, nz1
, sign
;
2602 CONST
char *s
, *s0
, *s1
;
2605 U aadj2
, adj
, rv
, rv0
;
2608 Bigint
*bb
, *bb1
, *bd
, *bd0
, *bs
, *delta
;
2609 #ifdef Avoid_Underflow
2615 #ifndef NO_STRTOD_BIGCOMP
2616 int req_bigcomp
= 0;
2618 #ifdef Honor_FLT_ROUNDS /*{*/
2619 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
2620 bc
.rounding
= Flt_Rounds
;
2623 switch(fegetround()) {
2624 case FE_TOWARDZERO
: bc
.rounding
= 0; break;
2625 case FE_UPWARD
: bc
.rounding
= 2; break;
2626 case FE_DOWNWARD
: bc
.rounding
= 3;
2634 sign
= nz0
= nz1
= nz
= bc
.dplen
= bc
.uflchk
= 0;
2636 for(s
= s00
;; s
++) switch(*s
) {
2659 #ifndef NO_HEX_FP /*{*/
2663 #ifdef Honor_FLT_ROUNDS
2664 gethex(&s
, &rv
, bc
.rounding
, sign
);
2666 gethex(&s
, &rv
, 1, sign
);
2672 while(*++s
== '0') ;
2679 for(nd
= nf
= 0; (c
= *s
) >= '0' && c
<= '9'; nd
++, s
++)
2687 bc
.dp0
= bc
.dp1
= s
- s0
;
2688 for(s1
= s
; s1
> s0
&& *--s1
== '0'; ) {
2692 s1
= localeconv()->decimal_point
;
2713 bc
.dplen
= bc
.dp1
- bc
.dp0
;
2715 for(; c
== '0'; c
= *++s
) {
2718 if (c
> '0' && c
<= '9') {
2720 bc
.dp1
= bc
.dp0
+ bc
.dplen
;
2728 for(; c
>= '0' && c
<= '9'; c
= *++s
) {
2733 for(i
= 1; i
< nz
; i
++)
2737 else if (nd
<= DBL_DIG
+ 1) {
2743 else if (nd
<= DBL_DIG
+ 1) {
2752 if (c
== 'e' || c
== 'E') {
2753 if (!nd
&& !nz
&& !nz0
) {
2764 if (c
>= '0' && c
<= '9') {
2768 if (c
> '0' && c
<= '9') {
2771 while((c
= *++s
) >= '0' && c
<= '9') {
2774 if (s
- s1
> 8 || L
> 19999)
2775 /* Avoid confusion from exponents
2776 * so large that e might overflow.
2779 e
= 19999; /* safe for 16 bit ints */
2799 /* Check for Nan and Infinity */
2804 if (match(&s
,"nf")) {
2806 if (!match(&s
,"inity")) {
2809 word0(&rv
) = 0x7ff00000;
2816 if (match(&s
, "an")) {
2817 word0(&rv
) = NAN_WORD0
;
2818 word1(&rv
) = NAN_WORD1
;
2820 if (*s
== '(') { /*)*/
2827 #endif /* INFNAN_CHECK */
2834 bc
.e0
= e1
= e
-= nf
;
2836 /* Now we have nd0 digits, starting at s0, followed by a
2837 * decimal point, followed by nd-nd0 digits. The number we're
2838 * after is the integer represented by those digits times
2844 k
= nd
< DBL_DIG
+ 1 ? nd
: DBL_DIG
+ 1;
2849 oldinexact
= get_inexact();
2852 dval(&rv
) = tens
[k
- 9] * dval(&rv
) + z
;
2856 #ifndef RND_PRODQUOT
2857 #ifndef Honor_FLT_ROUNDS
2865 #ifndef ROUND_BIASED_without_Round_Up
2867 if (e
<= Ten_pmax
) {
2869 goto vax_ovfl_check
;
2871 #ifdef Honor_FLT_ROUNDS
2872 /* round correctly FLT_ROUNDS = 2 or 3 */
2878 /* rv = */ rounded_product(dval(&rv
), tens
[e
]);
2883 if (e
<= Ten_pmax
+ i
) {
2884 /* A fancier test would sometimes let us do
2885 * this for larger i values.
2887 #ifdef Honor_FLT_ROUNDS
2888 /* round correctly FLT_ROUNDS = 2 or 3 */
2895 dval(&rv
) *= tens
[i
];
2897 /* VAX exponent range is so narrow we must
2898 * worry about overflow here...
2901 word0(&rv
) -= P
*Exp_msk1
;
2902 /* rv = */ rounded_product(dval(&rv
), tens
[e
]);
2903 if ((word0(&rv
) & Exp_mask
)
2904 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
)) {
2907 word0(&rv
) += P
*Exp_msk1
;
2909 /* rv = */ rounded_product(dval(&rv
), tens
[e
]);
2914 #ifndef Inaccurate_Divide
2915 else if (e
>= -Ten_pmax
) {
2916 #ifdef Honor_FLT_ROUNDS
2917 /* round correctly FLT_ROUNDS = 2 or 3 */
2923 /* rv = */ rounded_quotient(dval(&rv
), tens
[-e
]);
2927 #endif /* ROUND_BIASED_without_Round_Up */
2935 oldinexact
= get_inexact();
2938 #ifdef Avoid_Underflow
2941 #ifdef Honor_FLT_ROUNDS
2942 if (bc
.rounding
>= 2) {
2944 bc
.rounding
= bc
.rounding
== 2 ? 0 : 2;
2946 else if (bc
.rounding
!= 2) {
2951 #endif /*IEEE_Arith*/
2953 /* Get starting approximation = rv * 10**e1 */
2956 if ((i
= e1
& 15)) {
2957 dval(&rv
) *= tens
[i
];
2960 if (e1
> DBL_MAX_10_EXP
) {
2962 /* Can't trust HUGE_VAL */
2964 #ifdef Honor_FLT_ROUNDS
2965 switch(bc
.rounding
) {
2966 case 0: /* toward 0 */
2967 case 3: /* toward -infinity */
2972 word0(&rv
) = Exp_mask
;
2975 #else /*Honor_FLT_ROUNDS*/
2976 word0(&rv
) = Exp_mask
;
2978 #endif /*Honor_FLT_ROUNDS*/
2980 /* set overflow bit */
2982 dval(&rv0
) *= dval(&rv0
);
2984 #else /*IEEE_Arith*/
2987 #endif /*IEEE_Arith*/
3002 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
3004 dval(&rv
) *= bigtens
[j
];
3006 /* The last multiplication could overflow. */
3007 word0(&rv
) -= P
*Exp_msk1
;
3008 dval(&rv
) *= bigtens
[j
];
3009 if ((z
= word0(&rv
) & Exp_mask
)
3010 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
)) {
3013 if (z
> Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
)) {
3014 /* set to largest number */
3015 /* (Can't trust DBL_MAX) */
3020 word0(&rv
) += P
*Exp_msk1
;
3026 if ((i
= e1
& 15)) {
3027 dval(&rv
) /= tens
[i
];
3030 if (e1
>= 1 << n_bigtens
) {
3033 #ifdef Avoid_Underflow
3034 if (e1
& Scale_Bit
) {
3037 for(j
= 0; e1
> 0; j
++, e1
>>= 1)
3039 dval(&rv
) *= tinytens
[j
];
3041 if (bc
.scale
&& (j
= 2*P
+ 1 - ((word0(&rv
) & Exp_mask
)
3042 >> Exp_shift
)) > 0) {
3043 /* scaled rv is denormal; clear j low bits */
3050 word0(&rv
) = (P
+2)*Exp_msk1
;
3053 word0(&rv
) &= 0xffffffff << (j
-32);
3057 word1(&rv
) &= 0xffffffff << j
;
3061 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
3063 dval(&rv
) *= tinytens
[j
];
3065 /* The last multiplication could underflow. */
3066 dval(&rv0
) = dval(&rv
);
3067 dval(&rv
) *= tinytens
[j
];
3069 dval(&rv
) = 2.*dval(&rv0
);
3070 dval(&rv
) *= tinytens
[j
];
3077 #ifndef Avoid_Underflow
3080 /* The refinement below will clean
3081 * this approximation up.
3088 /* Now the hard part -- adjusting rv to the correct value.*/
3090 /* Put digits into bd: true value = bd * 10^e */
3093 #ifndef NO_STRTOD_BIGCOMP
3094 bc
.nd0
= nd0
; /* Only needed if nd > strtod_diglim, but done here */
3095 /* to silence an erroneous warning about bc.nd0 */
3096 /* possibly not being initialized. */
3097 if (nd
> strtod_diglim
) {
3098 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */
3099 /* minimum number of decimal digits to distinguish double values */
3100 /* in IEEE arithmetic. */
3106 if (--j
< bc
.dp1
&& j
>= bc
.dp0
) {
3119 if (nd
< 9) { /* must recompute y */
3121 for(i
= 0; i
< nd0
; ++i
) {
3122 y
= 10*y
+ s0
[i
] - '0';
3124 for(j
= bc
.dp1
; i
< nd
; ++i
) {
3125 y
= 10*y
+ s0
[j
++] - '0';
3130 bd0
= s2b(s0
, nd0
, nd
, y
, bc
.dplen
);
3133 bd
= Balloc(bd0
->k
);
3135 bb
= d2b(&rv
, &bbe
, &bbbits
); /* rv = bb * 2^bbe */
3153 #ifdef Honor_FLT_ROUNDS
3154 if (bc
.rounding
!= 1) {
3158 #ifdef Avoid_Underflow
3162 i
= j
+ bbbits
- 1; /* logb(rv) */
3164 if (i
< Emin
) { /* denormal */
3171 Lsb1
= Lsb
<< (i
-32);
3177 #else /*Avoid_Underflow*/
3178 #ifdef Sudden_Underflow
3180 j
= 1 + 4*P
- 3 - bbbits
+ ((bbe
+ bbbits
- 1) & 3);
3184 #else /*Sudden_Underflow*/
3186 i
= j
+ bbbits
- 1; /* logb(rv) */
3187 if (i
< Emin
) { /* denormal */
3193 #endif /*Sudden_Underflow*/
3194 #endif /*Avoid_Underflow*/
3197 #ifdef Avoid_Underflow
3200 i
= bb2
< bd2
? bb2
: bd2
;
3210 bs
= pow5mult(bs
, bb5
);
3216 bb
= lshift(bb
, bb2
);
3219 bd
= pow5mult(bd
, bd5
);
3222 bd
= lshift(bd
, bd2
);
3225 bs
= lshift(bs
, bs2
);
3227 delta
= diff(bb
, bd
);
3228 bc
.dsign
= delta
->sign
;
3231 #ifndef NO_STRTOD_BIGCOMP /*{*/
3232 if (bc
.nd
> nd
&& i
<= 0) {
3234 /* Must use bigcomp(). */
3238 #ifdef Honor_FLT_ROUNDS
3239 if (bc
.rounding
!= 1) {
3247 i
= -1; /* Discarded digits make delta smaller. */
3250 #ifdef Honor_FLT_ROUNDS /*{*/
3251 if (bc
.rounding
!= 1) {
3253 /* Error is less than an ulp */
3254 if (!delta
->x
[0] && delta
->wds
<= 1) {
3267 else if (!bc
.dsign
) {
3270 && !(word0(&rv
) & Frac_mask
)) {
3271 y
= word0(&rv
) & Exp_mask
;
3272 #ifdef Avoid_Underflow
3273 if (!bc
.scale
|| y
> 2*P
*Exp_msk1
)
3278 delta
= lshift(delta
,Log2P
);
3279 if (cmp(delta
, bs
) <= 0) {
3285 #ifdef Avoid_Underflow /*{*/
3286 if (bc
.scale
&& (y
= word0(&rv
) & Exp_mask
)
3288 word0(&adj
) += (2*P
+1)*Exp_msk1
- y
;
3291 #ifdef Sudden_Underflow
3292 if ((word0(&rv
) & Exp_mask
) <=
3294 word0(&rv
) += P
*Exp_msk1
;
3295 dval(&rv
) += adj
.d
*ulp(dval(&rv
));
3296 word0(&rv
) -= P
*Exp_msk1
;
3299 #endif /*Sudden_Underflow*/
3300 #endif /*Avoid_Underflow}*/
3301 dval(&rv
) += adj
.d
*ulp(&rv
);
3305 adj
.d
= ratio(delta
, bs
);
3309 if (adj
.d
<= 0x7ffffffe) {
3310 /* adj = rounding ? ceil(adj) : floor(adj); */
3313 if (!((bc
.rounding
>>1) ^ bc
.dsign
)) {
3319 #ifdef Avoid_Underflow /*{*/
3320 if (bc
.scale
&& (y
= word0(&rv
) & Exp_mask
) <= 2*P
*Exp_msk1
) {
3321 word0(&adj
) += (2*P
+1)*Exp_msk1
- y
;
3324 #ifdef Sudden_Underflow
3325 if ((word0(&rv
) & Exp_mask
) <= P
*Exp_msk1
) {
3326 word0(&rv
) += P
*Exp_msk1
;
3327 adj
.d
*= ulp(dval(&rv
));
3334 word0(&rv
) -= P
*Exp_msk1
;
3337 #endif /*Sudden_Underflow*/
3338 #endif /*Avoid_Underflow}*/
3341 if (word0(&rv
) == Big0
&& word1(&rv
) == Big1
) {
3351 #endif /*}Honor_FLT_ROUNDS*/
3354 /* Error is less than half an ulp -- check for
3355 * special case of mantissa a power of two.
3357 if (bc
.dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
3358 #ifdef IEEE_Arith /*{*/
3359 #ifdef Avoid_Underflow
3360 || (word0(&rv
) & Exp_mask
) <= (2*P
+1)*Exp_msk1
3362 || (word0(&rv
) & Exp_mask
) <= Exp_msk1
3367 if (!delta
->x
[0] && delta
->wds
<= 1) {
3373 if (!delta
->x
[0] && delta
->wds
<= 1) {
3380 delta
= lshift(delta
,Log2P
);
3381 if (cmp(delta
, bs
) > 0) {
3387 /* exactly half-way between */
3389 if ((word0(&rv
) & Bndry_mask1
) == Bndry_mask1
3391 #ifdef Avoid_Underflow
3392 (bc
.scale
&& (y
= word0(&rv
) & Exp_mask
) <= 2*P
*Exp_msk1
)
3393 ? (0xffffffff & (0xffffffff << (2*P
+1-(y
>>Exp_shift
)))) :
3396 /*boundary case -- increment exponent*/
3397 if (word0(&rv
) == Big0
&& word1(&rv
) == Big1
) {
3400 word0(&rv
) = (word0(&rv
) & Exp_mask
)
3407 #ifdef Avoid_Underflow
3413 else if (!(word0(&rv
) & Bndry_mask
) && !word1(&rv
)) {
3415 /* boundary case -- decrement exponent */
3416 #ifdef Sudden_Underflow /*{{*/
3417 L
= word0(&rv
) & Exp_mask
;
3421 #ifdef Avoid_Underflow
3422 if (L
<= (bc
.scale
? (2*P
+1)*Exp_msk1
: Exp_msk1
))
3425 #endif /*Avoid_Underflow*/
3435 #else /*Sudden_Underflow}{*/
3436 #ifdef Avoid_Underflow
3438 L
= word0(&rv
) & Exp_mask
;
3439 if (L
<= (2*P
+1)*Exp_msk1
) {
3440 if (L
> (P
+2)*Exp_msk1
)
3441 /* round even ==> */
3446 /* rv = smallest denormal */
3454 #endif /*Avoid_Underflow*/
3455 L
= (word0(&rv
) & Exp_mask
) - Exp_msk1
;
3456 #endif /*Sudden_Underflow}}*/
3457 word0(&rv
) = L
| Bndry_mask1
;
3458 word1(&rv
) = 0xffffffff;
3462 #ifndef NO_STRTOD_BIGCOMP
3470 #ifndef ROUND_BIASED
3471 #ifdef Avoid_Underflow
3473 if (!(word0(&rv
) & Lsb1
)) {
3477 else if (!(word1(&rv
) & Lsb
)) {
3481 if (!(word1(&rv
) & LSB
)) {
3487 #ifdef Avoid_Underflow
3488 dval(&rv
) += sulp(&rv
, &bc
);
3490 dval(&rv
) += ulp(&rv
);
3492 #ifndef ROUND_BIASED
3494 #ifdef Avoid_Underflow
3495 dval(&rv
) -= sulp(&rv
, &bc
);
3497 dval(&rv
) -= ulp(&rv
);
3499 #ifndef Sudden_Underflow
3509 #ifdef Avoid_Underflow
3510 bc
.dsign
= 1 - bc
.dsign
;
3515 if ((aadj
= ratio(delta
, bs
)) <= 2.) {
3519 else if (word1(&rv
) || word0(&rv
) & Bndry_mask
) {
3520 #ifndef Sudden_Underflow
3521 if (word1(&rv
) == Tiny1
&& !word0(&rv
)) {
3533 /* special case -- power of FLT_RADIX to be */
3534 /* rounded down... */
3536 if (aadj
< 2./FLT_RADIX
) {
3537 aadj
= 1./FLT_RADIX
;
3547 aadj1
= bc
.dsign
? aadj
: -aadj
;
3548 #ifdef Check_FLT_ROUNDS
3549 switch(bc
.rounding
) {
3550 case 2: /* towards +infinity */
3553 case 0: /* towards 0 */
3554 case 3: /* towards -infinity */
3558 if (Flt_Rounds
== 0) {
3561 #endif /*Check_FLT_ROUNDS*/
3563 y
= word0(&rv
) & Exp_mask
;
3565 /* Check for overflow */
3567 if (y
== Exp_msk1
*(DBL_MAX_EXP
+Bias
-1)) {
3568 dval(&rv0
) = dval(&rv
);
3569 word0(&rv
) -= P
*Exp_msk1
;
3570 adj
.d
= aadj1
* ulp(&rv
);
3572 if ((word0(&rv
) & Exp_mask
) >=
3573 Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
)) {
3574 if (word0(&rv0
) == Big0
&& word1(&rv0
) == Big1
) {
3582 word0(&rv
) += P
*Exp_msk1
;
3586 #ifdef Avoid_Underflow
3587 if (bc
.scale
&& y
<= 2*P
*Exp_msk1
) {
3588 if (aadj
<= 0x7fffffff) {
3589 if ((z
= aadj
) <= 0) {
3593 aadj1
= bc
.dsign
? aadj
: -aadj
;
3595 dval(&aadj2
) = aadj1
;
3596 word0(&aadj2
) += (2*P
+1)*Exp_msk1
- y
;
3597 aadj1
= dval(&aadj2
);
3598 adj
.d
= aadj1
* ulp(&rv
);
3601 #ifdef NO_STRTOD_BIGCOMP
3613 adj
.d
= aadj1
* ulp(&rv
);
3617 #ifdef Sudden_Underflow
3618 if ((word0(&rv
) & Exp_mask
) <= P
*Exp_msk1
) {
3619 dval(&rv0
) = dval(&rv
);
3620 word0(&rv
) += P
*Exp_msk1
;
3621 adj
.d
= aadj1
* ulp(&rv
);
3624 if ((word0(&rv
) & Exp_mask
) < P
*Exp_msk1
)
3626 if ((word0(&rv
) & Exp_mask
) <= P
*Exp_msk1
)
3629 if (word0(&rv0
) == Tiny0
3630 && word1(&rv0
) == Tiny1
) {
3642 word0(&rv
) -= P
*Exp_msk1
;
3646 adj
.d
= aadj1
* ulp(&rv
);
3649 #else /*Sudden_Underflow*/
3650 /* Compute adj so that the IEEE rounding rules will
3651 * correctly round rv + adj in some half-way cases.
3652 * If rv * ulp(rv) is denormalized (i.e.,
3653 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
3654 * trouble from bits lost to denormalization;
3655 * example: 1.2e-307 .
3657 if (y
<= (P
-1)*Exp_msk1
&& aadj
> 1.) {
3658 aadj1
= (double)(int)(aadj
+ 0.5);
3663 adj
.d
= aadj1
* ulp(&rv
);
3665 #endif /*Sudden_Underflow*/
3666 #endif /*Avoid_Underflow*/
3668 z
= word0(&rv
) & Exp_mask
;
3671 #ifdef Avoid_Underflow
3675 /* Can we stop now? */
3678 /* The tolerances below are conservative. */
3679 if (bc
.dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
) {
3680 if (aadj
< .4999999 || aadj
> .5000001) {
3684 else if (aadj
< .4999999/FLT_RADIX
) {
3701 #ifndef NO_STRTOD_BIGCOMP
3705 bigcomp(&rv
, s0
, &bc
);
3706 y
= word0(&rv
) & Exp_mask
;
3707 if (y
== Exp_mask
) {
3710 if (y
== 0 && rv
.d
== 0.) {
3718 word0(&rv0
) = Exp_1
+ (70 << Exp_shift
);
3723 else if (!oldinexact
) {
3727 #ifdef Avoid_Underflow
3729 word0(&rv0
) = Exp_1
- 2*P
*Exp_msk1
;
3731 dval(&rv
) *= dval(&rv0
);
3733 /* try to avoid the bug of testing an 8087 register value */
3735 if (!(word0(&rv
) & Exp_mask
))
3737 if (word0(&rv
) == 0 && word1(&rv
) == 0)
3742 #endif /* Avoid_Underflow */
3744 if (bc
.inexact
&& !(word0(&rv
) & Exp_mask
)) {
3745 /* set underflow bit */
3746 dval(&rv0
) = 1e-300;
3747 dval(&rv0
) *= dval(&rv0
);
3754 return sign
? -dval(&rv
) : dval(&rv
);
3757 #ifndef MULTIPLE_THREADS
3758 static char *dtoa_result
;
3772 sizeof(Bigint
) - sizeof(ULong
) - sizeof(int) + j
<= i
;
3776 r
= (int*)Balloc(k
);
3779 #ifndef MULTIPLE_THREADS
3787 nrv_alloc(s
, rve
, n
) char *s
, **rve
; int n
;
3789 nrv_alloc(const char *s
, char **rve
, int n
)
3794 t
= rv
= rv_alloc(n
);
3795 while((*t
= *s
++)) {
3804 /* freedtoa(s) must be used to free values s returned by dtoa
3805 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
3806 * but for consistency with earlier versions of dtoa, it is optional
3807 * when MULTIPLE_THREADS is not defined.
3812 freedtoa(s
) char *s
;
3817 Bigint
*b
= (Bigint
*)((int *)s
- 1);
3818 b
->maxwds
= 1 << (b
->k
= *(int*)b
);
3820 #ifndef MULTIPLE_THREADS
3821 if (s
== dtoa_result
) {
3827 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
3829 * Inspired by "How to Print Floating-Point Numbers Accurately" by
3830 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
3833 * 1. Rather than iterating, we use a simple numeric overestimate
3834 * to determine k = floor(log10(d)). We scale relevant
3835 * quantities using O(log2(k)) rather than O(k) multiplications.
3836 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
3837 * try to generate digits strictly left to right. Instead, we
3838 * compute with fewer bits and propagate the carry if necessary
3839 * when rounding the final digit up. This is often faster.
3840 * 3. Under the assumption that input will be rounded nearest,
3841 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
3842 * That is, we allow equality in stopping tests when the
3843 * round-nearest rule will give the same floating-point value
3844 * as would satisfaction of the stopping test with strict
3846 * 4. We remove common factors of powers of 2 from relevant
3848 * 5. When converting floating-point integers less than 1e16,
3849 * we use floating-point arithmetic rather than resorting
3850 * to multiple-precision integers.
3851 * 6. When asked to produce fewer than 15 digits, we first try
3852 * to get by with floating-point arithmetic; we resort to
3853 * multiple-precision integer arithmetic only if we cannot
3854 * guarantee that the floating-point calculation has given
3855 * the correctly rounded result. For k requested digits and
3856 * "uniformly" distributed input, the probability is
3857 * something like 10^(k-15) that we must resort to the Long
3864 (dd
, mode
, ndigits
, decpt
, sign
, rve
)
3865 double dd
; int mode
, ndigits
, *decpt
, *sign
; char **rve
;
3867 (double dd
, int mode
, int ndigits
, int *decpt
, int *sign
, char **rve
)
3870 /* Arguments ndigits, decpt, sign are similar to those
3871 of ecvt and fcvt; trailing zeros are suppressed from
3872 the returned string. If not null, *rve is set to point
3873 to the end of the return value. If d is +-Infinity or NaN,
3874 then *decpt is set to 9999.
3877 0 ==> shortest string that yields d when read in
3878 and rounded to nearest.
3879 1 ==> like 0, but with Steele & White stopping rule;
3880 e.g. with IEEE P754 arithmetic , mode 0 gives
3881 1e23 whereas mode 1 gives 9.999999999999999e22.
3882 2 ==> max(1,ndigits) significant digits. This gives a
3883 return value similar to that of ecvt, except
3884 that trailing zeros are suppressed.
3885 3 ==> through ndigits past the decimal point. This
3886 gives a return value similar to that from fcvt,
3887 except that trailing zeros are suppressed, and
3888 ndigits can be negative.
3889 4,5 ==> similar to 2 and 3, respectively, but (in
3890 round-nearest mode) with the tests of mode 0 to
3891 possibly return a shorter string that rounds to d.
3892 With IEEE arithmetic and compilation with
3893 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
3894 as modes 2 and 3 when FLT_ROUNDS != 1.
3895 6-9 ==> Debugging modes similar to mode - 4: don't try
3896 fast floating-point estimate (if applicable).
3898 Values of mode other than 0-9 are treated as mode 0.
3900 Sufficient space is allocated to the return value
3901 to hold the suppressed trailing zeros.
3904 int bbits
, b2
, b5
, be
, dig
, i
, ieps
, ilim
, ilim0
, ilim1
,
3905 j
, j1
, k
, k0
, k_check
, leftright
, m2
, m5
, s2
, s5
,
3906 spec_case
, try_quick
;
3908 #ifndef Sudden_Underflow
3912 Bigint
*b
, *b1
, *delta
, *mlo
, *mhi
, *S
;
3916 #ifndef No_leftright
3922 int inexact
, oldinexact
;
3924 #ifdef Honor_FLT_ROUNDS /*{*/
3926 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
3927 Rounding
= Flt_Rounds
;
3930 switch(fegetround()) {
3931 case FE_TOWARDZERO
: Rounding
= 0; break;
3932 case FE_UPWARD
: Rounding
= 2; break;
3933 case FE_DOWNWARD
: Rounding
= 3;
3938 #ifndef MULTIPLE_THREADS
3940 freedtoa(dtoa_result
);
3946 if (word0(&u
) & Sign_bit
) {
3947 /* set sign for everything, including 0's and NaNs */
3949 word0(&u
) &= ~Sign_bit
; /* clear sign bit */
3955 #if defined(IEEE_Arith) + defined(VAX)
3957 if ((word0(&u
) & Exp_mask
) == Exp_mask
)
3959 if (word0(&u
) == 0x8000)
3962 /* Infinity or NaN */
3965 if (!word1(&u
) && !(word0(&u
) & 0xfffff)) {
3966 return nrv_alloc("Infinity", rve
, 8);
3969 return nrv_alloc("NaN", rve
, 3);
3973 dval(&u
) += 0; /* normalize */
3977 return nrv_alloc("0", rve
, 1);
3981 try_quick
= oldinexact
= get_inexact();
3984 #ifdef Honor_FLT_ROUNDS
3985 if (Rounding
>= 2) {
3987 Rounding
= Rounding
== 2 ? 0 : 2;
3989 else if (Rounding
!= 2) {
3995 b
= d2b(&u
, &be
, &bbits
);
3996 #ifdef Sudden_Underflow
3997 i
= (int)(word0(&u
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
));
3999 if ((i
= (int)(word0(&u
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
)))) {
4001 dval(&d2
) = dval(&u
);
4002 word0(&d2
) &= Frac_mask1
;
4003 word0(&d2
) |= Exp_11
;
4005 if (j
= 11 - hi0bits(word0(&d2
) & Frac_mask
)) {
4006 dval(&d2
) /= 1 << j
;
4010 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
4011 * log10(x) = log(x) / log(10)
4012 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
4013 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
4015 * This suggests computing an approximation k to log10(d) by
4017 * k = (i - Bias)*0.301029995663981
4018 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
4020 * We want k to be too large rather than too small.
4021 * The error in the first-order Taylor series approximation
4022 * is in our favor, so we just round up the constant enough
4023 * to compensate for any error in the multiplication of
4024 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
4025 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
4026 * adding 1e-13 to the constant term more than suffices.
4027 * Hence we adjust the constant term to 0.1760912590558.
4028 * (We could get a more accurate k by invoking log10,
4029 * but this is probably not worthwhile.)
4037 #ifndef Sudden_Underflow
4041 /* d is denormalized */
4043 i
= bbits
+ be
+ (Bias
+ (P
-1) - 1);
4044 x
= i
> 32 ? word0(&u
) << (64 - i
) | word1(&u
) >> (i
- 32)
4045 : word1(&u
) << (32 - i
);
4047 word0(&d2
) -= 31*Exp_msk1
; /* adjust exponent */
4048 i
-= (Bias
+ (P
-1) - 1) + 1;
4052 ds
= (dval(&d2
)-1.5)*0.289529654602168 + 0.1760912590558 + i
*0.301029995663981;
4054 if (ds
< 0. && ds
!= k
) {
4055 k
--; /* want k = floor(ds) */
4058 if (k
>= 0 && k
<= Ten_pmax
) {
4059 if (dval(&u
) < tens
[k
]) {
4083 if (mode
< 0 || mode
> 9) {
4088 #ifdef Check_FLT_ROUNDS
4089 try_quick
= Rounding
== 1;
4093 #endif /*SET_INEXACT*/
4100 ilim
= ilim1
= -1; /* Values for cases 0 and 1; done here to */
4101 /* silence erroneous "gcc -Wall" warning. */
4115 ilim
= ilim1
= i
= ndigits
;
4121 i
= ndigits
+ k
+ 1;
4128 s
= s0
= rv_alloc(i
);
4130 #ifdef Honor_FLT_ROUNDS
4131 if (mode
> 1 && Rounding
!= 1) {
4136 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
) {
4138 /* Try to get by with floating-point arithmetic. */
4141 dval(&d2
) = dval(&u
);
4144 ieps
= 2; /* conservative */
4149 /* prevent overflows */
4151 dval(&u
) /= bigtens
[n_bigtens
-1];
4154 for(; j
; j
>>= 1, i
++)
4161 else if ((j1
= -k
)) {
4162 dval(&u
) *= tens
[j1
& 0xf];
4163 for(j
= j1
>> 4; j
; j
>>= 1, i
++)
4166 dval(&u
) *= bigtens
[i
];
4169 if (k_check
&& dval(&u
) < 1. && ilim
> 0) {
4178 dval(&eps
) = ieps
*dval(&u
) + 7.;
4179 word0(&eps
) -= (P
-1)*Exp_msk1
;
4183 if (dval(&u
) > dval(&eps
)) {
4186 if (dval(&u
) < -dval(&eps
)) {
4191 #ifndef No_leftright
4193 /* Use Steele & White method of only
4194 * generating digits needed.
4196 dval(&eps
) = 0.5/tens
[ilim
-1] - dval(&eps
);
4198 if (k0
< 0 && j1
>= 307) {
4199 eps1
.d
= 1.01e256
; /* 1.01 allows roundoff in the next few lines */
4200 word0(&eps1
) -= Exp_msk1
* (Bias
+P
-1);
4201 dval(&eps1
) *= tens
[j1
& 0xf];
4202 for(i
= 0, j
= (j1
-256) >> 4; j
; j
>>= 1, i
++)
4204 dval(&eps1
) *= bigtens
[i
];
4206 if (eps
.d
< eps1
.d
) {
4214 *s
++ = '0' + (int)L
;
4215 if (1. - dval(&u
) < dval(&eps
)) {
4218 if (dval(&u
) < dval(&eps
)) {
4230 /* Generate ilim digits, then fix them up. */
4231 dval(&eps
) *= tens
[ilim
-1];
4232 for(i
= 1;; i
++, dval(&u
) *= 10.) {
4233 L
= (Long
)(dval(&u
));
4234 if (!(dval(&u
) -= L
)) {
4237 *s
++ = '0' + (int)L
;
4239 if (dval(&u
) > 0.5 + dval(&eps
)) {
4242 else if (dval(&u
) < 0.5 - dval(&eps
)) {
4250 #ifndef No_leftright
4255 dval(&u
) = dval(&d2
);
4260 /* Do we have a "small" integer? */
4262 if (be
>= 0 && k
<= Int_max
) {
4265 if (ndigits
< 0 && ilim
<= 0) {
4267 if (ilim
< 0 || dval(&u
) <= 5*ds
) {
4272 for(i
= 1;; i
++, dval(&u
) *= 10.) {
4273 L
= (Long
)(dval(&u
) / ds
);
4275 #ifdef Check_FLT_ROUNDS
4276 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
4282 *s
++ = '0' + (int)L
;
4290 #ifdef Honor_FLT_ROUNDS
4294 case 2: goto bump_up
;
4297 dval(&u
) += dval(&u
);
4301 if (dval(&u
) > ds
|| (dval(&u
) == ds
&& L
& 1))
4324 #ifndef Sudden_Underflow
4325 denorm
? be
+ (Bias
+ (P
-1) - 1 + 1) :
4328 1 + 4*P
- 3 - bbits
+ ((bbits
+ be
- 1) & 3);
4336 if (m2
> 0 && s2
> 0) {
4337 i
= m2
< s2
? m2
: s2
;
4345 mhi
= pow5mult(mhi
, m5
);
4350 if ((j
= b5
- m5
)) {
4355 b
= pow5mult(b
, b5
);
4360 S
= pow5mult(S
, s5
);
4363 /* Check for special case that d is a normalized power of 2. */
4366 if ((mode
< 2 || leftright
)
4367 #ifdef Honor_FLT_ROUNDS
4371 if (!word1(&u
) && !(word0(&u
) & Bndry_mask
)
4372 #ifndef Sudden_Underflow
4373 && word0(&u
) & (Exp_mask
& ~Exp_msk1
)
4376 /* The special case */
4383 /* Arrange for convenient computation of quotients:
4384 * shift left if necessary so divisor has 4 leading 0 bits.
4386 * Perhaps we should just compute leading 28 bits of S once
4387 * and for all and pass them and a shift to quorem, so it
4388 * can do shifts and ors to compute the numerator for q.
4403 b
= multadd(b
, 10, 0); /* we botched the k estimate */
4405 mhi
= multadd(mhi
, 10, 0);
4410 if (ilim
<= 0 && (mode
== 3 || mode
== 5)) {
4411 if (ilim
< 0 || cmp(b
,S
= multadd(S
,5,0)) <= 0) {
4412 /* no digits, fcvt style */
4424 mhi
= lshift(mhi
, m2
);
4427 /* Compute mlo -- check for special case
4428 * that d is a normalized power of 2.
4433 mhi
= Balloc(mhi
->k
);
4435 mhi
= lshift(mhi
, Log2P
);
4439 dig
= quorem(b
,S
) + '0';
4440 /* Do we yet have the shortest decimal string
4441 * that will round to d?
4444 delta
= diff(S
, mhi
);
4445 j1
= delta
->sign
? 1 : cmp(b
, delta
);
4447 #ifndef ROUND_BIASED
4448 if (j1
== 0 && mode
!= 1 && !(word1(&u
) & 1)
4449 #ifdef Honor_FLT_ROUNDS
4460 else if (!b
->x
[0] && b
->wds
<= 1) {
4468 if (j
< 0 || (j
== 0 && mode
!= 1
4469 #ifndef ROUND_BIASED
4473 if (!b
->x
[0] && b
->wds
<= 1) {
4479 #ifdef Honor_FLT_ROUNDS
4482 case 0: goto accept_dig
;
4483 case 2: goto keep_dig
;
4485 #endif /*Honor_FLT_ROUNDS*/
4492 if ((j1
> 0 || (j1
== 0 && dig
& 1))
4502 #ifdef Honor_FLT_ROUNDS
4507 if (dig
== '9') { /* possible if i == 1 */
4515 #ifdef Honor_FLT_ROUNDS
4522 b
= multadd(b
, 10, 0);
4524 mlo
= mhi
= multadd(mhi
, 10, 0);
4527 mlo
= multadd(mlo
, 10, 0);
4528 mhi
= multadd(mhi
, 10, 0);
4534 *s
++ = dig
= quorem(b
,S
) + '0';
4535 if (!b
->x
[0] && b
->wds
<= 1) {
4544 b
= multadd(b
, 10, 0);
4547 /* Round off last digit */
4549 #ifdef Honor_FLT_ROUNDS
4551 case 0: goto trimzeros
;
4552 case 2: goto roundoff
;
4560 if (j
> 0 || (j
== 0 && dig
& 1))
4573 #ifdef Honor_FLT_ROUNDS
4582 if (mlo
&& mlo
!= mhi
) {
4591 word0(&u
) = Exp_1
+ (70 << Exp_shift
);
4596 else if (!oldinexact
) {