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 __BIG_ENDIAN__ || defined _BIG_ENDIAN
178 #elif defined __LITTLE_ENDIAN__ || defined _LITTLE_ENDIAN
182 #if defined(TARGET_X86) || defined(mips) && defined(MIPSEL) || defined (__arm__) || defined(__aarch64__)
186 #elif defined(TARGET_AMD64) || defined(__alpha__)
190 #elif defined(__ia64)
198 #elif defined(__hppa)
203 #warning byte order unknown, assuming big endian
210 #define ULong guint32
214 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
227 extern char *MALLOC();
229 extern void *MALLOC(size_t);
232 #define MALLOC malloc
235 #define Omit_Private_Memory
236 #ifndef Omit_Private_Memory
238 #define PRIVATE_MEM 2304
240 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
241 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
245 #undef Avoid_Underflow
259 #define DBL_MAX_10_EXP 308
260 #define DBL_MAX_EXP 1024
262 #endif /*IEEE_Arith*/
266 #define DBL_MAX_10_EXP 75
267 #define DBL_MAX_EXP 63
269 #define DBL_MAX 7.2370055773322621e+75
274 #define DBL_MAX_10_EXP 38
275 #define DBL_MAX_EXP 127
277 #define DBL_MAX 1.7014118346046923e+38
281 #define LONG_MAX 2147483647
284 #else /* ifndef Bad_float_h */
286 #endif /* Bad_float_h */
298 #define CONST /* blank */
304 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
305 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
308 typedef union { double d; ULong L[2]; } U;
313 #define word0(x) ((ULong *)&x)[1]
314 #define word1(x) ((ULong *)&x)[0]
316 #define word0(x) ((ULong *)&x)[0]
317 #define word1(x) ((ULong *)&x)[1]
321 #define word0(x) ((U*)&x)->L[1]
322 #define word1(x) ((U*)&x)->L[0]
324 #define word0(x) ((U*)&x)->L[0]
325 #define word1(x) ((U*)&x)->L[1]
327 #define dval(x) ((U*)&x)->d
330 /* The following definition of Storeinc is appropriate for MIPS processors.
331 * An alternative that might be better on some machines is
332 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
334 #if defined(IEEE_8087) + defined(VAX)
335 #define Storeinc(a,b,c) do { (((unsigned short *)a)[1] = (unsigned short)b, \
336 ((unsigned short *)a)[0] = (unsigned short)c, a++) } while (0)
338 #define Storeinc(a,b,c) do { (((unsigned short *)a)[0] = (unsigned short)b, \
339 ((unsigned short *)a)[1] = (unsigned short)c, a++) } while (0)
342 /* #define P DBL_MANT_DIG */
343 /* Ten_pmax = floor(P*log(2)/log(5)) */
344 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
345 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
346 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
350 #define Exp_shift1 20
351 #define Exp_msk1 0x100000
352 #define Exp_msk11 0x100000
353 #define Exp_mask 0x7ff00000
357 #define Exp_1 0x3ff00000
358 #define Exp_11 0x3ff00000
360 #define Frac_mask 0xfffff
361 #define Frac_mask1 0xfffff
364 #define Bndry_mask 0xfffff
365 #define Bndry_mask1 0xfffff
367 #define Sign_bit 0x80000000
373 #ifndef NO_IEEE_Scale
374 #define Avoid_Underflow
375 #ifdef Flush_Denorm /* debugging option */
376 #undef Sudden_Underflow
382 #define Flt_Rounds FLT_ROUNDS
386 #endif /*Flt_Rounds*/
388 #ifdef Honor_FLT_ROUNDS
389 #define Rounding rounding
390 #undef Check_FLT_ROUNDS
391 #define Check_FLT_ROUNDS
393 #define Rounding Flt_Rounds
396 #else /* ifndef IEEE_Arith */
397 #undef Check_FLT_ROUNDS
398 #undef Honor_FLT_ROUNDS
400 #undef Sudden_Underflow
401 #define Sudden_Underflow
406 #define Exp_shift1 24
407 #define Exp_msk1 0x1000000
408 #define Exp_msk11 0x1000000
409 #define Exp_mask 0x7f000000
412 #define Exp_1 0x41000000
413 #define Exp_11 0x41000000
414 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
415 #define Frac_mask 0xffffff
416 #define Frac_mask1 0xffffff
419 #define Bndry_mask 0xefffff
420 #define Bndry_mask1 0xffffff
422 #define Sign_bit 0x80000000
424 #define Tiny0 0x100000
433 #define Exp_msk1 0x80
434 #define Exp_msk11 0x800000
435 #define Exp_mask 0x7f80
438 #define Exp_1 0x40800000
439 #define Exp_11 0x4080
441 #define Frac_mask 0x7fffff
442 #define Frac_mask1 0xffff007f
445 #define Bndry_mask 0xffff007f
446 #define Bndry_mask1 0xffff007f
448 #define Sign_bit 0x8000
454 #endif /* IBM, VAX */
455 #endif /* IEEE_Arith */
462 #define rounded_product(a,b) a = rnd_prod(a, b)
463 #define rounded_quotient(a,b) a = rnd_quot(a, b)
465 extern double rnd_prod(), rnd_quot();
467 extern double rnd_prod(double, double), rnd_quot(double, double);
470 #define rounded_product(a,b) a *= b
471 #define rounded_quotient(a,b) a /= b
474 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
475 #define Big1 0xffffffff
482 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
484 #define FFFFFFFF 0xffffffffUL
491 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
492 * This makes some inner loops simpler and sometimes saves work
493 * during multiplications, but it often seems to make things slightly
494 * slower. Hence the default is now to store 32 bits per Long.
497 #else /* long long available */
499 #define Llong long long
502 #define ULLong unsigned Llong
504 #endif /* NO_LONG_LONG */
506 #ifndef MULTIPLE_THREADS
507 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
508 #define FREE_DTOA_LOCK(n) /*nothing*/
514 extern "C" double strtod(const char *s00, char **se);
515 extern "C" char *dtoa(double d, int mode, int ndigits,
516 int *decpt, int *sign, char **rve);
522 int k, maxwds, sign, wds;
526 typedef struct Bigint Bigint;
528 static Bigint *freelist[Kmax+1];
540 #ifndef Omit_Private_Memory
544 ACQUIRE_DTOA_LOCK(0);
545 if ((rv = freelist[k])) {
546 freelist[k] = rv->next;
550 #ifdef Omit_Private_Memory
551 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
553 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
555 if (pmem_next - private_mem + len <= PRIVATE_mem) {
556 rv = (Bigint*)pmem_next;
560 rv = (Bigint*)MALLOC(len*sizeof(double));
566 rv->sign = rv->wds = 0;
578 #ifdef Omit_Private_Memory
582 ACQUIRE_DTOA_LOCK(0);
583 v->next = freelist[v->k];
590 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
591 y->wds*sizeof(Long) + 2*sizeof(int))
596 (b, m, a) Bigint *b; int m, a;
598 (Bigint *b, int m, int a) /* multiply by m and add a */
619 y = *x * (ULLong)m + carry;
625 y = (xi & 0xffff) * m + carry;
626 z = (xi >> 16) * m + (y >> 16);
628 *x++ = (z << 16) + (y & 0xffff);
638 if (wds >= b->maxwds) {
653 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
655 (CONST char *s, int nd0, int nd, ULong y9)
663 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
670 b->x[0] = y9 & 0xffff;
671 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
677 do b = multadd(b, 10, *s++ - '0');
684 b = multadd(b, 10, *s++ - '0');
691 (x) register ULong x;
698 if (!(x & 0xffff0000)) {
702 if (!(x & 0xff000000)) {
706 if (!(x & 0xf0000000)) {
710 if (!(x & 0xc0000000)) {
714 if (!(x & 0x80000000)) {
716 if (!(x & 0x40000000))
731 register ULong x = *y;
789 (a, b) Bigint *a, *b;
791 (Bigint *a, Bigint *b)
796 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
807 if (a->wds < b->wds) {
819 for(x = c->x, xa = x + wc; x < xa; x++)
827 for(; xb < xbe; xc0++) {
833 z = *x++ * (ULLong)y + *xc + carry;
835 *xc++ = z & FFFFFFFF;
843 for(; xb < xbe; xb++, xc0++) {
844 if (y = *xb & 0xffff) {
849 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
851 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
864 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
867 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
875 for(; xb < xbe; xc0++) {
881 z = *x++ * y + *xc + carry;
891 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
901 (b, k) Bigint *b; int k;
906 Bigint *b1, *p5, *p51;
908 static int p05[3] = { 5, 25, 125 };
911 b = multadd(b, p05[i-1], 0);
917 #ifdef MULTIPLE_THREADS
918 ACQUIRE_DTOA_LOCK(1);
937 if (!(p51 = p5->next)) {
938 #ifdef MULTIPLE_THREADS
939 ACQUIRE_DTOA_LOCK(1);
940 if (!(p51 = p5->next)) {
941 p51 = p5->next = mult(p5,p5);
946 p51 = p5->next = mult(p5,p5);
958 (b, k) Bigint *b; int k;
965 ULong *x, *x1, *xe, z;
974 for(i = b->maxwds; n1 > i; i <<= 1)
978 for(i = 0; i < n; i++)
999 *x1++ = *x << k & 0xffff | z;
1018 (a, b) Bigint *a, *b;
1020 (Bigint *a, Bigint *b)
1023 ULong *xa, *xa0, *xb, *xb0;
1029 if (i > 1 && !a->x[i-1])
1030 Bug("cmp called with a->x[a->wds-1] == 0");
1031 if (j > 1 && !b->x[j-1])
1032 Bug("cmp called with b->x[b->wds-1] == 0");
1042 return *xa < *xb ? -1 : 1;
1052 (a, b) Bigint *a, *b;
1054 (Bigint *a, Bigint *b)
1059 ULong *xa, *xae, *xb, *xbe, *xc;
1096 y = (ULLong)*xa++ - *xb++ - borrow;
1097 borrow = y >> 32 & (ULong)1;
1098 *xc++ = y & FFFFFFFF;
1103 borrow = y >> 32 & (ULong)1;
1104 *xc++ = y & FFFFFFFF;
1109 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1110 borrow = (y & 0x10000) >> 16;
1111 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1112 borrow = (z & 0x10000) >> 16;
1117 y = (*xa & 0xffff) - borrow;
1118 borrow = (y & 0x10000) >> 16;
1119 z = (*xa++ >> 16) - borrow;
1120 borrow = (z & 0x10000) >> 16;
1125 y = *xa++ - *xb++ - borrow;
1126 borrow = (y & 0x10000) >> 16;
1132 borrow = (y & 0x10000) >> 16;
1154 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1155 #ifndef Avoid_Underflow
1156 #ifndef Sudden_Underflow
1165 #ifndef Avoid_Underflow
1166 #ifndef Sudden_Underflow
1169 L = -L >> Exp_shift;
1170 if (L < Exp_shift) {
1171 word0(a) = 0x80000 >> L;
1177 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
1188 (a, e) Bigint *a; int *e;
1193 ULong *xa, *xa0, w, y, z;
1207 if (!y) Bug("zero y in b2d");
1213 d0 = Exp_1 | (y >> (Ebits - k));
1214 w = xa > xa0 ? *--xa : 0;
1215 d1 = y << ((32-Ebits) + k) | (w >> (Ebits - k));
1218 z = xa > xa0 ? *--xa : 0;
1220 d0 = Exp_1 | y << k | (z >> (32 - k));
1221 y = xa > xa0 ? *--xa : 0;
1222 d1 = z << k | (y >> (32 - k));
1229 if (k < Ebits + 16) {
1230 z = xa > xa0 ? *--xa : 0;
1231 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1232 w = xa > xa0 ? *--xa : 0;
1233 y = xa > xa0 ? *--xa : 0;
1234 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1237 z = xa > xa0 ? *--xa : 0;
1238 w = xa > xa0 ? *--xa : 0;
1240 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1241 y = xa > xa0 ? *--xa : 0;
1242 d1 = w << k + 16 | y << k;
1246 word0(d) = d0 >> 16 | d0 << 16;
1247 word1(d) = d1 >> 16 | d1 << 16;
1258 (d, e, bits) double d; int *e, *bits;
1260 (double d, int *e, int *bits)
1266 #ifndef Sudden_Underflow
1271 d0 = word0(d) >> 16 | word0(d) << 16;
1272 d1 = word1(d) >> 16 | word1(d) << 16;
1286 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1287 #ifdef Sudden_Underflow
1288 de = (int)(d0 >> Exp_shift);
1293 if ((de = (int)(d0 >> Exp_shift)))
1298 if ((k = lo0bits(&y))) {
1299 x[0] = y | (z << (32 - k));
1304 #ifndef Sudden_Underflow
1307 b->wds = (x[1] = z) ? 2 : 1;
1312 Bug("Zero passed to d2b");
1316 #ifndef Sudden_Underflow
1324 if (k = lo0bits(&y))
1326 x[0] = y | z << 32 - k & 0xffff;
1327 x[1] = z >> k - 16 & 0xffff;
1333 x[1] = y >> 16 | z << 16 - k & 0xffff;
1334 x[2] = z >> k & 0xffff;
1349 Bug("Zero passed to d2b");
1367 #ifndef Sudden_Underflow
1371 *e = (de - Bias - (P-1) << 2) + k;
1372 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1374 *e = de - Bias - (P-1) + k;
1377 #ifndef Sudden_Underflow
1380 *e = de - Bias - (P-1) + 1 + k;
1382 *bits = 32*i - hi0bits(x[i-1]);
1384 *bits = (i+2)*16 - hi0bits(x[i]);
1396 (a, b) Bigint *a, *b;
1398 (Bigint *a, Bigint *b)
1404 dval(da) = b2d(a, &ka);
1405 dval(db) = b2d(b, &kb);
1407 k = ka - kb + 32*(a->wds - b->wds);
1409 k = ka - kb + 16*(a->wds - b->wds);
1413 word0(da) += (k >> 2)*Exp_msk1;
1419 word0(db) += (k >> 2)*Exp_msk1;
1425 word0(da) += k*Exp_msk1;
1428 word0(db) += k*Exp_msk1;
1431 return dval(da) / dval(db);
1436 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1437 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1446 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1447 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1448 #ifdef Avoid_Underflow
1449 9007199254740992.*9007199254740992.e-256
1450 /* = 2^106 * 1e-53 */
1455 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1456 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1457 #define Scale_Bit 0x10
1461 bigtens[] = { 1e16, 1e32, 1e64 };
1462 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1465 bigtens[] = { 1e16, 1e32 };
1466 static CONST double tinytens[] = { 1e-16, 1e-32 };
1478 #define NAN_WORD0 0x7ff80000
1488 (sp, t) char **sp, *t;
1490 (CONST char **sp, char *t)
1494 CONST char *s = *sp;
1497 if ((c = *++s) >= 'A' && c <= 'Z')
1510 (rvp, sp) double *rvp; CONST char **sp;
1512 (double *rvp, CONST char **sp)
1517 int havedig, udx0, xshift;
1520 havedig = xshift = 0;
1523 while(c = *(CONST unsigned char*)++s) {
1524 if (c >= '0' && c <= '9')
1526 else if (c >= 'a' && c <= 'f')
1528 else if (c >= 'A' && c <= 'F')
1530 else if (c <= ' ') {
1531 if (udx0 && havedig) {
1537 else if (/*(*/ c == ')' && havedig) {
1542 return; /* invalid form: don't change *sp */
1550 x[0] = (x[0] << 4) | (x[1] >> 28);
1551 x[1] = (x[1] << 4) | c;
1553 if ((x[0] &= 0xfffff) || x[1]) {
1554 word0(*rvp) = Exp_mask | x[0];
1558 #endif /*No_Hex_NaN*/
1559 #endif /* INFNAN_CHECK */
1562 * LOCKING: This is not thread-safe, since the locking macros are defined as no-ops,
1563 * the caller should lock.
1569 (s00, se) CONST char *s00; char **se;
1571 (CONST char *s00, char **se)
1574 #ifdef Avoid_Underflow
1577 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1578 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1579 CONST char *s, *s0, *s1;
1580 double aadj, aadj1, adj, rv, rv0;
1583 Bigint *bb = NULL, *bb1, *bd = NULL, *bd0, *bs = NULL, *delta = NULL;
1585 int inexact, oldinexact;
1587 #ifdef Honor_FLT_ROUNDS
1594 sign = nz0 = nz = 0;
1596 for(s = s00;;s++) switch(*s) {
1619 while(*++s == '0') ;
1625 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1632 s1 = localeconv()->decimal_point;
1653 for(; c == '0'; c = *++s)
1655 if (c > '0' && c <= '9') {
1663 for(; c >= '0' && c <= '9'; c = *++s) {
1668 for(i = 1; i < nz; i++)
1671 else if (nd <= DBL_DIG + 1)
1675 else if (nd <= DBL_DIG + 1)
1683 if (c == 'e' || c == 'E') {
1684 if (!nd && !nz && !nz0) {
1695 if (c >= '0' && c <= '9') {
1698 if (c > '0' && c <= '9') {
1701 while((c = *++s) >= '0' && c <= '9')
1703 if (s - s1 > 8 || L > 19999)
1704 /* Avoid confusion from exponents
1705 * so large that e might overflow.
1707 e = 19999; /* safe for 16 bit ints */
1722 /* Check for Nan and Infinity */
1726 if (match(&s,"nf")) {
1728 if (!match(&s,"inity"))
1730 word0(rv) = 0x7ff00000;
1737 if (match(&s, "an")) {
1738 word0(rv) = NAN_WORD0;
1739 word1(rv) = NAN_WORD1;
1741 if (*s == '(') /*)*/
1747 #endif /* INFNAN_CHECK */
1756 /* Now we have nd0 digits, starting at s0, followed by a
1757 * decimal point, followed by nd-nd0 digits. The number we're
1758 * after is the integer represented by those digits times
1763 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1768 oldinexact = get_inexact();
1770 dval(rv) = tens[k - 9] * dval(rv) + z;
1774 #ifndef RND_PRODQUOT
1775 #ifndef Honor_FLT_ROUNDS
1783 if (e <= Ten_pmax) {
1785 goto vax_ovfl_check;
1787 #ifdef Honor_FLT_ROUNDS
1788 /* round correctly FLT_ROUNDS = 2 or 3 */
1794 /* rv = */ rounded_product(dval(rv), tens[e]);
1799 if (e <= Ten_pmax + i) {
1800 /* A fancier test would sometimes let us do
1801 * this for larger i values.
1803 #ifdef Honor_FLT_ROUNDS
1804 /* round correctly FLT_ROUNDS = 2 or 3 */
1811 dval(rv) *= tens[i];
1813 /* VAX exponent range is so narrow we must
1814 * worry about overflow here...
1817 word0(rv) -= P*Exp_msk1;
1818 /* rv = */ rounded_product(dval(rv), tens[e]);
1819 if ((word0(rv) & Exp_mask)
1820 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1822 word0(rv) += P*Exp_msk1;
1824 /* rv = */ rounded_product(dval(rv), tens[e]);
1829 #ifndef Inaccurate_Divide
1830 else if (e >= -Ten_pmax) {
1831 #ifdef Honor_FLT_ROUNDS
1832 /* round correctly FLT_ROUNDS = 2 or 3 */
1838 /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1849 oldinexact = get_inexact();
1851 #ifdef Avoid_Underflow
1854 #ifdef Honor_FLT_ROUNDS
1855 if ((rounding = Flt_Rounds) >= 2) {
1857 rounding = rounding == 2 ? 0 : 2;
1863 #endif /*IEEE_Arith*/
1865 /* Get starting approximation = rv * 10**e1 */
1869 dval(rv) *= tens[i];
1871 if (e1 > DBL_MAX_10_EXP) {
1876 /* Can't trust HUGE_VAL */
1878 #ifdef Honor_FLT_ROUNDS
1880 case 0: /* toward 0 */
1881 case 3: /* toward -infinity */
1886 word0(rv) = Exp_mask;
1889 #else /*Honor_FLT_ROUNDS*/
1890 word0(rv) = Exp_mask;
1892 #endif /*Honor_FLT_ROUNDS*/
1894 /* set overflow bit */
1896 dval(rv0) *= dval(rv0);
1898 #else /*IEEE_Arith*/
1901 #endif /*IEEE_Arith*/
1907 for(j = 0; e1 > 1; j++, e1 >>= 1)
1909 dval(rv) *= bigtens[j];
1910 /* The last multiplication could overflow. */
1911 word0(rv) -= P*Exp_msk1;
1912 dval(rv) *= bigtens[j];
1913 if ((z = word0(rv) & Exp_mask)
1914 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1916 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1917 /* set to largest number */
1918 /* (Can't trust DBL_MAX) */
1923 word0(rv) += P*Exp_msk1;
1929 dval(rv) /= tens[i];
1931 if (e1 >= 1 << n_bigtens)
1933 #ifdef Avoid_Underflow
1936 for(j = 0; e1 > 0; j++, e1 >>= 1)
1938 dval(rv) *= tinytens[j];
1939 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
1940 >> Exp_shift)) > 0) {
1941 /* scaled rv is denormal; zap j low bits */
1945 word0(rv) = (P+2)*Exp_msk1;
1947 word0(rv) &= 0xffffffff << (j-32);
1950 word1(rv) &= 0xffffffff << j;
1953 for(j = 0; e1 > 1; j++, e1 >>= 1)
1955 dval(rv) *= tinytens[j];
1956 /* The last multiplication could underflow. */
1957 dval(rv0) = dval(rv);
1958 dval(rv) *= tinytens[j];
1960 dval(rv) = 2.*dval(rv0);
1961 dval(rv) *= tinytens[j];
1973 #ifndef Avoid_Underflow
1976 /* The refinement below will clean
1977 * this approximation up.
1984 /* Now the hard part -- adjusting rv to the correct value.*/
1986 /* Put digits into bd: true value = bd * 10^e */
1988 bd0 = s2b(s0, nd0, nd, y);
1991 bd = Balloc(bd0->k);
1993 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
2009 #ifdef Honor_FLT_ROUNDS
2013 #ifdef Avoid_Underflow
2015 i = j + bbbits - 1; /* logb(rv) */
2016 if (i < Emin) /* denormal */
2020 #else /*Avoid_Underflow*/
2021 #ifdef Sudden_Underflow
2023 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
2027 #else /*Sudden_Underflow*/
2029 i = j + bbbits - 1; /* logb(rv) */
2030 if (i < Emin) /* denormal */
2034 #endif /*Sudden_Underflow*/
2035 #endif /*Avoid_Underflow*/
2038 #ifdef Avoid_Underflow
2041 i = bb2 < bd2 ? bb2 : bd2;
2050 bs = pow5mult(bs, bb5);
2056 bb = lshift(bb, bb2);
2058 bd = pow5mult(bd, bd5);
2060 bd = lshift(bd, bd2);
2062 bs = lshift(bs, bs2);
2063 delta = diff(bb, bd);
2064 dsign = delta->sign;
2067 #ifdef Honor_FLT_ROUNDS
2068 if (rounding != 1) {
2070 /* Error is less than an ulp */
2071 if (!delta->x[0] && delta->wds <= 1) {
2087 && !(word0(rv) & Frac_mask)) {
2088 y = word0(rv) & Exp_mask;
2089 #ifdef Avoid_Underflow
2090 if (!scale || y > 2*P*Exp_msk1)
2095 delta = lshift(delta,Log2P);
2096 if (cmp(delta, bs) <= 0)
2101 #ifdef Avoid_Underflow
2102 if (scale && (y = word0(rv) & Exp_mask)
2104 word0(adj) += (2*P+1)*Exp_msk1 - y;
2106 #ifdef Sudden_Underflow
2107 if ((word0(rv) & Exp_mask) <=
2109 word0(rv) += P*Exp_msk1;
2110 dval(rv) += adj*ulp(dval(rv));
2111 word0(rv) -= P*Exp_msk1;
2114 #endif /*Sudden_Underflow*/
2115 #endif /*Avoid_Underflow*/
2116 dval(rv) += adj*ulp(dval(rv));
2120 adj = ratio(delta, bs);
2123 if (adj <= 0x7ffffffe) {
2124 /* adj = rounding ? ceil(adj) : floor(adj); */
2127 if (!((rounding>>1) ^ dsign))
2132 #ifdef Avoid_Underflow
2133 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2134 word0(adj) += (2*P+1)*Exp_msk1 - y;
2136 #ifdef Sudden_Underflow
2137 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2138 word0(rv) += P*Exp_msk1;
2139 adj *= ulp(dval(rv));
2144 word0(rv) -= P*Exp_msk1;
2147 #endif /*Sudden_Underflow*/
2148 #endif /*Avoid_Underflow*/
2149 adj *= ulp(dval(rv));
2156 #endif /*Honor_FLT_ROUNDS*/
2159 /* Error is less than half an ulp -- check for
2160 * special case of mantissa a power of two.
2162 if (dsign || word1(rv) || word0(rv) & Bndry_mask
2164 #ifdef Avoid_Underflow
2165 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2167 || (word0(rv) & Exp_mask) <= Exp_msk1
2172 if (!delta->x[0] && delta->wds <= 1)
2177 if (!delta->x[0] && delta->wds <= 1) {
2184 delta = lshift(delta,Log2P);
2185 if (cmp(delta, bs) > 0)
2190 /* exactly half-way between */
2192 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2194 #ifdef Avoid_Underflow
2195 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2196 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2199 /*boundary case -- increment exponent*/
2200 word0(rv) = (word0(rv) & Exp_mask)
2207 #ifdef Avoid_Underflow
2213 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2215 /* boundary case -- decrement exponent */
2216 #ifdef Sudden_Underflow /*{{*/
2217 L = word0(rv) & Exp_mask;
2221 #ifdef Avoid_Underflow
2222 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
2225 #endif /*Avoid_Underflow*/
2229 #else /*Sudden_Underflow}{*/
2230 #ifdef Avoid_Underflow
2232 L = word0(rv) & Exp_mask;
2233 if (L <= (2*P+1)*Exp_msk1) {
2234 if (L > (P+2)*Exp_msk1)
2235 /* round even ==> */
2238 /* rv = smallest denormal */
2242 #endif /*Avoid_Underflow*/
2243 L = (word0(rv) & Exp_mask) - Exp_msk1;
2244 #endif /*Sudden_Underflow}}*/
2245 word0(rv) = L | Bndry_mask1;
2246 word1(rv) = 0xffffffff;
2253 #ifndef ROUND_BIASED
2254 if (!(word1(rv) & LSB))
2258 dval(rv) += ulp(dval(rv));
2259 #ifndef ROUND_BIASED
2261 dval(rv) -= ulp(dval(rv));
2262 #ifndef Sudden_Underflow
2267 #ifdef Avoid_Underflow
2273 if ((aadj = ratio(delta, bs)) <= 2.) {
2276 else if (word1(rv) || word0(rv) & Bndry_mask) {
2277 #ifndef Sudden_Underflow
2278 if (word1(rv) == Tiny1 && !word0(rv))
2285 /* special case -- power of FLT_RADIX to be */
2286 /* rounded down... */
2288 if (aadj < 2./FLT_RADIX)
2289 aadj = 1./FLT_RADIX;
2297 aadj1 = dsign ? aadj : -aadj;
2298 #ifdef Check_FLT_ROUNDS
2300 case 2: /* towards +infinity */
2303 case 0: /* towards 0 */
2304 case 3: /* towards -infinity */
2308 if (Flt_Rounds == 0)
2310 #endif /*Check_FLT_ROUNDS*/
2312 y = word0(rv) & Exp_mask;
2314 /* Check for overflow */
2316 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2317 dval(rv0) = dval(rv);
2318 word0(rv) -= P*Exp_msk1;
2319 adj = aadj1 * ulp(dval(rv));
2321 if ((word0(rv) & Exp_mask) >=
2322 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2323 if (word0(rv0) == Big0 && word1(rv0) == Big1)
2330 word0(rv) += P*Exp_msk1;
2333 #ifdef Avoid_Underflow
2334 if (scale && y <= 2*P*Exp_msk1) {
2335 if (aadj <= 0x7fffffff) {
2336 if ((z = aadj) <= 0)
2339 aadj1 = dsign ? aadj : -aadj;
2341 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
2343 adj = aadj1 * ulp(dval(rv));
2346 #ifdef Sudden_Underflow
2347 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2348 dval(rv0) = dval(rv);
2349 word0(rv) += P*Exp_msk1;
2350 adj = aadj1 * ulp(dval(rv));
2353 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
2355 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2358 if (word0(rv0) == Tiny0
2359 && word1(rv0) == Tiny1)
2366 word0(rv) -= P*Exp_msk1;
2369 adj = aadj1 * ulp(dval(rv));
2372 #else /*Sudden_Underflow*/
2373 /* Compute adj so that the IEEE rounding rules will
2374 * correctly round rv + adj in some half-way cases.
2375 * If rv * ulp(rv) is denormalized (i.e.,
2376 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2377 * trouble from bits lost to denormalization;
2378 * example: 1.2e-307 .
2380 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
2381 aadj1 = (double)(int)(aadj + 0.5);
2385 adj = aadj1 * ulp(dval(rv));
2387 #endif /*Sudden_Underflow*/
2388 #endif /*Avoid_Underflow*/
2390 z = word0(rv) & Exp_mask;
2392 #ifdef Avoid_Underflow
2396 /* Can we stop now? */
2399 /* The tolerances below are conservative. */
2400 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2401 if (aadj < .4999999 || aadj > .5000001)
2404 else if (aadj < .4999999/FLT_RADIX)
2417 word0(rv0) = Exp_1 + (70 << Exp_shift);
2422 else if (!oldinexact)
2425 #ifdef Avoid_Underflow
2427 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
2429 dval(rv) *= dval(rv0);
2431 /* try to avoid the bug of testing an 8087 register value */
2432 if (word0(rv) == 0 && word1(rv) == 0)
2436 #endif /* Avoid_Underflow */
2438 if (inexact && !(word0(rv) & Exp_mask)) {
2439 /* set underflow bit */
2441 dval(rv0) *= dval(rv0);
2453 return sign ? -dval(rv) : dval(rv);
2460 (b, S) Bigint *b, *S;
2462 (Bigint *b, Bigint *S)
2466 ULong *bx, *bxe, q, *sx, *sxe;
2468 ULLong borrow, carry, y, ys;
2470 ULong borrow, carry, y, ys;
2478 /*debug*/ if (b->wds > n)
2479 /*debug*/ Bug("oversize b in quorem");
2487 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
2489 /*debug*/ if (q > 9)
2490 /*debug*/ Bug("oversized quotient in quorem");
2497 ys = *sx++ * (ULLong)q + carry;
2499 y = *bx - (ys & FFFFFFFF) - borrow;
2500 borrow = y >> 32 & (ULong)1;
2501 *bx++ = y & FFFFFFFF;
2505 ys = (si & 0xffff) * q + carry;
2506 zs = (si >> 16) * q + (ys >> 16);
2508 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2509 borrow = (y & 0x10000) >> 16;
2510 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2511 borrow = (z & 0x10000) >> 16;
2514 ys = *sx++ * q + carry;
2516 y = *bx - (ys & 0xffff) - borrow;
2517 borrow = (y & 0x10000) >> 16;
2525 while(--bxe > bx && !*bxe)
2530 if (cmp(b, S) >= 0) {
2540 y = *bx - (ys & FFFFFFFF) - borrow;
2541 borrow = y >> 32 & (ULong)1;
2542 *bx++ = y & FFFFFFFF;
2546 ys = (si & 0xffff) + carry;
2547 zs = (si >> 16) + (ys >> 16);
2549 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2550 borrow = (y & 0x10000) >> 16;
2551 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2552 borrow = (z & 0x10000) >> 16;
2557 y = *bx - (ys & 0xffff) - borrow;
2558 borrow = (y & 0x10000) >> 16;
2567 while(--bxe > bx && !*bxe)
2576 #ifndef MULTIPLE_THREADS
2577 static char *dtoa_result;
2592 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
2595 r = (int*)Balloc(k);
2598 #ifndef MULTIPLE_THREADS
2606 nrv_alloc(s, rve, n) char *s, **rve; int n;
2608 nrv_alloc(char *s, char **rve, int n)
2613 t = rv = rv_alloc(n);
2614 while((*t = *s++)) t++;
2620 /* freedtoa(s) must be used to free values s returned by dtoa
2621 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2622 * but for consistency with earlier versions of dtoa, it is optional
2623 * when MULTIPLE_THREADS is not defined.
2626 static void freedtoa (char *s);
2630 freedtoa(s) char *s;
2635 Bigint *b = (Bigint *)((int *)s - 1);
2636 b->maxwds = 1 << (b->k = *(int*)b);
2638 #ifndef MULTIPLE_THREADS
2639 if (s == dtoa_result)
2644 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2646 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2647 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2650 * 1. Rather than iterating, we use a simple numeric overestimate
2651 * to determine k = floor(log10(d)). We scale relevant
2652 * quantities using O(log2(k)) rather than O(k) multiplications.
2653 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2654 * try to generate digits strictly left to right. Instead, we
2655 * compute with fewer bits and propagate the carry if necessary
2656 * when rounding the final digit up. This is often faster.
2657 * 3. Under the assumption that input will be rounded nearest,
2658 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2659 * That is, we allow equality in stopping tests when the
2660 * round-nearest rule will give the same floating-point value
2661 * as would satisfaction of the stopping test with strict
2663 * 4. We remove common factors of powers of 2 from relevant
2665 * 5. When converting floating-point integers less than 1e16,
2666 * we use floating-point arithmetic rather than resorting
2667 * to multiple-precision integers.
2668 * 6. When asked to produce fewer than 15 digits, we first try
2669 * to get by with floating-point arithmetic; we resort to
2670 * multiple-precision integer arithmetic only if we cannot
2671 * guarantee that the floating-point calculation has given
2672 * the correctly rounded result. For k requested digits and
2673 * "uniformly" distributed input, the probability is
2674 * something like 10^(k-15) that we must resort to the Long
2681 (d, mode, ndigits, decpt, sign, rve)
2682 double d; int mode, ndigits, *decpt, *sign; char **rve;
2684 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
2687 /* Arguments ndigits, decpt, sign are similar to those
2688 of ecvt and fcvt; trailing zeros are suppressed from
2689 the returned string. If not null, *rve is set to point
2690 to the end of the return value. If d is +-Infinity or NaN,
2691 then *decpt is set to 9999.
2694 0 ==> shortest string that yields d when read in
2695 and rounded to nearest.
2696 1 ==> like 0, but with Steele & White stopping rule;
2697 e.g. with IEEE P754 arithmetic , mode 0 gives
2698 1e23 whereas mode 1 gives 9.999999999999999e22.
2699 2 ==> max(1,ndigits) significant digits. This gives a
2700 return value similar to that of ecvt, except
2701 that trailing zeros are suppressed.
2702 3 ==> through ndigits past the decimal point. This
2703 gives a return value similar to that from fcvt,
2704 except that trailing zeros are suppressed, and
2705 ndigits can be negative.
2706 4,5 ==> similar to 2 and 3, respectively, but (in
2707 round-nearest mode) with the tests of mode 0 to
2708 possibly return a shorter string that rounds to d.
2709 With IEEE arithmetic and compilation with
2710 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2711 as modes 2 and 3 when FLT_ROUNDS != 1.
2712 6-9 ==> Debugging modes similar to mode - 4: don't try
2713 fast floating-point estimate (if applicable).
2715 Values of mode other than 0-9 are treated as mode 0.
2717 Sufficient space is allocated to the return value
2718 to hold the suppressed trailing zeros.
2721 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
2722 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2723 spec_case, try_quick;
2725 #ifndef Sudden_Underflow
2729 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
2732 #ifdef Honor_FLT_ROUNDS
2736 int inexact, oldinexact;
2739 #ifndef MULTIPLE_THREADS
2741 freedtoa(dtoa_result);
2746 if (word0(d) & Sign_bit) {
2747 /* set sign for everything, including 0's and NaNs */
2749 word0(d) &= ~Sign_bit; /* clear sign bit */
2754 #if defined(IEEE_Arith) + defined(VAX)
2756 if ((word0(d) & Exp_mask) == Exp_mask)
2758 if (word0(d) == 0x8000)
2761 /* Infinity or NaN */
2764 if (!word1(d) && !(word0(d) & 0xfffff))
2765 return nrv_alloc("Infinity", rve, 8);
2767 return nrv_alloc("NaN", rve, 3);
2771 dval(d) += 0; /* normalize */
2775 return nrv_alloc("0", rve, 1);
2779 try_quick = oldinexact = get_inexact();
2782 #ifdef Honor_FLT_ROUNDS
2783 if ((rounding = Flt_Rounds) >= 2) {
2785 rounding = rounding == 2 ? 0 : 2;
2792 b = d2b(dval(d), &be, &bbits);
2793 #ifdef Sudden_Underflow
2794 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
2796 if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
2799 word0(d2) &= Frac_mask1;
2800 word0(d2) |= Exp_11;
2802 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
2806 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2807 * log10(x) = log(x) / log(10)
2808 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2809 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2811 * This suggests computing an approximation k to log10(d) by
2813 * k = (i - Bias)*0.301029995663981
2814 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2816 * We want k to be too large rather than too small.
2817 * The error in the first-order Taylor series approximation
2818 * is in our favor, so we just round up the constant enough
2819 * to compensate for any error in the multiplication of
2820 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2821 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2822 * adding 1e-13 to the constant term more than suffices.
2823 * Hence we adjust the constant term to 0.1760912590558.
2824 * (We could get a more accurate k by invoking log10,
2825 * but this is probably not worthwhile.)
2833 #ifndef Sudden_Underflow
2837 /* d is denormalized */
2839 i = bbits + be + (Bias + (P-1) - 1);
2840 x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32
2841 : word1(d) << 32 - i;
2843 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2844 i -= (Bias + (P-1) - 1) + 1;
2848 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2850 if (ds < 0. && ds != k)
2851 k--; /* want k = floor(ds) */
2853 if (k >= 0 && k <= Ten_pmax) {
2854 if (dval(d) < tens[k])
2877 if (mode < 0 || mode > 9)
2881 #ifdef Check_FLT_ROUNDS
2882 try_quick = Rounding == 1;
2886 #endif /*SET_INEXACT*/
2906 ilim = ilim1 = i = ndigits;
2912 i = ndigits + k + 1;
2918 s = s0 = rv_alloc(i);
2920 #ifdef Honor_FLT_ROUNDS
2921 if (mode > 1 && rounding != 1)
2925 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2927 /* Try to get by with floating-point arithmetic. */
2933 ieps = 2; /* conservative */
2938 /* prevent overflows */
2940 dval(d) /= bigtens[n_bigtens-1];
2943 for(; j; j >>= 1, i++)
2951 dval(d) *= tens[j1 & 0xf];
2952 for(j = j1 >> 4; j; j >>= 1, i++)
2955 dval(d) *= bigtens[i];
2958 if (k_check && dval(d) < 1. && ilim > 0) {
2966 dval(eps) = ieps*dval(d) + 7.;
2967 word0(eps) -= (P-1)*Exp_msk1;
2971 if (dval(d) > dval(eps))
2973 if (dval(d) < -dval(eps))
2977 #ifndef No_leftright
2979 /* Use Steele & White method of only
2980 * generating digits needed.
2982 dval(eps) = 0.5/tens[ilim-1] - dval(eps);
2986 *s++ = '0' + (int)L;
2987 if (dval(d) < dval(eps))
2989 if (1. - dval(d) < dval(eps))
2999 /* Generate ilim digits, then fix them up. */
3000 dval(eps) *= tens[ilim-1];
3001 for(i = 1;; i++, dval(d) *= 10.) {
3002 L = (Long)(dval(d));
3003 if (!(dval(d) -= L))
3005 *s++ = '0' + (int)L;
3007 if (dval(d) > 0.5 + dval(eps))
3009 else if (dval(d) < 0.5 - dval(eps)) {
3017 #ifndef No_leftright
3027 /* Do we have a "small" integer? */
3029 if (be >= 0 && k <= Int_max) {
3032 if (ndigits < 0 && ilim <= 0) {
3034 if (ilim < 0 || dval(d) <= 5*ds)
3038 for(i = 1;; i++, dval(d) *= 10.) {
3039 L = (Long)(dval(d) / ds);
3041 #ifdef Check_FLT_ROUNDS
3042 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
3048 *s++ = '0' + (int)L;
3056 #ifdef Honor_FLT_ROUNDS
3060 case 2: goto bump_up;
3064 if (dval(d) > ds || dval(d) == ds && L & 1) {
3085 #ifndef Sudden_Underflow
3086 denorm ? be + (Bias + (P-1) - 1 + 1) :
3089 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
3097 if (m2 > 0 && s2 > 0) {
3098 i = m2 < s2 ? m2 : s2;
3106 mhi = pow5mult(mhi, m5);
3115 b = pow5mult(b, b5);
3119 S = pow5mult(S, s5);
3121 /* Check for special case that d is a normalized power of 2. */
3124 if ((mode < 2 || leftright)
3125 #ifdef Honor_FLT_ROUNDS
3129 if (!word1(d) && !(word0(d) & Bndry_mask)
3130 #ifndef Sudden_Underflow
3131 && word0(d) & (Exp_mask & ~Exp_msk1)
3134 /* The special case */
3141 /* Arrange for convenient computation of quotients:
3142 * shift left if necessary so divisor has 4 leading 0 bits.
3144 * Perhaps we should just compute leading 28 bits of S once
3145 * and for all and pass them and a shift to quorem, so it
3146 * can do shifts and ors to compute the numerator for q.
3149 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
3152 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
3174 b = multadd(b, 10, 0); /* we botched the k estimate */
3176 mhi = multadd(mhi, 10, 0);
3180 if (ilim <= 0 && (mode == 3 || mode == 5)) {
3181 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
3182 /* no digits, fcvt style */
3194 mhi = lshift(mhi, m2);
3196 /* Compute mlo -- check for special case
3197 * that d is a normalized power of 2.
3202 mhi = Balloc(mhi->k);
3204 mhi = lshift(mhi, Log2P);
3208 dig = quorem(b,S) + '0';
3209 /* Do we yet have the shortest decimal string
3210 * that will round to d?
3213 delta = diff(S, mhi);
3214 j1 = delta->sign ? 1 : cmp(b, delta);
3216 #ifndef ROUND_BIASED
3217 if (j1 == 0 && mode != 1 && !(word1(d) & 1)
3218 #ifdef Honor_FLT_ROUNDS
3227 else if (!b->x[0] && b->wds <= 1)
3234 if (j < 0 || j == 0 && mode != 1
3235 #ifndef ROUND_BIASED
3239 if (!b->x[0] && b->wds <= 1) {
3245 #ifdef Honor_FLT_ROUNDS
3248 case 0: goto accept_dig;
3249 case 2: goto keep_dig;
3251 #endif /*Honor_FLT_ROUNDS*/
3255 if ((j1 > 0 || j1 == 0 && dig & 1)
3264 #ifdef Honor_FLT_ROUNDS
3268 if (dig == '9') { /* possible if i == 1 */
3276 #ifdef Honor_FLT_ROUNDS
3282 b = multadd(b, 10, 0);
3284 mlo = mhi = multadd(mhi, 10, 0);
3286 mlo = multadd(mlo, 10, 0);
3287 mhi = multadd(mhi, 10, 0);
3293 *s++ = dig = quorem(b,S) + '0';
3294 if (!b->x[0] && b->wds <= 1) {
3302 b = multadd(b, 10, 0);
3305 /* Round off last digit */
3307 #ifdef Honor_FLT_ROUNDS
3309 case 0: goto trimzeros;
3310 case 2: goto roundoff;
3315 if (j > 0 || j == 0 && dig & 1) {
3333 if (mlo && mlo != mhi)
3341 word0(d) = Exp_1 + (70 << Exp_shift);
3346 else if (!oldinexact)