1 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
2 /****************************************************************
4 * The author of this software is David M. Gay.
6 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
8 * Permission to use, copy, modify, and distribute this software for any
9 * purpose without fee is hereby granted, provided that this entire notice
10 * is included in all copies of any software which is or includes a copy
11 * or modification of this software and in all copies of the supporting
12 * documentation for such software.
14 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
15 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
16 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
17 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
19 ***************************************************************/
21 /* Please send bug reports to David M. Gay (dmg at acm dot org,
22 * with " at " changed at "@" and " dot " changed to "."). */
24 /* On a machine with IEEE extended-precision registers, it is
25 * necessary to specify double-precision (53-bit) rounding precision
26 * before invoking strtod or dtoa. If the machine uses (the equivalent
27 * of) Intel 80x87 arithmetic, the call
28 * _control87(PC_53, MCW_PC);
29 * does this with many compilers. Whether this or another call is
30 * appropriate depends on the compiler; for this to work, it may be
31 * necessary to #include "float.h" or another system-dependent header
35 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
37 * This strtod returns a nearest machine number to the input decimal
38 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
39 * broken by the IEEE round-even rule. Otherwise ties are broken by
40 * biased rounding (add half and chop).
42 * Inspired loosely by William D. Clinger's paper "How to Read Floating
43 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
47 * 1. We only require IEEE, IBM, or VAX double-precision
48 * arithmetic (not IEEE double-extended).
49 * 2. We get by with floating-point arithmetic in a case that
50 * Clinger missed -- when we're computing d * 10^n
51 * for a small integer d and the integer n is not too
52 * much larger than 22 (the maximum integer k for which
53 * we can represent 10^k exactly), we may be able to
54 * compute (d*10^k) * 10^(e-k) with just one roundoff.
55 * 3. Rather than a bit-at-a-time adjustment of the binary
56 * result in the hard case, we use floating-point
57 * arithmetic to determine the adjustment to within
58 * one bit; only in really hard cases do we need to
59 * compute a second residual.
60 * 4. Because of 3., we don't need a large table of powers of 10
61 * for ten-to-e (just some small tables, e.g. of 10^k
66 * #define IEEE_8087 for IEEE-arithmetic machines where the least
67 * significant byte has the lowest address.
68 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
69 * significant byte has the lowest address.
70 * #define Long int on machines with 32-bit ints and 64-bit longs.
71 * #define IBM for IBM mainframe-style floating-point arithmetic.
72 * #define VAX for VAX-style floating-point arithmetic (D_floating).
73 * #define No_leftright to omit left-right logic in fast floating-point
74 * computation of dtoa.
75 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
76 * and strtod and dtoa should round accordingly.
77 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
78 * and Honor_FLT_ROUNDS is not #defined.
79 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
80 * that use extended-precision instructions to compute rounded
81 * products and quotients) with IBM.
82 * #define ROUND_BIASED for IEEE-format with biased rounding.
83 * #define Inaccurate_Divide for IEEE-format with correctly rounded
84 * products but inaccurate quotients, e.g., for Intel i860.
85 * #define NO_LONG_LONG on machines that do not have a "long long"
86 * integer type (of >= 64 bits). On such machines, you can
87 * #define Just_16 to store 16 bits per 32-bit Long when doing
88 * high-precision integer arithmetic. Whether this speeds things
89 * up or slows things down depends on the machine and the number
90 * being converted. If long long is available and the name is
91 * something other than "long long", #define Llong to be the name,
92 * and if "unsigned Llong" does not work as an unsigned version of
93 * Llong, #define #ULLong to be the corresponding unsigned type.
94 * #define KR_headers for old-style C function headers.
95 * #define Bad_float_h if your system lacks a float.h or if it does not
96 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
97 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
98 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
99 * if memory is available and otherwise does something you deem
100 * appropriate. If MALLOC is undefined, malloc will be invoked
101 * directly -- and assumed always to succeed. Similarly, if you
102 * want something other than the system's free() to be called to
103 * recycle memory acquired from MALLOC, #define FREE to be the
104 * name of the alternate routine. (Unless you #define
105 * NO_GLOBAL_STATE and call destroydtoa, FREE or free is only
106 * called in pathological cases, e.g., in a dtoa call after a dtoa
107 * return in mode 3 with thousands of digits requested.)
108 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
109 * memory allocations from a private pool of memory when possible.
110 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
111 * unless #defined to be a different length. This default length
112 * suffices to get rid of MALLOC calls except for unusual cases,
113 * such as decimal-to-binary conversion of a very long string of
114 * digits. The longest string dtoa can return is about 751 bytes
115 * long. For conversions by strtod of strings of 800 digits and
116 * all dtoa conversions in single-threaded executions with 8-byte
117 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
118 * pointers, PRIVATE_MEM >= 7112 appears adequate.
119 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
120 * multiple threads. In this case, you must provide (or suitably
121 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
122 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
123 * in pow5mult, ensures lazy evaluation of only one copy of high
124 * powers of 5; omitting this lock would introduce a small
125 * probability of wasting memory, but would otherwise be harmless.)
126 * You must also invoke freedtoa(s) to free the value s returned by
127 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
128 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
129 * avoids underflows on inputs whose result does not underflow.
130 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
131 * floating-point numbers and flushes underflows to zero rather
132 * than implementing gradual underflow, then you must also #define
134 * #define USE_LOCALE to use the current locale's decimal_point value.
135 * #define SET_INEXACT if IEEE arithmetic is being used and extra
136 * computation should be done to set the inexact flag when the
137 * result is inexact and avoid setting inexact when the result
138 * is exact. In this case, dtoa.c must be compiled in
139 * an environment, perhaps provided by #include "dtoa.c" in a
140 * suitable wrapper, that defines two functions,
141 * int get_inexact(void);
142 * void clear_inexact(void);
143 * such that get_inexact() returns a nonzero value if the
144 * inexact bit is already set, and clear_inexact() sets the
145 * inexact bit to 0. When SET_INEXACT is #defined, strtod
146 * also does extra computations to set the underflow and overflow
147 * flags when appropriate (i.e., when the result is tiny and
148 * inexact or when it is a numeric value rounded to +-infinity).
149 * #define NO_ERRNO if strtod should not assign errno = ERANGE when
150 * the result overflows to +-Infinity or underflows to 0.
151 * #define NO_GLOBAL_STATE to avoid defining any non-const global or
152 * static variables. Instead the necessary state is stored in an
153 * opaque struct, DtoaState, a pointer to which must be passed to
154 * every entry point. Two new functions are added to the API:
155 * DtoaState *newdtoa(void);
156 * void destroydtoa(DtoaState *);
163 typedef unsigned Long ULong;
168 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
180 extern char *MALLOC();
182 extern void *MALLOC(size_t);
185 #define MALLOC malloc
192 #ifndef Omit_Private_Memory
194 #define PRIVATE_MEM 2304
196 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
200 #undef Avoid_Underflow
214 #define DBL_MAX_10_EXP 308
215 #define DBL_MAX_EXP 1024
217 #endif /*IEEE_Arith*/
221 #define DBL_MAX_10_EXP 75
222 #define DBL_MAX_EXP 63
224 #define DBL_MAX 7.2370055773322621e+75
229 #define DBL_MAX_10_EXP 38
230 #define DBL_MAX_EXP 127
232 #define DBL_MAX 1.7014118346046923e+38
236 #define LONG_MAX 2147483647
239 #else /* ifndef Bad_float_h */
241 #endif /* Bad_float_h */
253 #define CONST /* blank */
259 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
260 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
263 typedef union { double d; ULong L[2]; } U;
265 #define dval(x) ((x).d)
267 #define word0(x) ((x).L[1])
268 #define word1(x) ((x).L[0])
270 #define word0(x) ((x).L[0])
271 #define word1(x) ((x).L[1])
274 /* The following definition of Storeinc is appropriate for MIPS processors.
275 * An alternative that might be better on some machines is
276 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
278 #if defined(IEEE_8087) + defined(VAX)
279 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
280 ((unsigned short *)a)[0] = (unsigned short)c, a++)
282 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
283 ((unsigned short *)a)[1] = (unsigned short)c, a++)
286 /* #define P DBL_MANT_DIG */
287 /* Ten_pmax = floor(P*log(2)/log(5)) */
288 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
289 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
290 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
294 #define Exp_shift1 20
295 #define Exp_msk1 0x100000
296 #define Exp_msk11 0x100000
297 #define Exp_mask 0x7ff00000
301 #define Exp_1 0x3ff00000
302 #define Exp_11 0x3ff00000
304 #define Frac_mask 0xfffff
305 #define Frac_mask1 0xfffff
308 #define Bndry_mask 0xfffff
309 #define Bndry_mask1 0xfffff
311 #define Sign_bit 0x80000000
317 #ifndef NO_IEEE_Scale
318 #define Avoid_Underflow
319 #ifdef Flush_Denorm /* debugging option */
320 #undef Sudden_Underflow
326 #define Flt_Rounds FLT_ROUNDS
330 #endif /*Flt_Rounds*/
332 #ifdef Honor_FLT_ROUNDS
333 #define Rounding rounding
334 #undef Check_FLT_ROUNDS
335 #define Check_FLT_ROUNDS
337 #define Rounding Flt_Rounds
340 #else /* ifndef IEEE_Arith */
341 #undef Check_FLT_ROUNDS
342 #undef Honor_FLT_ROUNDS
344 #undef Sudden_Underflow
345 #define Sudden_Underflow
350 #define Exp_shift1 24
351 #define Exp_msk1 0x1000000
352 #define Exp_msk11 0x1000000
353 #define Exp_mask 0x7f000000
356 #define Exp_1 0x41000000
357 #define Exp_11 0x41000000
358 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
359 #define Frac_mask 0xffffff
360 #define Frac_mask1 0xffffff
363 #define Bndry_mask 0xefffff
364 #define Bndry_mask1 0xffffff
366 #define Sign_bit 0x80000000
368 #define Tiny0 0x100000
377 #define Exp_msk1 0x80
378 #define Exp_msk11 0x800000
379 #define Exp_mask 0x7f80
382 #define Exp_1 0x40800000
383 #define Exp_11 0x4080
385 #define Frac_mask 0x7fffff
386 #define Frac_mask1 0xffff007f
389 #define Bndry_mask 0xffff007f
390 #define Bndry_mask1 0xffff007f
392 #define Sign_bit 0x8000
398 #endif /* IBM, VAX */
399 #endif /* IEEE_Arith */
406 #define rounded_product(a,b) a = rnd_prod(a, b)
407 #define rounded_quotient(a,b) a = rnd_quot(a, b)
409 extern double rnd_prod(), rnd_quot();
411 extern double rnd_prod(double, double), rnd_quot(double, double);
414 #define rounded_product(a,b) a *= b
415 #define rounded_quotient(a,b) a /= b
418 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
419 #define Big1 0xffffffff
426 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
428 #define FFFFFFFF 0xffffffffUL
435 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
436 * This makes some inner loops simpler and sometimes saves work
437 * during multiplications, but it often seems to make things slightly
438 * slower. Hence the default is now to store 32 bits per Long.
441 #else /* long long available */
443 #define Llong long long
446 #define ULLong unsigned Llong
448 #endif /* NO_LONG_LONG */
450 #ifndef MULTIPLE_THREADS
451 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
452 #define FREE_DTOA_LOCK(n) /*nothing*/
460 int k, maxwds, sign, wds;
464 typedef struct Bigint Bigint;
466 #ifdef NO_GLOBAL_STATE
467 #ifdef MULTIPLE_THREADS
468 #error "cannot have both NO_GLOBAL_STATE and MULTIPLE_THREADS"
472 #define DECLARE_GLOBAL_STATE /* nothing */
474 #define DECLARE_GLOBAL_STATE static
477 DECLARE_GLOBAL_STATE Bigint *freelist[Kmax+1];
478 DECLARE_GLOBAL_STATE Bigint *p5s;
479 #ifndef Omit_Private_Memory
480 DECLARE_GLOBAL_STATE double private_mem[PRIVATE_mem];
481 DECLARE_GLOBAL_STATE double *pmem_next
482 #ifndef NO_GLOBAL_STATE
487 #ifdef NO_GLOBAL_STATE
489 typedef struct DtoaState DtoaState;
491 #define STATE_PARAM state,
492 #define STATE_PARAM_DECL DtoaState *state;
494 #define STATE_PARAM DtoaState *state,
496 #define PASS_STATE state,
497 #define GET_STATE(field) (state->field)
502 DtoaState *state = (DtoaState *) MALLOC(sizeof(DtoaState));
504 memset(state, 0, sizeof(DtoaState));
505 #ifndef Omit_Private_Memory
506 state->pmem_next = state->private_mem;
515 (state) STATE_PARAM_DECL
523 for (i = 0; i <= Kmax; i++) {
524 for (v = GET_STATE(freelist)[i]; v; v = next) {
526 #ifndef Omit_Private_Memory
527 if ((double*)v < GET_STATE(private_mem) ||
528 (double*)v >= GET_STATE(private_mem) + PRIVATE_mem)
537 #define STATE_PARAM /* nothing */
538 #define STATE_PARAM_DECL /* nothing */
539 #define PASS_STATE /* nothing */
540 #define GET_STATE(name) name
546 (STATE_PARAM k) STATE_PARAM_DECL int k;
553 #ifndef Omit_Private_Memory
557 ACQUIRE_DTOA_LOCK(0);
558 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
559 /* but this case seems very unlikely. */
560 if (k <= Kmax && (rv = GET_STATE(freelist)[k]))
561 GET_STATE(freelist)[k] = rv->next;
564 #ifdef Omit_Private_Memory
565 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
567 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
569 if (k <= Kmax && GET_STATE(pmem_next) - GET_STATE(private_mem) + len <= PRIVATE_mem) {
570 rv = (Bigint*)GET_STATE(pmem_next);
571 GET_STATE(pmem_next) += len;
574 rv = (Bigint*)MALLOC(len*sizeof(double));
580 rv->sign = rv->wds = 0;
587 (STATE_PARAM v) STATE_PARAM_DECL Bigint *v;
589 (STATE_PARAM Bigint *v)
596 ACQUIRE_DTOA_LOCK(0);
597 v->next = GET_STATE(freelist)[v->k];
598 GET_STATE(freelist)[v->k] = v;
604 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
605 y->wds*sizeof(Long) + 2*sizeof(int))
610 (STATE_PARAM b, m, a) STATE_PARAM_DECL Bigint *b; int m, a;
612 (STATE_PARAM Bigint *b, int m, int a) /* multiply by m and add a */
633 y = *x * (ULLong)m + carry;
635 *x++ = (ULong) y & FFFFFFFF;
639 y = (xi & 0xffff) * m + carry;
640 z = (xi >> 16) * m + (y >> 16);
642 *x++ = (z << 16) + (y & 0xffff);
652 if (wds >= b->maxwds) {
653 b1 = Balloc(PASS_STATE b->k+1);
658 b->x[wds++] = (ULong) carry;
667 (STATE_PARAM s, nd0, nd, y9) STATE_PARAM_DECL CONST char *s; int nd0, nd; ULong y9;
669 (STATE_PARAM CONST char *s, int nd0, int nd, ULong y9)
677 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
679 b = Balloc(PASS_STATE k);
683 b = Balloc(PASS_STATE k+1);
684 b->x[0] = y9 & 0xffff;
685 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
691 do b = multadd(PASS_STATE b, 10, *s++ - '0');
698 b = multadd(PASS_STATE b, 10, *s++ - '0');
705 (x) register ULong x;
712 if (!(x & 0xffff0000)) {
716 if (!(x & 0xff000000)) {
720 if (!(x & 0xf0000000)) {
724 if (!(x & 0xc0000000)) {
728 if (!(x & 0x80000000)) {
730 if (!(x & 0x40000000))
745 register ULong x = *y;
787 (STATE_PARAM i) STATE_PARAM_DECL int i;
794 b = Balloc(PASS_STATE 1);
803 (STATE_PARAM a, b) STATE_PARAM_DECL Bigint *a, *b;
805 (STATE_PARAM Bigint *a, Bigint *b)
810 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
821 if (a->wds < b->wds) {
832 c = Balloc(PASS_STATE k);
833 for(x = c->x, xa = x + wc; x < xa; x++)
841 for(; xb < xbe; xc0++) {
847 z = *x++ * (ULLong)y + *xc + carry;
849 *xc++ = (ULong) z & FFFFFFFF;
857 for(; xb < xbe; xb++, xc0++) {
858 if (y = *xb & 0xffff) {
863 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
865 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
878 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
881 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
889 for(; xb < xbe; xc0++) {
895 z = *x++ * y + *xc + carry;
905 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
913 (STATE_PARAM b, k) STATE_PARAM_DECL Bigint *b; int k;
915 (STATE_PARAM Bigint *b, int k)
918 Bigint *b1, *p5, *p51;
920 static CONST int p05[3] = { 5, 25, 125 };
923 b = multadd(PASS_STATE b, p05[i-1], 0);
927 if (!(p5 = GET_STATE(p5s))) {
929 #ifdef MULTIPLE_THREADS
930 ACQUIRE_DTOA_LOCK(1);
937 p5 = GET_STATE(p5s) = i2b(PASS_STATE 625);
943 b1 = mult(PASS_STATE b, p5);
949 if (!(p51 = p5->next)) {
950 #ifdef MULTIPLE_THREADS
951 ACQUIRE_DTOA_LOCK(1);
952 if (!(p51 = p5->next)) {
953 p51 = p5->next = mult(p5,p5);
958 p51 = p5->next = mult(PASS_STATE p5,p5);
970 (STATE_PARAM b, k) STATE_PARAM_DECL Bigint *b; int k;
972 (STATE_PARAM Bigint *b, int k)
977 ULong *x, *x1, *xe, z;
986 for(i = b->maxwds; n1 > i; i <<= 1)
988 b1 = Balloc(PASS_STATE k1);
990 for(i = 0; i < n; i++)
1011 *x1++ = *x << k & 0xffff | z;
1023 Bfree(PASS_STATE b);
1030 (a, b) Bigint *a, *b;
1032 (Bigint *a, Bigint *b)
1035 ULong *xa, *xa0, *xb, *xb0;
1041 if (i > 1 && !a->x[i-1])
1042 Bug("cmp called with a->x[a->wds-1] == 0");
1043 if (j > 1 && !b->x[j-1])
1044 Bug("cmp called with b->x[b->wds-1] == 0");
1054 return *xa < *xb ? -1 : 1;
1064 (STATE_PARAM a, b) STATE_PARAM_DECL Bigint *a, *b;
1066 (STATE_PARAM Bigint *a, Bigint *b)
1071 ULong *xa, *xae, *xb, *xbe, *xc;
1083 c = Balloc(PASS_STATE 0);
1096 c = Balloc(PASS_STATE a->k);
1108 y = (ULLong)*xa++ - *xb++ - borrow;
1109 borrow = y >> 32 & (ULong)1;
1110 *xc++ = (ULong) y & FFFFFFFF;
1115 borrow = y >> 32 & (ULong)1;
1116 *xc++ = (ULong) y & FFFFFFFF;
1121 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1122 borrow = (y & 0x10000) >> 16;
1123 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1124 borrow = (z & 0x10000) >> 16;
1129 y = (*xa & 0xffff) - borrow;
1130 borrow = (y & 0x10000) >> 16;
1131 z = (*xa++ >> 16) - borrow;
1132 borrow = (z & 0x10000) >> 16;
1137 y = *xa++ - *xb++ - borrow;
1138 borrow = (y & 0x10000) >> 16;
1144 borrow = (y & 0x10000) >> 16;
1166 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1167 #ifndef Avoid_Underflow
1168 #ifndef Sudden_Underflow
1177 #ifndef Avoid_Underflow
1178 #ifndef Sudden_Underflow
1181 L = -L >> Exp_shift;
1182 if (L < Exp_shift) {
1183 word0(a) = 0x80000 >> L;
1189 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
1200 (a, e) Bigint *a; int *e;
1205 ULong *xa, *xa0, w, y, z;
1219 if (!y) Bug("zero y in b2d");
1225 d0 = Exp_1 | y >> (Ebits - k);
1226 w = xa > xa0 ? *--xa : 0;
1227 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1230 z = xa > xa0 ? *--xa : 0;
1232 d0 = Exp_1 | y << k | z >> (32 - k);
1233 y = xa > xa0 ? *--xa : 0;
1234 d1 = z << k | y >> (32 - k);
1241 if (k < Ebits + 16) {
1242 z = xa > xa0 ? *--xa : 0;
1243 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1244 w = xa > xa0 ? *--xa : 0;
1245 y = xa > xa0 ? *--xa : 0;
1246 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1249 z = xa > xa0 ? *--xa : 0;
1250 w = xa > xa0 ? *--xa : 0;
1252 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1253 y = xa > xa0 ? *--xa : 0;
1254 d1 = w << k + 16 | y << k;
1258 word0(d) = d0 >> 16 | d0 << 16;
1259 word1(d) = d1 >> 16 | d1 << 16;
1270 (STATE_PARAM d, e, bits) STATE_PARAM_DECL U d; int *e, *bits;
1272 (STATE_PARAM U d, int *e, int *bits)
1278 #ifndef Sudden_Underflow
1283 d0 = word0(d) >> 16 | word0(d) << 16;
1284 d1 = word1(d) >> 16 | word1(d) << 16;
1291 b = Balloc(PASS_STATE 1);
1293 b = Balloc(PASS_STATE 2);
1298 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1299 #ifdef Sudden_Underflow
1300 de = (int)(d0 >> Exp_shift);
1305 if ((de = (int)(d0 >> Exp_shift)))
1310 if ((k = lo0bits(&y))) {
1311 x[0] = y | z << (32 - k);
1316 #ifndef Sudden_Underflow
1319 b->wds = (x[1] = z) ? 2 : 1;
1324 #ifndef Sudden_Underflow
1332 if (k = lo0bits(&y))
1334 x[0] = y | z << 32 - k & 0xffff;
1335 x[1] = z >> k - 16 & 0xffff;
1341 x[1] = y >> 16 | z << 16 - k & 0xffff;
1342 x[2] = z >> k & 0xffff;
1357 Bug("Zero passed to d2b");
1375 #ifndef Sudden_Underflow
1379 *e = (de - Bias - (P-1) << 2) + k;
1380 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1382 *e = de - Bias - (P-1) + k;
1385 #ifndef Sudden_Underflow
1388 *e = de - Bias - (P-1) + 1 + k;
1390 *bits = 32*i - hi0bits(x[i-1]);
1392 *bits = (i+2)*16 - hi0bits(x[i]);
1404 (a, b) Bigint *a, *b;
1406 (Bigint *a, Bigint *b)
1412 dval(da) = b2d(a, &ka);
1413 dval(db) = b2d(b, &kb);
1415 k = ka - kb + 32*(a->wds - b->wds);
1417 k = ka - kb + 16*(a->wds - b->wds);
1421 word0(da) += (k >> 2)*Exp_msk1;
1427 word0(db) += (k >> 2)*Exp_msk1;
1433 word0(da) += k*Exp_msk1;
1436 word0(db) += k*Exp_msk1;
1439 return dval(da) / dval(db);
1444 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1445 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1454 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1455 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1456 #ifdef Avoid_Underflow
1457 9007199254740992.*9007199254740992.e-256
1458 /* = 2^106 * 1e-53 */
1463 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1464 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1465 #define Scale_Bit 0x10
1469 bigtens[] = { 1e16, 1e32, 1e64 };
1470 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1473 bigtens[] = { 1e16, 1e32 };
1474 static CONST double tinytens[] = { 1e-16, 1e-32 };
1482 (STATE_PARAM s00, se) STATE_PARAM_DECL CONST char *s00; char **se;
1484 (STATE_PARAM CONST char *s00, char **se)
1487 #ifdef Avoid_Underflow
1490 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1491 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1492 CONST char *s, *s0, *s1;
1497 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1499 int inexact, oldinexact;
1501 #ifdef Honor_FLT_ROUNDS
1509 delta = bb = bd = bs = 0;
1512 sign = nz0 = nz = 0;
1514 for(s = s00;;s++) switch(*s) {
1537 while(*++s == '0') ;
1543 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1550 s1 = localeconv()->decimal_point;
1571 for(; c == '0'; c = *++s)
1573 if (c > '0' && c <= '9') {
1581 for(; c >= '0' && c <= '9'; c = *++s) {
1586 for(i = 1; i < nz; i++)
1589 else if (nd <= DBL_DIG + 1)
1593 else if (nd <= DBL_DIG + 1)
1601 if (c == 'e' || c == 'E') {
1602 if (!nd && !nz && !nz0) {
1613 if (c >= '0' && c <= '9') {
1616 if (c > '0' && c <= '9') {
1619 while((c = *++s) >= '0' && c <= '9')
1621 if (s - s1 > 8 || L > 19999)
1622 /* Avoid confusion from exponents
1623 * so large that e might overflow.
1625 e = 19999; /* safe for 16 bit ints */
1647 /* Now we have nd0 digits, starting at s0, followed by a
1648 * decimal point, followed by nd-nd0 digits. The number we're
1649 * after is the integer represented by those digits times
1654 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1659 oldinexact = get_inexact();
1661 dval(rv) = tens[k - 9] * dval(rv) + z;
1665 #ifndef RND_PRODQUOT
1666 #ifndef Honor_FLT_ROUNDS
1674 if (e <= Ten_pmax) {
1676 goto vax_ovfl_check;
1678 #ifdef Honor_FLT_ROUNDS
1679 /* round correctly FLT_ROUNDS = 2 or 3 */
1685 /* rv = */ rounded_product(dval(rv), tens[e]);
1690 if (e <= Ten_pmax + i) {
1691 /* A fancier test would sometimes let us do
1692 * this for larger i values.
1694 #ifdef Honor_FLT_ROUNDS
1695 /* round correctly FLT_ROUNDS = 2 or 3 */
1702 dval(rv) *= tens[i];
1704 /* VAX exponent range is so narrow we must
1705 * worry about overflow here...
1708 word0(rv) -= P*Exp_msk1;
1709 /* rv = */ rounded_product(dval(rv), tens[e]);
1710 if ((word0(rv) & Exp_mask)
1711 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1713 word0(rv) += P*Exp_msk1;
1715 /* rv = */ rounded_product(dval(rv), tens[e]);
1720 #ifndef Inaccurate_Divide
1721 else if (e >= -Ten_pmax) {
1722 #ifdef Honor_FLT_ROUNDS
1723 /* round correctly FLT_ROUNDS = 2 or 3 */
1729 /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1740 oldinexact = get_inexact();
1742 #ifdef Avoid_Underflow
1745 #ifdef Honor_FLT_ROUNDS
1746 if ((rounding = Flt_Rounds) >= 2) {
1748 rounding = rounding == 2 ? 0 : 2;
1754 #endif /*IEEE_Arith*/
1756 /* Get starting approximation = rv * 10**e1 */
1760 dval(rv) *= tens[i];
1762 if (e1 > DBL_MAX_10_EXP) {
1767 /* Can't trust HUGE_VAL */
1769 #ifdef Honor_FLT_ROUNDS
1771 case 0: /* toward 0 */
1772 case 3: /* toward -infinity */
1777 word0(rv) = Exp_mask;
1780 #else /*Honor_FLT_ROUNDS*/
1781 word0(rv) = Exp_mask;
1783 #endif /*Honor_FLT_ROUNDS*/
1785 /* set overflow bit */
1787 dval(rv0) *= dval(rv0);
1789 #else /*IEEE_Arith*/
1792 #endif /*IEEE_Arith*/
1798 for(j = 0; e1 > 1; j++, e1 >>= 1)
1800 dval(rv) *= bigtens[j];
1801 /* The last multiplication could overflow. */
1802 word0(rv) -= P*Exp_msk1;
1803 dval(rv) *= bigtens[j];
1804 if ((z = word0(rv) & Exp_mask)
1805 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1807 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1808 /* set to largest number */
1809 /* (Can't trust DBL_MAX) */
1814 word0(rv) += P*Exp_msk1;
1820 dval(rv) /= tens[i];
1822 if (e1 >= 1 << n_bigtens)
1824 #ifdef Avoid_Underflow
1827 for(j = 0; e1 > 0; j++, e1 >>= 1)
1829 dval(rv) *= tinytens[j];
1830 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
1831 >> Exp_shift)) > 0) {
1832 /* scaled rv is denormal; zap j low bits */
1836 word0(rv) = (P+2)*Exp_msk1;
1838 word0(rv) &= 0xffffffff << (j-32);
1841 word1(rv) &= 0xffffffff << j;
1844 for(j = 0; e1 > 1; j++, e1 >>= 1)
1846 dval(rv) *= tinytens[j];
1847 /* The last multiplication could underflow. */
1848 dval(rv0) = dval(rv);
1849 dval(rv) *= tinytens[j];
1851 dval(rv) = 2.*dval(rv0);
1852 dval(rv) *= tinytens[j];
1864 #ifndef Avoid_Underflow
1867 /* The refinement below will clean
1868 * this approximation up.
1875 /* Now the hard part -- adjusting rv to the correct value.*/
1877 /* Put digits into bd: true value = bd * 10^e */
1879 bd0 = s2b(PASS_STATE s0, nd0, nd, y);
1882 bd = Balloc(PASS_STATE bd0->k);
1884 bb = d2b(PASS_STATE rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1885 bs = i2b(PASS_STATE 1);
1900 #ifdef Honor_FLT_ROUNDS
1904 #ifdef Avoid_Underflow
1906 i = j + bbbits - 1; /* logb(rv) */
1907 if (i < Emin) /* denormal */
1911 #else /*Avoid_Underflow*/
1912 #ifdef Sudden_Underflow
1914 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1918 #else /*Sudden_Underflow*/
1920 i = j + bbbits - 1; /* logb(rv) */
1921 if (i < Emin) /* denormal */
1925 #endif /*Sudden_Underflow*/
1926 #endif /*Avoid_Underflow*/
1929 #ifdef Avoid_Underflow
1932 i = bb2 < bd2 ? bb2 : bd2;
1941 bs = pow5mult(PASS_STATE bs, bb5);
1942 bb1 = mult(PASS_STATE bs, bb);
1943 Bfree(PASS_STATE bb);
1947 bb = lshift(PASS_STATE bb, bb2);
1949 bd = pow5mult(PASS_STATE bd, bd5);
1951 bd = lshift(PASS_STATE bd, bd2);
1953 bs = lshift(PASS_STATE bs, bs2);
1954 delta = diff(PASS_STATE bb, bd);
1955 dsign = delta->sign;
1958 #ifdef Honor_FLT_ROUNDS
1959 if (rounding != 1) {
1961 /* Error is less than an ulp */
1962 if (!delta->x[0] && delta->wds <= 1) {
1978 && !(word0(rv) & Frac_mask)) {
1979 y = word0(rv) & Exp_mask;
1980 #ifdef Avoid_Underflow
1981 if (!scale || y > 2*P*Exp_msk1)
1986 delta = lshift(PASS_STATE delta,Log2P);
1987 if (cmp(delta, bs) <= 0)
1992 #ifdef Avoid_Underflow
1993 if (scale && (y = word0(rv) & Exp_mask)
1995 word0(adj) += (2*P+1)*Exp_msk1 - y;
1997 #ifdef Sudden_Underflow
1998 if ((word0(rv) & Exp_mask) <=
2000 word0(rv) += P*Exp_msk1;
2001 dval(rv) += adj*ulp(rv);
2002 word0(rv) -= P*Exp_msk1;
2005 #endif /*Sudden_Underflow*/
2006 #endif /*Avoid_Underflow*/
2007 dval(rv) += adj*ulp(rv);
2011 adj = ratio(delta, bs);
2014 if (adj <= 0x7ffffffe) {
2015 /* adj = rounding ? ceil(adj) : floor(adj); */
2018 if (!((rounding>>1) ^ dsign))
2023 #ifdef Avoid_Underflow
2024 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2025 word0(adj) += (2*P+1)*Exp_msk1 - y;
2027 #ifdef Sudden_Underflow
2028 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2029 word0(rv) += P*Exp_msk1;
2035 word0(rv) -= P*Exp_msk1;
2038 #endif /*Sudden_Underflow*/
2039 #endif /*Avoid_Underflow*/
2047 #endif /*Honor_FLT_ROUNDS*/
2050 /* Error is less than half an ulp -- check for
2051 * special case of mantissa a power of two.
2053 if (dsign || word1(rv) || word0(rv) & Bndry_mask
2055 #ifdef Avoid_Underflow
2056 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2058 || (word0(rv) & Exp_mask) <= Exp_msk1
2063 if (!delta->x[0] && delta->wds <= 1)
2068 if (!delta->x[0] && delta->wds <= 1) {
2075 delta = lshift(PASS_STATE delta,Log2P);
2076 if (cmp(delta, bs) > 0)
2081 /* exactly half-way between */
2083 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2085 #ifdef Avoid_Underflow
2086 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2087 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2090 /*boundary case -- increment exponent*/
2091 word0(rv) = (word0(rv) & Exp_mask)
2098 #ifdef Avoid_Underflow
2104 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2106 /* boundary case -- decrement exponent */
2107 #ifdef Sudden_Underflow /*{{*/
2108 L = word0(rv) & Exp_mask;
2112 #ifdef Avoid_Underflow
2113 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
2116 #endif /*Avoid_Underflow*/
2120 #else /*Sudden_Underflow}{*/
2121 #ifdef Avoid_Underflow
2123 L = word0(rv) & Exp_mask;
2124 if (L <= (2*P+1)*Exp_msk1) {
2125 if (L > (P+2)*Exp_msk1)
2126 /* round even ==> */
2129 /* rv = smallest denormal */
2133 #endif /*Avoid_Underflow*/
2134 L = (word0(rv) & Exp_mask) - Exp_msk1;
2135 #endif /*Sudden_Underflow}}*/
2136 word0(rv) = L | Bndry_mask1;
2137 word1(rv) = 0xffffffff;
2144 #ifndef ROUND_BIASED
2145 if (!(word1(rv) & LSB))
2149 dval(rv) += ulp(rv);
2150 #ifndef ROUND_BIASED
2152 dval(rv) -= ulp(rv);
2153 #ifndef Sudden_Underflow
2158 #ifdef Avoid_Underflow
2164 if ((aadj = ratio(delta, bs)) <= 2.) {
2166 aadj = dval(aadj1) = 1.;
2167 else if (word1(rv) || word0(rv) & Bndry_mask) {
2168 #ifndef Sudden_Underflow
2169 if (word1(rv) == Tiny1 && !word0(rv))
2176 /* special case -- power of FLT_RADIX to be */
2177 /* rounded down... */
2179 if (aadj < 2./FLT_RADIX)
2180 aadj = 1./FLT_RADIX;
2183 dval(aadj1) = -aadj;
2188 dval(aadj1) = dsign ? aadj : -aadj;
2189 #ifdef Check_FLT_ROUNDS
2191 case 2: /* towards +infinity */
2194 case 0: /* towards 0 */
2195 case 3: /* towards -infinity */
2199 if (Flt_Rounds == 0)
2201 #endif /*Check_FLT_ROUNDS*/
2203 y = word0(rv) & Exp_mask;
2205 /* Check for overflow */
2207 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2208 dval(rv0) = dval(rv);
2209 word0(rv) -= P*Exp_msk1;
2210 adj = dval(aadj1) * ulp(rv);
2212 if ((word0(rv) & Exp_mask) >=
2213 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2214 if (word0(rv0) == Big0 && word1(rv0) == Big1)
2221 word0(rv) += P*Exp_msk1;
2224 #ifdef Avoid_Underflow
2225 if (scale && y <= 2*P*Exp_msk1) {
2226 if (aadj <= 0x7fffffff) {
2227 if ((z = (ULong) aadj) <= 0)
2230 dval(aadj1) = dsign ? aadj : -aadj;
2232 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
2234 adj = dval(aadj1) * ulp(rv);
2237 #ifdef Sudden_Underflow
2238 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2239 dval(rv0) = dval(rv);
2240 word0(rv) += P*Exp_msk1;
2241 adj = dval(aadj1) * ulp(rv);
2244 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
2246 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2249 if (word0(rv0) == Tiny0
2250 && word1(rv0) == Tiny1)
2257 word0(rv) -= P*Exp_msk1;
2260 adj = dval(aadj1) * ulp(rv);
2263 #else /*Sudden_Underflow*/
2264 /* Compute adj so that the IEEE rounding rules will
2265 * correctly round rv + adj in some half-way cases.
2266 * If rv * ulp(rv) is denormalized (i.e.,
2267 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2268 * trouble from bits lost to denormalization;
2269 * example: 1.2e-307 .
2271 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
2272 dval(aadj1) = (double)(int)(aadj + 0.5);
2274 dval(aadj1) = -dval(aadj1);
2276 adj = dval(aadj1) * ulp(rv);
2278 #endif /*Sudden_Underflow*/
2279 #endif /*Avoid_Underflow*/
2281 z = word0(rv) & Exp_mask;
2283 #ifdef Avoid_Underflow
2287 /* Can we stop now? */
2290 /* The tolerances below are conservative. */
2291 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2292 if (aadj < .4999999 || aadj > .5000001)
2295 else if (aadj < .4999999/FLT_RADIX)
2300 Bfree(PASS_STATE bb);
2301 Bfree(PASS_STATE bd);
2302 Bfree(PASS_STATE bs);
2303 Bfree(PASS_STATE delta);
2308 word0(rv0) = Exp_1 + (70 << Exp_shift);
2313 else if (!oldinexact)
2316 #ifdef Avoid_Underflow
2318 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
2320 dval(rv) *= dval(rv0);
2322 /* try to avoid the bug of testing an 8087 register value */
2323 if (word0(rv) == 0 && word1(rv) == 0)
2327 #endif /* Avoid_Underflow */
2329 if (inexact && !(word0(rv) & Exp_mask)) {
2330 /* set underflow bit */
2332 dval(rv0) *= dval(rv0);
2336 Bfree(PASS_STATE bb);
2337 Bfree(PASS_STATE bd);
2338 Bfree(PASS_STATE bs);
2339 Bfree(PASS_STATE bd0);
2340 Bfree(PASS_STATE delta);
2344 return sign ? -dval(rv) : dval(rv);
2350 (b, S) Bigint *b, *S;
2352 (Bigint *b, Bigint *S)
2356 ULong *bx, *bxe, q, *sx, *sxe;
2358 ULLong borrow, carry, y, ys;
2360 ULong borrow, carry, y, ys;
2368 /*debug*/ if (b->wds > n)
2369 /*debug*/ Bug("oversize b in quorem");
2377 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
2379 /*debug*/ if (q > 9)
2380 /*debug*/ Bug("oversized quotient in quorem");
2387 ys = *sx++ * (ULLong)q + carry;
2389 y = *bx - (ys & FFFFFFFF) - borrow;
2390 borrow = y >> 32 & (ULong)1;
2391 *bx++ = (ULong) y & FFFFFFFF;
2395 ys = (si & 0xffff) * q + carry;
2396 zs = (si >> 16) * q + (ys >> 16);
2398 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2399 borrow = (y & 0x10000) >> 16;
2400 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2401 borrow = (z & 0x10000) >> 16;
2404 ys = *sx++ * q + carry;
2406 y = *bx - (ys & 0xffff) - borrow;
2407 borrow = (y & 0x10000) >> 16;
2415 while(--bxe > bx && !*bxe)
2420 if (cmp(b, S) >= 0) {
2430 y = *bx - (ys & FFFFFFFF) - borrow;
2431 borrow = y >> 32 & (ULong)1;
2432 *bx++ = (ULong) y & FFFFFFFF;
2436 ys = (si & 0xffff) + carry;
2437 zs = (si >> 16) + (ys >> 16);
2439 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2440 borrow = (y & 0x10000) >> 16;
2441 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2442 borrow = (z & 0x10000) >> 16;
2447 y = *bx - (ys & 0xffff) - borrow;
2448 borrow = (y & 0x10000) >> 16;
2457 while(--bxe > bx && !*bxe)
2465 #if !defined(MULTIPLE_THREADS) && !defined(NO_GLOBAL_STATE)
2466 #define USE_DTOA_RESULT 1
2467 static char *dtoa_result;
2472 rv_alloc(STATE_PARAM i) STATE_PARAM_DECL int i;
2474 rv_alloc(STATE_PARAM int i)
2481 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned) i;
2484 r = (int*)Balloc(PASS_STATE k);
2487 #ifdef USE_DTOA_RESULT
2495 nrv_alloc(STATE_PARAM s, rve, n) STATE_PARAM_DECL char *s, **rve; int n;
2497 nrv_alloc(STATE_PARAM CONST char *s, char **rve, int n)
2502 t = rv = rv_alloc(PASS_STATE n);
2503 while((*t = *s++)) t++;
2509 /* freedtoa(s) must be used to free values s returned by dtoa
2510 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2511 * but for consistency with earlier versions of dtoa, it is optional
2512 * when MULTIPLE_THREADS is not defined.
2517 freedtoa(STATE_PARAM s) STATE_PARAM_DECL char *s;
2519 freedtoa(STATE_PARAM char *s)
2522 Bigint *b = (Bigint *)((int *)s - 1);
2523 b->maxwds = 1 << (b->k = *(int*)b);
2524 Bfree(PASS_STATE b);
2525 #ifdef USE_DTOA_RESULT
2526 if (s == dtoa_result)
2531 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2533 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2534 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2537 * 1. Rather than iterating, we use a simple numeric overestimate
2538 * to determine k = floor(log10(d)). We scale relevant
2539 * quantities using O(log2(k)) rather than O(k) multiplications.
2540 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2541 * try to generate digits strictly left to right. Instead, we
2542 * compute with fewer bits and propagate the carry if necessary
2543 * when rounding the final digit up. This is often faster.
2544 * 3. Under the assumption that input will be rounded nearest,
2545 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2546 * That is, we allow equality in stopping tests when the
2547 * round-nearest rule will give the same floating-point value
2548 * as would satisfaction of the stopping test with strict
2550 * 4. We remove common factors of powers of 2 from relevant
2552 * 5. When converting floating-point integers less than 1e16,
2553 * we use floating-point arithmetic rather than resorting
2554 * to multiple-precision integers.
2555 * 6. When asked to produce fewer than 15 digits, we first try
2556 * to get by with floating-point arithmetic; we resort to
2557 * multiple-precision integer arithmetic only if we cannot
2558 * guarantee that the floating-point calculation has given
2559 * the correctly rounded result. For k requested digits and
2560 * "uniformly" distributed input, the probability is
2561 * something like 10^(k-15) that we must resort to the Long
2568 (STATE_PARAM d, mode, ndigits, decpt, sign, rve)
2569 STATE_PARAM_DECL U d; int mode, ndigits, *decpt, *sign; char **rve;
2571 (STATE_PARAM U d, int mode, int ndigits, int *decpt, int *sign, char **rve)
2574 /* Arguments ndigits, decpt, sign are similar to those
2575 of ecvt and fcvt; trailing zeros are suppressed from
2576 the returned string. If not null, *rve is set to point
2577 to the end of the return value. If d is +-Infinity or NaN,
2578 then *decpt is set to 9999.
2581 0 ==> shortest string that yields d when read in
2582 and rounded to nearest.
2583 1 ==> like 0, but with Steele & White stopping rule;
2584 e.g. with IEEE P754 arithmetic , mode 0 gives
2585 1e23 whereas mode 1 gives 9.999999999999999e22.
2586 2 ==> max(1,ndigits) significant digits. This gives a
2587 return value similar to that of ecvt, except
2588 that trailing zeros are suppressed.
2589 3 ==> through ndigits past the decimal point. This
2590 gives a return value similar to that from fcvt,
2591 except that trailing zeros are suppressed, and
2592 ndigits can be negative.
2593 4,5 ==> similar to 2 and 3, respectively, but (in
2594 round-nearest mode) with the tests of mode 0 to
2595 possibly return a shorter string that rounds to d.
2596 With IEEE arithmetic and compilation with
2597 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2598 as modes 2 and 3 when FLT_ROUNDS != 1.
2599 6-9 ==> Debugging modes similar to mode - 4: don't try
2600 fast floating-point estimate (if applicable).
2602 Values of mode other than 0-9 are treated as mode 0.
2604 Sufficient space is allocated to the return value
2605 to hold the suppressed trailing zeros.
2608 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
2609 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2610 spec_case, try_quick;
2612 #ifndef Sudden_Underflow
2616 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
2620 #ifdef Honor_FLT_ROUNDS
2624 int inexact, oldinexact;
2632 #ifdef USE_DTOA_RESULT
2634 freedtoa(PASS_STATE dtoa_result);
2639 if (word0(d) & Sign_bit) {
2640 /* set sign for everything, including 0's and NaNs */
2642 word0(d) &= ~Sign_bit; /* clear sign bit */
2647 #if defined(IEEE_Arith) + defined(VAX)
2649 if ((word0(d) & Exp_mask) == Exp_mask)
2651 if (word0(d) == 0x8000)
2654 /* Infinity or NaN */
2657 if (!word1(d) && !(word0(d) & 0xfffff))
2658 return nrv_alloc(PASS_STATE "Infinity", rve, 8);
2660 return nrv_alloc(PASS_STATE "NaN", rve, 3);
2664 dval(d) += 0; /* normalize */
2668 return nrv_alloc(PASS_STATE "0", rve, 1);
2672 try_quick = oldinexact = get_inexact();
2675 #ifdef Honor_FLT_ROUNDS
2676 if ((rounding = Flt_Rounds) >= 2) {
2678 rounding = rounding == 2 ? 0 : 2;
2685 b = d2b(PASS_STATE d, &be, &bbits);
2686 #ifdef Sudden_Underflow
2687 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
2689 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
2692 word0(d2) &= Frac_mask1;
2693 word0(d2) |= Exp_11;
2695 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
2699 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2700 * log10(x) = log(x) / log(10)
2701 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2702 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2704 * This suggests computing an approximation k to log10(d) by
2706 * k = (i - Bias)*0.301029995663981
2707 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2709 * We want k to be too large rather than too small.
2710 * The error in the first-order Taylor series approximation
2711 * is in our favor, so we just round up the constant enough
2712 * to compensate for any error in the multiplication of
2713 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2714 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2715 * adding 1e-13 to the constant term more than suffices.
2716 * Hence we adjust the constant term to 0.1760912590558.
2717 * (We could get a more accurate k by invoking log10,
2718 * but this is probably not worthwhile.)
2726 #ifndef Sudden_Underflow
2730 /* d is denormalized */
2732 i = bbits + be + (Bias + (P-1) - 1);
2733 x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
2734 : word1(d) << (32 - i);
2736 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2737 i -= (Bias + (P-1) - 1) + 1;
2741 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2743 if (ds < 0. && ds != k)
2744 k--; /* want k = floor(ds) */
2746 if (k >= 0 && k <= Ten_pmax) {
2747 if (dval(d) < tens[k])
2770 if (mode < 0 || mode > 9)
2774 #ifdef Check_FLT_ROUNDS
2775 try_quick = Rounding == 1;
2779 #endif /*SET_INEXACT*/
2799 ilim = ilim1 = i = ndigits;
2805 i = ndigits + k + 1;
2811 s = s0 = rv_alloc(PASS_STATE i);
2813 #ifdef Honor_FLT_ROUNDS
2814 if (mode > 1 && rounding != 1)
2818 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2820 /* Try to get by with floating-point arithmetic. */
2826 ieps = 2; /* conservative */
2831 /* prevent overflows */
2833 dval(d) /= bigtens[n_bigtens-1];
2836 for(; j; j >>= 1, i++)
2843 else if ((j1 = -k)) {
2844 dval(d) *= tens[j1 & 0xf];
2845 for(j = j1 >> 4; j; j >>= 1, i++)
2848 dval(d) *= bigtens[i];
2851 if (k_check && dval(d) < 1. && ilim > 0) {
2859 dval(eps) = ieps*dval(d) + 7.;
2860 word0(eps) -= (P-1)*Exp_msk1;
2864 if (dval(d) > dval(eps))
2866 if (dval(d) < -dval(eps))
2870 #ifndef No_leftright
2872 /* Use Steele & White method of only
2873 * generating digits needed.
2875 dval(eps) = 0.5/tens[ilim-1] - dval(eps);
2877 L = (ULong) dval(d);
2879 *s++ = '0' + (int)L;
2880 if (dval(d) < dval(eps))
2882 if (1. - dval(d) < dval(eps))
2892 /* Generate ilim digits, then fix them up. */
2893 dval(eps) *= tens[ilim-1];
2894 for(i = 1;; i++, dval(d) *= 10.) {
2895 L = (Long)(dval(d));
2896 if (!(dval(d) -= L))
2898 *s++ = '0' + (int)L;
2900 if (dval(d) > 0.5 + dval(eps))
2902 else if (dval(d) < 0.5 - dval(eps)) {
2910 #ifndef No_leftright
2920 /* Do we have a "small" integer? */
2922 if (be >= 0 && k <= Int_max) {
2925 if (ndigits < 0 && ilim <= 0) {
2927 if (ilim < 0 || dval(d) < 5*ds)
2931 for(i = 1;; i++, dval(d) *= 10.) {
2932 L = (Long)(dval(d) / ds);
2934 #ifdef Check_FLT_ROUNDS
2935 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2941 *s++ = '0' + (int)L;
2949 #ifdef Honor_FLT_ROUNDS
2953 case 2: goto bump_up;
2957 if (dval(d) > ds || (dval(d) == ds && L & 1)) {
2978 #ifndef Sudden_Underflow
2979 denorm ? be + (Bias + (P-1) - 1 + 1) :
2982 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2988 mhi = i2b(PASS_STATE 1);
2990 if (m2 > 0 && s2 > 0) {
2991 i = m2 < s2 ? m2 : s2;
2999 mhi = pow5mult(PASS_STATE mhi, m5);
3000 b1 = mult(PASS_STATE mhi, b);
3001 Bfree(PASS_STATE b);
3005 b = pow5mult(PASS_STATE b, j);
3008 b = pow5mult(PASS_STATE b, b5);
3010 S = i2b(PASS_STATE 1);
3012 S = pow5mult(PASS_STATE S, s5);
3014 /* Check for special case that d is a normalized power of 2. */
3017 if ((mode < 2 || leftright)
3018 #ifdef Honor_FLT_ROUNDS
3022 if (!word1(d) && !(word0(d) & Bndry_mask)
3023 #ifndef Sudden_Underflow
3024 && word0(d) & (Exp_mask & ~Exp_msk1)
3027 /* The special case */
3034 /* Arrange for convenient computation of quotients:
3035 * shift left if necessary so divisor has 4 leading 0 bits.
3037 * Perhaps we should just compute leading 28 bits of S once
3038 * and for all and pass them and a shift to quorem, so it
3039 * can do shifts and ors to compute the numerator for q.
3042 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
3045 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
3061 b = lshift(PASS_STATE b, b2);
3063 S = lshift(PASS_STATE S, s2);
3067 b = multadd(PASS_STATE b, 10, 0); /* we botched the k estimate */
3069 mhi = multadd(PASS_STATE mhi, 10, 0);
3073 if (ilim <= 0 && (mode == 3 || mode == 5)) {
3074 if (ilim < 0 || cmp(b,S = multadd(PASS_STATE S,5,0)) < 0) {
3075 /* no digits, fcvt style */
3077 /* MOZILLA CHANGE: Always return a non-empty string. */
3089 mhi = lshift(PASS_STATE mhi, m2);
3091 /* Compute mlo -- check for special case
3092 * that d is a normalized power of 2.
3097 mhi = Balloc(PASS_STATE mhi->k);
3099 mhi = lshift(PASS_STATE mhi, Log2P);
3103 dig = quorem(b,S) + '0';
3104 /* Do we yet have the shortest decimal string
3105 * that will round to d?
3108 delta = diff(PASS_STATE S, mhi);
3109 j1 = delta->sign ? 1 : cmp(b, delta);
3110 Bfree(PASS_STATE delta);
3111 #ifndef ROUND_BIASED
3112 if (j1 == 0 && mode != 1 && !(word1(d) & 1)
3113 #ifdef Honor_FLT_ROUNDS
3122 else if (!b->x[0] && b->wds <= 1)
3129 if (j < 0 || (j == 0 && mode != 1
3130 #ifndef ROUND_BIASED
3134 if (!b->x[0] && b->wds <= 1) {
3140 #ifdef Honor_FLT_ROUNDS
3143 case 0: goto accept_dig;
3144 case 2: goto keep_dig;
3146 #endif /*Honor_FLT_ROUNDS*/
3148 b = lshift(PASS_STATE b, 1);
3150 if ((j1 > 0 || (j1 == 0 && dig & 1))
3159 #ifdef Honor_FLT_ROUNDS
3163 if (dig == '9') { /* possible if i == 1 */
3171 #ifdef Honor_FLT_ROUNDS
3177 b = multadd(PASS_STATE b, 10, 0);
3179 mlo = mhi = multadd(PASS_STATE mhi, 10, 0);
3181 mlo = multadd(PASS_STATE mlo, 10, 0);
3182 mhi = multadd(PASS_STATE mhi, 10, 0);
3188 *s++ = dig = quorem(b,S) + '0';
3189 if (!b->x[0] && b->wds <= 1) {
3197 b = multadd(PASS_STATE b, 10, 0);
3200 /* Round off last digit */
3202 #ifdef Honor_FLT_ROUNDS
3204 case 0: goto trimzeros;
3205 case 2: goto roundoff;
3208 b = lshift(PASS_STATE b, 1);
3210 if (j >= 0) { /* ECMA compatible rounding needed by Spidermonkey */
3221 #ifdef Honor_FLT_ROUNDS
3228 Bfree(PASS_STATE S);
3230 if (mlo && mlo != mhi)
3231 Bfree(PASS_STATE mlo);
3232 Bfree(PASS_STATE mhi);
3238 word0(d) = Exp_1 + (70 << Exp_shift);
3243 else if (!oldinexact)
3246 Bfree(PASS_STATE b);