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.
74 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
75 * and strtod and dtoa should round accordingly. Unless Trust_FLT_ROUNDS
76 * is also #defined, fegetround() will be queried for the rounding mode.
77 * Note that both FLT_ROUNDS and fegetround() are specified by the C99
78 * standard (and are specified to be consistent, with fesetround()
79 * affecting the value of FLT_ROUNDS), but that some (Linux) systems
80 * do not work correctly in this regard, so using fegetround() is more
81 * portable than using FLT_FOUNDS directly.
82 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
83 * and Honor_FLT_ROUNDS is not #defined.
84 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
85 * that use extended-precision instructions to compute rounded
86 * products and quotients) with IBM.
87 * #define ROUND_BIASED for IEEE-format with biased rounding.
88 * #define Inaccurate_Divide for IEEE-format with correctly rounded
89 * products but inaccurate quotients, e.g., for Intel i860.
90 * #define NO_LONG_LONG on machines that do not have a "long long"
91 * integer type (of >= 64 bits). On such machines, you can
92 * #define Just_16 to store 16 bits per 32-bit Long when doing
93 * high-precision integer arithmetic. Whether this speeds things
94 * up or slows things down depends on the machine and the number
95 * being converted. If long long is available and the name is
96 * something other than "long long", #define Llong to be the name,
97 * and if "unsigned Llong" does not work as an unsigned version of
98 * Llong, #define #ULLong to be the corresponding unsigned type.
99 * #define KR_headers for old-style C function headers.
100 * #define Bad_float_h if your system lacks a float.h or if it does not
101 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
102 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
103 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
104 * if memory is available and otherwise does something you deem
105 * appropriate. If MALLOC is undefined, malloc will be invoked
106 * directly -- and assumed always to succeed. Similarly, if you
107 * want something other than the system's free() to be called to
108 * recycle memory acquired from MALLOC, #define FREE to be the
109 * name of the alternate routine. (FREE or free is only called in
110 * pathological cases, e.g., in a dtoa call after a dtoa return in
111 * mode 3 with thousands of digits requested.)
112 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
113 * memory allocations from a private pool of memory when possible.
114 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
115 * unless #defined to be a different length. This default length
116 * suffices to get rid of MALLOC calls except for unusual cases,
117 * such as decimal-to-binary conversion of a very long string of
118 * digits. The longest string dtoa can return is about 751 bytes
119 * long. For conversions by strtod of strings of 800 digits and
120 * all dtoa conversions in single-threaded executions with 8-byte
121 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
122 * pointers, PRIVATE_MEM >= 7112 appears adequate.
123 * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
124 * #defined automatically on IEEE systems. On such systems,
125 * when INFNAN_CHECK is #defined, strtod checks
126 * for Infinity and NaN (case insensitively). On some systems
127 * (e.g., some HP systems), it may be necessary to #define NAN_WORD0
128 * appropriately -- to the most significant word of a quiet NaN.
129 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
130 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
131 * strtod also accepts (case insensitively) strings of the form
132 * NaN(x), where x is a string of hexadecimal digits and spaces;
133 * if there is only one string of hexadecimal digits, it is taken
134 * for the 52 fraction bits of the resulting NaN; if there are two
135 * or more strings of hex digits, the first is for the high 20 bits,
136 * the second and subsequent for the low 32 bits, with intervening
137 * white space ignored; but if this results in none of the 52
138 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
139 * and NAN_WORD1 are used instead.
140 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
141 * multiple threads. In this case, you must provide (or suitably
142 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
143 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
144 * in pow5mult, ensures lazy evaluation of only one copy of high
145 * powers of 5; omitting this lock would introduce a small
146 * probability of wasting memory, but would otherwise be harmless.)
147 * You must also invoke freedtoa(s) to free the value s returned by
148 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
149 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
150 * avoids underflows on inputs whose result does not underflow.
151 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
152 * floating-point numbers and flushes underflows to zero rather
153 * than implementing gradual underflow, then you must also #define
155 * #define USE_LOCALE to use the current locale's decimal_point value.
156 * #define SET_INEXACT if IEEE arithmetic is being used and extra
157 * computation should be done to set the inexact flag when the
158 * result is inexact and avoid setting inexact when the result
159 * is exact. In this case, dtoa.c must be compiled in
160 * an environment, perhaps provided by #include "dtoa.c" in a
161 * suitable wrapper, that defines two functions,
162 * int get_inexact(void);
163 * void clear_inexact(void);
164 * such that get_inexact() returns a nonzero value if the
165 * inexact bit is already set, and clear_inexact() sets the
166 * inexact bit to 0. When SET_INEXACT is #defined, strtod
167 * also does extra computations to set the underflow and overflow
168 * flags when appropriate (i.e., when the result is tiny and
169 * inexact or when it is a numeric value rounded to +-infinity).
170 * #define NO_ERRNO if strtod should not assign errno = ERANGE when
171 * the result overflows to +-Infinity or underflows to 0.
172 * #define NO_HEX_FP to omit recognition of hexadecimal floating-point
174 * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now)
175 * to disable logic for "fast" testing of very long input strings
176 * to strtod. This testing proceeds by initially truncating the
177 * input string, then if necessary comparing the whole string with
178 * a decimal expansion to decide close cases. This logic is only
179 * used for input more than STRTOD_DIGLIM digits long (default 40).
182 #if defined _MSC_VER && _MSC_VER == 1800
183 // TODO(scottmg): VS2013 RC ICEs on a bunch of functions in this file.
184 // This should be removed after RTM. See http://crbug.com/288948.
185 #pragma optimize("", off)
199 typedef unsigned Long ULong;
204 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
214 #ifdef Honor_FLT_ROUNDS
215 #ifndef Trust_FLT_ROUNDS
222 extern char *MALLOC();
224 extern void *MALLOC(size_t);
227 #define MALLOC malloc
230 #ifndef Omit_Private_Memory
232 #define PRIVATE_MEM 2304
234 #define PRIVATE_mem ((unsigned)((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)))
235 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
239 #undef Avoid_Underflow
248 #ifndef NO_INFNAN_CHECK
254 #define NO_STRTOD_BIGCOMP
263 #define DBL_MAX_10_EXP 308
264 #define DBL_MAX_EXP 1024
266 #endif /*IEEE_Arith*/
270 #define DBL_MAX_10_EXP 75
271 #define DBL_MAX_EXP 63
273 #define DBL_MAX 7.2370055773322621e+75
278 #define DBL_MAX_10_EXP 38
279 #define DBL_MAX_EXP 127
281 #define DBL_MAX 1.7014118346046923e+38
285 #define LONG_MAX 2147483647
288 #else /* ifndef Bad_float_h */
290 #endif /* Bad_float_h */
300 #define CONST /* blank */
306 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
307 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
310 typedef union { double d; ULong L[2]; } U;
313 #define word0(x) (x)->L[1]
314 #define word1(x) (x)->L[0]
316 #define word0(x) (x)->L[0]
317 #define word1(x) (x)->L[1]
319 #define dval(x) (x)->d
321 #ifndef STRTOD_DIGLIM
322 #define STRTOD_DIGLIM 40
326 extern int strtod_diglim;
328 #define strtod_diglim STRTOD_DIGLIM
331 /* The following definition of Storeinc is appropriate for MIPS processors.
332 * An alternative that might be better on some machines is
333 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
335 #if defined(IEEE_8087) + defined(VAX)
336 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
337 ((unsigned short *)a)[0] = (unsigned short)c, a++)
339 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
340 ((unsigned short *)a)[1] = (unsigned short)c, a++)
343 /* #define P DBL_MANT_DIG */
344 /* Ten_pmax = floor(P*log(2)/log(5)) */
345 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
346 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
347 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
351 #define Exp_shift1 20
352 #define Exp_msk1 0x100000
353 #define Exp_msk11 0x100000
354 #define Exp_mask 0x7ff00000
360 #define Exp_1 0x3ff00000
361 #define Exp_11 0x3ff00000
363 #define Frac_mask 0xfffff
364 #define Frac_mask1 0xfffff
367 #define Bndry_mask 0xfffff
368 #define Bndry_mask1 0xfffff
370 #define Sign_bit 0x80000000
376 #ifndef NO_IEEE_Scale
377 #define Avoid_Underflow
378 #ifdef Flush_Denorm /* debugging option */
379 #undef Sudden_Underflow
385 #define Flt_Rounds FLT_ROUNDS
389 #endif /*Flt_Rounds*/
391 #ifdef Honor_FLT_ROUNDS
392 #undef Check_FLT_ROUNDS
393 #define Check_FLT_ROUNDS
395 #define Rounding Flt_Rounds
398 #else /* ifndef IEEE_Arith */
399 #undef Check_FLT_ROUNDS
400 #undef Honor_FLT_ROUNDS
402 #undef Sudden_Underflow
403 #define Sudden_Underflow
408 #define Exp_shift1 24
409 #define Exp_msk1 0x1000000
410 #define Exp_msk11 0x1000000
411 #define Exp_mask 0x7f000000
417 #define Exp_1 0x41000000
418 #define Exp_11 0x41000000
419 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
420 #define Frac_mask 0xffffff
421 #define Frac_mask1 0xffffff
424 #define Bndry_mask 0xefffff
425 #define Bndry_mask1 0xffffff
427 #define Sign_bit 0x80000000
429 #define Tiny0 0x100000
438 #define Exp_msk1 0x80
439 #define Exp_msk11 0x800000
440 #define Exp_mask 0x7f80
446 #define Exp_1 0x40800000
447 #define Exp_11 0x4080
449 #define Frac_mask 0x7fffff
450 #define Frac_mask1 0xffff007f
453 #define Bndry_mask 0xffff007f
454 #define Bndry_mask1 0xffff007f
456 #define Sign_bit 0x8000
462 #endif /* IBM, VAX */
463 #endif /* IEEE_Arith */
470 #define rounded_product(a,b) a = rnd_prod(a, b)
471 #define rounded_quotient(a,b) a = rnd_quot(a, b)
473 extern double rnd_prod(), rnd_quot();
475 extern double rnd_prod(double, double), rnd_quot(double, double);
478 #define rounded_product(a,b) a *= b
479 #define rounded_quotient(a,b) a /= b
482 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
483 #define Big1 0xffffffff
489 typedef struct BCinfo BCinfo;
491 BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflchk; };
494 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
496 #define FFFFFFFF 0xffffffffUL
503 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
504 * This makes some inner loops simpler and sometimes saves work
505 * during multiplications, but it often seems to make things slightly
506 * slower. Hence the default is now to store 32 bits per Long.
509 #else /* long long available */
511 #define Llong long long
514 #define ULLong unsigned Llong
516 #endif /* NO_LONG_LONG */
518 #ifndef MULTIPLE_THREADS
519 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
520 #define FREE_DTOA_LOCK(n) /*nothing*/
525 double strtod(const char *s00, char **se);
526 char *dtoa(double d, int mode, int ndigits,
527 int *decpt, int *sign, char **rve);
532 int k, maxwds, sign, wds;
536 typedef struct Bigint Bigint;
538 static Bigint *freelist[Kmax+1];
550 #ifndef Omit_Private_Memory
554 ACQUIRE_DTOA_LOCK(0);
555 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
556 /* but this case seems very unlikely. */
557 if (k <= Kmax && (rv = freelist[k]))
558 freelist[k] = rv->next;
561 #ifdef Omit_Private_Memory
562 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
564 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
566 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
567 rv = (Bigint*)pmem_next;
571 rv = (Bigint*)MALLOC(len*sizeof(double));
577 rv->sign = rv->wds = 0;
597 ACQUIRE_DTOA_LOCK(0);
598 v->next = freelist[v->k];
605 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
606 y->wds*sizeof(Long) + 2*sizeof(int))
611 (b, m, a) Bigint *b; int m, a;
613 (Bigint *b, int m, int a) /* multiply by m and add a */
634 y = *x * (ULLong)m + carry;
640 y = (xi & 0xffff) * m + carry;
641 z = (xi >> 16) * m + (y >> 16);
643 *x++ = (z << 16) + (y & 0xffff);
653 if (wds >= b->maxwds) {
668 (s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9;
670 (CONST char *s, int nd0, int nd, ULong y9, int dplen)
678 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
685 b->x[0] = y9 & 0xffff;
686 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
692 do b = multadd(b, 10, *s++ - '0');
699 b = multadd(b, 10, *s++ - '0');
713 if (!(x & 0xffff0000)) {
717 if (!(x & 0xff000000)) {
721 if (!(x & 0xf0000000)) {
725 if (!(x & 0xc0000000)) {
729 if (!(x & 0x80000000)) {
731 if (!(x & 0x40000000))
804 (a, b) Bigint *a, *b;
806 (Bigint *a, Bigint *b)
811 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
822 if (a->wds < b->wds) {
834 for(x = c->x, xa = x + wc; x < xa; x++)
842 for(; xb < xbe; xc0++) {
848 z = *x++ * (ULLong)y + *xc + carry;
850 *xc++ = z & FFFFFFFF;
858 for(; xb < xbe; xb++, xc0++) {
859 if (y = *xb & 0xffff) {
864 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
866 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
879 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
882 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
890 for(; xb < xbe; xc0++) {
896 z = *x++ * y + *xc + carry;
906 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
916 (b, k) Bigint *b; int k;
921 Bigint *b1, *p5, *p51;
923 static int p05[3] = { 5, 25, 125 };
926 b = multadd(b, p05[i-1], 0);
932 #ifdef MULTIPLE_THREADS
933 ACQUIRE_DTOA_LOCK(1);
952 if (!(p51 = p5->next)) {
953 #ifdef MULTIPLE_THREADS
954 ACQUIRE_DTOA_LOCK(1);
955 if (!(p51 = p5->next)) {
956 p51 = p5->next = mult(p5,p5);
961 p51 = p5->next = mult(p5,p5);
973 (b, k) Bigint *b; int k;
980 ULong *x, *x1, *xe, z;
989 for(i = b->maxwds; n1 > i; i <<= 1)
993 for(i = 0; i < n; i++)
1002 *x1++ = *x << k | z;
1014 *x1++ = *x << k & 0xffff | z;
1033 (a, b) Bigint *a, *b;
1035 (Bigint *a, Bigint *b)
1038 ULong *xa, *xa0, *xb, *xb0;
1044 if (i > 1 && !a->x[i-1])
1045 Bug("cmp called with a->x[a->wds-1] == 0");
1046 if (j > 1 && !b->x[j-1])
1047 Bug("cmp called with b->x[b->wds-1] == 0");
1057 return *xa < *xb ? -1 : 1;
1067 (a, b) Bigint *a, *b;
1069 (Bigint *a, Bigint *b)
1074 ULong *xa, *xae, *xb, *xbe, *xc;
1111 y = (ULLong)*xa++ - *xb++ - borrow;
1112 borrow = y >> 32 & (ULong)1;
1113 *xc++ = y & FFFFFFFF;
1118 borrow = y >> 32 & (ULong)1;
1119 *xc++ = y & FFFFFFFF;
1124 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1125 borrow = (y & 0x10000) >> 16;
1126 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1127 borrow = (z & 0x10000) >> 16;
1132 y = (*xa & 0xffff) - borrow;
1133 borrow = (y & 0x10000) >> 16;
1134 z = (*xa++ >> 16) - borrow;
1135 borrow = (z & 0x10000) >> 16;
1140 y = *xa++ - *xb++ - borrow;
1141 borrow = (y & 0x10000) >> 16;
1147 borrow = (y & 0x10000) >> 16;
1169 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1170 #ifndef Avoid_Underflow
1171 #ifndef Sudden_Underflow
1180 #ifndef Avoid_Underflow
1181 #ifndef Sudden_Underflow
1184 L = -L >> Exp_shift;
1185 if (L < Exp_shift) {
1186 word0(&u) = 0x80000 >> L;
1192 word1(&u) = L >= 31 ? 1 : 1 << 31 - L;
1203 (a, e) Bigint *a; int *e;
1208 ULong *xa, *xa0, w, y, z;
1214 #define d0 word0(&d)
1215 #define d1 word1(&d)
1222 if (!y) Bug("zero y in b2d");
1228 d0 = Exp_1 | y >> (Ebits - k);
1229 w = xa > xa0 ? *--xa : 0;
1230 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1233 z = xa > xa0 ? *--xa : 0;
1235 d0 = Exp_1 | y << k | z >> (32 - k);
1236 y = xa > xa0 ? *--xa : 0;
1237 d1 = z << k | y >> (32 - k);
1244 if (k < Ebits + 16) {
1245 z = xa > xa0 ? *--xa : 0;
1246 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1247 w = xa > xa0 ? *--xa : 0;
1248 y = xa > xa0 ? *--xa : 0;
1249 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1252 z = xa > xa0 ? *--xa : 0;
1253 w = xa > xa0 ? *--xa : 0;
1255 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1256 y = xa > xa0 ? *--xa : 0;
1257 d1 = w << k + 16 | y << k;
1261 word0(&d) = d0 >> 16 | d0 << 16;
1262 word1(&d) = d1 >> 16 | d1 << 16;
1273 (d, e, bits) U *d; int *e, *bits;
1275 (U *d, int *e, int *bits)
1281 #ifndef Sudden_Underflow
1286 d0 = word0(d) >> 16 | word0(d) << 16;
1287 d1 = word1(d) >> 16 | word1(d) << 16;
1301 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1302 #ifdef Sudden_Underflow
1303 de = (int)(d0 >> Exp_shift);
1308 if ((de = (int)(d0 >> Exp_shift)))
1313 if ((k = lo0bits(&y))) {
1314 x[0] = y | z << (32 - k);
1319 #ifndef Sudden_Underflow
1322 b->wds = (x[1] = z) ? 2 : 1;
1327 #ifndef Sudden_Underflow
1335 if (k = lo0bits(&y))
1337 x[0] = y | z << 32 - k & 0xffff;
1338 x[1] = z >> k - 16 & 0xffff;
1344 x[1] = y >> 16 | z << 16 - k & 0xffff;
1345 x[2] = z >> k & 0xffff;
1360 Bug("Zero passed to d2b");
1378 #ifndef Sudden_Underflow
1382 *e = (de - Bias - (P-1) << 2) + k;
1383 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1385 *e = de - Bias - (P-1) + k;
1388 #ifndef Sudden_Underflow
1391 *e = de - Bias - (P-1) + 1 + k;
1393 *bits = 32*i - hi0bits(x[i-1]);
1395 *bits = (i+2)*16 - hi0bits(x[i]);
1407 (a, b) Bigint *a, *b;
1409 (Bigint *a, Bigint *b)
1415 dval(&da) = b2d(a, &ka);
1416 dval(&db) = b2d(b, &kb);
1418 k = ka - kb + 32*(a->wds - b->wds);
1420 k = ka - kb + 16*(a->wds - b->wds);
1424 word0(&da) += (k >> 2)*Exp_msk1;
1426 dval(&da) *= 1 << k;
1430 word0(&db) += (k >> 2)*Exp_msk1;
1432 dval(&db) *= 1 << k;
1436 word0(&da) += k*Exp_msk1;
1439 word0(&db) += k*Exp_msk1;
1442 return dval(&da) / dval(&db);
1447 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1448 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1457 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1458 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1459 #ifdef Avoid_Underflow
1460 9007199254740992.*9007199254740992.e-256
1461 /* = 2^106 * 1e-256 */
1466 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1467 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1468 #define Scale_Bit 0x10
1472 bigtens[] = { 1e16, 1e32, 1e64 };
1473 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1476 bigtens[] = { 1e16, 1e32 };
1477 static CONST double tinytens[] = { 1e-16, 1e-32 };
1495 #ifdef Need_Hexdig /*{*/
1496 static unsigned char hexdig[256];
1500 htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc;
1502 htinit(unsigned char *h, unsigned char *s, int inc)
1506 for(i = 0; (j = s[i]) !=0; i++)
1517 #define USC (unsigned char *)
1518 htinit(hexdig, USC "0123456789", 0x10);
1519 htinit(hexdig, USC "abcdef", 0x10 + 10);
1520 htinit(hexdig, USC "ABCDEF", 0x10 + 10);
1522 #endif /* } Need_Hexdig */
1527 #define NAN_WORD0 0x7ff80000
1537 (sp, t) char **sp, *t;
1539 (CONST char **sp, CONST char *t)
1543 CONST char *s = *sp;
1546 if ((c = *++s) >= 'A' && c <= 'Z')
1559 (rvp, sp) U *rvp; CONST char **sp;
1561 (U *rvp, CONST char **sp)
1566 int c1, havedig, udx0, xshift;
1571 havedig = xshift = 0;
1574 /* allow optional initial 0x or 0X */
1575 while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
1577 if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
1579 while((c = *(CONST unsigned char*)++s)) {
1580 if ((c1 = hexdig[c]))
1582 else if (c <= ' ') {
1583 if (udx0 && havedig) {
1589 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
1590 else if (/*(*/ c == ')' && havedig) {
1595 return; /* invalid form: don't change *sp */
1599 if (/*(*/ c == ')') {
1603 } while((c = *++s));
1614 x[0] = (x[0] << 4) | (x[1] >> 28);
1615 x[1] = (x[1] << 4) | c;
1617 if ((x[0] &= 0xfffff) || x[1]) {
1618 word0(rvp) = Exp_mask | x[0];
1622 #endif /*No_Hex_NaN*/
1623 #endif /* INFNAN_CHECK */
1634 #ifndef NO_HEX_FP /*{*/
1638 rshift(b, k) Bigint *b; int k;
1640 rshift(Bigint *b, int k)
1643 ULong *x, *x1, *xe, y;
1655 *x1++ = (y | (*x << n)) & 0xffffffff;
1665 if ((b->wds = x1 - b->x) == 0)
1671 any_on(b, k) Bigint *b; int k;
1673 any_on(Bigint *b, int k)
1677 ULong *x, *x0, x1, x2;
1684 else if (n < nwds && (k &= kmask)) {
1699 enum { /* rounding values: same as FLT_ROUNDS */
1708 increment(b) Bigint *b;
1710 increment(Bigint *b)
1719 if (*x < (ULong)0xffffffffL) {
1726 if (b->wds >= b->maxwds) {
1727 b1 = Balloc(b->k+1);
1739 gethex(sp, rvp, rounding, sign)
1740 CONST char **sp; U *rvp; int rounding, sign;
1742 gethex( CONST char **sp, U *rvp, int rounding, int sign)
1746 CONST unsigned char *decpt, *s0, *s, *s1;
1748 ULong L, lostbits, *x;
1749 int big, denorm, esign, havedig, k, n, nbits, up, zret;
1754 #ifdef IEEE_Arith /*{{*/
1755 emax = 0x7fe - Bias - P + 1,
1760 emax = 0x7ff - Bias - P + 1
1763 emax = 0x7f - Bias - P
1769 #ifdef NO_LOCALE_CACHE
1770 const unsigned char *decimalpoint = (unsigned char*)
1771 localeconv()->decimal_point;
1773 const unsigned char *decimalpoint;
1774 static unsigned char *decimalpoint_cache;
1775 if (!(s0 = decimalpoint_cache)) {
1776 s0 = (unsigned char*)localeconv()->decimal_point;
1777 if ((decimalpoint_cache = (unsigned char*)
1778 MALLOC(strlen((CONST char*)s0) + 1))) {
1779 strcpy((char*)decimalpoint_cache, (CONST char*)s0);
1780 s0 = decimalpoint_cache;
1790 s0 = *(CONST unsigned char **)sp + 2;
1791 while(s0[havedig] == '0')
1803 for(i = 0; decimalpoint[i]; ++i) {
1804 if (s[i] != decimalpoint[i])
1825 if (*s == *decimalpoint && !decpt) {
1826 for(i = 1; decimalpoint[i]; ++i) {
1827 if (s[i] != decimalpoint[i])
1832 if (*s == '.' && !decpt) {
1839 e = -(((Long)(s-decpt)) << 2);
1853 if ((n = hexdig[*s]) == 0 || n > 0x19) {
1858 while((n = hexdig[*++s]) !=0 && n <= 0x19) {
1859 if (e1 & 0xf8000000)
1861 e1 = 10*e1 + n - 0x10;
1869 *sp = (char*)s0 - 1;
1895 #endif /* IEEE_Arith */
1915 for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
1922 for(i = 0; decimalpoint[i+1]; ++i);
1926 if (*--s1 == decimalpoint[i]) {
1939 L |= (hexdig[*s1] & 0x0f) << n;
1943 b->wds = n = x - b->x;
1944 n = ULbits*n - hi0bits(L);
1953 if (x[k>>kshift] & 1 << (k & kmask)) {
1955 if (k > 0 && any_on(b,k))
1962 else if (n < nbits) {
1975 word0(rvp) = Exp_mask;
1984 #ifdef IEEE_Arith /*{*/
1987 if (n == nbits && (n < 2 || any_on(b,n-1)))
1998 #endif /* } IEEE_Arith */
2012 lostbits = any_on(b,k);
2013 if (x[k>>kshift] & 1 << (k & kmask))
2026 && (lostbits & 1) | (x[0] & 1))
2041 if (nbits == Nbits - 1
2042 && x[nbits >> kshift] & 1 << (nbits & kmask))
2043 denorm = 0; /* not currently used */
2047 || ((n = nbits & kmask) !=0
2048 && hi0bits(x[k-1]) < 32-n)) {
2057 word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0;
2059 word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20);
2060 word1(rvp) = b->x[0];
2064 k = b->x[0] & ((1 << j) - 1);
2078 if (k & j && ((k & (j-1)) | lostbits))
2084 word0(rvp) = b->x[1] | ((e + 65 + 13) << 24);
2085 word1(rvp) = b->x[0];
2088 /* The next two lines ignore swap of low- and high-order 2 bytes. */
2089 /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */
2090 /* word1(rvp) = b->x[0]; */
2091 word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->x[1] << 16);
2092 word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16);
2096 #endif /*}!NO_HEX_FP*/
2100 dshift(b, p2) Bigint *b; int p2;
2102 dshift(Bigint *b, int p2)
2105 int rv = hi0bits(b->x[b->wds-1]) - 4;
2114 (b, S) Bigint *b, *S;
2116 (Bigint *b, Bigint *S)
2120 ULong *bx, *bxe, q, *sx, *sxe;
2122 ULLong borrow, carry, y, ys;
2124 ULong borrow, carry, y, ys;
2132 /*debug*/ if (b->wds > n)
2133 /*debug*/ Bug("oversize b in quorem");
2141 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
2143 /*debug*/ if (q > 9)
2144 /*debug*/ Bug("oversized quotient in quorem");
2151 ys = *sx++ * (ULLong)q + carry;
2153 y = *bx - (ys & FFFFFFFF) - borrow;
2154 borrow = y >> 32 & (ULong)1;
2155 *bx++ = y & FFFFFFFF;
2159 ys = (si & 0xffff) * q + carry;
2160 zs = (si >> 16) * q + (ys >> 16);
2162 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2163 borrow = (y & 0x10000) >> 16;
2164 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2165 borrow = (z & 0x10000) >> 16;
2168 ys = *sx++ * q + carry;
2170 y = *bx - (ys & 0xffff) - borrow;
2171 borrow = (y & 0x10000) >> 16;
2179 while(--bxe > bx && !*bxe)
2184 if (cmp(b, S) >= 0) {
2194 y = *bx - (ys & FFFFFFFF) - borrow;
2195 borrow = y >> 32 & (ULong)1;
2196 *bx++ = y & FFFFFFFF;
2200 ys = (si & 0xffff) + carry;
2201 zs = (si >> 16) + (ys >> 16);
2203 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2204 borrow = (y & 0x10000) >> 16;
2205 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2206 borrow = (z & 0x10000) >> 16;
2211 y = *bx - (ys & 0xffff) - borrow;
2212 borrow = (y & 0x10000) >> 16;
2221 while(--bxe > bx && !*bxe)
2229 #ifndef NO_STRTOD_BIGCOMP
2235 U *rv; CONST char *s0; BCinfo *bc;
2237 (U *rv, CONST char *s0, BCinfo *bc)
2241 int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase;
2246 p5 = nd + bc->e0 - 1;
2248 #ifndef Sudden_Underflow
2249 if (rv->d == 0.) { /* special case: value near underflow-to-zero */
2250 /* threshold was rounded to zero */
2254 #ifdef Avoid_Underflow
2255 word0(rv) = (P+2) << Exp_shift;
2260 #ifdef Honor_FLT_ROUNDS
2261 if (bc->rounding == 1)
2272 b = d2b(rv, &p2, &bbits);
2273 #ifdef Avoid_Underflow
2276 /* floor(log2(rv)) == bbits - 1 + p2 */
2277 /* Check for denormal case. */
2279 if (i > (j = P - Emin - 1 + p2)) {
2280 #ifdef Sudden_Underflow
2285 #ifdef Avoid_Underflow
2286 word0(rv) = (1 + bc->scale) << Exp_shift;
2288 word0(rv) = Exp_msk1;
2295 #ifdef Honor_FLT_ROUNDS
2296 if (bc->rounding != 1) {
2308 #ifndef Sudden_Underflow
2313 /* Arrange for convenient computation of quotients:
2314 * shift left if necessary so divisor has 4 leading 0 bits.
2317 d = pow5mult(d, p5);
2319 b = pow5mult(b, -p5);
2334 /* Now b/d = exactly half-way between the two floating-point values */
2335 /* on either side of the input string. Compute first digit of b/d. */
2337 if (!(dig = quorem(b,d))) {
2338 b = multadd(b, 10, 0); /* very unlikely */
2342 /* Compare b/d with s0 */
2344 for(i = 0; i < nd0; ) {
2345 if ((dd = s0[i++] - '0' - dig))
2347 if (!b->x[0] && b->wds == 1) {
2352 b = multadd(b, 10, 0);
2355 for(j = bc->dp1; i++ < nd;) {
2356 if ((dd = s0[j++] - '0' - dig))
2358 if (!b->x[0] && b->wds == 1) {
2363 b = multadd(b, 10, 0);
2366 if (b->x[0] || b->wds > 1)
2371 #ifdef Honor_FLT_ROUNDS
2372 if (bc->rounding != 1) {
2374 if (bc->rounding == 0) {
2382 if (bc->rounding == 0) {
2389 dval(rv) += 2.*ulp(rv);
2404 if (!dsign) /* does not happen for round-near */
2406 dval(rv) -= ulp(rv);
2411 dval(rv) += ulp(rv);
2415 /* Exact half-way case: apply round-even rule. */
2416 if (word1(rv) & 1) {
2423 #ifdef Honor_FLT_ROUNDS
2428 #endif /* NO_STRTOD_BIGCOMP */
2433 (s00, se) CONST char *s00; char **se;
2435 (CONST char *s00, char **se)
2438 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1;
2439 int esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
2440 CONST char *s, *s0, *s1;
2443 U aadj2, adj, rv, rv0;
2446 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
2450 #ifdef Honor_FLT_ROUNDS /*{*/
2451 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
2452 bc.rounding = Flt_Rounds;
2455 switch(fegetround()) {
2456 case FE_TOWARDZERO: bc.rounding = 0; break;
2457 case FE_UPWARD: bc.rounding = 2; break;
2458 case FE_DOWNWARD: bc.rounding = 3;
2466 sign = nz0 = nz = bc.dplen = bc.uflchk = 0;
2468 for(s = s00;;s++) switch(*s) {
2490 #ifndef NO_HEX_FP /*{*/
2494 #ifdef Honor_FLT_ROUNDS
2495 gethex(&s, &rv, bc.rounding, sign);
2497 gethex(&s, &rv, 1, sign);
2503 while(*++s == '0') ;
2509 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
2515 bc.dp0 = bc.dp1 = s - s0;
2517 s1 = localeconv()->decimal_point;
2538 bc.dplen = bc.dp1 - bc.dp0;
2540 for(; c == '0'; c = *++s)
2542 if (c > '0' && c <= '9') {
2550 for(; c >= '0' && c <= '9'; c = *++s) {
2555 for(i = 1; i < nz; i++)
2558 else if (nd <= DBL_DIG + 1)
2562 else if (nd <= DBL_DIG + 1)
2570 if (c == 'e' || c == 'E') {
2571 if (!nd && !nz && !nz0) {
2582 if (c >= '0' && c <= '9') {
2585 if (c > '0' && c <= '9') {
2588 while((c = *++s) >= '0' && c <= '9')
2590 if (s - s1 > 8 || L > 19999)
2591 /* Avoid confusion from exponents
2592 * so large that e might overflow.
2594 e = 19999; /* safe for 16 bit ints */
2609 /* Check for Nan and Infinity */
2614 if (match(&s,"nf")) {
2616 if (!match(&s,"inity"))
2618 word0(&rv) = 0x7ff00000;
2625 if (match(&s, "an")) {
2626 word0(&rv) = NAN_WORD0;
2627 word1(&rv) = NAN_WORD1;
2629 if (*s == '(') /*)*/
2635 #endif /* INFNAN_CHECK */
2642 bc.e0 = e1 = e -= nf;
2644 /* Now we have nd0 digits, starting at s0, followed by a
2645 * decimal point, followed by nd-nd0 digits. The number we're
2646 * after is the integer represented by those digits times
2651 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
2656 oldinexact = get_inexact();
2658 dval(&rv) = tens[k - 9] * dval(&rv) + z;
2662 #ifndef RND_PRODQUOT
2663 #ifndef Honor_FLT_ROUNDS
2671 if (e <= Ten_pmax) {
2673 goto vax_ovfl_check;
2675 #ifdef Honor_FLT_ROUNDS
2676 /* round correctly FLT_ROUNDS = 2 or 3 */
2682 /* rv = */ rounded_product(dval(&rv), tens[e]);
2687 if (e <= Ten_pmax + i) {
2688 /* A fancier test would sometimes let us do
2689 * this for larger i values.
2691 #ifdef Honor_FLT_ROUNDS
2692 /* round correctly FLT_ROUNDS = 2 or 3 */
2699 dval(&rv) *= tens[i];
2701 /* VAX exponent range is so narrow we must
2702 * worry about overflow here...
2705 word0(&rv) -= P*Exp_msk1;
2706 /* rv = */ rounded_product(dval(&rv), tens[e]);
2707 if ((word0(&rv) & Exp_mask)
2708 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
2710 word0(&rv) += P*Exp_msk1;
2712 /* rv = */ rounded_product(dval(&rv), tens[e]);
2717 #ifndef Inaccurate_Divide
2718 else if (e >= -Ten_pmax) {
2719 #ifdef Honor_FLT_ROUNDS
2720 /* round correctly FLT_ROUNDS = 2 or 3 */
2726 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
2737 oldinexact = get_inexact();
2739 #ifdef Avoid_Underflow
2742 #ifdef Honor_FLT_ROUNDS
2743 if (bc.rounding >= 2) {
2745 bc.rounding = bc.rounding == 2 ? 0 : 2;
2747 if (bc.rounding != 2)
2751 #endif /*IEEE_Arith*/
2753 /* Get starting approximation = rv * 10**e1 */
2757 dval(&rv) *= tens[i];
2759 if (e1 > DBL_MAX_10_EXP) {
2764 /* Can't trust HUGE_VAL */
2766 #ifdef Honor_FLT_ROUNDS
2767 switch(bc.rounding) {
2768 case 0: /* toward 0 */
2769 case 3: /* toward -infinity */
2774 word0(&rv) = Exp_mask;
2777 #else /*Honor_FLT_ROUNDS*/
2778 word0(&rv) = Exp_mask;
2780 #endif /*Honor_FLT_ROUNDS*/
2782 /* set overflow bit */
2784 dval(&rv0) *= dval(&rv0);
2786 #else /*IEEE_Arith*/
2789 #endif /*IEEE_Arith*/
2793 for(j = 0; e1 > 1; j++, e1 >>= 1)
2795 dval(&rv) *= bigtens[j];
2796 /* The last multiplication could overflow. */
2797 word0(&rv) -= P*Exp_msk1;
2798 dval(&rv) *= bigtens[j];
2799 if ((z = word0(&rv) & Exp_mask)
2800 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
2802 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
2803 /* set to largest number */
2804 /* (Can't trust DBL_MAX) */
2809 word0(&rv) += P*Exp_msk1;
2815 dval(&rv) /= tens[i];
2817 if (e1 >= 1 << n_bigtens)
2819 #ifdef Avoid_Underflow
2822 for(j = 0; e1 > 0; j++, e1 >>= 1)
2824 dval(&rv) *= tinytens[j];
2825 if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
2826 >> Exp_shift)) > 0) {
2827 /* scaled rv is denormal; clear j low bits */
2831 word0(&rv) = (P+2)*Exp_msk1;
2833 word0(&rv) &= 0xffffffff << (j-32);
2836 word1(&rv) &= 0xffffffff << j;
2839 for(j = 0; e1 > 1; j++, e1 >>= 1)
2841 dval(&rv) *= tinytens[j];
2842 /* The last multiplication could underflow. */
2843 dval(&rv0) = dval(&rv);
2844 dval(&rv) *= tinytens[j];
2846 dval(&rv) = 2.*dval(&rv0);
2847 dval(&rv) *= tinytens[j];
2857 #ifndef Avoid_Underflow
2860 /* The refinement below will clean
2861 * this approximation up.
2868 /* Now the hard part -- adjusting rv to the correct value.*/
2870 /* Put digits into bd: true value = bd * 10^e */
2873 #ifndef NO_STRTOD_BIGCOMP
2874 bc.nd0 = nd0; /* Only needed if nd > strtod_diglim, but done here */
2875 /* to silence an erroneous warning about bc.nd0 */
2876 /* possibly not being initialized. */
2877 if (nd > strtod_diglim) {
2878 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */
2879 /* minimum number of decimal digits to distinguish double values */
2880 /* in IEEE arithmetic. */
2885 if (--j <= bc.dp1 && j >= bc.dp0)
2895 if (nd < 9) { /* must recompute y */
2897 for(i = 0; i < nd0; ++i)
2898 y = 10*y + s0[i] - '0';
2899 for(j = bc.dp1; i < nd; ++i)
2900 y = 10*y + s0[j++] - '0';
2904 bd0 = s2b(s0, nd0, nd, y, bc.dplen);
2907 bd = Balloc(bd0->k);
2909 bb = d2b(&rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
2925 #ifdef Honor_FLT_ROUNDS
2926 if (bc.rounding != 1)
2929 #ifdef Avoid_Underflow
2931 i = j + bbbits - 1; /* logb(rv) */
2932 if (i < Emin) /* denormal */
2936 #else /*Avoid_Underflow*/
2937 #ifdef Sudden_Underflow
2939 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
2943 #else /*Sudden_Underflow*/
2945 i = j + bbbits - 1; /* logb(rv) */
2946 if (i < Emin) /* denormal */
2950 #endif /*Sudden_Underflow*/
2951 #endif /*Avoid_Underflow*/
2954 #ifdef Avoid_Underflow
2957 i = bb2 < bd2 ? bb2 : bd2;
2966 bs = pow5mult(bs, bb5);
2972 bb = lshift(bb, bb2);
2974 bd = pow5mult(bd, bd5);
2976 bd = lshift(bd, bd2);
2978 bs = lshift(bs, bs2);
2979 delta = diff(bb, bd);
2980 bc.dsign = delta->sign;
2983 #ifndef NO_STRTOD_BIGCOMP
2984 if (bc.nd > nd && i <= 0) {
2986 break; /* Must use bigcomp(). */
2987 #ifdef Honor_FLT_ROUNDS
2988 if (bc.rounding != 1) {
2996 i = -1; /* Discarded digits make delta smaller. */
3000 #ifdef Honor_FLT_ROUNDS
3001 if (bc.rounding != 1) {
3003 /* Error is less than an ulp */
3004 if (!delta->x[0] && delta->wds <= 1) {
3017 else if (!bc.dsign) {
3020 && !(word0(&rv) & Frac_mask)) {
3021 y = word0(&rv) & Exp_mask;
3022 #ifdef Avoid_Underflow
3023 if (!bc.scale || y > 2*P*Exp_msk1)
3028 delta = lshift(delta,Log2P);
3029 if (cmp(delta, bs) <= 0)
3034 #ifdef Avoid_Underflow
3035 if (bc.scale && (y = word0(&rv) & Exp_mask)
3037 word0(&adj) += (2*P+1)*Exp_msk1 - y;
3039 #ifdef Sudden_Underflow
3040 if ((word0(&rv) & Exp_mask) <=
3042 word0(&rv) += P*Exp_msk1;
3043 dval(&rv) += adj.d*ulp(dval(&rv));
3044 word0(&rv) -= P*Exp_msk1;
3047 #endif /*Sudden_Underflow*/
3048 #endif /*Avoid_Underflow*/
3049 dval(&rv) += adj.d*ulp(&rv);
3053 adj.d = ratio(delta, bs);
3056 if (adj.d <= 0x7ffffffe) {
3057 /* adj = rounding ? ceil(adj) : floor(adj); */
3060 if (!((bc.rounding>>1) ^ bc.dsign))
3065 #ifdef Avoid_Underflow
3066 if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
3067 word0(&adj) += (2*P+1)*Exp_msk1 - y;
3069 #ifdef Sudden_Underflow
3070 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
3071 word0(&rv) += P*Exp_msk1;
3072 adj.d *= ulp(dval(&rv));
3077 word0(&rv) -= P*Exp_msk1;
3080 #endif /*Sudden_Underflow*/
3081 #endif /*Avoid_Underflow*/
3084 if (word0(&rv) == Big0 && word1(&rv) == Big1)
3092 #endif /*Honor_FLT_ROUNDS*/
3095 /* Error is less than half an ulp -- check for
3096 * special case of mantissa a power of two.
3098 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask
3100 #ifdef Avoid_Underflow
3101 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
3103 || (word0(&rv) & Exp_mask) <= Exp_msk1
3108 if (!delta->x[0] && delta->wds <= 1)
3113 if (!delta->x[0] && delta->wds <= 1) {
3120 delta = lshift(delta,Log2P);
3121 if (cmp(delta, bs) > 0)
3126 /* exactly half-way between */
3128 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
3130 #ifdef Avoid_Underflow
3131 (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
3132 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
3135 /*boundary case -- increment exponent*/
3136 word0(&rv) = (word0(&rv) & Exp_mask)
3143 #ifdef Avoid_Underflow
3149 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
3151 /* boundary case -- decrement exponent */
3152 #ifdef Sudden_Underflow /*{{*/
3153 L = word0(&rv) & Exp_mask;
3157 #ifdef Avoid_Underflow
3158 if (L <= (bc.scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
3161 #endif /*Avoid_Underflow*/
3171 #else /*Sudden_Underflow}{*/
3172 #ifdef Avoid_Underflow
3174 L = word0(&rv) & Exp_mask;
3175 if (L <= (2*P+1)*Exp_msk1) {
3176 if (L > (P+2)*Exp_msk1)
3177 /* round even ==> */
3180 /* rv = smallest denormal */
3188 #endif /*Avoid_Underflow*/
3189 L = (word0(&rv) & Exp_mask) - Exp_msk1;
3190 #endif /*Sudden_Underflow}}*/
3191 word0(&rv) = L | Bndry_mask1;
3192 word1(&rv) = 0xffffffff;
3199 #ifndef ROUND_BIASED
3200 if (!(word1(&rv) & LSB))
3204 dval(&rv) += ulp(&rv);
3205 #ifndef ROUND_BIASED
3207 dval(&rv) -= ulp(&rv);
3208 #ifndef Sudden_Underflow
3218 #ifdef Avoid_Underflow
3219 bc.dsign = 1 - bc.dsign;
3224 if ((aadj = ratio(delta, bs)) <= 2.) {
3227 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
3228 #ifndef Sudden_Underflow
3229 if (word1(&rv) == Tiny1 && !word0(&rv)) {
3241 /* special case -- power of FLT_RADIX to be */
3242 /* rounded down... */
3244 if (aadj < 2./FLT_RADIX)
3245 aadj = 1./FLT_RADIX;
3253 aadj1 = bc.dsign ? aadj : -aadj;
3254 #ifdef Check_FLT_ROUNDS
3255 switch(bc.rounding) {
3256 case 2: /* towards +infinity */
3259 case 0: /* towards 0 */
3260 case 3: /* towards -infinity */
3264 if (Flt_Rounds == 0)
3266 #endif /*Check_FLT_ROUNDS*/
3268 y = word0(&rv) & Exp_mask;
3270 /* Check for overflow */
3272 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
3273 dval(&rv0) = dval(&rv);
3274 word0(&rv) -= P*Exp_msk1;
3275 adj.d = aadj1 * ulp(&rv);
3277 if ((word0(&rv) & Exp_mask) >=
3278 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
3279 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
3286 word0(&rv) += P*Exp_msk1;
3289 #ifdef Avoid_Underflow
3290 if (bc.scale && y <= 2*P*Exp_msk1) {
3291 if (aadj <= 0x7fffffff) {
3292 if ((z = aadj) <= 0)
3295 aadj1 = bc.dsign ? aadj : -aadj;
3297 dval(&aadj2) = aadj1;
3298 word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
3299 aadj1 = dval(&aadj2);
3301 adj.d = aadj1 * ulp(&rv);
3304 #ifdef Sudden_Underflow
3305 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
3306 dval(&rv0) = dval(&rv);
3307 word0(&rv) += P*Exp_msk1;
3308 adj.d = aadj1 * ulp(&rv);
3311 if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
3313 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
3316 if (word0(&rv0) == Tiny0
3317 && word1(&rv0) == Tiny1) {
3329 word0(&rv) -= P*Exp_msk1;
3332 adj.d = aadj1 * ulp(&rv);
3335 #else /*Sudden_Underflow*/
3336 /* Compute adj so that the IEEE rounding rules will
3337 * correctly round rv + adj in some half-way cases.
3338 * If rv * ulp(rv) is denormalized (i.e.,
3339 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
3340 * trouble from bits lost to denormalization;
3341 * example: 1.2e-307 .
3343 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
3344 aadj1 = (double)(int)(aadj + 0.5);
3348 adj.d = aadj1 * ulp(&rv);
3350 #endif /*Sudden_Underflow*/
3351 #endif /*Avoid_Underflow*/
3353 z = word0(&rv) & Exp_mask;
3356 #ifdef Avoid_Underflow
3360 /* Can we stop now? */
3363 /* The tolerances below are conservative. */
3364 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
3365 if (aadj < .4999999 || aadj > .5000001)
3368 else if (aadj < .4999999/FLT_RADIX)
3384 #ifndef NO_STRTOD_BIGCOMP
3386 bigcomp(&rv, s0, &bc);
3391 word0(&rv0) = Exp_1 + (70 << Exp_shift);
3396 else if (!oldinexact)
3399 #ifdef Avoid_Underflow
3401 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
3403 dval(&rv) *= dval(&rv0);
3405 /* try to avoid the bug of testing an 8087 register value */
3407 if (!(word0(&rv) & Exp_mask))
3409 if (word0(&rv) == 0 && word1(&rv) == 0)
3414 #endif /* Avoid_Underflow */
3416 if (bc.inexact && !(word0(&rv) & Exp_mask)) {
3417 /* set underflow bit */
3418 dval(&rv0) = 1e-300;
3419 dval(&rv0) *= dval(&rv0);
3425 return sign ? -dval(&rv) : dval(&rv);
3428 #ifndef MULTIPLE_THREADS
3429 static char *dtoa_result;
3443 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (size_t)i;
3446 r = (int*)Balloc(k);
3449 #ifndef MULTIPLE_THREADS
3457 nrv_alloc(s, rve, n) char *s, **rve; int n;
3459 nrv_alloc(CONST char *s, char **rve, int n)
3464 t = rv = rv_alloc(n);
3465 while((*t = *s++)) t++;
3471 /* freedtoa(s) must be used to free values s returned by dtoa
3472 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
3473 * but for consistency with earlier versions of dtoa, it is optional
3474 * when MULTIPLE_THREADS is not defined.
3479 freedtoa(s) char *s;
3484 Bigint *b = (Bigint *)((int *)s - 1);
3485 b->maxwds = 1 << (b->k = *(int*)b);
3487 #ifndef MULTIPLE_THREADS
3488 if (s == dtoa_result)
3493 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
3495 * Inspired by "How to Print Floating-Point Numbers Accurately" by
3496 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
3499 * 1. Rather than iterating, we use a simple numeric overestimate
3500 * to determine k = floor(log10(d)). We scale relevant
3501 * quantities using O(log2(k)) rather than O(k) multiplications.
3502 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
3503 * try to generate digits strictly left to right. Instead, we
3504 * compute with fewer bits and propagate the carry if necessary
3505 * when rounding the final digit up. This is often faster.
3506 * 3. Under the assumption that input will be rounded nearest,
3507 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
3508 * That is, we allow equality in stopping tests when the
3509 * round-nearest rule will give the same floating-point value
3510 * as would satisfaction of the stopping test with strict
3512 * 4. We remove common factors of powers of 2 from relevant
3514 * 5. When converting floating-point integers less than 1e16,
3515 * we use floating-point arithmetic rather than resorting
3516 * to multiple-precision integers.
3517 * 6. When asked to produce fewer than 15 digits, we first try
3518 * to get by with floating-point arithmetic; we resort to
3519 * multiple-precision integer arithmetic only if we cannot
3520 * guarantee that the floating-point calculation has given
3521 * the correctly rounded result. For k requested digits and
3522 * "uniformly" distributed input, the probability is
3523 * something like 10^(k-15) that we must resort to the Long
3530 (dd, mode, ndigits, decpt, sign, rve)
3531 double dd; int mode, ndigits, *decpt, *sign; char **rve;
3533 (double dd, int mode, int ndigits, int *decpt, int *sign, char **rve)
3536 /* Arguments ndigits, decpt, sign are similar to those
3537 of ecvt and fcvt; trailing zeros are suppressed from
3538 the returned string. If not null, *rve is set to point
3539 to the end of the return value. If d is +-Infinity or NaN,
3540 then *decpt is set to 9999.
3543 0 ==> shortest string that yields d when read in
3544 and rounded to nearest.
3545 1 ==> like 0, but with Steele & White stopping rule;
3546 e.g. with IEEE P754 arithmetic , mode 0 gives
3547 1e23 whereas mode 1 gives 9.999999999999999e22.
3548 2 ==> max(1,ndigits) significant digits. This gives a
3549 return value similar to that of ecvt, except
3550 that trailing zeros are suppressed.
3551 3 ==> through ndigits past the decimal point. This
3552 gives a return value similar to that from fcvt,
3553 except that trailing zeros are suppressed, and
3554 ndigits can be negative.
3555 4,5 ==> similar to 2 and 3, respectively, but (in
3556 round-nearest mode) with the tests of mode 0 to
3557 possibly return a shorter string that rounds to d.
3558 With IEEE arithmetic and compilation with
3559 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
3560 as modes 2 and 3 when FLT_ROUNDS != 1.
3561 6-9 ==> Debugging modes similar to mode - 4: don't try
3562 fast floating-point estimate (if applicable).
3564 Values of mode other than 0-9 are treated as mode 0.
3566 Sufficient space is allocated to the return value
3567 to hold the suppressed trailing zeros.
3570 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
3571 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
3572 spec_case, try_quick;
3574 #ifndef Sudden_Underflow
3578 Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
3583 int inexact, oldinexact;
3585 #ifdef Honor_FLT_ROUNDS /*{*/
3587 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
3588 Rounding = Flt_Rounds;
3591 switch(fegetround()) {
3592 case FE_TOWARDZERO: Rounding = 0; break;
3593 case FE_UPWARD: Rounding = 2; break;
3594 case FE_DOWNWARD: Rounding = 3;
3599 #ifndef MULTIPLE_THREADS
3601 freedtoa(dtoa_result);
3607 if (word0(&u) & Sign_bit) {
3608 /* set sign for everything, including 0's and NaNs */
3610 word0(&u) &= ~Sign_bit; /* clear sign bit */
3615 #if defined(IEEE_Arith) + defined(VAX)
3617 if ((word0(&u) & Exp_mask) == Exp_mask)
3619 if (word0(&u) == 0x8000)
3622 /* Infinity or NaN */
3625 if (!word1(&u) && !(word0(&u) & 0xfffff))
3626 return nrv_alloc("Infinity", rve, 8);
3628 return nrv_alloc("NaN", rve, 3);
3632 dval(&u) += 0; /* normalize */
3636 return nrv_alloc("0", rve, 1);
3640 try_quick = oldinexact = get_inexact();
3643 #ifdef Honor_FLT_ROUNDS
3644 if (Rounding >= 2) {
3646 Rounding = Rounding == 2 ? 0 : 2;
3653 b = d2b(&u, &be, &bbits);
3654 #ifdef Sudden_Underflow
3655 i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
3657 if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
3659 dval(&d2) = dval(&u);
3660 word0(&d2) &= Frac_mask1;
3661 word0(&d2) |= Exp_11;
3663 if (j = 11 - hi0bits(word0(&d2) & Frac_mask))
3664 dval(&d2) /= 1 << j;
3667 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
3668 * log10(x) = log(x) / log(10)
3669 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
3670 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
3672 * This suggests computing an approximation k to log10(d) by
3674 * k = (i - Bias)*0.301029995663981
3675 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
3677 * We want k to be too large rather than too small.
3678 * The error in the first-order Taylor series approximation
3679 * is in our favor, so we just round up the constant enough
3680 * to compensate for any error in the multiplication of
3681 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
3682 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
3683 * adding 1e-13 to the constant term more than suffices.
3684 * Hence we adjust the constant term to 0.1760912590558.
3685 * (We could get a more accurate k by invoking log10,
3686 * but this is probably not worthwhile.)
3694 #ifndef Sudden_Underflow
3698 /* d is denormalized */
3700 i = bbits + be + (Bias + (P-1) - 1);
3701 x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
3702 : word1(&u) << (32 - i);
3704 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
3705 i -= (Bias + (P-1) - 1) + 1;
3709 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
3711 if (ds < 0. && ds != k)
3712 k--; /* want k = floor(ds) */
3714 if (k >= 0 && k <= Ten_pmax) {
3715 if (dval(&u) < tens[k])
3738 if (mode < 0 || mode > 9)
3742 #ifdef Check_FLT_ROUNDS
3743 try_quick = Rounding == 1;
3747 #endif /*SET_INEXACT*/
3754 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
3755 /* silence erroneous "gcc -Wall" warning. */
3768 ilim = ilim1 = i = ndigits;
3774 i = ndigits + k + 1;
3780 s = s0 = rv_alloc(i);
3782 #ifdef Honor_FLT_ROUNDS
3783 if (mode > 1 && Rounding != 1)
3787 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
3789 /* Try to get by with floating-point arithmetic. */
3792 dval(&d2) = dval(&u);
3795 ieps = 2; /* conservative */
3800 /* prevent overflows */
3802 dval(&u) /= bigtens[n_bigtens-1];
3805 for(; j; j >>= 1, i++)
3812 else if ((j1 = -k)) {
3813 dval(&u) *= tens[j1 & 0xf];
3814 for(j = j1 >> 4; j; j >>= 1, i++)
3817 dval(&u) *= bigtens[i];
3820 if (k_check && dval(&u) < 1. && ilim > 0) {
3828 dval(&eps) = ieps*dval(&u) + 7.;
3829 word0(&eps) -= (P-1)*Exp_msk1;
3833 if (dval(&u) > dval(&eps))
3835 if (dval(&u) < -dval(&eps))
3839 #ifndef No_leftright
3841 /* Use Steele & White method of only
3842 * generating digits needed.
3844 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
3848 *s++ = '0' + (int)L;
3849 if (dval(&u) < dval(&eps))
3851 if (1. - dval(&u) < dval(&eps))
3861 /* Generate ilim digits, then fix them up. */
3862 dval(&eps) *= tens[ilim-1];
3863 for(i = 1;; i++, dval(&u) *= 10.) {
3864 L = (Long)(dval(&u));
3865 if (!(dval(&u) -= L))
3867 *s++ = '0' + (int)L;
3869 if (dval(&u) > 0.5 + dval(&eps))
3871 else if (dval(&u) < 0.5 - dval(&eps)) {
3872 while(*--s == '0') {}
3879 #ifndef No_leftright
3884 dval(&u) = dval(&d2);
3889 /* Do we have a "small" integer? */
3891 if (be >= 0 && k <= Int_max) {
3894 if (ndigits < 0 && ilim <= 0) {
3896 if (ilim < 0 || dval(&u) <= 5*ds)
3900 for(i = 1; i <= k + 1; i++, dval(&u) *= 10.) {
3901 L = (Long)(dval(&u) / ds);
3903 #ifdef Check_FLT_ROUNDS
3904 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
3910 *s++ = '0' + (int)L;
3918 #ifdef Honor_FLT_ROUNDS
3922 case 2: goto bump_up;
3925 dval(&u) += dval(&u);
3926 if (dval(&u) > ds || (dval(&u) == ds && L & 1)) {
3947 #ifndef Sudden_Underflow
3948 denorm ? be + (Bias + (P-1) - 1 + 1) :
3951 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
3959 if (m2 > 0 && s2 > 0) {
3960 i = m2 < s2 ? m2 : s2;
3968 mhi = pow5mult(mhi, m5);
3977 b = pow5mult(b, b5);
3981 S = pow5mult(S, s5);
3983 /* Check for special case that d is a normalized power of 2. */
3986 if ((mode < 2 || leftright)
3987 #ifdef Honor_FLT_ROUNDS
3991 if (!word1(&u) && !(word0(&u) & Bndry_mask)
3992 #ifndef Sudden_Underflow
3993 && word0(&u) & (Exp_mask & ~Exp_msk1)
3996 /* The special case */
4003 /* Arrange for convenient computation of quotients:
4004 * shift left if necessary so divisor has 4 leading 0 bits.
4006 * Perhaps we should just compute leading 28 bits of S once
4007 * and for all and pass them and a shift to quorem, so it
4008 * can do shifts and ors to compute the numerator for q.
4011 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
4015 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
4030 b = multadd(b, 10, 0); /* we botched the k estimate */
4032 mhi = multadd(mhi, 10, 0);
4036 if (ilim <= 0 && (mode == 3 || mode == 5)) {
4037 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
4038 /* no digits, fcvt style */
4050 mhi = lshift(mhi, m2);
4052 /* Compute mlo -- check for special case
4053 * that d is a normalized power of 2.
4058 mhi = Balloc(mhi->k);
4060 mhi = lshift(mhi, Log2P);
4064 dig = quorem(b,S) + '0';
4065 /* Do we yet have the shortest decimal string
4066 * that will round to d?
4069 delta = diff(S, mhi);
4070 j1 = delta->sign ? 1 : cmp(b, delta);
4072 #ifndef ROUND_BIASED
4073 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
4074 #ifdef Honor_FLT_ROUNDS
4083 else if (!b->x[0] && b->wds <= 1)
4090 if (j < 0 || (j == 0 && mode != 1
4091 #ifndef ROUND_BIASED
4095 if (!b->x[0] && b->wds <= 1) {
4101 #ifdef Honor_FLT_ROUNDS
4104 case 0: goto accept_dig;
4105 case 2: goto keep_dig;
4107 #endif /*Honor_FLT_ROUNDS*/
4111 if ((j1 > 0 || (j1 == 0 && dig & 1))
4120 #ifdef Honor_FLT_ROUNDS
4124 if (dig == '9') { /* possible if i == 1 */
4132 #ifdef Honor_FLT_ROUNDS
4138 b = multadd(b, 10, 0);
4140 mlo = mhi = multadd(mhi, 10, 0);
4142 mlo = multadd(mlo, 10, 0);
4143 mhi = multadd(mhi, 10, 0);
4149 *s++ = dig = quorem(b,S) + '0';
4150 if (!b->x[0] && b->wds <= 1) {
4158 b = multadd(b, 10, 0);
4161 /* Round off last digit */
4163 #ifdef Honor_FLT_ROUNDS
4165 case 0: goto trimzeros;
4166 case 2: goto roundoff;
4171 if (j > 0 || (j == 0 && dig & 1)) {
4182 #ifdef Honor_FLT_ROUNDS
4185 while(*--s == '0') {}
4191 if (mlo && mlo != mhi)
4199 word0(&u) = Exp_1 + (70 << Exp_shift);
4204 else if (!oldinexact)
4215 } // namespace dmg_fp