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 ***************************************************************/
21 #define freedtoa __freedtoa
24 #define Omit_Private_Memory
25 #define MULTIPLE_THREADS 1
26 /* Lock 0 is not used because of USE_MALLOC, Lock 1 protects a lazy-initialized table */
27 #define ACQUIRE_DTOA_LOCK(n)
28 #define FREE_DTOA_LOCK(n)
30 /* Please send bug reports to David M. Gay (dmg at acm dot org,
31 * with " at " changed at "@" and " dot " changed to "."). */
33 /* On a machine with IEEE extended-precision registers, it is
34 * necessary to specify double-precision (53-bit) rounding precision
35 * before invoking strtod or dtoa. If the machine uses (the equivalent
36 * of) Intel 80x87 arithmetic, the call
37 * _control87(PC_53, MCW_PC);
38 * does this with many compilers. Whether this or another call is
39 * appropriate depends on the compiler; for this to work, it may be
40 * necessary to #include "float.h" or another system-dependent header
44 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
46 * This strtod returns a nearest machine number to the input decimal
47 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
48 * broken by the IEEE round-even rule. Otherwise ties are broken by
49 * biased rounding (add half and chop).
51 * Inspired loosely by William D. Clinger's paper "How to Read Floating
52 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
56 * 1. We only require IEEE, IBM, or VAX double-precision
57 * arithmetic (not IEEE double-extended).
58 * 2. We get by with floating-point arithmetic in a case that
59 * Clinger missed -- when we're computing d * 10^n
60 * for a small integer d and the integer n is not too
61 * much larger than 22 (the maximum integer k for which
62 * we can represent 10^k exactly), we may be able to
63 * compute (d*10^k) * 10^(e-k) with just one roundoff.
64 * 3. Rather than a bit-at-a-time adjustment of the binary
65 * result in the hard case, we use floating-point
66 * arithmetic to determine the adjustment to within
67 * one bit; only in really hard cases do we need to
68 * compute a second residual.
69 * 4. Because of 3., we don't need a large table of powers of 10
70 * for ten-to-e (just some small tables, e.g. of 10^k
75 * #define IEEE_8087 for IEEE-arithmetic machines where the least
76 * significant byte has the lowest address.
77 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
78 * significant byte has the lowest address.
79 * #define Long int on machines with 32-bit ints and 64-bit longs.
80 * #define IBM for IBM mainframe-style floating-point arithmetic.
81 * #define VAX for VAX-style floating-point arithmetic (D_floating).
82 * #define No_leftright to omit left-right logic in fast floating-point
83 * computation of dtoa.
84 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
85 * and strtod and dtoa should round accordingly.
86 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
87 * and Honor_FLT_ROUNDS is not #defined.
88 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
89 * that use extended-precision instructions to compute rounded
90 * products and quotients) with IBM.
91 * #define ROUND_BIASED for IEEE-format with biased rounding.
92 * #define Inaccurate_Divide for IEEE-format with correctly rounded
93 * products but inaccurate quotients, e.g., for Intel i860.
94 * #define NO_LONG_LONG on machines that do not have a "long long"
95 * integer type (of >= 64 bits). On such machines, you can
96 * #define Just_16 to store 16 bits per 32-bit Long when doing
97 * high-precision integer arithmetic. Whether this speeds things
98 * up or slows things down depends on the machine and the number
99 * being converted. If long long is available and the name is
100 * something other than "long long", #define Llong to be the name,
101 * and if "unsigned Llong" does not work as an unsigned version of
102 * Llong, #define #ULLong to be the corresponding unsigned type.
103 * #define KR_headers for old-style C function headers.
104 * #define Bad_float_h if your system lacks a float.h or if it does not
105 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
106 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
107 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
108 * if memory is available and otherwise does something you deem
109 * appropriate. If MALLOC is undefined, malloc will be invoked
110 * directly -- and assumed always to succeed.
111 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
112 * memory allocations from a private pool of memory when possible.
113 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
114 * unless #defined to be a different length. This default length
115 * suffices to get rid of MALLOC calls except for unusual cases,
116 * such as decimal-to-binary conversion of a very long string of
117 * digits. The longest string dtoa can return is about 751 bytes
118 * long. For conversions by strtod of strings of 800 digits and
119 * all dtoa conversions in single-threaded executions with 8-byte
120 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
121 * pointers, PRIVATE_MEM >= 7112 appears adequate.
122 * #define INFNAN_CHECK on IEEE systems to cause strtod to check for
123 * Infinity and NaN (case insensitively). On some systems (e.g.,
124 * some HP systems), it may be necessary to #define NAN_WORD0
125 * appropriately -- to the most significant word of a quiet NaN.
126 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
127 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
128 * strtod also accepts (case insensitively) strings of the form
129 * NaN(x), where x is a string of hexadecimal digits and spaces;
130 * if there is only one string of hexadecimal digits, it is taken
131 * for the 52 fraction bits of the resulting NaN; if there are two
132 * or more strings of hex digits, the first is for the high 20 bits,
133 * the second and subsequent for the low 32 bits, with intervening
134 * white space ignored; but if this results in none of the 52
135 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
136 * and NAN_WORD1 are used instead.
137 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
138 * multiple threads. In this case, you must provide (or suitably
139 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
140 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
141 * in pow5mult, ensures lazy evaluation of only one copy of high
142 * powers of 5; omitting this lock would introduce a small
143 * probability of wasting memory, but would otherwise be harmless.)
144 * You must also invoke freedtoa(s) to free the value s returned by
145 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
146 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
147 * avoids underflows on inputs whose result does not underflow.
148 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
149 * floating-point numbers and flushes underflows to zero rather
150 * than implementing gradual underflow, then you must also #define
152 * #define YES_ALIAS to permit aliasing certain double values with
153 * arrays of ULongs. This leads to slightly better code with
154 * some compilers and was always used prior to 19990916, but it
155 * is not strictly legal and can cause trouble with aggressively
156 * optimizing compilers (e.g., gcc 2.95.1 under -O2).
157 * #define USE_LOCALE to use the current locale's decimal_point value.
158 * #define SET_INEXACT if IEEE arithmetic is being used and extra
159 * computation should be done to set the inexact flag when the
160 * result is inexact and avoid setting inexact when the result
161 * is exact. In this case, dtoa.c must be compiled in
162 * an environment, perhaps provided by #include "dtoa.c" in a
163 * suitable wrapper, that defines two functions,
164 * int get_inexact(void);
165 * void clear_inexact(void);
166 * such that get_inexact() returns a nonzero value if the
167 * inexact bit is already set, and clear_inexact() sets the
168 * inexact bit to 0. When SET_INEXACT is #defined, strtod
169 * also does extra computations to set the underflow and overflow
170 * flags when appropriate (i.e., when the result is tiny and
171 * inexact or when it is a numeric value rounded to +-infinity).
172 * #define NO_ERRNO if strtod should not assign errno = ERANGE when
173 * the result overflows to +-Infinity or underflows to 0.
175 #if defined(i386) || defined(mips) && defined(MIPSEL) || defined (__arm__)
179 #elif defined(__x86_64__) || defined(__alpha__)
183 #elif defined(__ia64)
191 #elif defined(__hppa)
200 #define ULong guint32
204 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
217 extern char *MALLOC();
219 extern void *MALLOC(size_t);
222 #define MALLOC malloc
225 #define Omit_Private_Memory
226 #ifndef Omit_Private_Memory
228 #define PRIVATE_MEM 2304
230 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
231 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
235 #undef Avoid_Underflow
249 #define DBL_MAX_10_EXP 308
250 #define DBL_MAX_EXP 1024
252 #endif /*IEEE_Arith*/
256 #define DBL_MAX_10_EXP 75
257 #define DBL_MAX_EXP 63
259 #define DBL_MAX 7.2370055773322621e+75
264 #define DBL_MAX_10_EXP 38
265 #define DBL_MAX_EXP 127
267 #define DBL_MAX 1.7014118346046923e+38
271 #define LONG_MAX 2147483647
274 #else /* ifndef Bad_float_h */
276 #endif /* Bad_float_h */
288 #define CONST /* blank */
294 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
295 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
298 typedef union { double d; ULong L[2]; } U;
303 #define word0(x) ((ULong *)&x)[1]
304 #define word1(x) ((ULong *)&x)[0]
306 #define word0(x) ((ULong *)&x)[0]
307 #define word1(x) ((ULong *)&x)[1]
311 #define word0(x) ((U*)&x)->L[1]
312 #define word1(x) ((U*)&x)->L[0]
314 #define word0(x) ((U*)&x)->L[0]
315 #define word1(x) ((U*)&x)->L[1]
317 #define dval(x) ((U*)&x)->d
320 /* The following definition of Storeinc is appropriate for MIPS processors.
321 * An alternative that might be better on some machines is
322 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
324 #if defined(IEEE_8087) + defined(VAX)
325 #define Storeinc(a,b,c) do { (((unsigned short *)a)[1] = (unsigned short)b, \
326 ((unsigned short *)a)[0] = (unsigned short)c, a++) } while (0)
328 #define Storeinc(a,b,c) do { (((unsigned short *)a)[0] = (unsigned short)b, \
329 ((unsigned short *)a)[1] = (unsigned short)c, a++) } while (0)
332 /* #define P DBL_MANT_DIG */
333 /* Ten_pmax = floor(P*log(2)/log(5)) */
334 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
335 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
336 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
340 #define Exp_shift1 20
341 #define Exp_msk1 0x100000
342 #define Exp_msk11 0x100000
343 #define Exp_mask 0x7ff00000
347 #define Exp_1 0x3ff00000
348 #define Exp_11 0x3ff00000
350 #define Frac_mask 0xfffff
351 #define Frac_mask1 0xfffff
354 #define Bndry_mask 0xfffff
355 #define Bndry_mask1 0xfffff
357 #define Sign_bit 0x80000000
363 #ifndef NO_IEEE_Scale
364 #define Avoid_Underflow
365 #ifdef Flush_Denorm /* debugging option */
366 #undef Sudden_Underflow
372 #define Flt_Rounds FLT_ROUNDS
376 #endif /*Flt_Rounds*/
378 #ifdef Honor_FLT_ROUNDS
379 #define Rounding rounding
380 #undef Check_FLT_ROUNDS
381 #define Check_FLT_ROUNDS
383 #define Rounding Flt_Rounds
386 #else /* ifndef IEEE_Arith */
387 #undef Check_FLT_ROUNDS
388 #undef Honor_FLT_ROUNDS
390 #undef Sudden_Underflow
391 #define Sudden_Underflow
396 #define Exp_shift1 24
397 #define Exp_msk1 0x1000000
398 #define Exp_msk11 0x1000000
399 #define Exp_mask 0x7f000000
402 #define Exp_1 0x41000000
403 #define Exp_11 0x41000000
404 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
405 #define Frac_mask 0xffffff
406 #define Frac_mask1 0xffffff
409 #define Bndry_mask 0xefffff
410 #define Bndry_mask1 0xffffff
412 #define Sign_bit 0x80000000
414 #define Tiny0 0x100000
423 #define Exp_msk1 0x80
424 #define Exp_msk11 0x800000
425 #define Exp_mask 0x7f80
428 #define Exp_1 0x40800000
429 #define Exp_11 0x4080
431 #define Frac_mask 0x7fffff
432 #define Frac_mask1 0xffff007f
435 #define Bndry_mask 0xffff007f
436 #define Bndry_mask1 0xffff007f
438 #define Sign_bit 0x8000
444 #endif /* IBM, VAX */
445 #endif /* IEEE_Arith */
452 #define rounded_product(a,b) a = rnd_prod(a, b)
453 #define rounded_quotient(a,b) a = rnd_quot(a, b)
455 extern double rnd_prod(), rnd_quot();
457 extern double rnd_prod(double, double), rnd_quot(double, double);
460 #define rounded_product(a,b) a *= b
461 #define rounded_quotient(a,b) a /= b
464 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
465 #define Big1 0xffffffff
472 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
474 #define FFFFFFFF 0xffffffffUL
481 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
482 * This makes some inner loops simpler and sometimes saves work
483 * during multiplications, but it often seems to make things slightly
484 * slower. Hence the default is now to store 32 bits per Long.
487 #else /* long long available */
489 #define Llong long long
492 #define ULLong unsigned Llong
494 #endif /* NO_LONG_LONG */
496 #ifndef MULTIPLE_THREADS
497 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
498 #define FREE_DTOA_LOCK(n) /*nothing*/
504 extern "C" double strtod(const char *s00, char **se);
505 extern "C" char *dtoa(double d, int mode, int ndigits,
506 int *decpt, int *sign, char **rve);
512 int k, maxwds, sign, wds;
516 typedef struct Bigint Bigint;
518 static Bigint *freelist[Kmax+1];
530 #ifndef Omit_Private_Memory
534 ACQUIRE_DTOA_LOCK(0);
535 if ((rv = freelist[k])) {
536 freelist[k] = rv->next;
540 #ifdef Omit_Private_Memory
541 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
543 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
545 if (pmem_next - private_mem + len <= PRIVATE_mem) {
546 rv = (Bigint*)pmem_next;
550 rv = (Bigint*)MALLOC(len*sizeof(double));
556 rv->sign = rv->wds = 0;
568 #ifdef Omit_Private_Memory
572 ACQUIRE_DTOA_LOCK(0);
573 v->next = freelist[v->k];
580 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
581 y->wds*sizeof(Long) + 2*sizeof(int))
586 (b, m, a) Bigint *b; int m, a;
588 (Bigint *b, int m, int a) /* multiply by m and add a */
609 y = *x * (ULLong)m + carry;
615 y = (xi & 0xffff) * m + carry;
616 z = (xi >> 16) * m + (y >> 16);
618 *x++ = (z << 16) + (y & 0xffff);
628 if (wds >= b->maxwds) {
643 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
645 (CONST char *s, int nd0, int nd, ULong y9)
653 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
660 b->x[0] = y9 & 0xffff;
661 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
667 do b = multadd(b, 10, *s++ - '0');
674 b = multadd(b, 10, *s++ - '0');
681 (x) register ULong x;
688 if (!(x & 0xffff0000)) {
692 if (!(x & 0xff000000)) {
696 if (!(x & 0xf0000000)) {
700 if (!(x & 0xc0000000)) {
704 if (!(x & 0x80000000)) {
706 if (!(x & 0x40000000))
721 register ULong x = *y;
779 (a, b) Bigint *a, *b;
781 (Bigint *a, Bigint *b)
786 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
797 if (a->wds < b->wds) {
809 for(x = c->x, xa = x + wc; x < xa; x++)
817 for(; xb < xbe; xc0++) {
823 z = *x++ * (ULLong)y + *xc + carry;
825 *xc++ = z & FFFFFFFF;
833 for(; xb < xbe; xb++, xc0++) {
834 if (y = *xb & 0xffff) {
839 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
841 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
854 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
857 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
865 for(; xb < xbe; xc0++) {
871 z = *x++ * y + *xc + carry;
881 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
891 (b, k) Bigint *b; int k;
896 Bigint *b1, *p5, *p51;
898 static int p05[3] = { 5, 25, 125 };
901 b = multadd(b, p05[i-1], 0);
907 #ifdef MULTIPLE_THREADS
908 ACQUIRE_DTOA_LOCK(1);
927 if (!(p51 = p5->next)) {
928 #ifdef MULTIPLE_THREADS
929 ACQUIRE_DTOA_LOCK(1);
930 if (!(p51 = p5->next)) {
931 p51 = p5->next = mult(p5,p5);
936 p51 = p5->next = mult(p5,p5);
948 (b, k) Bigint *b; int k;
955 ULong *x, *x1, *xe, z;
964 for(i = b->maxwds; n1 > i; i <<= 1)
968 for(i = 0; i < n; i++)
989 *x1++ = *x << k & 0xffff | z;
1008 (a, b) Bigint *a, *b;
1010 (Bigint *a, Bigint *b)
1013 ULong *xa, *xa0, *xb, *xb0;
1019 if (i > 1 && !a->x[i-1])
1020 Bug("cmp called with a->x[a->wds-1] == 0");
1021 if (j > 1 && !b->x[j-1])
1022 Bug("cmp called with b->x[b->wds-1] == 0");
1032 return *xa < *xb ? -1 : 1;
1042 (a, b) Bigint *a, *b;
1044 (Bigint *a, Bigint *b)
1049 ULong *xa, *xae, *xb, *xbe, *xc;
1086 y = (ULLong)*xa++ - *xb++ - borrow;
1087 borrow = y >> 32 & (ULong)1;
1088 *xc++ = y & FFFFFFFF;
1093 borrow = y >> 32 & (ULong)1;
1094 *xc++ = y & FFFFFFFF;
1099 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1100 borrow = (y & 0x10000) >> 16;
1101 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1102 borrow = (z & 0x10000) >> 16;
1107 y = (*xa & 0xffff) - borrow;
1108 borrow = (y & 0x10000) >> 16;
1109 z = (*xa++ >> 16) - borrow;
1110 borrow = (z & 0x10000) >> 16;
1115 y = *xa++ - *xb++ - borrow;
1116 borrow = (y & 0x10000) >> 16;
1122 borrow = (y & 0x10000) >> 16;
1144 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1145 #ifndef Avoid_Underflow
1146 #ifndef Sudden_Underflow
1155 #ifndef Avoid_Underflow
1156 #ifndef Sudden_Underflow
1159 L = -L >> Exp_shift;
1160 if (L < Exp_shift) {
1161 word0(a) = 0x80000 >> L;
1167 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
1178 (a, e) Bigint *a; int *e;
1183 ULong *xa, *xa0, w, y, z;
1197 if (!y) Bug("zero y in b2d");
1203 d0 = Exp_1 | (y >> (Ebits - k));
1204 w = xa > xa0 ? *--xa : 0;
1205 d1 = y << ((32-Ebits) + k) | (w >> (Ebits - k));
1208 z = xa > xa0 ? *--xa : 0;
1210 d0 = Exp_1 | y << k | (z >> (32 - k));
1211 y = xa > xa0 ? *--xa : 0;
1212 d1 = z << k | (y >> (32 - k));
1219 if (k < Ebits + 16) {
1220 z = xa > xa0 ? *--xa : 0;
1221 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1222 w = xa > xa0 ? *--xa : 0;
1223 y = xa > xa0 ? *--xa : 0;
1224 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1227 z = xa > xa0 ? *--xa : 0;
1228 w = xa > xa0 ? *--xa : 0;
1230 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1231 y = xa > xa0 ? *--xa : 0;
1232 d1 = w << k + 16 | y << k;
1236 word0(d) = d0 >> 16 | d0 << 16;
1237 word1(d) = d1 >> 16 | d1 << 16;
1248 (d, e, bits) double d; int *e, *bits;
1250 (double d, int *e, int *bits)
1256 #ifndef Sudden_Underflow
1261 d0 = word0(d) >> 16 | word0(d) << 16;
1262 d1 = word1(d) >> 16 | word1(d) << 16;
1276 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1277 #ifdef Sudden_Underflow
1278 de = (int)(d0 >> Exp_shift);
1283 if ((de = (int)(d0 >> Exp_shift)))
1288 if ((k = lo0bits(&y))) {
1289 x[0] = y | (z << (32 - k));
1294 #ifndef Sudden_Underflow
1297 b->wds = (x[1] = z) ? 2 : 1;
1302 Bug("Zero passed to d2b");
1306 #ifndef Sudden_Underflow
1314 if (k = lo0bits(&y))
1316 x[0] = y | z << 32 - k & 0xffff;
1317 x[1] = z >> k - 16 & 0xffff;
1323 x[1] = y >> 16 | z << 16 - k & 0xffff;
1324 x[2] = z >> k & 0xffff;
1339 Bug("Zero passed to d2b");
1357 #ifndef Sudden_Underflow
1361 *e = (de - Bias - (P-1) << 2) + k;
1362 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1364 *e = de - Bias - (P-1) + k;
1367 #ifndef Sudden_Underflow
1370 *e = de - Bias - (P-1) + 1 + k;
1372 *bits = 32*i - hi0bits(x[i-1]);
1374 *bits = (i+2)*16 - hi0bits(x[i]);
1386 (a, b) Bigint *a, *b;
1388 (Bigint *a, Bigint *b)
1394 dval(da) = b2d(a, &ka);
1395 dval(db) = b2d(b, &kb);
1397 k = ka - kb + 32*(a->wds - b->wds);
1399 k = ka - kb + 16*(a->wds - b->wds);
1403 word0(da) += (k >> 2)*Exp_msk1;
1409 word0(db) += (k >> 2)*Exp_msk1;
1415 word0(da) += k*Exp_msk1;
1418 word0(db) += k*Exp_msk1;
1421 return dval(da) / dval(db);
1426 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1427 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1436 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1437 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1438 #ifdef Avoid_Underflow
1439 9007199254740992.*9007199254740992.e-256
1440 /* = 2^106 * 1e-53 */
1445 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1446 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1447 #define Scale_Bit 0x10
1451 bigtens[] = { 1e16, 1e32, 1e64 };
1452 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1455 bigtens[] = { 1e16, 1e32 };
1456 static CONST double tinytens[] = { 1e-16, 1e-32 };
1468 #define NAN_WORD0 0x7ff80000
1478 (sp, t) char **sp, *t;
1480 (CONST char **sp, char *t)
1484 CONST char *s = *sp;
1487 if ((c = *++s) >= 'A' && c <= 'Z')
1500 (rvp, sp) double *rvp; CONST char **sp;
1502 (double *rvp, CONST char **sp)
1507 int havedig, udx0, xshift;
1510 havedig = xshift = 0;
1513 while(c = *(CONST unsigned char*)++s) {
1514 if (c >= '0' && c <= '9')
1516 else if (c >= 'a' && c <= 'f')
1518 else if (c >= 'A' && c <= 'F')
1520 else if (c <= ' ') {
1521 if (udx0 && havedig) {
1527 else if (/*(*/ c == ')' && havedig) {
1532 return; /* invalid form: don't change *sp */
1540 x[0] = (x[0] << 4) | (x[1] >> 28);
1541 x[1] = (x[1] << 4) | c;
1543 if ((x[0] &= 0xfffff) || x[1]) {
1544 word0(*rvp) = Exp_mask | x[0];
1548 #endif /*No_Hex_NaN*/
1549 #endif /* INFNAN_CHECK */
1552 * LOCKING: This is not thread-safe, since the locking macros are defined as no-ops,
1553 * the caller should lock.
1559 (s00, se) CONST char *s00; char **se;
1561 (CONST char *s00, char **se)
1564 #ifdef Avoid_Underflow
1567 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1568 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1569 CONST char *s, *s0, *s1;
1570 double aadj, aadj1, adj, rv, rv0;
1573 Bigint *bb = NULL, *bb1, *bd = NULL, *bd0, *bs = NULL, *delta = NULL;
1575 int inexact, oldinexact;
1577 #ifdef Honor_FLT_ROUNDS
1584 sign = nz0 = nz = 0;
1586 for(s = s00;;s++) switch(*s) {
1609 while(*++s == '0') ;
1615 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1622 s1 = localeconv()->decimal_point;
1643 for(; c == '0'; c = *++s)
1645 if (c > '0' && c <= '9') {
1653 for(; c >= '0' && c <= '9'; c = *++s) {
1658 for(i = 1; i < nz; i++)
1661 else if (nd <= DBL_DIG + 1)
1665 else if (nd <= DBL_DIG + 1)
1673 if (c == 'e' || c == 'E') {
1674 if (!nd && !nz && !nz0) {
1685 if (c >= '0' && c <= '9') {
1688 if (c > '0' && c <= '9') {
1691 while((c = *++s) >= '0' && c <= '9')
1693 if (s - s1 > 8 || L > 19999)
1694 /* Avoid confusion from exponents
1695 * so large that e might overflow.
1697 e = 19999; /* safe for 16 bit ints */
1712 /* Check for Nan and Infinity */
1716 if (match(&s,"nf")) {
1718 if (!match(&s,"inity"))
1720 word0(rv) = 0x7ff00000;
1727 if (match(&s, "an")) {
1728 word0(rv) = NAN_WORD0;
1729 word1(rv) = NAN_WORD1;
1731 if (*s == '(') /*)*/
1737 #endif /* INFNAN_CHECK */
1746 /* Now we have nd0 digits, starting at s0, followed by a
1747 * decimal point, followed by nd-nd0 digits. The number we're
1748 * after is the integer represented by those digits times
1753 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1758 oldinexact = get_inexact();
1760 dval(rv) = tens[k - 9] * dval(rv) + z;
1764 #ifndef RND_PRODQUOT
1765 #ifndef Honor_FLT_ROUNDS
1773 if (e <= Ten_pmax) {
1775 goto vax_ovfl_check;
1777 #ifdef Honor_FLT_ROUNDS
1778 /* round correctly FLT_ROUNDS = 2 or 3 */
1784 /* rv = */ rounded_product(dval(rv), tens[e]);
1789 if (e <= Ten_pmax + i) {
1790 /* A fancier test would sometimes let us do
1791 * this for larger i values.
1793 #ifdef Honor_FLT_ROUNDS
1794 /* round correctly FLT_ROUNDS = 2 or 3 */
1801 dval(rv) *= tens[i];
1803 /* VAX exponent range is so narrow we must
1804 * worry about overflow here...
1807 word0(rv) -= P*Exp_msk1;
1808 /* rv = */ rounded_product(dval(rv), tens[e]);
1809 if ((word0(rv) & Exp_mask)
1810 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1812 word0(rv) += P*Exp_msk1;
1814 /* rv = */ rounded_product(dval(rv), tens[e]);
1819 #ifndef Inaccurate_Divide
1820 else if (e >= -Ten_pmax) {
1821 #ifdef Honor_FLT_ROUNDS
1822 /* round correctly FLT_ROUNDS = 2 or 3 */
1828 /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1839 oldinexact = get_inexact();
1841 #ifdef Avoid_Underflow
1844 #ifdef Honor_FLT_ROUNDS
1845 if ((rounding = Flt_Rounds) >= 2) {
1847 rounding = rounding == 2 ? 0 : 2;
1853 #endif /*IEEE_Arith*/
1855 /* Get starting approximation = rv * 10**e1 */
1859 dval(rv) *= tens[i];
1861 if (e1 > DBL_MAX_10_EXP) {
1866 /* Can't trust HUGE_VAL */
1868 #ifdef Honor_FLT_ROUNDS
1870 case 0: /* toward 0 */
1871 case 3: /* toward -infinity */
1876 word0(rv) = Exp_mask;
1879 #else /*Honor_FLT_ROUNDS*/
1880 word0(rv) = Exp_mask;
1882 #endif /*Honor_FLT_ROUNDS*/
1884 /* set overflow bit */
1886 dval(rv0) *= dval(rv0);
1888 #else /*IEEE_Arith*/
1891 #endif /*IEEE_Arith*/
1897 for(j = 0; e1 > 1; j++, e1 >>= 1)
1899 dval(rv) *= bigtens[j];
1900 /* The last multiplication could overflow. */
1901 word0(rv) -= P*Exp_msk1;
1902 dval(rv) *= bigtens[j];
1903 if ((z = word0(rv) & Exp_mask)
1904 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1906 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1907 /* set to largest number */
1908 /* (Can't trust DBL_MAX) */
1913 word0(rv) += P*Exp_msk1;
1919 dval(rv) /= tens[i];
1921 if (e1 >= 1 << n_bigtens)
1923 #ifdef Avoid_Underflow
1926 for(j = 0; e1 > 0; j++, e1 >>= 1)
1928 dval(rv) *= tinytens[j];
1929 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
1930 >> Exp_shift)) > 0) {
1931 /* scaled rv is denormal; zap j low bits */
1935 word0(rv) = (P+2)*Exp_msk1;
1937 word0(rv) &= 0xffffffff << (j-32);
1940 word1(rv) &= 0xffffffff << j;
1943 for(j = 0; e1 > 1; j++, e1 >>= 1)
1945 dval(rv) *= tinytens[j];
1946 /* The last multiplication could underflow. */
1947 dval(rv0) = dval(rv);
1948 dval(rv) *= tinytens[j];
1950 dval(rv) = 2.*dval(rv0);
1951 dval(rv) *= tinytens[j];
1963 #ifndef Avoid_Underflow
1966 /* The refinement below will clean
1967 * this approximation up.
1974 /* Now the hard part -- adjusting rv to the correct value.*/
1976 /* Put digits into bd: true value = bd * 10^e */
1978 bd0 = s2b(s0, nd0, nd, y);
1981 bd = Balloc(bd0->k);
1983 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
1999 #ifdef Honor_FLT_ROUNDS
2003 #ifdef Avoid_Underflow
2005 i = j + bbbits - 1; /* logb(rv) */
2006 if (i < Emin) /* denormal */
2010 #else /*Avoid_Underflow*/
2011 #ifdef Sudden_Underflow
2013 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
2017 #else /*Sudden_Underflow*/
2019 i = j + bbbits - 1; /* logb(rv) */
2020 if (i < Emin) /* denormal */
2024 #endif /*Sudden_Underflow*/
2025 #endif /*Avoid_Underflow*/
2028 #ifdef Avoid_Underflow
2031 i = bb2 < bd2 ? bb2 : bd2;
2040 bs = pow5mult(bs, bb5);
2046 bb = lshift(bb, bb2);
2048 bd = pow5mult(bd, bd5);
2050 bd = lshift(bd, bd2);
2052 bs = lshift(bs, bs2);
2053 delta = diff(bb, bd);
2054 dsign = delta->sign;
2057 #ifdef Honor_FLT_ROUNDS
2058 if (rounding != 1) {
2060 /* Error is less than an ulp */
2061 if (!delta->x[0] && delta->wds <= 1) {
2077 && !(word0(rv) & Frac_mask)) {
2078 y = word0(rv) & Exp_mask;
2079 #ifdef Avoid_Underflow
2080 if (!scale || y > 2*P*Exp_msk1)
2085 delta = lshift(delta,Log2P);
2086 if (cmp(delta, bs) <= 0)
2091 #ifdef Avoid_Underflow
2092 if (scale && (y = word0(rv) & Exp_mask)
2094 word0(adj) += (2*P+1)*Exp_msk1 - y;
2096 #ifdef Sudden_Underflow
2097 if ((word0(rv) & Exp_mask) <=
2099 word0(rv) += P*Exp_msk1;
2100 dval(rv) += adj*ulp(dval(rv));
2101 word0(rv) -= P*Exp_msk1;
2104 #endif /*Sudden_Underflow*/
2105 #endif /*Avoid_Underflow*/
2106 dval(rv) += adj*ulp(dval(rv));
2110 adj = ratio(delta, bs);
2113 if (adj <= 0x7ffffffe) {
2114 /* adj = rounding ? ceil(adj) : floor(adj); */
2117 if (!((rounding>>1) ^ dsign))
2122 #ifdef Avoid_Underflow
2123 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2124 word0(adj) += (2*P+1)*Exp_msk1 - y;
2126 #ifdef Sudden_Underflow
2127 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2128 word0(rv) += P*Exp_msk1;
2129 adj *= ulp(dval(rv));
2134 word0(rv) -= P*Exp_msk1;
2137 #endif /*Sudden_Underflow*/
2138 #endif /*Avoid_Underflow*/
2139 adj *= ulp(dval(rv));
2146 #endif /*Honor_FLT_ROUNDS*/
2149 /* Error is less than half an ulp -- check for
2150 * special case of mantissa a power of two.
2152 if (dsign || word1(rv) || word0(rv) & Bndry_mask
2154 #ifdef Avoid_Underflow
2155 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2157 || (word0(rv) & Exp_mask) <= Exp_msk1
2162 if (!delta->x[0] && delta->wds <= 1)
2167 if (!delta->x[0] && delta->wds <= 1) {
2174 delta = lshift(delta,Log2P);
2175 if (cmp(delta, bs) > 0)
2180 /* exactly half-way between */
2182 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2184 #ifdef Avoid_Underflow
2185 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2186 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2189 /*boundary case -- increment exponent*/
2190 word0(rv) = (word0(rv) & Exp_mask)
2197 #ifdef Avoid_Underflow
2203 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2205 /* boundary case -- decrement exponent */
2206 #ifdef Sudden_Underflow /*{{*/
2207 L = word0(rv) & Exp_mask;
2211 #ifdef Avoid_Underflow
2212 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
2215 #endif /*Avoid_Underflow*/
2219 #else /*Sudden_Underflow}{*/
2220 #ifdef Avoid_Underflow
2222 L = word0(rv) & Exp_mask;
2223 if (L <= (2*P+1)*Exp_msk1) {
2224 if (L > (P+2)*Exp_msk1)
2225 /* round even ==> */
2228 /* rv = smallest denormal */
2232 #endif /*Avoid_Underflow*/
2233 L = (word0(rv) & Exp_mask) - Exp_msk1;
2234 #endif /*Sudden_Underflow}}*/
2235 word0(rv) = L | Bndry_mask1;
2236 word1(rv) = 0xffffffff;
2243 #ifndef ROUND_BIASED
2244 if (!(word1(rv) & LSB))
2248 dval(rv) += ulp(dval(rv));
2249 #ifndef ROUND_BIASED
2251 dval(rv) -= ulp(dval(rv));
2252 #ifndef Sudden_Underflow
2257 #ifdef Avoid_Underflow
2263 if ((aadj = ratio(delta, bs)) <= 2.) {
2266 else if (word1(rv) || word0(rv) & Bndry_mask) {
2267 #ifndef Sudden_Underflow
2268 if (word1(rv) == Tiny1 && !word0(rv))
2275 /* special case -- power of FLT_RADIX to be */
2276 /* rounded down... */
2278 if (aadj < 2./FLT_RADIX)
2279 aadj = 1./FLT_RADIX;
2287 aadj1 = dsign ? aadj : -aadj;
2288 #ifdef Check_FLT_ROUNDS
2290 case 2: /* towards +infinity */
2293 case 0: /* towards 0 */
2294 case 3: /* towards -infinity */
2298 if (Flt_Rounds == 0)
2300 #endif /*Check_FLT_ROUNDS*/
2302 y = word0(rv) & Exp_mask;
2304 /* Check for overflow */
2306 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2307 dval(rv0) = dval(rv);
2308 word0(rv) -= P*Exp_msk1;
2309 adj = aadj1 * ulp(dval(rv));
2311 if ((word0(rv) & Exp_mask) >=
2312 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2313 if (word0(rv0) == Big0 && word1(rv0) == Big1)
2320 word0(rv) += P*Exp_msk1;
2323 #ifdef Avoid_Underflow
2324 if (scale && y <= 2*P*Exp_msk1) {
2325 if (aadj <= 0x7fffffff) {
2326 if ((z = aadj) <= 0)
2329 aadj1 = dsign ? aadj : -aadj;
2331 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
2333 adj = aadj1 * ulp(dval(rv));
2336 #ifdef Sudden_Underflow
2337 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2338 dval(rv0) = dval(rv);
2339 word0(rv) += P*Exp_msk1;
2340 adj = aadj1 * ulp(dval(rv));
2343 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
2345 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2348 if (word0(rv0) == Tiny0
2349 && word1(rv0) == Tiny1)
2356 word0(rv) -= P*Exp_msk1;
2359 adj = aadj1 * ulp(dval(rv));
2362 #else /*Sudden_Underflow*/
2363 /* Compute adj so that the IEEE rounding rules will
2364 * correctly round rv + adj in some half-way cases.
2365 * If rv * ulp(rv) is denormalized (i.e.,
2366 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2367 * trouble from bits lost to denormalization;
2368 * example: 1.2e-307 .
2370 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
2371 aadj1 = (double)(int)(aadj + 0.5);
2375 adj = aadj1 * ulp(dval(rv));
2377 #endif /*Sudden_Underflow*/
2378 #endif /*Avoid_Underflow*/
2380 z = word0(rv) & Exp_mask;
2382 #ifdef Avoid_Underflow
2386 /* Can we stop now? */
2389 /* The tolerances below are conservative. */
2390 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2391 if (aadj < .4999999 || aadj > .5000001)
2394 else if (aadj < .4999999/FLT_RADIX)
2407 word0(rv0) = Exp_1 + (70 << Exp_shift);
2412 else if (!oldinexact)
2415 #ifdef Avoid_Underflow
2417 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
2419 dval(rv) *= dval(rv0);
2421 /* try to avoid the bug of testing an 8087 register value */
2422 if (word0(rv) == 0 && word1(rv) == 0)
2426 #endif /* Avoid_Underflow */
2428 if (inexact && !(word0(rv) & Exp_mask)) {
2429 /* set underflow bit */
2431 dval(rv0) *= dval(rv0);
2443 return sign ? -dval(rv) : dval(rv);
2449 (b, S) Bigint *b, *S;
2451 (Bigint *b, Bigint *S)
2455 ULong *bx, *bxe, q, *sx, *sxe;
2457 ULLong borrow, carry, y, ys;
2459 ULong borrow, carry, y, ys;
2467 /*debug*/ if (b->wds > n)
2468 /*debug*/ Bug("oversize b in quorem");
2476 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
2478 /*debug*/ if (q > 9)
2479 /*debug*/ Bug("oversized quotient in quorem");
2486 ys = *sx++ * (ULLong)q + carry;
2488 y = *bx - (ys & FFFFFFFF) - borrow;
2489 borrow = y >> 32 & (ULong)1;
2490 *bx++ = y & FFFFFFFF;
2494 ys = (si & 0xffff) * q + carry;
2495 zs = (si >> 16) * q + (ys >> 16);
2497 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2498 borrow = (y & 0x10000) >> 16;
2499 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2500 borrow = (z & 0x10000) >> 16;
2503 ys = *sx++ * q + carry;
2505 y = *bx - (ys & 0xffff) - borrow;
2506 borrow = (y & 0x10000) >> 16;
2514 while(--bxe > bx && !*bxe)
2519 if (cmp(b, S) >= 0) {
2529 y = *bx - (ys & FFFFFFFF) - borrow;
2530 borrow = y >> 32 & (ULong)1;
2531 *bx++ = y & FFFFFFFF;
2535 ys = (si & 0xffff) + carry;
2536 zs = (si >> 16) + (ys >> 16);
2538 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2539 borrow = (y & 0x10000) >> 16;
2540 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2541 borrow = (z & 0x10000) >> 16;
2546 y = *bx - (ys & 0xffff) - borrow;
2547 borrow = (y & 0x10000) >> 16;
2556 while(--bxe > bx && !*bxe)
2564 #ifndef MULTIPLE_THREADS
2565 static char *dtoa_result;
2579 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
2582 r = (int*)Balloc(k);
2585 #ifndef MULTIPLE_THREADS
2593 nrv_alloc(s, rve, n) char *s, **rve; int n;
2595 nrv_alloc(char *s, char **rve, int n)
2600 t = rv = rv_alloc(n);
2601 while((*t = *s++)) t++;
2607 /* freedtoa(s) must be used to free values s returned by dtoa
2608 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2609 * but for consistency with earlier versions of dtoa, it is optional
2610 * when MULTIPLE_THREADS is not defined.
2613 static void freedtoa (char *s);
2617 freedtoa(s) char *s;
2622 Bigint *b = (Bigint *)((int *)s - 1);
2623 b->maxwds = 1 << (b->k = *(int*)b);
2625 #ifndef MULTIPLE_THREADS
2626 if (s == dtoa_result)
2632 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2634 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2635 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2638 * 1. Rather than iterating, we use a simple numeric overestimate
2639 * to determine k = floor(log10(d)). We scale relevant
2640 * quantities using O(log2(k)) rather than O(k) multiplications.
2641 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2642 * try to generate digits strictly left to right. Instead, we
2643 * compute with fewer bits and propagate the carry if necessary
2644 * when rounding the final digit up. This is often faster.
2645 * 3. Under the assumption that input will be rounded nearest,
2646 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2647 * That is, we allow equality in stopping tests when the
2648 * round-nearest rule will give the same floating-point value
2649 * as would satisfaction of the stopping test with strict
2651 * 4. We remove common factors of powers of 2 from relevant
2653 * 5. When converting floating-point integers less than 1e16,
2654 * we use floating-point arithmetic rather than resorting
2655 * to multiple-precision integers.
2656 * 6. When asked to produce fewer than 15 digits, we first try
2657 * to get by with floating-point arithmetic; we resort to
2658 * multiple-precision integer arithmetic only if we cannot
2659 * guarantee that the floating-point calculation has given
2660 * the correctly rounded result. For k requested digits and
2661 * "uniformly" distributed input, the probability is
2662 * something like 10^(k-15) that we must resort to the Long
2669 (d, mode, ndigits, decpt, sign, rve)
2670 double d; int mode, ndigits, *decpt, *sign; char **rve;
2672 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
2675 /* Arguments ndigits, decpt, sign are similar to those
2676 of ecvt and fcvt; trailing zeros are suppressed from
2677 the returned string. If not null, *rve is set to point
2678 to the end of the return value. If d is +-Infinity or NaN,
2679 then *decpt is set to 9999.
2682 0 ==> shortest string that yields d when read in
2683 and rounded to nearest.
2684 1 ==> like 0, but with Steele & White stopping rule;
2685 e.g. with IEEE P754 arithmetic , mode 0 gives
2686 1e23 whereas mode 1 gives 9.999999999999999e22.
2687 2 ==> max(1,ndigits) significant digits. This gives a
2688 return value similar to that of ecvt, except
2689 that trailing zeros are suppressed.
2690 3 ==> through ndigits past the decimal point. This
2691 gives a return value similar to that from fcvt,
2692 except that trailing zeros are suppressed, and
2693 ndigits can be negative.
2694 4,5 ==> similar to 2 and 3, respectively, but (in
2695 round-nearest mode) with the tests of mode 0 to
2696 possibly return a shorter string that rounds to d.
2697 With IEEE arithmetic and compilation with
2698 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2699 as modes 2 and 3 when FLT_ROUNDS != 1.
2700 6-9 ==> Debugging modes similar to mode - 4: don't try
2701 fast floating-point estimate (if applicable).
2703 Values of mode other than 0-9 are treated as mode 0.
2705 Sufficient space is allocated to the return value
2706 to hold the suppressed trailing zeros.
2709 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
2710 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2711 spec_case, try_quick;
2713 #ifndef Sudden_Underflow
2717 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
2720 #ifdef Honor_FLT_ROUNDS
2724 int inexact, oldinexact;
2727 #ifndef MULTIPLE_THREADS
2729 freedtoa(dtoa_result);
2734 if (word0(d) & Sign_bit) {
2735 /* set sign for everything, including 0's and NaNs */
2737 word0(d) &= ~Sign_bit; /* clear sign bit */
2742 #if defined(IEEE_Arith) + defined(VAX)
2744 if ((word0(d) & Exp_mask) == Exp_mask)
2746 if (word0(d) == 0x8000)
2749 /* Infinity or NaN */
2752 if (!word1(d) && !(word0(d) & 0xfffff))
2753 return nrv_alloc("Infinity", rve, 8);
2755 return nrv_alloc("NaN", rve, 3);
2759 dval(d) += 0; /* normalize */
2763 return nrv_alloc("0", rve, 1);
2767 try_quick = oldinexact = get_inexact();
2770 #ifdef Honor_FLT_ROUNDS
2771 if ((rounding = Flt_Rounds) >= 2) {
2773 rounding = rounding == 2 ? 0 : 2;
2780 b = d2b(dval(d), &be, &bbits);
2781 #ifdef Sudden_Underflow
2782 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
2784 if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
2787 word0(d2) &= Frac_mask1;
2788 word0(d2) |= Exp_11;
2790 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
2794 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2795 * log10(x) = log(x) / log(10)
2796 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2797 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2799 * This suggests computing an approximation k to log10(d) by
2801 * k = (i - Bias)*0.301029995663981
2802 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2804 * We want k to be too large rather than too small.
2805 * The error in the first-order Taylor series approximation
2806 * is in our favor, so we just round up the constant enough
2807 * to compensate for any error in the multiplication of
2808 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2809 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2810 * adding 1e-13 to the constant term more than suffices.
2811 * Hence we adjust the constant term to 0.1760912590558.
2812 * (We could get a more accurate k by invoking log10,
2813 * but this is probably not worthwhile.)
2821 #ifndef Sudden_Underflow
2825 /* d is denormalized */
2827 i = bbits + be + (Bias + (P-1) - 1);
2828 x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32
2829 : word1(d) << 32 - i;
2831 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2832 i -= (Bias + (P-1) - 1) + 1;
2836 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2838 if (ds < 0. && ds != k)
2839 k--; /* want k = floor(ds) */
2841 if (k >= 0 && k <= Ten_pmax) {
2842 if (dval(d) < tens[k])
2865 if (mode < 0 || mode > 9)
2869 #ifdef Check_FLT_ROUNDS
2870 try_quick = Rounding == 1;
2874 #endif /*SET_INEXACT*/
2894 ilim = ilim1 = i = ndigits;
2900 i = ndigits + k + 1;
2906 s = s0 = rv_alloc(i);
2908 #ifdef Honor_FLT_ROUNDS
2909 if (mode > 1 && rounding != 1)
2913 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2915 /* Try to get by with floating-point arithmetic. */
2921 ieps = 2; /* conservative */
2926 /* prevent overflows */
2928 dval(d) /= bigtens[n_bigtens-1];
2931 for(; j; j >>= 1, i++)
2939 dval(d) *= tens[j1 & 0xf];
2940 for(j = j1 >> 4; j; j >>= 1, i++)
2943 dval(d) *= bigtens[i];
2946 if (k_check && dval(d) < 1. && ilim > 0) {
2954 dval(eps) = ieps*dval(d) + 7.;
2955 word0(eps) -= (P-1)*Exp_msk1;
2959 if (dval(d) > dval(eps))
2961 if (dval(d) < -dval(eps))
2965 #ifndef No_leftright
2967 /* Use Steele & White method of only
2968 * generating digits needed.
2970 dval(eps) = 0.5/tens[ilim-1] - dval(eps);
2974 *s++ = '0' + (int)L;
2975 if (dval(d) < dval(eps))
2977 if (1. - dval(d) < dval(eps))
2987 /* Generate ilim digits, then fix them up. */
2988 dval(eps) *= tens[ilim-1];
2989 for(i = 1;; i++, dval(d) *= 10.) {
2990 L = (Long)(dval(d));
2991 if (!(dval(d) -= L))
2993 *s++ = '0' + (int)L;
2995 if (dval(d) > 0.5 + dval(eps))
2997 else if (dval(d) < 0.5 - dval(eps)) {
3005 #ifndef No_leftright
3015 /* Do we have a "small" integer? */
3017 if (be >= 0 && k <= Int_max) {
3020 if (ndigits < 0 && ilim <= 0) {
3022 if (ilim < 0 || dval(d) <= 5*ds)
3026 for(i = 1;; i++, dval(d) *= 10.) {
3027 L = (Long)(dval(d) / ds);
3029 #ifdef Check_FLT_ROUNDS
3030 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
3036 *s++ = '0' + (int)L;
3044 #ifdef Honor_FLT_ROUNDS
3048 case 2: goto bump_up;
3052 if (dval(d) > ds || dval(d) == ds && L & 1) {
3073 #ifndef Sudden_Underflow
3074 denorm ? be + (Bias + (P-1) - 1 + 1) :
3077 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
3085 if (m2 > 0 && s2 > 0) {
3086 i = m2 < s2 ? m2 : s2;
3094 mhi = pow5mult(mhi, m5);
3103 b = pow5mult(b, b5);
3107 S = pow5mult(S, s5);
3109 /* Check for special case that d is a normalized power of 2. */
3112 if ((mode < 2 || leftright)
3113 #ifdef Honor_FLT_ROUNDS
3117 if (!word1(d) && !(word0(d) & Bndry_mask)
3118 #ifndef Sudden_Underflow
3119 && word0(d) & (Exp_mask & ~Exp_msk1)
3122 /* The special case */
3129 /* Arrange for convenient computation of quotients:
3130 * shift left if necessary so divisor has 4 leading 0 bits.
3132 * Perhaps we should just compute leading 28 bits of S once
3133 * and for all and pass them and a shift to quorem, so it
3134 * can do shifts and ors to compute the numerator for q.
3137 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
3140 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
3162 b = multadd(b, 10, 0); /* we botched the k estimate */
3164 mhi = multadd(mhi, 10, 0);
3168 if (ilim <= 0 && (mode == 3 || mode == 5)) {
3169 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
3170 /* no digits, fcvt style */
3182 mhi = lshift(mhi, m2);
3184 /* Compute mlo -- check for special case
3185 * that d is a normalized power of 2.
3190 mhi = Balloc(mhi->k);
3192 mhi = lshift(mhi, Log2P);
3196 dig = quorem(b,S) + '0';
3197 /* Do we yet have the shortest decimal string
3198 * that will round to d?
3201 delta = diff(S, mhi);
3202 j1 = delta->sign ? 1 : cmp(b, delta);
3204 #ifndef ROUND_BIASED
3205 if (j1 == 0 && mode != 1 && !(word1(d) & 1)
3206 #ifdef Honor_FLT_ROUNDS
3215 else if (!b->x[0] && b->wds <= 1)
3222 if (j < 0 || j == 0 && mode != 1
3223 #ifndef ROUND_BIASED
3227 if (!b->x[0] && b->wds <= 1) {
3233 #ifdef Honor_FLT_ROUNDS
3236 case 0: goto accept_dig;
3237 case 2: goto keep_dig;
3239 #endif /*Honor_FLT_ROUNDS*/
3243 if ((j1 > 0 || j1 == 0 && dig & 1)
3252 #ifdef Honor_FLT_ROUNDS
3256 if (dig == '9') { /* possible if i == 1 */
3264 #ifdef Honor_FLT_ROUNDS
3270 b = multadd(b, 10, 0);
3272 mlo = mhi = multadd(mhi, 10, 0);
3274 mlo = multadd(mlo, 10, 0);
3275 mhi = multadd(mhi, 10, 0);
3281 *s++ = dig = quorem(b,S) + '0';
3282 if (!b->x[0] && b->wds <= 1) {
3290 b = multadd(b, 10, 0);
3293 /* Round off last digit */
3295 #ifdef Honor_FLT_ROUNDS
3297 case 0: goto trimzeros;
3298 case 2: goto roundoff;
3303 if (j > 0 || j == 0 && dig & 1) {
3321 if (mlo && mlo != mhi)
3329 word0(d) = Exp_1 + (70 << Exp_shift);
3334 else if (!oldinexact)