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 ***************************************************************/
22 #define freedtoa __freedtoa
25 #define Omit_Private_Memory
26 #define MULTIPLE_THREADS 1
27 /* Lock 0 is not used because of USE_MALLOC, Lock 1 protects a lazy-initialized table */
28 #define ACQUIRE_DTOA_LOCK(n)
29 #define FREE_DTOA_LOCK(n)
31 /* Please send bug reports to David M. Gay (dmg at acm dot org,
32 * with " at " changed at "@" and " dot " changed to "."). */
34 /* On a machine with IEEE extended-precision registers, it is
35 * necessary to specify double-precision (53-bit) rounding precision
36 * before invoking strtod or dtoa. If the machine uses (the equivalent
37 * of) Intel 80x87 arithmetic, the call
38 * _control87(PC_53, MCW_PC);
39 * does this with many compilers. Whether this or another call is
40 * appropriate depends on the compiler; for this to work, it may be
41 * necessary to #include "float.h" or another system-dependent header
45 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
47 * This strtod returns a nearest machine number to the input decimal
48 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
49 * broken by the IEEE round-even rule. Otherwise ties are broken by
50 * biased rounding (add half and chop).
52 * Inspired loosely by William D. Clinger's paper "How to Read Floating
53 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
57 * 1. We only require IEEE, IBM, or VAX double-precision
58 * arithmetic (not IEEE double-extended).
59 * 2. We get by with floating-point arithmetic in a case that
60 * Clinger missed -- when we're computing d * 10^n
61 * for a small integer d and the integer n is not too
62 * much larger than 22 (the maximum integer k for which
63 * we can represent 10^k exactly), we may be able to
64 * compute (d*10^k) * 10^(e-k) with just one roundoff.
65 * 3. Rather than a bit-at-a-time adjustment of the binary
66 * result in the hard case, we use floating-point
67 * arithmetic to determine the adjustment to within
68 * one bit; only in really hard cases do we need to
69 * compute a second residual.
70 * 4. Because of 3., we don't need a large table of powers of 10
71 * for ten-to-e (just some small tables, e.g. of 10^k
76 * #define IEEE_8087 for IEEE-arithmetic machines where the least
77 * significant byte has the lowest address.
78 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
79 * significant byte has the lowest address.
80 * #define Long int on machines with 32-bit ints and 64-bit longs.
81 * #define IBM for IBM mainframe-style floating-point arithmetic.
82 * #define VAX for VAX-style floating-point arithmetic (D_floating).
83 * #define No_leftright to omit left-right logic in fast floating-point
84 * computation of dtoa.
85 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
86 * and strtod and dtoa should round accordingly.
87 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
88 * and Honor_FLT_ROUNDS is not #defined.
89 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
90 * that use extended-precision instructions to compute rounded
91 * products and quotients) with IBM.
92 * #define ROUND_BIASED for IEEE-format with biased rounding.
93 * #define Inaccurate_Divide for IEEE-format with correctly rounded
94 * products but inaccurate quotients, e.g., for Intel i860.
95 * #define NO_LONG_LONG on machines that do not have a "long long"
96 * integer type (of >= 64 bits). On such machines, you can
97 * #define Just_16 to store 16 bits per 32-bit Long when doing
98 * high-precision integer arithmetic. Whether this speeds things
99 * up or slows things down depends on the machine and the number
100 * being converted. If long long is available and the name is
101 * something other than "long long", #define Llong to be the name,
102 * and if "unsigned Llong" does not work as an unsigned version of
103 * Llong, #define #ULLong to be the corresponding unsigned type.
104 * #define KR_headers for old-style C function headers.
105 * #define Bad_float_h if your system lacks a float.h or if it does not
106 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
107 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
108 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
109 * if memory is available and otherwise does something you deem
110 * appropriate. If MALLOC is undefined, malloc will be invoked
111 * directly -- and assumed always to succeed.
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 INFNAN_CHECK on IEEE systems to cause strtod to check for
124 * Infinity and NaN (case insensitively). On some systems (e.g.,
125 * some HP systems), it may be necessary to #define NAN_WORD0
126 * appropriately -- to the most significant word of a quiet NaN.
127 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
128 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
129 * strtod also accepts (case insensitively) strings of the form
130 * NaN(x), where x is a string of hexadecimal digits and spaces;
131 * if there is only one string of hexadecimal digits, it is taken
132 * for the 52 fraction bits of the resulting NaN; if there are two
133 * or more strings of hex digits, the first is for the high 20 bits,
134 * the second and subsequent for the low 32 bits, with intervening
135 * white space ignored; but if this results in none of the 52
136 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
137 * and NAN_WORD1 are used instead.
138 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
139 * multiple threads. In this case, you must provide (or suitably
140 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
141 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
142 * in pow5mult, ensures lazy evaluation of only one copy of high
143 * powers of 5; omitting this lock would introduce a small
144 * probability of wasting memory, but would otherwise be harmless.)
145 * You must also invoke freedtoa(s) to free the value s returned by
146 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
147 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
148 * avoids underflows on inputs whose result does not underflow.
149 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
150 * floating-point numbers and flushes underflows to zero rather
151 * than implementing gradual underflow, then you must also #define
153 * #define YES_ALIAS to permit aliasing certain double values with
154 * arrays of ULongs. This leads to slightly better code with
155 * some compilers and was always used prior to 19990916, but it
156 * is not strictly legal and can cause trouble with aggressively
157 * optimizing compilers (e.g., gcc 2.95.1 under -O2).
158 * #define USE_LOCALE to use the current locale's decimal_point value.
159 * #define SET_INEXACT if IEEE arithmetic is being used and extra
160 * computation should be done to set the inexact flag when the
161 * result is inexact and avoid setting inexact when the result
162 * is exact. In this case, dtoa.c must be compiled in
163 * an environment, perhaps provided by #include "dtoa.c" in a
164 * suitable wrapper, that defines two functions,
165 * int get_inexact(void);
166 * void clear_inexact(void);
167 * such that get_inexact() returns a nonzero value if the
168 * inexact bit is already set, and clear_inexact() sets the
169 * inexact bit to 0. When SET_INEXACT is #defined, strtod
170 * also does extra computations to set the underflow and overflow
171 * flags when appropriate (i.e., when the result is tiny and
172 * inexact or when it is a numeric value rounded to +-infinity).
173 * #define NO_ERRNO if strtod should not assign errno = ERANGE when
174 * the result overflows to +-Infinity or underflows to 0.
176 #if defined(TARGET_X86) || defined(mips) && defined(MIPSEL) || defined (__arm__)
180 #elif defined(TARGET_AMD64) || defined(__alpha__)
184 #elif defined(__ia64)
192 #elif defined(__hppa)
201 #define ULong guint32
205 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
218 extern char *MALLOC();
220 extern void *MALLOC(size_t);
223 #define MALLOC malloc
226 #define Omit_Private_Memory
227 #ifndef Omit_Private_Memory
229 #define PRIVATE_MEM 2304
231 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
232 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
236 #undef Avoid_Underflow
250 #define DBL_MAX_10_EXP 308
251 #define DBL_MAX_EXP 1024
253 #endif /*IEEE_Arith*/
257 #define DBL_MAX_10_EXP 75
258 #define DBL_MAX_EXP 63
260 #define DBL_MAX 7.2370055773322621e+75
265 #define DBL_MAX_10_EXP 38
266 #define DBL_MAX_EXP 127
268 #define DBL_MAX 1.7014118346046923e+38
272 #define LONG_MAX 2147483647
275 #else /* ifndef Bad_float_h */
277 #endif /* Bad_float_h */
289 #define CONST /* blank */
295 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
296 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
299 typedef union { double d; ULong L[2]; } U;
304 #define word0(x) ((ULong *)&x)[1]
305 #define word1(x) ((ULong *)&x)[0]
307 #define word0(x) ((ULong *)&x)[0]
308 #define word1(x) ((ULong *)&x)[1]
312 #define word0(x) ((U*)&x)->L[1]
313 #define word1(x) ((U*)&x)->L[0]
315 #define word0(x) ((U*)&x)->L[0]
316 #define word1(x) ((U*)&x)->L[1]
318 #define dval(x) ((U*)&x)->d
321 /* The following definition of Storeinc is appropriate for MIPS processors.
322 * An alternative that might be better on some machines is
323 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
325 #if defined(IEEE_8087) + defined(VAX)
326 #define Storeinc(a,b,c) do { (((unsigned short *)a)[1] = (unsigned short)b, \
327 ((unsigned short *)a)[0] = (unsigned short)c, a++) } while (0)
329 #define Storeinc(a,b,c) do { (((unsigned short *)a)[0] = (unsigned short)b, \
330 ((unsigned short *)a)[1] = (unsigned short)c, a++) } while (0)
333 /* #define P DBL_MANT_DIG */
334 /* Ten_pmax = floor(P*log(2)/log(5)) */
335 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
336 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
337 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
341 #define Exp_shift1 20
342 #define Exp_msk1 0x100000
343 #define Exp_msk11 0x100000
344 #define Exp_mask 0x7ff00000
348 #define Exp_1 0x3ff00000
349 #define Exp_11 0x3ff00000
351 #define Frac_mask 0xfffff
352 #define Frac_mask1 0xfffff
355 #define Bndry_mask 0xfffff
356 #define Bndry_mask1 0xfffff
358 #define Sign_bit 0x80000000
364 #ifndef NO_IEEE_Scale
365 #define Avoid_Underflow
366 #ifdef Flush_Denorm /* debugging option */
367 #undef Sudden_Underflow
373 #define Flt_Rounds FLT_ROUNDS
377 #endif /*Flt_Rounds*/
379 #ifdef Honor_FLT_ROUNDS
380 #define Rounding rounding
381 #undef Check_FLT_ROUNDS
382 #define Check_FLT_ROUNDS
384 #define Rounding Flt_Rounds
387 #else /* ifndef IEEE_Arith */
388 #undef Check_FLT_ROUNDS
389 #undef Honor_FLT_ROUNDS
391 #undef Sudden_Underflow
392 #define Sudden_Underflow
397 #define Exp_shift1 24
398 #define Exp_msk1 0x1000000
399 #define Exp_msk11 0x1000000
400 #define Exp_mask 0x7f000000
403 #define Exp_1 0x41000000
404 #define Exp_11 0x41000000
405 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
406 #define Frac_mask 0xffffff
407 #define Frac_mask1 0xffffff
410 #define Bndry_mask 0xefffff
411 #define Bndry_mask1 0xffffff
413 #define Sign_bit 0x80000000
415 #define Tiny0 0x100000
424 #define Exp_msk1 0x80
425 #define Exp_msk11 0x800000
426 #define Exp_mask 0x7f80
429 #define Exp_1 0x40800000
430 #define Exp_11 0x4080
432 #define Frac_mask 0x7fffff
433 #define Frac_mask1 0xffff007f
436 #define Bndry_mask 0xffff007f
437 #define Bndry_mask1 0xffff007f
439 #define Sign_bit 0x8000
445 #endif /* IBM, VAX */
446 #endif /* IEEE_Arith */
453 #define rounded_product(a,b) a = rnd_prod(a, b)
454 #define rounded_quotient(a,b) a = rnd_quot(a, b)
456 extern double rnd_prod(), rnd_quot();
458 extern double rnd_prod(double, double), rnd_quot(double, double);
461 #define rounded_product(a,b) a *= b
462 #define rounded_quotient(a,b) a /= b
465 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
466 #define Big1 0xffffffff
473 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
475 #define FFFFFFFF 0xffffffffUL
482 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
483 * This makes some inner loops simpler and sometimes saves work
484 * during multiplications, but it often seems to make things slightly
485 * slower. Hence the default is now to store 32 bits per Long.
488 #else /* long long available */
490 #define Llong long long
493 #define ULLong unsigned Llong
495 #endif /* NO_LONG_LONG */
497 #ifndef MULTIPLE_THREADS
498 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
499 #define FREE_DTOA_LOCK(n) /*nothing*/
505 extern "C" double strtod(const char *s00, char **se);
506 extern "C" char *dtoa(double d, int mode, int ndigits,
507 int *decpt, int *sign, char **rve);
513 int k, maxwds, sign, wds;
517 typedef struct Bigint Bigint;
519 static Bigint *freelist[Kmax+1];
531 #ifndef Omit_Private_Memory
535 ACQUIRE_DTOA_LOCK(0);
536 if ((rv = freelist[k])) {
537 freelist[k] = rv->next;
541 #ifdef Omit_Private_Memory
542 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
544 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
546 if (pmem_next - private_mem + len <= PRIVATE_mem) {
547 rv = (Bigint*)pmem_next;
551 rv = (Bigint*)MALLOC(len*sizeof(double));
557 rv->sign = rv->wds = 0;
569 #ifdef Omit_Private_Memory
573 ACQUIRE_DTOA_LOCK(0);
574 v->next = freelist[v->k];
581 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
582 y->wds*sizeof(Long) + 2*sizeof(int))
587 (b, m, a) Bigint *b; int m, a;
589 (Bigint *b, int m, int a) /* multiply by m and add a */
610 y = *x * (ULLong)m + carry;
616 y = (xi & 0xffff) * m + carry;
617 z = (xi >> 16) * m + (y >> 16);
619 *x++ = (z << 16) + (y & 0xffff);
629 if (wds >= b->maxwds) {
644 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
646 (CONST char *s, int nd0, int nd, ULong y9)
654 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
661 b->x[0] = y9 & 0xffff;
662 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
668 do b = multadd(b, 10, *s++ - '0');
675 b = multadd(b, 10, *s++ - '0');
682 (x) register ULong x;
689 if (!(x & 0xffff0000)) {
693 if (!(x & 0xff000000)) {
697 if (!(x & 0xf0000000)) {
701 if (!(x & 0xc0000000)) {
705 if (!(x & 0x80000000)) {
707 if (!(x & 0x40000000))
722 register ULong x = *y;
780 (a, b) Bigint *a, *b;
782 (Bigint *a, Bigint *b)
787 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
798 if (a->wds < b->wds) {
810 for(x = c->x, xa = x + wc; x < xa; x++)
818 for(; xb < xbe; xc0++) {
824 z = *x++ * (ULLong)y + *xc + carry;
826 *xc++ = z & FFFFFFFF;
834 for(; xb < xbe; xb++, xc0++) {
835 if (y = *xb & 0xffff) {
840 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
842 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
855 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
858 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
866 for(; xb < xbe; xc0++) {
872 z = *x++ * y + *xc + carry;
882 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
892 (b, k) Bigint *b; int k;
897 Bigint *b1, *p5, *p51;
899 static int p05[3] = { 5, 25, 125 };
902 b = multadd(b, p05[i-1], 0);
908 #ifdef MULTIPLE_THREADS
909 ACQUIRE_DTOA_LOCK(1);
928 if (!(p51 = p5->next)) {
929 #ifdef MULTIPLE_THREADS
930 ACQUIRE_DTOA_LOCK(1);
931 if (!(p51 = p5->next)) {
932 p51 = p5->next = mult(p5,p5);
937 p51 = p5->next = mult(p5,p5);
949 (b, k) Bigint *b; int k;
956 ULong *x, *x1, *xe, z;
965 for(i = b->maxwds; n1 > i; i <<= 1)
969 for(i = 0; i < n; i++)
990 *x1++ = *x << k & 0xffff | z;
1009 (a, b) Bigint *a, *b;
1011 (Bigint *a, Bigint *b)
1014 ULong *xa, *xa0, *xb, *xb0;
1020 if (i > 1 && !a->x[i-1])
1021 Bug("cmp called with a->x[a->wds-1] == 0");
1022 if (j > 1 && !b->x[j-1])
1023 Bug("cmp called with b->x[b->wds-1] == 0");
1033 return *xa < *xb ? -1 : 1;
1043 (a, b) Bigint *a, *b;
1045 (Bigint *a, Bigint *b)
1050 ULong *xa, *xae, *xb, *xbe, *xc;
1087 y = (ULLong)*xa++ - *xb++ - borrow;
1088 borrow = y >> 32 & (ULong)1;
1089 *xc++ = y & FFFFFFFF;
1094 borrow = y >> 32 & (ULong)1;
1095 *xc++ = y & FFFFFFFF;
1100 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1101 borrow = (y & 0x10000) >> 16;
1102 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1103 borrow = (z & 0x10000) >> 16;
1108 y = (*xa & 0xffff) - borrow;
1109 borrow = (y & 0x10000) >> 16;
1110 z = (*xa++ >> 16) - borrow;
1111 borrow = (z & 0x10000) >> 16;
1116 y = *xa++ - *xb++ - borrow;
1117 borrow = (y & 0x10000) >> 16;
1123 borrow = (y & 0x10000) >> 16;
1145 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1146 #ifndef Avoid_Underflow
1147 #ifndef Sudden_Underflow
1156 #ifndef Avoid_Underflow
1157 #ifndef Sudden_Underflow
1160 L = -L >> Exp_shift;
1161 if (L < Exp_shift) {
1162 word0(a) = 0x80000 >> L;
1168 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
1179 (a, e) Bigint *a; int *e;
1184 ULong *xa, *xa0, w, y, z;
1198 if (!y) Bug("zero y in b2d");
1204 d0 = Exp_1 | (y >> (Ebits - k));
1205 w = xa > xa0 ? *--xa : 0;
1206 d1 = y << ((32-Ebits) + k) | (w >> (Ebits - k));
1209 z = xa > xa0 ? *--xa : 0;
1211 d0 = Exp_1 | y << k | (z >> (32 - k));
1212 y = xa > xa0 ? *--xa : 0;
1213 d1 = z << k | (y >> (32 - k));
1220 if (k < Ebits + 16) {
1221 z = xa > xa0 ? *--xa : 0;
1222 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1223 w = xa > xa0 ? *--xa : 0;
1224 y = xa > xa0 ? *--xa : 0;
1225 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1228 z = xa > xa0 ? *--xa : 0;
1229 w = xa > xa0 ? *--xa : 0;
1231 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1232 y = xa > xa0 ? *--xa : 0;
1233 d1 = w << k + 16 | y << k;
1237 word0(d) = d0 >> 16 | d0 << 16;
1238 word1(d) = d1 >> 16 | d1 << 16;
1249 (d, e, bits) double d; int *e, *bits;
1251 (double d, int *e, int *bits)
1257 #ifndef Sudden_Underflow
1262 d0 = word0(d) >> 16 | word0(d) << 16;
1263 d1 = word1(d) >> 16 | word1(d) << 16;
1277 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1278 #ifdef Sudden_Underflow
1279 de = (int)(d0 >> Exp_shift);
1284 if ((de = (int)(d0 >> Exp_shift)))
1289 if ((k = lo0bits(&y))) {
1290 x[0] = y | (z << (32 - k));
1295 #ifndef Sudden_Underflow
1298 b->wds = (x[1] = z) ? 2 : 1;
1303 Bug("Zero passed to d2b");
1307 #ifndef Sudden_Underflow
1315 if (k = lo0bits(&y))
1317 x[0] = y | z << 32 - k & 0xffff;
1318 x[1] = z >> k - 16 & 0xffff;
1324 x[1] = y >> 16 | z << 16 - k & 0xffff;
1325 x[2] = z >> k & 0xffff;
1340 Bug("Zero passed to d2b");
1358 #ifndef Sudden_Underflow
1362 *e = (de - Bias - (P-1) << 2) + k;
1363 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1365 *e = de - Bias - (P-1) + k;
1368 #ifndef Sudden_Underflow
1371 *e = de - Bias - (P-1) + 1 + k;
1373 *bits = 32*i - hi0bits(x[i-1]);
1375 *bits = (i+2)*16 - hi0bits(x[i]);
1387 (a, b) Bigint *a, *b;
1389 (Bigint *a, Bigint *b)
1395 dval(da) = b2d(a, &ka);
1396 dval(db) = b2d(b, &kb);
1398 k = ka - kb + 32*(a->wds - b->wds);
1400 k = ka - kb + 16*(a->wds - b->wds);
1404 word0(da) += (k >> 2)*Exp_msk1;
1410 word0(db) += (k >> 2)*Exp_msk1;
1416 word0(da) += k*Exp_msk1;
1419 word0(db) += k*Exp_msk1;
1422 return dval(da) / dval(db);
1427 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1428 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1437 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1438 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1439 #ifdef Avoid_Underflow
1440 9007199254740992.*9007199254740992.e-256
1441 /* = 2^106 * 1e-53 */
1446 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1447 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1448 #define Scale_Bit 0x10
1452 bigtens[] = { 1e16, 1e32, 1e64 };
1453 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1456 bigtens[] = { 1e16, 1e32 };
1457 static CONST double tinytens[] = { 1e-16, 1e-32 };
1469 #define NAN_WORD0 0x7ff80000
1479 (sp, t) char **sp, *t;
1481 (CONST char **sp, char *t)
1485 CONST char *s = *sp;
1488 if ((c = *++s) >= 'A' && c <= 'Z')
1501 (rvp, sp) double *rvp; CONST char **sp;
1503 (double *rvp, CONST char **sp)
1508 int havedig, udx0, xshift;
1511 havedig = xshift = 0;
1514 while(c = *(CONST unsigned char*)++s) {
1515 if (c >= '0' && c <= '9')
1517 else if (c >= 'a' && c <= 'f')
1519 else if (c >= 'A' && c <= 'F')
1521 else if (c <= ' ') {
1522 if (udx0 && havedig) {
1528 else if (/*(*/ c == ')' && havedig) {
1533 return; /* invalid form: don't change *sp */
1541 x[0] = (x[0] << 4) | (x[1] >> 28);
1542 x[1] = (x[1] << 4) | c;
1544 if ((x[0] &= 0xfffff) || x[1]) {
1545 word0(*rvp) = Exp_mask | x[0];
1549 #endif /*No_Hex_NaN*/
1550 #endif /* INFNAN_CHECK */
1553 * LOCKING: This is not thread-safe, since the locking macros are defined as no-ops,
1554 * the caller should lock.
1560 (s00, se) CONST char *s00; char **se;
1562 (CONST char *s00, char **se)
1565 #ifdef Avoid_Underflow
1568 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1569 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1570 CONST char *s, *s0, *s1;
1571 double aadj, aadj1, adj, rv, rv0;
1574 Bigint *bb = NULL, *bb1, *bd = NULL, *bd0, *bs = NULL, *delta = NULL;
1576 int inexact, oldinexact;
1578 #ifdef Honor_FLT_ROUNDS
1585 sign = nz0 = nz = 0;
1587 for(s = s00;;s++) switch(*s) {
1610 while(*++s == '0') ;
1616 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1623 s1 = localeconv()->decimal_point;
1644 for(; c == '0'; c = *++s)
1646 if (c > '0' && c <= '9') {
1654 for(; c >= '0' && c <= '9'; c = *++s) {
1659 for(i = 1; i < nz; i++)
1662 else if (nd <= DBL_DIG + 1)
1666 else if (nd <= DBL_DIG + 1)
1674 if (c == 'e' || c == 'E') {
1675 if (!nd && !nz && !nz0) {
1686 if (c >= '0' && c <= '9') {
1689 if (c > '0' && c <= '9') {
1692 while((c = *++s) >= '0' && c <= '9')
1694 if (s - s1 > 8 || L > 19999)
1695 /* Avoid confusion from exponents
1696 * so large that e might overflow.
1698 e = 19999; /* safe for 16 bit ints */
1713 /* Check for Nan and Infinity */
1717 if (match(&s,"nf")) {
1719 if (!match(&s,"inity"))
1721 word0(rv) = 0x7ff00000;
1728 if (match(&s, "an")) {
1729 word0(rv) = NAN_WORD0;
1730 word1(rv) = NAN_WORD1;
1732 if (*s == '(') /*)*/
1738 #endif /* INFNAN_CHECK */
1747 /* Now we have nd0 digits, starting at s0, followed by a
1748 * decimal point, followed by nd-nd0 digits. The number we're
1749 * after is the integer represented by those digits times
1754 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1759 oldinexact = get_inexact();
1761 dval(rv) = tens[k - 9] * dval(rv) + z;
1765 #ifndef RND_PRODQUOT
1766 #ifndef Honor_FLT_ROUNDS
1774 if (e <= Ten_pmax) {
1776 goto vax_ovfl_check;
1778 #ifdef Honor_FLT_ROUNDS
1779 /* round correctly FLT_ROUNDS = 2 or 3 */
1785 /* rv = */ rounded_product(dval(rv), tens[e]);
1790 if (e <= Ten_pmax + i) {
1791 /* A fancier test would sometimes let us do
1792 * this for larger i values.
1794 #ifdef Honor_FLT_ROUNDS
1795 /* round correctly FLT_ROUNDS = 2 or 3 */
1802 dval(rv) *= tens[i];
1804 /* VAX exponent range is so narrow we must
1805 * worry about overflow here...
1808 word0(rv) -= P*Exp_msk1;
1809 /* rv = */ rounded_product(dval(rv), tens[e]);
1810 if ((word0(rv) & Exp_mask)
1811 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1813 word0(rv) += P*Exp_msk1;
1815 /* rv = */ rounded_product(dval(rv), tens[e]);
1820 #ifndef Inaccurate_Divide
1821 else if (e >= -Ten_pmax) {
1822 #ifdef Honor_FLT_ROUNDS
1823 /* round correctly FLT_ROUNDS = 2 or 3 */
1829 /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1840 oldinexact = get_inexact();
1842 #ifdef Avoid_Underflow
1845 #ifdef Honor_FLT_ROUNDS
1846 if ((rounding = Flt_Rounds) >= 2) {
1848 rounding = rounding == 2 ? 0 : 2;
1854 #endif /*IEEE_Arith*/
1856 /* Get starting approximation = rv * 10**e1 */
1860 dval(rv) *= tens[i];
1862 if (e1 > DBL_MAX_10_EXP) {
1867 /* Can't trust HUGE_VAL */
1869 #ifdef Honor_FLT_ROUNDS
1871 case 0: /* toward 0 */
1872 case 3: /* toward -infinity */
1877 word0(rv) = Exp_mask;
1880 #else /*Honor_FLT_ROUNDS*/
1881 word0(rv) = Exp_mask;
1883 #endif /*Honor_FLT_ROUNDS*/
1885 /* set overflow bit */
1887 dval(rv0) *= dval(rv0);
1889 #else /*IEEE_Arith*/
1892 #endif /*IEEE_Arith*/
1898 for(j = 0; e1 > 1; j++, e1 >>= 1)
1900 dval(rv) *= bigtens[j];
1901 /* The last multiplication could overflow. */
1902 word0(rv) -= P*Exp_msk1;
1903 dval(rv) *= bigtens[j];
1904 if ((z = word0(rv) & Exp_mask)
1905 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1907 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1908 /* set to largest number */
1909 /* (Can't trust DBL_MAX) */
1914 word0(rv) += P*Exp_msk1;
1920 dval(rv) /= tens[i];
1922 if (e1 >= 1 << n_bigtens)
1924 #ifdef Avoid_Underflow
1927 for(j = 0; e1 > 0; j++, e1 >>= 1)
1929 dval(rv) *= tinytens[j];
1930 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
1931 >> Exp_shift)) > 0) {
1932 /* scaled rv is denormal; zap j low bits */
1936 word0(rv) = (P+2)*Exp_msk1;
1938 word0(rv) &= 0xffffffff << (j-32);
1941 word1(rv) &= 0xffffffff << j;
1944 for(j = 0; e1 > 1; j++, e1 >>= 1)
1946 dval(rv) *= tinytens[j];
1947 /* The last multiplication could underflow. */
1948 dval(rv0) = dval(rv);
1949 dval(rv) *= tinytens[j];
1951 dval(rv) = 2.*dval(rv0);
1952 dval(rv) *= tinytens[j];
1964 #ifndef Avoid_Underflow
1967 /* The refinement below will clean
1968 * this approximation up.
1975 /* Now the hard part -- adjusting rv to the correct value.*/
1977 /* Put digits into bd: true value = bd * 10^e */
1979 bd0 = s2b(s0, nd0, nd, y);
1982 bd = Balloc(bd0->k);
1984 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
2000 #ifdef Honor_FLT_ROUNDS
2004 #ifdef Avoid_Underflow
2006 i = j + bbbits - 1; /* logb(rv) */
2007 if (i < Emin) /* denormal */
2011 #else /*Avoid_Underflow*/
2012 #ifdef Sudden_Underflow
2014 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
2018 #else /*Sudden_Underflow*/
2020 i = j + bbbits - 1; /* logb(rv) */
2021 if (i < Emin) /* denormal */
2025 #endif /*Sudden_Underflow*/
2026 #endif /*Avoid_Underflow*/
2029 #ifdef Avoid_Underflow
2032 i = bb2 < bd2 ? bb2 : bd2;
2041 bs = pow5mult(bs, bb5);
2047 bb = lshift(bb, bb2);
2049 bd = pow5mult(bd, bd5);
2051 bd = lshift(bd, bd2);
2053 bs = lshift(bs, bs2);
2054 delta = diff(bb, bd);
2055 dsign = delta->sign;
2058 #ifdef Honor_FLT_ROUNDS
2059 if (rounding != 1) {
2061 /* Error is less than an ulp */
2062 if (!delta->x[0] && delta->wds <= 1) {
2078 && !(word0(rv) & Frac_mask)) {
2079 y = word0(rv) & Exp_mask;
2080 #ifdef Avoid_Underflow
2081 if (!scale || y > 2*P*Exp_msk1)
2086 delta = lshift(delta,Log2P);
2087 if (cmp(delta, bs) <= 0)
2092 #ifdef Avoid_Underflow
2093 if (scale && (y = word0(rv) & Exp_mask)
2095 word0(adj) += (2*P+1)*Exp_msk1 - y;
2097 #ifdef Sudden_Underflow
2098 if ((word0(rv) & Exp_mask) <=
2100 word0(rv) += P*Exp_msk1;
2101 dval(rv) += adj*ulp(dval(rv));
2102 word0(rv) -= P*Exp_msk1;
2105 #endif /*Sudden_Underflow*/
2106 #endif /*Avoid_Underflow*/
2107 dval(rv) += adj*ulp(dval(rv));
2111 adj = ratio(delta, bs);
2114 if (adj <= 0x7ffffffe) {
2115 /* adj = rounding ? ceil(adj) : floor(adj); */
2118 if (!((rounding>>1) ^ dsign))
2123 #ifdef Avoid_Underflow
2124 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2125 word0(adj) += (2*P+1)*Exp_msk1 - y;
2127 #ifdef Sudden_Underflow
2128 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2129 word0(rv) += P*Exp_msk1;
2130 adj *= ulp(dval(rv));
2135 word0(rv) -= P*Exp_msk1;
2138 #endif /*Sudden_Underflow*/
2139 #endif /*Avoid_Underflow*/
2140 adj *= ulp(dval(rv));
2147 #endif /*Honor_FLT_ROUNDS*/
2150 /* Error is less than half an ulp -- check for
2151 * special case of mantissa a power of two.
2153 if (dsign || word1(rv) || word0(rv) & Bndry_mask
2155 #ifdef Avoid_Underflow
2156 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2158 || (word0(rv) & Exp_mask) <= Exp_msk1
2163 if (!delta->x[0] && delta->wds <= 1)
2168 if (!delta->x[0] && delta->wds <= 1) {
2175 delta = lshift(delta,Log2P);
2176 if (cmp(delta, bs) > 0)
2181 /* exactly half-way between */
2183 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2185 #ifdef Avoid_Underflow
2186 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2187 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2190 /*boundary case -- increment exponent*/
2191 word0(rv) = (word0(rv) & Exp_mask)
2198 #ifdef Avoid_Underflow
2204 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2206 /* boundary case -- decrement exponent */
2207 #ifdef Sudden_Underflow /*{{*/
2208 L = word0(rv) & Exp_mask;
2212 #ifdef Avoid_Underflow
2213 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
2216 #endif /*Avoid_Underflow*/
2220 #else /*Sudden_Underflow}{*/
2221 #ifdef Avoid_Underflow
2223 L = word0(rv) & Exp_mask;
2224 if (L <= (2*P+1)*Exp_msk1) {
2225 if (L > (P+2)*Exp_msk1)
2226 /* round even ==> */
2229 /* rv = smallest denormal */
2233 #endif /*Avoid_Underflow*/
2234 L = (word0(rv) & Exp_mask) - Exp_msk1;
2235 #endif /*Sudden_Underflow}}*/
2236 word0(rv) = L | Bndry_mask1;
2237 word1(rv) = 0xffffffff;
2244 #ifndef ROUND_BIASED
2245 if (!(word1(rv) & LSB))
2249 dval(rv) += ulp(dval(rv));
2250 #ifndef ROUND_BIASED
2252 dval(rv) -= ulp(dval(rv));
2253 #ifndef Sudden_Underflow
2258 #ifdef Avoid_Underflow
2264 if ((aadj = ratio(delta, bs)) <= 2.) {
2267 else if (word1(rv) || word0(rv) & Bndry_mask) {
2268 #ifndef Sudden_Underflow
2269 if (word1(rv) == Tiny1 && !word0(rv))
2276 /* special case -- power of FLT_RADIX to be */
2277 /* rounded down... */
2279 if (aadj < 2./FLT_RADIX)
2280 aadj = 1./FLT_RADIX;
2288 aadj1 = dsign ? aadj : -aadj;
2289 #ifdef Check_FLT_ROUNDS
2291 case 2: /* towards +infinity */
2294 case 0: /* towards 0 */
2295 case 3: /* towards -infinity */
2299 if (Flt_Rounds == 0)
2301 #endif /*Check_FLT_ROUNDS*/
2303 y = word0(rv) & Exp_mask;
2305 /* Check for overflow */
2307 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2308 dval(rv0) = dval(rv);
2309 word0(rv) -= P*Exp_msk1;
2310 adj = aadj1 * ulp(dval(rv));
2312 if ((word0(rv) & Exp_mask) >=
2313 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2314 if (word0(rv0) == Big0 && word1(rv0) == Big1)
2321 word0(rv) += P*Exp_msk1;
2324 #ifdef Avoid_Underflow
2325 if (scale && y <= 2*P*Exp_msk1) {
2326 if (aadj <= 0x7fffffff) {
2327 if ((z = aadj) <= 0)
2330 aadj1 = dsign ? aadj : -aadj;
2332 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
2334 adj = aadj1 * ulp(dval(rv));
2337 #ifdef Sudden_Underflow
2338 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2339 dval(rv0) = dval(rv);
2340 word0(rv) += P*Exp_msk1;
2341 adj = aadj1 * ulp(dval(rv));
2344 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
2346 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2349 if (word0(rv0) == Tiny0
2350 && word1(rv0) == Tiny1)
2357 word0(rv) -= P*Exp_msk1;
2360 adj = aadj1 * ulp(dval(rv));
2363 #else /*Sudden_Underflow*/
2364 /* Compute adj so that the IEEE rounding rules will
2365 * correctly round rv + adj in some half-way cases.
2366 * If rv * ulp(rv) is denormalized (i.e.,
2367 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2368 * trouble from bits lost to denormalization;
2369 * example: 1.2e-307 .
2371 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
2372 aadj1 = (double)(int)(aadj + 0.5);
2376 adj = aadj1 * ulp(dval(rv));
2378 #endif /*Sudden_Underflow*/
2379 #endif /*Avoid_Underflow*/
2381 z = word0(rv) & Exp_mask;
2383 #ifdef Avoid_Underflow
2387 /* Can we stop now? */
2390 /* The tolerances below are conservative. */
2391 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2392 if (aadj < .4999999 || aadj > .5000001)
2395 else if (aadj < .4999999/FLT_RADIX)
2408 word0(rv0) = Exp_1 + (70 << Exp_shift);
2413 else if (!oldinexact)
2416 #ifdef Avoid_Underflow
2418 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
2420 dval(rv) *= dval(rv0);
2422 /* try to avoid the bug of testing an 8087 register value */
2423 if (word0(rv) == 0 && word1(rv) == 0)
2427 #endif /* Avoid_Underflow */
2429 if (inexact && !(word0(rv) & Exp_mask)) {
2430 /* set underflow bit */
2432 dval(rv0) *= dval(rv0);
2444 return sign ? -dval(rv) : dval(rv);
2450 (b, S) Bigint *b, *S;
2452 (Bigint *b, Bigint *S)
2456 ULong *bx, *bxe, q, *sx, *sxe;
2458 ULLong borrow, carry, y, ys;
2460 ULong borrow, carry, y, ys;
2468 /*debug*/ if (b->wds > n)
2469 /*debug*/ Bug("oversize b in quorem");
2477 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
2479 /*debug*/ if (q > 9)
2480 /*debug*/ Bug("oversized quotient in quorem");
2487 ys = *sx++ * (ULLong)q + carry;
2489 y = *bx - (ys & FFFFFFFF) - borrow;
2490 borrow = y >> 32 & (ULong)1;
2491 *bx++ = y & FFFFFFFF;
2495 ys = (si & 0xffff) * q + carry;
2496 zs = (si >> 16) * q + (ys >> 16);
2498 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2499 borrow = (y & 0x10000) >> 16;
2500 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2501 borrow = (z & 0x10000) >> 16;
2504 ys = *sx++ * q + carry;
2506 y = *bx - (ys & 0xffff) - borrow;
2507 borrow = (y & 0x10000) >> 16;
2515 while(--bxe > bx && !*bxe)
2520 if (cmp(b, S) >= 0) {
2530 y = *bx - (ys & FFFFFFFF) - borrow;
2531 borrow = y >> 32 & (ULong)1;
2532 *bx++ = y & FFFFFFFF;
2536 ys = (si & 0xffff) + carry;
2537 zs = (si >> 16) + (ys >> 16);
2539 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2540 borrow = (y & 0x10000) >> 16;
2541 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2542 borrow = (z & 0x10000) >> 16;
2547 y = *bx - (ys & 0xffff) - borrow;
2548 borrow = (y & 0x10000) >> 16;
2557 while(--bxe > bx && !*bxe)
2565 #ifndef MULTIPLE_THREADS
2566 static char *dtoa_result;
2580 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
2583 r = (int*)Balloc(k);
2586 #ifndef MULTIPLE_THREADS
2594 nrv_alloc(s, rve, n) char *s, **rve; int n;
2596 nrv_alloc(char *s, char **rve, int n)
2601 t = rv = rv_alloc(n);
2602 while((*t = *s++)) t++;
2608 /* freedtoa(s) must be used to free values s returned by dtoa
2609 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2610 * but for consistency with earlier versions of dtoa, it is optional
2611 * when MULTIPLE_THREADS is not defined.
2614 static void freedtoa (char *s);
2618 freedtoa(s) char *s;
2623 Bigint *b = (Bigint *)((int *)s - 1);
2624 b->maxwds = 1 << (b->k = *(int*)b);
2626 #ifndef MULTIPLE_THREADS
2627 if (s == dtoa_result)
2633 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2635 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2636 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2639 * 1. Rather than iterating, we use a simple numeric overestimate
2640 * to determine k = floor(log10(d)). We scale relevant
2641 * quantities using O(log2(k)) rather than O(k) multiplications.
2642 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2643 * try to generate digits strictly left to right. Instead, we
2644 * compute with fewer bits and propagate the carry if necessary
2645 * when rounding the final digit up. This is often faster.
2646 * 3. Under the assumption that input will be rounded nearest,
2647 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2648 * That is, we allow equality in stopping tests when the
2649 * round-nearest rule will give the same floating-point value
2650 * as would satisfaction of the stopping test with strict
2652 * 4. We remove common factors of powers of 2 from relevant
2654 * 5. When converting floating-point integers less than 1e16,
2655 * we use floating-point arithmetic rather than resorting
2656 * to multiple-precision integers.
2657 * 6. When asked to produce fewer than 15 digits, we first try
2658 * to get by with floating-point arithmetic; we resort to
2659 * multiple-precision integer arithmetic only if we cannot
2660 * guarantee that the floating-point calculation has given
2661 * the correctly rounded result. For k requested digits and
2662 * "uniformly" distributed input, the probability is
2663 * something like 10^(k-15) that we must resort to the Long
2670 (d, mode, ndigits, decpt, sign, rve)
2671 double d; int mode, ndigits, *decpt, *sign; char **rve;
2673 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
2676 /* Arguments ndigits, decpt, sign are similar to those
2677 of ecvt and fcvt; trailing zeros are suppressed from
2678 the returned string. If not null, *rve is set to point
2679 to the end of the return value. If d is +-Infinity or NaN,
2680 then *decpt is set to 9999.
2683 0 ==> shortest string that yields d when read in
2684 and rounded to nearest.
2685 1 ==> like 0, but with Steele & White stopping rule;
2686 e.g. with IEEE P754 arithmetic , mode 0 gives
2687 1e23 whereas mode 1 gives 9.999999999999999e22.
2688 2 ==> max(1,ndigits) significant digits. This gives a
2689 return value similar to that of ecvt, except
2690 that trailing zeros are suppressed.
2691 3 ==> through ndigits past the decimal point. This
2692 gives a return value similar to that from fcvt,
2693 except that trailing zeros are suppressed, and
2694 ndigits can be negative.
2695 4,5 ==> similar to 2 and 3, respectively, but (in
2696 round-nearest mode) with the tests of mode 0 to
2697 possibly return a shorter string that rounds to d.
2698 With IEEE arithmetic and compilation with
2699 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2700 as modes 2 and 3 when FLT_ROUNDS != 1.
2701 6-9 ==> Debugging modes similar to mode - 4: don't try
2702 fast floating-point estimate (if applicable).
2704 Values of mode other than 0-9 are treated as mode 0.
2706 Sufficient space is allocated to the return value
2707 to hold the suppressed trailing zeros.
2710 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
2711 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2712 spec_case, try_quick;
2714 #ifndef Sudden_Underflow
2718 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
2721 #ifdef Honor_FLT_ROUNDS
2725 int inexact, oldinexact;
2728 #ifndef MULTIPLE_THREADS
2730 freedtoa(dtoa_result);
2735 if (word0(d) & Sign_bit) {
2736 /* set sign for everything, including 0's and NaNs */
2738 word0(d) &= ~Sign_bit; /* clear sign bit */
2743 #if defined(IEEE_Arith) + defined(VAX)
2745 if ((word0(d) & Exp_mask) == Exp_mask)
2747 if (word0(d) == 0x8000)
2750 /* Infinity or NaN */
2753 if (!word1(d) && !(word0(d) & 0xfffff))
2754 return nrv_alloc("Infinity", rve, 8);
2756 return nrv_alloc("NaN", rve, 3);
2760 dval(d) += 0; /* normalize */
2764 return nrv_alloc("0", rve, 1);
2768 try_quick = oldinexact = get_inexact();
2771 #ifdef Honor_FLT_ROUNDS
2772 if ((rounding = Flt_Rounds) >= 2) {
2774 rounding = rounding == 2 ? 0 : 2;
2781 b = d2b(dval(d), &be, &bbits);
2782 #ifdef Sudden_Underflow
2783 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
2785 if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
2788 word0(d2) &= Frac_mask1;
2789 word0(d2) |= Exp_11;
2791 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
2795 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2796 * log10(x) = log(x) / log(10)
2797 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2798 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2800 * This suggests computing an approximation k to log10(d) by
2802 * k = (i - Bias)*0.301029995663981
2803 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2805 * We want k to be too large rather than too small.
2806 * The error in the first-order Taylor series approximation
2807 * is in our favor, so we just round up the constant enough
2808 * to compensate for any error in the multiplication of
2809 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2810 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2811 * adding 1e-13 to the constant term more than suffices.
2812 * Hence we adjust the constant term to 0.1760912590558.
2813 * (We could get a more accurate k by invoking log10,
2814 * but this is probably not worthwhile.)
2822 #ifndef Sudden_Underflow
2826 /* d is denormalized */
2828 i = bbits + be + (Bias + (P-1) - 1);
2829 x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32
2830 : word1(d) << 32 - i;
2832 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2833 i -= (Bias + (P-1) - 1) + 1;
2837 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2839 if (ds < 0. && ds != k)
2840 k--; /* want k = floor(ds) */
2842 if (k >= 0 && k <= Ten_pmax) {
2843 if (dval(d) < tens[k])
2866 if (mode < 0 || mode > 9)
2870 #ifdef Check_FLT_ROUNDS
2871 try_quick = Rounding == 1;
2875 #endif /*SET_INEXACT*/
2895 ilim = ilim1 = i = ndigits;
2901 i = ndigits + k + 1;
2907 s = s0 = rv_alloc(i);
2909 #ifdef Honor_FLT_ROUNDS
2910 if (mode > 1 && rounding != 1)
2914 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2916 /* Try to get by with floating-point arithmetic. */
2922 ieps = 2; /* conservative */
2927 /* prevent overflows */
2929 dval(d) /= bigtens[n_bigtens-1];
2932 for(; j; j >>= 1, i++)
2940 dval(d) *= tens[j1 & 0xf];
2941 for(j = j1 >> 4; j; j >>= 1, i++)
2944 dval(d) *= bigtens[i];
2947 if (k_check && dval(d) < 1. && ilim > 0) {
2955 dval(eps) = ieps*dval(d) + 7.;
2956 word0(eps) -= (P-1)*Exp_msk1;
2960 if (dval(d) > dval(eps))
2962 if (dval(d) < -dval(eps))
2966 #ifndef No_leftright
2968 /* Use Steele & White method of only
2969 * generating digits needed.
2971 dval(eps) = 0.5/tens[ilim-1] - dval(eps);
2975 *s++ = '0' + (int)L;
2976 if (dval(d) < dval(eps))
2978 if (1. - dval(d) < dval(eps))
2988 /* Generate ilim digits, then fix them up. */
2989 dval(eps) *= tens[ilim-1];
2990 for(i = 1;; i++, dval(d) *= 10.) {
2991 L = (Long)(dval(d));
2992 if (!(dval(d) -= L))
2994 *s++ = '0' + (int)L;
2996 if (dval(d) > 0.5 + dval(eps))
2998 else if (dval(d) < 0.5 - dval(eps)) {
3006 #ifndef No_leftright
3016 /* Do we have a "small" integer? */
3018 if (be >= 0 && k <= Int_max) {
3021 if (ndigits < 0 && ilim <= 0) {
3023 if (ilim < 0 || dval(d) <= 5*ds)
3027 for(i = 1;; i++, dval(d) *= 10.) {
3028 L = (Long)(dval(d) / ds);
3030 #ifdef Check_FLT_ROUNDS
3031 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
3037 *s++ = '0' + (int)L;
3045 #ifdef Honor_FLT_ROUNDS
3049 case 2: goto bump_up;
3053 if (dval(d) > ds || dval(d) == ds && L & 1) {
3074 #ifndef Sudden_Underflow
3075 denorm ? be + (Bias + (P-1) - 1 + 1) :
3078 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
3086 if (m2 > 0 && s2 > 0) {
3087 i = m2 < s2 ? m2 : s2;
3095 mhi = pow5mult(mhi, m5);
3104 b = pow5mult(b, b5);
3108 S = pow5mult(S, s5);
3110 /* Check for special case that d is a normalized power of 2. */
3113 if ((mode < 2 || leftright)
3114 #ifdef Honor_FLT_ROUNDS
3118 if (!word1(d) && !(word0(d) & Bndry_mask)
3119 #ifndef Sudden_Underflow
3120 && word0(d) & (Exp_mask & ~Exp_msk1)
3123 /* The special case */
3130 /* Arrange for convenient computation of quotients:
3131 * shift left if necessary so divisor has 4 leading 0 bits.
3133 * Perhaps we should just compute leading 28 bits of S once
3134 * and for all and pass them and a shift to quorem, so it
3135 * can do shifts and ors to compute the numerator for q.
3138 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
3141 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
3163 b = multadd(b, 10, 0); /* we botched the k estimate */
3165 mhi = multadd(mhi, 10, 0);
3169 if (ilim <= 0 && (mode == 3 || mode == 5)) {
3170 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
3171 /* no digits, fcvt style */
3183 mhi = lshift(mhi, m2);
3185 /* Compute mlo -- check for special case
3186 * that d is a normalized power of 2.
3191 mhi = Balloc(mhi->k);
3193 mhi = lshift(mhi, Log2P);
3197 dig = quorem(b,S) + '0';
3198 /* Do we yet have the shortest decimal string
3199 * that will round to d?
3202 delta = diff(S, mhi);
3203 j1 = delta->sign ? 1 : cmp(b, delta);
3205 #ifndef ROUND_BIASED
3206 if (j1 == 0 && mode != 1 && !(word1(d) & 1)
3207 #ifdef Honor_FLT_ROUNDS
3216 else if (!b->x[0] && b->wds <= 1)
3223 if (j < 0 || j == 0 && mode != 1
3224 #ifndef ROUND_BIASED
3228 if (!b->x[0] && b->wds <= 1) {
3234 #ifdef Honor_FLT_ROUNDS
3237 case 0: goto accept_dig;
3238 case 2: goto keep_dig;
3240 #endif /*Honor_FLT_ROUNDS*/
3244 if ((j1 > 0 || j1 == 0 && dig & 1)
3253 #ifdef Honor_FLT_ROUNDS
3257 if (dig == '9') { /* possible if i == 1 */
3265 #ifdef Honor_FLT_ROUNDS
3271 b = multadd(b, 10, 0);
3273 mlo = mhi = multadd(mhi, 10, 0);
3275 mlo = multadd(mlo, 10, 0);
3276 mhi = multadd(mhi, 10, 0);
3282 *s++ = dig = quorem(b,S) + '0';
3283 if (!b->x[0] && b->wds <= 1) {
3291 b = multadd(b, 10, 0);
3294 /* Round off last digit */
3296 #ifdef Honor_FLT_ROUNDS
3298 case 0: goto trimzeros;
3299 case 2: goto roundoff;
3304 if (j > 0 || j == 0 && dig & 1) {
3322 if (mlo && mlo != mhi)
3330 word0(d) = Exp_1 + (70 << Exp_shift);
3335 else if (!oldinexact)