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 NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
120 * #defined automatically on IEEE systems. On such systems,
121 * when INFNAN_CHECK is #defined, strtod checks
122 * for Infinity and NaN (case insensitively). On some systems
123 * (e.g., some HP systems), it may be necessary to #define NAN_WORD0
124 * appropriately -- to the most significant word of a quiet NaN.
125 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
126 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
127 * strtod also accepts (case insensitively) strings of the form
128 * NaN(x), where x is a string of hexadecimal digits and spaces;
129 * if there is only one string of hexadecimal digits, it is taken
130 * for the 52 fraction bits of the resulting NaN; if there are two
131 * or more strings of hex digits, the first is for the high 20 bits,
132 * the second and subsequent for the low 32 bits, with intervening
133 * white space ignored; but if this results in none of the 52
134 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
135 * and NAN_WORD1 are used instead.
136 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
137 * multiple threads. In this case, you must provide (or suitably
138 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
139 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
140 * in pow5mult, ensures lazy evaluation of only one copy of high
141 * powers of 5; omitting this lock would introduce a small
142 * probability of wasting memory, but would otherwise be harmless.)
143 * You must also invoke freedtoa(s) to free the value s returned by
144 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
145 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
146 * avoids underflows on inputs whose result does not underflow.
147 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
148 * floating-point numbers and flushes underflows to zero rather
149 * than implementing gradual underflow, then you must also #define
151 * #define USE_LOCALE to use the current locale's decimal_point value.
152 * #define SET_INEXACT if IEEE arithmetic is being used and extra
153 * computation should be done to set the inexact flag when the
154 * result is inexact and avoid setting inexact when the result
155 * is exact. In this case, dtoa.c must be compiled in
156 * an environment, perhaps provided by #include "dtoa.c" in a
157 * suitable wrapper, that defines two functions,
158 * int get_inexact(void);
159 * void clear_inexact(void);
160 * such that get_inexact() returns a nonzero value if the
161 * inexact bit is already set, and clear_inexact() sets the
162 * inexact bit to 0. When SET_INEXACT is #defined, strtod
163 * also does extra computations to set the underflow and overflow
164 * flags when appropriate (i.e., when the result is tiny and
165 * inexact or when it is a numeric value rounded to +-infinity).
166 * #define NO_ERRNO if strtod should not assign errno = ERANGE when
167 * the result overflows to +-Infinity or underflows to 0.
168 * #define NO_GLOBAL_STATE to avoid defining any non-const global or
169 * static variables. Instead the necessary state is stored in an
170 * opaque struct, DtoaState, a pointer to which must be passed to
171 * every entry point. Two new functions are added to the API:
172 * DtoaState *newdtoa(void);
173 * void destroydtoa(DtoaState *);
180 typedef unsigned Long ULong
;
185 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
197 extern char *MALLOC();
199 extern void *MALLOC(size_t);
202 #define MALLOC malloc
209 #ifndef Omit_Private_Memory
211 #define PRIVATE_MEM 2304
213 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
217 #undef Avoid_Underflow
226 #ifndef NO_INFNAN_CHECK
240 #define DBL_MAX_10_EXP 308
241 #define DBL_MAX_EXP 1024
243 #endif /*IEEE_Arith*/
247 #define DBL_MAX_10_EXP 75
248 #define DBL_MAX_EXP 63
250 #define DBL_MAX 7.2370055773322621e+75
255 #define DBL_MAX_10_EXP 38
256 #define DBL_MAX_EXP 127
258 #define DBL_MAX 1.7014118346046923e+38
262 #define LONG_MAX 2147483647
265 #else /* ifndef Bad_float_h */
267 #endif /* Bad_float_h */
279 #define CONST /* blank */
285 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
286 Exactly one of IEEE_8087
, IEEE_MC68k
, VAX
, or IBM should be defined
.
289 typedef union { double d
; ULong L
[2]; } U
;
291 #define dval(x) ((x).d)
293 #define word0(x) ((x).L[1])
294 #define word1(x) ((x).L[0])
296 #define word0(x) ((x).L[0])
297 #define word1(x) ((x).L[1])
300 /* The following definition of Storeinc is appropriate for MIPS processors.
301 * An alternative that might be better on some machines is
302 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
304 #if defined(IEEE_8087) + defined(VAX)
305 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
306 ((unsigned short *)a)[0] = (unsigned short)c, a++)
308 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
309 ((unsigned short *)a)[1] = (unsigned short)c, a++)
312 /* #define P DBL_MANT_DIG */
313 /* Ten_pmax = floor(P*log(2)/log(5)) */
314 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
315 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
316 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
320 #define Exp_shift1 20
321 #define Exp_msk1 0x100000
322 #define Exp_msk11 0x100000
323 #define Exp_mask 0x7ff00000
327 #define Exp_1 0x3ff00000
328 #define Exp_11 0x3ff00000
330 #define Frac_mask 0xfffff
331 #define Frac_mask1 0xfffff
334 #define Bndry_mask 0xfffff
335 #define Bndry_mask1 0xfffff
337 #define Sign_bit 0x80000000
343 #ifndef NO_IEEE_Scale
344 #define Avoid_Underflow
345 #ifdef Flush_Denorm /* debugging option */
346 #undef Sudden_Underflow
352 #define Flt_Rounds FLT_ROUNDS
356 #endif /*Flt_Rounds*/
358 #ifdef Honor_FLT_ROUNDS
359 #define Rounding rounding
360 #undef Check_FLT_ROUNDS
361 #define Check_FLT_ROUNDS
363 #define Rounding Flt_Rounds
366 #else /* ifndef IEEE_Arith */
367 #undef Check_FLT_ROUNDS
368 #undef Honor_FLT_ROUNDS
370 #undef Sudden_Underflow
371 #define Sudden_Underflow
376 #define Exp_shift1 24
377 #define Exp_msk1 0x1000000
378 #define Exp_msk11 0x1000000
379 #define Exp_mask 0x7f000000
382 #define Exp_1 0x41000000
383 #define Exp_11 0x41000000
384 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
385 #define Frac_mask 0xffffff
386 #define Frac_mask1 0xffffff
389 #define Bndry_mask 0xefffff
390 #define Bndry_mask1 0xffffff
392 #define Sign_bit 0x80000000
394 #define Tiny0 0x100000
403 #define Exp_msk1 0x80
404 #define Exp_msk11 0x800000
405 #define Exp_mask 0x7f80
408 #define Exp_1 0x40800000
409 #define Exp_11 0x4080
411 #define Frac_mask 0x7fffff
412 #define Frac_mask1 0xffff007f
415 #define Bndry_mask 0xffff007f
416 #define Bndry_mask1 0xffff007f
418 #define Sign_bit 0x8000
424 #endif /* IBM, VAX */
425 #endif /* IEEE_Arith */
432 #define rounded_product(a,b) a = rnd_prod(a, b)
433 #define rounded_quotient(a,b) a = rnd_quot(a, b)
435 extern double rnd_prod(), rnd_quot();
437 extern double rnd_prod(double, double), rnd_quot(double, double);
440 #define rounded_product(a,b) a *= b
441 #define rounded_quotient(a,b) a /= b
444 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
445 #define Big1 0xffffffff
452 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
454 #define FFFFFFFF 0xffffffffUL
461 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
462 * This makes some inner loops simpler and sometimes saves work
463 * during multiplications, but it often seems to make things slightly
464 * slower. Hence the default is now to store 32 bits per Long.
467 #else /* long long available */
469 #define Llong long long
472 #define ULLong unsigned Llong
474 #endif /* NO_LONG_LONG */
476 #ifndef MULTIPLE_THREADS
477 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
478 #define FREE_DTOA_LOCK(n) /*nothing*/
486 int k
, maxwds
, sign
, wds
;
490 typedef struct Bigint Bigint
;
492 #ifdef NO_GLOBAL_STATE
493 #ifdef MULTIPLE_THREADS
494 #error "cannot have both NO_GLOBAL_STATE and MULTIPLE_THREADS"
498 #define DECLARE_GLOBAL_STATE /* nothing */
500 #define DECLARE_GLOBAL_STATE static
503 DECLARE_GLOBAL_STATE Bigint
*freelist
[Kmax
+1];
504 DECLARE_GLOBAL_STATE Bigint
*p5s
;
505 #ifndef Omit_Private_Memory
506 DECLARE_GLOBAL_STATE
double private_mem
[PRIVATE_mem
];
507 DECLARE_GLOBAL_STATE
double *pmem_next
508 #ifndef NO_GLOBAL_STATE
513 #ifdef NO_GLOBAL_STATE
515 typedef struct DtoaState DtoaState
;
517 #define STATE_PARAM state,
518 #define STATE_PARAM_DECL DtoaState *state;
520 #define STATE_PARAM DtoaState *state,
522 #define PASS_STATE state,
523 #define GET_STATE(field) (state->field)
528 DtoaState
*state
= (DtoaState
*) MALLOC(sizeof(DtoaState
));
530 memset(state
, 0, sizeof(DtoaState
));
531 state
->pmem_next
= state
->private_mem
;
539 (state
) STATE_PARAM_DECL
547 for (i
= 0; i
<= Kmax
; i
++) {
548 for (v
= GET_STATE(freelist
)[i
]; v
; v
= next
) {
550 #ifndef Omit_Private_Memory
551 if ((double*)v
< GET_STATE(private_mem
) ||
552 (double*)v
>= GET_STATE(private_mem
) + PRIVATE_mem
)
561 #define STATE_PARAM /* nothing */
562 #define STATE_PARAM_DECL /* nothing */
563 #define PASS_STATE /* nothing */
564 #define GET_STATE(name) name
570 (STATE_PARAM k
) STATE_PARAM_DECL
int k
;
577 #ifndef Omit_Private_Memory
581 ACQUIRE_DTOA_LOCK(0);
582 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
583 /* but this case seems very unlikely. */
584 if (k
<= Kmax
&& (rv
= GET_STATE(freelist
)[k
]))
585 GET_STATE(freelist
)[k
] = rv
->next
;
588 #ifdef Omit_Private_Memory
589 rv
= (Bigint
*)MALLOC(sizeof(Bigint
) + (x
-1)*sizeof(ULong
));
591 len
= (sizeof(Bigint
) + (x
-1)*sizeof(ULong
) + sizeof(double) - 1)
593 if (k
<= Kmax
&& GET_STATE(pmem_next
) - GET_STATE(private_mem
) + len
<= PRIVATE_mem
) {
594 rv
= (Bigint
*)GET_STATE(pmem_next
);
595 GET_STATE(pmem_next
) += len
;
598 rv
= (Bigint
*)MALLOC(len
*sizeof(double));
604 rv
->sign
= rv
->wds
= 0;
611 (STATE_PARAM v
) STATE_PARAM_DECL Bigint
*v
;
613 (STATE_PARAM Bigint
*v
)
620 ACQUIRE_DTOA_LOCK(0);
621 v
->next
= GET_STATE(freelist
)[v
->k
];
622 GET_STATE(freelist
)[v
->k
] = v
;
628 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
629 y->wds*sizeof(Long) + 2*sizeof(int))
634 (STATE_PARAM b
, m
, a
) STATE_PARAM_DECL Bigint
*b
; int m
, a
;
636 (STATE_PARAM Bigint
*b
, int m
, int a
) /* multiply by m and add a */
657 y
= *x
* (ULLong
)m
+ carry
;
659 *x
++ = (ULong
) y
& FFFFFFFF
;
663 y
= (xi
& 0xffff) * m
+ carry
;
664 z
= (xi
>> 16) * m
+ (y
>> 16);
666 *x
++ = (z
<< 16) + (y
& 0xffff);
676 if (wds
>= b
->maxwds
) {
677 b1
= Balloc(PASS_STATE b
->k
+1);
682 b
->x
[wds
++] = (ULong
) carry
;
691 (STATE_PARAM s
, nd0
, nd
, y9
) STATE_PARAM_DECL CONST
char *s
; int nd0
, nd
; ULong y9
;
693 (STATE_PARAM CONST
char *s
, int nd0
, int nd
, ULong y9
)
701 for(k
= 0, y
= 1; x
> y
; y
<<= 1, k
++) ;
703 b
= Balloc(PASS_STATE k
);
707 b
= Balloc(PASS_STATE k
+1);
708 b
->x
[0] = y9
& 0xffff;
709 b
->wds
= (b
->x
[1] = y9
>> 16) ? 2 : 1;
715 do b
= multadd(PASS_STATE b
, 10, *s
++ - '0');
722 b
= multadd(PASS_STATE b
, 10, *s
++ - '0');
729 (x
) register ULong x
;
736 if (!(x
& 0xffff0000)) {
740 if (!(x
& 0xff000000)) {
744 if (!(x
& 0xf0000000)) {
748 if (!(x
& 0xc0000000)) {
752 if (!(x
& 0x80000000)) {
754 if (!(x
& 0x40000000))
769 register ULong x
= *y
;
811 (STATE_PARAM i
) STATE_PARAM_DECL
int i
;
818 b
= Balloc(PASS_STATE
1);
827 (STATE_PARAM a
, b
) STATE_PARAM_DECL Bigint
*a
, *b
;
829 (STATE_PARAM Bigint
*a
, Bigint
*b
)
834 ULong
*x
, *xa
, *xae
, *xb
, *xbe
, *xc
, *xc0
;
845 if (a
->wds
< b
->wds
) {
856 c
= Balloc(PASS_STATE k
);
857 for(x
= c
->x
, xa
= x
+ wc
; x
< xa
; x
++)
865 for(; xb
< xbe
; xc0
++) {
871 z
= *x
++ * (ULLong
)y
+ *xc
+ carry
;
873 *xc
++ = (ULong
) z
& FFFFFFFF
;
881 for(; xb
< xbe
; xb
++, xc0
++) {
882 if (y
= *xb
& 0xffff) {
887 z
= (*x
& 0xffff) * y
+ (*xc
& 0xffff) + carry
;
889 z2
= (*x
++ >> 16) * y
+ (*xc
>> 16) + carry
;
902 z
= (*x
& 0xffff) * y
+ (*xc
>> 16) + carry
;
905 z2
= (*x
++ >> 16) * y
+ (*xc
& 0xffff) + carry
;
913 for(; xb
< xbe
; xc0
++) {
919 z
= *x
++ * y
+ *xc
+ carry
;
929 for(xc0
= c
->x
, xc
= xc0
+ wc
; wc
> 0 && !*--xc
; --wc
) ;
937 (STATE_PARAM b
, k
) STATE_PARAM_DECL Bigint
*b
; int k
;
939 (STATE_PARAM Bigint
*b
, int k
)
942 Bigint
*b1
, *p5
, *p51
;
944 static CONST
int p05
[3] = { 5, 25, 125 };
947 b
= multadd(PASS_STATE b
, p05
[i
-1], 0);
951 if (!(p5
= GET_STATE(p5s
))) {
953 #ifdef MULTIPLE_THREADS
954 ACQUIRE_DTOA_LOCK(1);
961 p5
= GET_STATE(p5s
) = i2b(PASS_STATE
625);
967 b1
= mult(PASS_STATE b
, p5
);
973 if (!(p51
= p5
->next
)) {
974 #ifdef MULTIPLE_THREADS
975 ACQUIRE_DTOA_LOCK(1);
976 if (!(p51
= p5
->next
)) {
977 p51
= p5
->next
= mult(p5
,p5
);
982 p51
= p5
->next
= mult(PASS_STATE p5
,p5
);
994 (STATE_PARAM b
, k
) STATE_PARAM_DECL Bigint
*b
; int k
;
996 (STATE_PARAM Bigint
*b
, int k
)
1001 ULong
*x
, *x1
, *xe
, z
;
1009 n1
= n
+ b
->wds
+ 1;
1010 for(i
= b
->maxwds
; n1
> i
; i
<<= 1)
1012 b1
= Balloc(PASS_STATE k1
);
1014 for(i
= 0; i
< n
; i
++)
1023 *x1
++ = *x
<< k
| z
;
1035 *x1
++ = *x
<< k
& 0xffff | z
;
1047 Bfree(PASS_STATE b
);
1054 (a
, b
) Bigint
*a
, *b
;
1056 (Bigint
*a
, Bigint
*b
)
1059 ULong
*xa
, *xa0
, *xb
, *xb0
;
1065 if (i
> 1 && !a
->x
[i
-1])
1066 Bug("cmp called with a->x[a->wds-1] == 0");
1067 if (j
> 1 && !b
->x
[j
-1])
1068 Bug("cmp called with b->x[b->wds-1] == 0");
1078 return *xa
< *xb
? -1 : 1;
1088 (STATE_PARAM a
, b
) STATE_PARAM_DECL Bigint
*a
, *b
;
1090 (STATE_PARAM Bigint
*a
, Bigint
*b
)
1095 ULong
*xa
, *xae
, *xb
, *xbe
, *xc
;
1107 c
= Balloc(PASS_STATE
0);
1120 c
= Balloc(PASS_STATE a
->k
);
1132 y
= (ULLong
)*xa
++ - *xb
++ - borrow
;
1133 borrow
= y
>> 32 & (ULong
)1;
1134 *xc
++ = (ULong
) y
& FFFFFFFF
;
1139 borrow
= y
>> 32 & (ULong
)1;
1140 *xc
++ = (ULong
) y
& FFFFFFFF
;
1145 y
= (*xa
& 0xffff) - (*xb
& 0xffff) - borrow
;
1146 borrow
= (y
& 0x10000) >> 16;
1147 z
= (*xa
++ >> 16) - (*xb
++ >> 16) - borrow
;
1148 borrow
= (z
& 0x10000) >> 16;
1153 y
= (*xa
& 0xffff) - borrow
;
1154 borrow
= (y
& 0x10000) >> 16;
1155 z
= (*xa
++ >> 16) - borrow
;
1156 borrow
= (z
& 0x10000) >> 16;
1161 y
= *xa
++ - *xb
++ - borrow
;
1162 borrow
= (y
& 0x10000) >> 16;
1168 borrow
= (y
& 0x10000) >> 16;
1190 L
= (word0(x
) & Exp_mask
) - (P
-1)*Exp_msk1
;
1191 #ifndef Avoid_Underflow
1192 #ifndef Sudden_Underflow
1201 #ifndef Avoid_Underflow
1202 #ifndef Sudden_Underflow
1205 L
= -L
>> Exp_shift
;
1206 if (L
< Exp_shift
) {
1207 word0(a
) = 0x80000 >> L
;
1213 word1(a
) = L
>= 31 ? 1 : 1 << 31 - L
;
1224 (a
, e
) Bigint
*a
; int *e
;
1229 ULong
*xa
, *xa0
, w
, y
, z
;
1243 if (!y
) Bug("zero y in b2d");
1249 d0
= Exp_1
| y
>> (Ebits
- k
);
1250 w
= xa
> xa0
? *--xa
: 0;
1251 d1
= y
<< ((32-Ebits
) + k
) | w
>> (Ebits
- k
);
1254 z
= xa
> xa0
? *--xa
: 0;
1256 d0
= Exp_1
| y
<< k
| z
>> (32 - k
);
1257 y
= xa
> xa0
? *--xa
: 0;
1258 d1
= z
<< k
| y
>> (32 - k
);
1265 if (k
< Ebits
+ 16) {
1266 z
= xa
> xa0
? *--xa
: 0;
1267 d0
= Exp_1
| y
<< k
- Ebits
| z
>> Ebits
+ 16 - k
;
1268 w
= xa
> xa0
? *--xa
: 0;
1269 y
= xa
> xa0
? *--xa
: 0;
1270 d1
= z
<< k
+ 16 - Ebits
| w
<< k
- Ebits
| y
>> 16 + Ebits
- k
;
1273 z
= xa
> xa0
? *--xa
: 0;
1274 w
= xa
> xa0
? *--xa
: 0;
1276 d0
= Exp_1
| y
<< k
+ 16 | z
<< k
| w
>> 16 - k
;
1277 y
= xa
> xa0
? *--xa
: 0;
1278 d1
= w
<< k
+ 16 | y
<< k
;
1282 word0(d
) = d0
>> 16 | d0
<< 16;
1283 word1(d
) = d1
>> 16 | d1
<< 16;
1294 (STATE_PARAM d
, e
, bits
) STATE_PARAM_DECL U d
; int *e
, *bits
;
1296 (STATE_PARAM U d
, int *e
, int *bits
)
1302 #ifndef Sudden_Underflow
1307 d0
= word0(d
) >> 16 | word0(d
) << 16;
1308 d1
= word1(d
) >> 16 | word1(d
) << 16;
1315 b
= Balloc(PASS_STATE
1);
1317 b
= Balloc(PASS_STATE
2);
1322 d0
&= 0x7fffffff; /* clear sign bit, which we ignore */
1323 #ifdef Sudden_Underflow
1324 de
= (int)(d0
>> Exp_shift
);
1329 if ((de
= (int)(d0
>> Exp_shift
)))
1334 if ((k
= lo0bits(&y
))) {
1335 x
[0] = y
| z
<< (32 - k
);
1340 #ifndef Sudden_Underflow
1343 b
->wds
= (x
[1] = z
) ? 2 : 1;
1348 #ifndef Sudden_Underflow
1356 if (k
= lo0bits(&y
))
1358 x
[0] = y
| z
<< 32 - k
& 0xffff;
1359 x
[1] = z
>> k
- 16 & 0xffff;
1365 x
[1] = y
>> 16 | z
<< 16 - k
& 0xffff;
1366 x
[2] = z
>> k
& 0xffff;
1381 Bug("Zero passed to d2b");
1399 #ifndef Sudden_Underflow
1403 *e
= (de
- Bias
- (P
-1) << 2) + k
;
1404 *bits
= 4*P
+ 8 - k
- hi0bits(word0(d
) & Frac_mask
);
1406 *e
= de
- Bias
- (P
-1) + k
;
1409 #ifndef Sudden_Underflow
1412 *e
= de
- Bias
- (P
-1) + 1 + k
;
1414 *bits
= 32*i
- hi0bits(x
[i
-1]);
1416 *bits
= (i
+2)*16 - hi0bits(x
[i
]);
1428 (a
, b
) Bigint
*a
, *b
;
1430 (Bigint
*a
, Bigint
*b
)
1436 dval(da
) = b2d(a
, &ka
);
1437 dval(db
) = b2d(b
, &kb
);
1439 k
= ka
- kb
+ 32*(a
->wds
- b
->wds
);
1441 k
= ka
- kb
+ 16*(a
->wds
- b
->wds
);
1445 word0(da
) += (k
>> 2)*Exp_msk1
;
1451 word0(db
) += (k
>> 2)*Exp_msk1
;
1457 word0(da
) += k
*Exp_msk1
;
1460 word0(db
) += k
*Exp_msk1
;
1463 return dval(da
) / dval(db
);
1468 1e0
, 1e1
, 1e2
, 1e3
, 1e4
, 1e5
, 1e6
, 1e7
, 1e8
, 1e9
,
1469 1e10
, 1e11
, 1e12
, 1e13
, 1e14
, 1e15
, 1e16
, 1e17
, 1e18
, 1e19
,
1478 bigtens
[] = { 1e16
, 1e32
, 1e64
, 1e128
, 1e256
};
1479 static CONST
double tinytens
[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1480 #ifdef Avoid_Underflow
1481 9007199254740992.*9007199254740992.e
-256
1482 /* = 2^106 * 1e-53 */
1487 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1488 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1489 #define Scale_Bit 0x10
1493 bigtens
[] = { 1e16
, 1e32
, 1e64
};
1494 static CONST
double tinytens
[] = { 1e-16, 1e-32, 1e-64 };
1497 bigtens
[] = { 1e16
, 1e32
};
1498 static CONST
double tinytens
[] = { 1e-16, 1e-32 };
1506 #define NAN_WORD0 0x7ff80000
1516 (sp
, t
) char **sp
, *t
;
1518 (CONST
char **sp
, CONST
char *t
)
1522 CONST
char *s
= *sp
;
1525 if ((c
= *++s
) >= 'A' && c
<= 'Z')
1538 (rvp
, sp
) U
*rvp
; CONST
char **sp
;
1540 (U
*rvp
, CONST
char **sp
)
1545 int havedig
, udx0
, xshift
;
1548 havedig
= xshift
= 0;
1551 /* allow optional initial 0x or 0X */
1552 while((c
= *(CONST
unsigned char*)(s
+1)) && c
<= ' ')
1554 if (s
[1] == '0' && (s
[2] == 'x' || s
[2] == 'X'))
1556 while((c
= *(CONST
unsigned char*)++s
)) {
1557 if (c
>= '0' && c
<= '9')
1559 else if (c
>= 'a' && c
<= 'f')
1561 else if (c
>= 'A' && c
<= 'F')
1563 else if (c
<= ' ') {
1564 if (udx0
&& havedig
) {
1570 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
1571 else if (/*(*/ c
== ')' && havedig
) {
1576 return; /* invalid form: don't change *sp */
1580 if (/*(*/ c
== ')') {
1584 } while((c
= *++s
));
1595 x
[0] = (x
[0] << 4) | (x
[1] >> 28);
1596 x
[1] = (x
[1] << 4) | c
;
1598 if ((x
[0] &= 0xfffff) || x
[1]) {
1599 word0(*rvp
) = Exp_mask
| x
[0];
1603 #endif /*No_Hex_NaN*/
1604 #endif /* INFNAN_CHECK */
1609 (STATE_PARAM s00
, se
) STATE_PARAM_DECL CONST
char *s00
; char **se
;
1611 (STATE_PARAM CONST
char *s00
, char **se
)
1614 #ifdef Avoid_Underflow
1617 int bb2
, bb5
, bbe
, bd2
, bd5
, bbbits
, bs2
, c
, dsign
,
1618 e
, e1
, esign
, i
, j
, k
, nd
, nd0
, nf
, nz
, nz0
, sign
;
1619 CONST
char *s
, *s0
, *s1
;
1624 Bigint
*bb
, *bb1
, *bd
, *bd0
, *bs
, *delta
;
1626 int inexact
, oldinexact
;
1628 #ifdef Honor_FLT_ROUNDS
1636 delta
= bb
= bd
= bs
= 0;
1639 sign
= nz0
= nz
= 0;
1641 for(s
= s00
;;s
++) switch(*s
) {
1664 while(*++s
== '0') ;
1670 for(nd
= nf
= 0; (c
= *s
) >= '0' && c
<= '9'; nd
++, s
++)
1677 s1
= localeconv()->decimal_point
;
1698 for(; c
== '0'; c
= *++s
)
1700 if (c
> '0' && c
<= '9') {
1708 for(; c
>= '0' && c
<= '9'; c
= *++s
) {
1713 for(i
= 1; i
< nz
; i
++)
1716 else if (nd
<= DBL_DIG
+ 1)
1720 else if (nd
<= DBL_DIG
+ 1)
1728 if (c
== 'e' || c
== 'E') {
1729 if (!nd
&& !nz
&& !nz0
) {
1740 if (c
>= '0' && c
<= '9') {
1743 if (c
> '0' && c
<= '9') {
1746 while((c
= *++s
) >= '0' && c
<= '9')
1748 if (s
- s1
> 8 || L
> 19999)
1749 /* Avoid confusion from exponents
1750 * so large that e might overflow.
1752 e
= 19999; /* safe for 16 bit ints */
1767 /* Check for Nan and Infinity */
1771 if (match(&s
,"nf")) {
1773 if (!match(&s
,"inity"))
1775 word0(rv
) = 0x7ff00000;
1782 if (match(&s
, "an")) {
1783 word0(rv
) = NAN_WORD0
;
1784 word1(rv
) = NAN_WORD1
;
1786 if (*s
== '(') /*)*/
1792 #endif /* INFNAN_CHECK */
1801 /* Now we have nd0 digits, starting at s0, followed by a
1802 * decimal point, followed by nd-nd0 digits. The number we're
1803 * after is the integer represented by those digits times
1808 k
= nd
< DBL_DIG
+ 1 ? nd
: DBL_DIG
+ 1;
1813 oldinexact
= get_inexact();
1815 dval(rv
) = tens
[k
- 9] * dval(rv
) + z
;
1819 #ifndef RND_PRODQUOT
1820 #ifndef Honor_FLT_ROUNDS
1828 if (e
<= Ten_pmax
) {
1830 goto vax_ovfl_check
;
1832 #ifdef Honor_FLT_ROUNDS
1833 /* round correctly FLT_ROUNDS = 2 or 3 */
1839 /* rv = */ rounded_product(dval(rv
), tens
[e
]);
1844 if (e
<= Ten_pmax
+ i
) {
1845 /* A fancier test would sometimes let us do
1846 * this for larger i values.
1848 #ifdef Honor_FLT_ROUNDS
1849 /* round correctly FLT_ROUNDS = 2 or 3 */
1856 dval(rv
) *= tens
[i
];
1858 /* VAX exponent range is so narrow we must
1859 * worry about overflow here...
1862 word0(rv
) -= P
*Exp_msk1
;
1863 /* rv = */ rounded_product(dval(rv
), tens
[e
]);
1864 if ((word0(rv
) & Exp_mask
)
1865 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
))
1867 word0(rv
) += P
*Exp_msk1
;
1869 /* rv = */ rounded_product(dval(rv
), tens
[e
]);
1874 #ifndef Inaccurate_Divide
1875 else if (e
>= -Ten_pmax
) {
1876 #ifdef Honor_FLT_ROUNDS
1877 /* round correctly FLT_ROUNDS = 2 or 3 */
1883 /* rv = */ rounded_quotient(dval(rv
), tens
[-e
]);
1894 oldinexact
= get_inexact();
1896 #ifdef Avoid_Underflow
1899 #ifdef Honor_FLT_ROUNDS
1900 if ((rounding
= Flt_Rounds
) >= 2) {
1902 rounding
= rounding
== 2 ? 0 : 2;
1908 #endif /*IEEE_Arith*/
1910 /* Get starting approximation = rv * 10**e1 */
1914 dval(rv
) *= tens
[i
];
1916 if (e1
> DBL_MAX_10_EXP
) {
1921 /* Can't trust HUGE_VAL */
1923 #ifdef Honor_FLT_ROUNDS
1925 case 0: /* toward 0 */
1926 case 3: /* toward -infinity */
1931 word0(rv
) = Exp_mask
;
1934 #else /*Honor_FLT_ROUNDS*/
1935 word0(rv
) = Exp_mask
;
1937 #endif /*Honor_FLT_ROUNDS*/
1939 /* set overflow bit */
1941 dval(rv0
) *= dval(rv0
);
1943 #else /*IEEE_Arith*/
1946 #endif /*IEEE_Arith*/
1952 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
1954 dval(rv
) *= bigtens
[j
];
1955 /* The last multiplication could overflow. */
1956 word0(rv
) -= P
*Exp_msk1
;
1957 dval(rv
) *= bigtens
[j
];
1958 if ((z
= word0(rv
) & Exp_mask
)
1959 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
))
1961 if (z
> Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
)) {
1962 /* set to largest number */
1963 /* (Can't trust DBL_MAX) */
1968 word0(rv
) += P
*Exp_msk1
;
1974 dval(rv
) /= tens
[i
];
1976 if (e1
>= 1 << n_bigtens
)
1978 #ifdef Avoid_Underflow
1981 for(j
= 0; e1
> 0; j
++, e1
>>= 1)
1983 dval(rv
) *= tinytens
[j
];
1984 if (scale
&& (j
= 2*P
+ 1 - ((word0(rv
) & Exp_mask
)
1985 >> Exp_shift
)) > 0) {
1986 /* scaled rv is denormal; zap j low bits */
1990 word0(rv
) = (P
+2)*Exp_msk1
;
1992 word0(rv
) &= 0xffffffff << (j
-32);
1995 word1(rv
) &= 0xffffffff << j
;
1998 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
2000 dval(rv
) *= tinytens
[j
];
2001 /* The last multiplication could underflow. */
2002 dval(rv0
) = dval(rv
);
2003 dval(rv
) *= tinytens
[j
];
2005 dval(rv
) = 2.*dval(rv0
);
2006 dval(rv
) *= tinytens
[j
];
2018 #ifndef Avoid_Underflow
2021 /* The refinement below will clean
2022 * this approximation up.
2029 /* Now the hard part -- adjusting rv to the correct value.*/
2031 /* Put digits into bd: true value = bd * 10^e */
2033 bd0
= s2b(PASS_STATE s0
, nd0
, nd
, y
);
2036 bd
= Balloc(PASS_STATE bd0
->k
);
2038 bb
= d2b(PASS_STATE rv
, &bbe
, &bbbits
); /* rv = bb * 2^bbe */
2039 bs
= i2b(PASS_STATE
1);
2054 #ifdef Honor_FLT_ROUNDS
2058 #ifdef Avoid_Underflow
2060 i
= j
+ bbbits
- 1; /* logb(rv) */
2061 if (i
< Emin
) /* denormal */
2065 #else /*Avoid_Underflow*/
2066 #ifdef Sudden_Underflow
2068 j
= 1 + 4*P
- 3 - bbbits
+ ((bbe
+ bbbits
- 1) & 3);
2072 #else /*Sudden_Underflow*/
2074 i
= j
+ bbbits
- 1; /* logb(rv) */
2075 if (i
< Emin
) /* denormal */
2079 #endif /*Sudden_Underflow*/
2080 #endif /*Avoid_Underflow*/
2083 #ifdef Avoid_Underflow
2086 i
= bb2
< bd2
? bb2
: bd2
;
2095 bs
= pow5mult(PASS_STATE bs
, bb5
);
2096 bb1
= mult(PASS_STATE bs
, bb
);
2097 Bfree(PASS_STATE bb
);
2101 bb
= lshift(PASS_STATE bb
, bb2
);
2103 bd
= pow5mult(PASS_STATE bd
, bd5
);
2105 bd
= lshift(PASS_STATE bd
, bd2
);
2107 bs
= lshift(PASS_STATE bs
, bs2
);
2108 delta
= diff(PASS_STATE bb
, bd
);
2109 dsign
= delta
->sign
;
2112 #ifdef Honor_FLT_ROUNDS
2113 if (rounding
!= 1) {
2115 /* Error is less than an ulp */
2116 if (!delta
->x
[0] && delta
->wds
<= 1) {
2132 && !(word0(rv
) & Frac_mask
)) {
2133 y
= word0(rv
) & Exp_mask
;
2134 #ifdef Avoid_Underflow
2135 if (!scale
|| y
> 2*P
*Exp_msk1
)
2140 delta
= lshift(PASS_STATE delta
,Log2P
);
2141 if (cmp(delta
, bs
) <= 0)
2146 #ifdef Avoid_Underflow
2147 if (scale
&& (y
= word0(rv
) & Exp_mask
)
2149 word0(adj
) += (2*P
+1)*Exp_msk1
- y
;
2151 #ifdef Sudden_Underflow
2152 if ((word0(rv
) & Exp_mask
) <=
2154 word0(rv
) += P
*Exp_msk1
;
2155 dval(rv
) += adj
*ulp(rv
);
2156 word0(rv
) -= P
*Exp_msk1
;
2159 #endif /*Sudden_Underflow*/
2160 #endif /*Avoid_Underflow*/
2161 dval(rv
) += adj
*ulp(rv
);
2165 adj
= ratio(delta
, bs
);
2168 if (adj
<= 0x7ffffffe) {
2169 /* adj = rounding ? ceil(adj) : floor(adj); */
2172 if (!((rounding
>>1) ^ dsign
))
2177 #ifdef Avoid_Underflow
2178 if (scale
&& (y
= word0(rv
) & Exp_mask
) <= 2*P
*Exp_msk1
)
2179 word0(adj
) += (2*P
+1)*Exp_msk1
- y
;
2181 #ifdef Sudden_Underflow
2182 if ((word0(rv
) & Exp_mask
) <= P
*Exp_msk1
) {
2183 word0(rv
) += P
*Exp_msk1
;
2189 word0(rv
) -= P
*Exp_msk1
;
2192 #endif /*Sudden_Underflow*/
2193 #endif /*Avoid_Underflow*/
2201 #endif /*Honor_FLT_ROUNDS*/
2204 /* Error is less than half an ulp -- check for
2205 * special case of mantissa a power of two.
2207 if (dsign
|| word1(rv
) || word0(rv
) & Bndry_mask
2209 #ifdef Avoid_Underflow
2210 || (word0(rv
) & Exp_mask
) <= (2*P
+1)*Exp_msk1
2212 || (word0(rv
) & Exp_mask
) <= Exp_msk1
2217 if (!delta
->x
[0] && delta
->wds
<= 1)
2222 if (!delta
->x
[0] && delta
->wds
<= 1) {
2229 delta
= lshift(PASS_STATE delta
,Log2P
);
2230 if (cmp(delta
, bs
) > 0)
2235 /* exactly half-way between */
2237 if ((word0(rv
) & Bndry_mask1
) == Bndry_mask1
2239 #ifdef Avoid_Underflow
2240 (scale
&& (y
= word0(rv
) & Exp_mask
) <= 2*P
*Exp_msk1
)
2241 ? (0xffffffff & (0xffffffff << (2*P
+1-(y
>>Exp_shift
)))) :
2244 /*boundary case -- increment exponent*/
2245 word0(rv
) = (word0(rv
) & Exp_mask
)
2252 #ifdef Avoid_Underflow
2258 else if (!(word0(rv
) & Bndry_mask
) && !word1(rv
)) {
2260 /* boundary case -- decrement exponent */
2261 #ifdef Sudden_Underflow /*{{*/
2262 L
= word0(rv
) & Exp_mask
;
2266 #ifdef Avoid_Underflow
2267 if (L
<= (scale
? (2*P
+1)*Exp_msk1
: Exp_msk1
))
2270 #endif /*Avoid_Underflow*/
2274 #else /*Sudden_Underflow}{*/
2275 #ifdef Avoid_Underflow
2277 L
= word0(rv
) & Exp_mask
;
2278 if (L
<= (2*P
+1)*Exp_msk1
) {
2279 if (L
> (P
+2)*Exp_msk1
)
2280 /* round even ==> */
2283 /* rv = smallest denormal */
2287 #endif /*Avoid_Underflow*/
2288 L
= (word0(rv
) & Exp_mask
) - Exp_msk1
;
2289 #endif /*Sudden_Underflow}}*/
2290 word0(rv
) = L
| Bndry_mask1
;
2291 word1(rv
) = 0xffffffff;
2298 #ifndef ROUND_BIASED
2299 if (!(word1(rv
) & LSB
))
2303 dval(rv
) += ulp(rv
);
2304 #ifndef ROUND_BIASED
2306 dval(rv
) -= ulp(rv
);
2307 #ifndef Sudden_Underflow
2312 #ifdef Avoid_Underflow
2318 if ((aadj
= ratio(delta
, bs
)) <= 2.) {
2320 aadj
= dval(aadj1
) = 1.;
2321 else if (word1(rv
) || word0(rv
) & Bndry_mask
) {
2322 #ifndef Sudden_Underflow
2323 if (word1(rv
) == Tiny1
&& !word0(rv
))
2330 /* special case -- power of FLT_RADIX to be */
2331 /* rounded down... */
2333 if (aadj
< 2./FLT_RADIX
)
2334 aadj
= 1./FLT_RADIX
;
2337 dval(aadj1
) = -aadj
;
2342 dval(aadj1
) = dsign
? aadj
: -aadj
;
2343 #ifdef Check_FLT_ROUNDS
2345 case 2: /* towards +infinity */
2348 case 0: /* towards 0 */
2349 case 3: /* towards -infinity */
2353 if (Flt_Rounds
== 0)
2355 #endif /*Check_FLT_ROUNDS*/
2357 y
= word0(rv
) & Exp_mask
;
2359 /* Check for overflow */
2361 if (y
== Exp_msk1
*(DBL_MAX_EXP
+Bias
-1)) {
2362 dval(rv0
) = dval(rv
);
2363 word0(rv
) -= P
*Exp_msk1
;
2364 adj
= dval(aadj1
) * ulp(rv
);
2366 if ((word0(rv
) & Exp_mask
) >=
2367 Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
)) {
2368 if (word0(rv0
) == Big0
&& word1(rv0
) == Big1
)
2375 word0(rv
) += P
*Exp_msk1
;
2378 #ifdef Avoid_Underflow
2379 if (scale
&& y
<= 2*P
*Exp_msk1
) {
2380 if (aadj
<= 0x7fffffff) {
2381 if ((z
= (ULong
) aadj
) <= 0)
2384 dval(aadj1
) = dsign
? aadj
: -aadj
;
2386 word0(aadj1
) += (2*P
+1)*Exp_msk1
- y
;
2388 adj
= dval(aadj1
) * ulp(rv
);
2391 #ifdef Sudden_Underflow
2392 if ((word0(rv
) & Exp_mask
) <= P
*Exp_msk1
) {
2393 dval(rv0
) = dval(rv
);
2394 word0(rv
) += P
*Exp_msk1
;
2395 adj
= dval(aadj1
) * ulp(rv
);
2398 if ((word0(rv
) & Exp_mask
) < P
*Exp_msk1
)
2400 if ((word0(rv
) & Exp_mask
) <= P
*Exp_msk1
)
2403 if (word0(rv0
) == Tiny0
2404 && word1(rv0
) == Tiny1
)
2411 word0(rv
) -= P
*Exp_msk1
;
2414 adj
= dval(aadj1
) * ulp(rv
);
2417 #else /*Sudden_Underflow*/
2418 /* Compute adj so that the IEEE rounding rules will
2419 * correctly round rv + adj in some half-way cases.
2420 * If rv * ulp(rv) is denormalized (i.e.,
2421 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2422 * trouble from bits lost to denormalization;
2423 * example: 1.2e-307 .
2425 if (y
<= (P
-1)*Exp_msk1
&& aadj
> 1.) {
2426 dval(aadj1
) = (double)(int)(aadj
+ 0.5);
2428 dval(aadj1
) = -dval(aadj1
);
2430 adj
= dval(aadj1
) * ulp(rv
);
2432 #endif /*Sudden_Underflow*/
2433 #endif /*Avoid_Underflow*/
2435 z
= word0(rv
) & Exp_mask
;
2437 #ifdef Avoid_Underflow
2441 /* Can we stop now? */
2444 /* The tolerances below are conservative. */
2445 if (dsign
|| word1(rv
) || word0(rv
) & Bndry_mask
) {
2446 if (aadj
< .4999999 || aadj
> .5000001)
2449 else if (aadj
< .4999999/FLT_RADIX
)
2454 Bfree(PASS_STATE bb
);
2455 Bfree(PASS_STATE bd
);
2456 Bfree(PASS_STATE bs
);
2457 Bfree(PASS_STATE delta
);
2462 word0(rv0
) = Exp_1
+ (70 << Exp_shift
);
2467 else if (!oldinexact
)
2470 #ifdef Avoid_Underflow
2472 word0(rv0
) = Exp_1
- 2*P
*Exp_msk1
;
2474 dval(rv
) *= dval(rv0
);
2476 /* try to avoid the bug of testing an 8087 register value */
2477 if (word0(rv
) == 0 && word1(rv
) == 0)
2481 #endif /* Avoid_Underflow */
2483 if (inexact
&& !(word0(rv
) & Exp_mask
)) {
2484 /* set underflow bit */
2486 dval(rv0
) *= dval(rv0
);
2490 Bfree(PASS_STATE bb
);
2491 Bfree(PASS_STATE bd
);
2492 Bfree(PASS_STATE bs
);
2493 Bfree(PASS_STATE bd0
);
2494 Bfree(PASS_STATE delta
);
2498 return sign
? -dval(rv
) : dval(rv
);
2504 (b
, S
) Bigint
*b
, *S
;
2506 (Bigint
*b
, Bigint
*S
)
2510 ULong
*bx
, *bxe
, q
, *sx
, *sxe
;
2512 ULLong borrow
, carry
, y
, ys
;
2514 ULong borrow
, carry
, y
, ys
;
2522 /*debug*/ if (b
->wds
> n
)
2523 /*debug*/ Bug("oversize b in quorem");
2531 q
= *bxe
/ (*sxe
+ 1); /* ensure q <= true quotient */
2533 /*debug*/ if (q
> 9)
2534 /*debug*/ Bug("oversized quotient in quorem");
2541 ys
= *sx
++ * (ULLong
)q
+ carry
;
2543 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
2544 borrow
= y
>> 32 & (ULong
)1;
2545 *bx
++ = (ULong
) y
& FFFFFFFF
;
2549 ys
= (si
& 0xffff) * q
+ carry
;
2550 zs
= (si
>> 16) * q
+ (ys
>> 16);
2552 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
2553 borrow
= (y
& 0x10000) >> 16;
2554 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
2555 borrow
= (z
& 0x10000) >> 16;
2558 ys
= *sx
++ * q
+ carry
;
2560 y
= *bx
- (ys
& 0xffff) - borrow
;
2561 borrow
= (y
& 0x10000) >> 16;
2569 while(--bxe
> bx
&& !*bxe
)
2574 if (cmp(b
, S
) >= 0) {
2584 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
2585 borrow
= y
>> 32 & (ULong
)1;
2586 *bx
++ = (ULong
) y
& FFFFFFFF
;
2590 ys
= (si
& 0xffff) + carry
;
2591 zs
= (si
>> 16) + (ys
>> 16);
2593 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
2594 borrow
= (y
& 0x10000) >> 16;
2595 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
2596 borrow
= (z
& 0x10000) >> 16;
2601 y
= *bx
- (ys
& 0xffff) - borrow
;
2602 borrow
= (y
& 0x10000) >> 16;
2611 while(--bxe
> bx
&& !*bxe
)
2619 #if !defined(MULTIPLE_THREADS) && !defined(NO_GLOBAL_STATE)
2620 #define USE_DTOA_RESULT 1
2621 static char *dtoa_result
;
2626 rv_alloc(STATE_PARAM i
) STATE_PARAM_DECL
int i
;
2628 rv_alloc(STATE_PARAM
int i
)
2635 sizeof(Bigint
) - sizeof(ULong
) - sizeof(int) + j
<= (unsigned) i
;
2638 r
= (int*)Balloc(PASS_STATE k
);
2641 #ifdef USE_DTOA_RESULT
2649 nrv_alloc(STATE_PARAM s
, rve
, n
) STATE_PARAM_DECL
char *s
, **rve
; int n
;
2651 nrv_alloc(STATE_PARAM CONST
char *s
, char **rve
, int n
)
2656 t
= rv
= rv_alloc(PASS_STATE n
);
2657 while((*t
= *s
++)) t
++;
2663 /* freedtoa(s) must be used to free values s returned by dtoa
2664 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2665 * but for consistency with earlier versions of dtoa, it is optional
2666 * when MULTIPLE_THREADS is not defined.
2671 freedtoa(STATE_PARAM s
) STATE_PARAM_DECL
char *s
;
2673 freedtoa(STATE_PARAM
char *s
)
2676 Bigint
*b
= (Bigint
*)((int *)s
- 1);
2677 b
->maxwds
= 1 << (b
->k
= *(int*)b
);
2678 Bfree(PASS_STATE b
);
2679 #ifdef USE_DTOA_RESULT
2680 if (s
== dtoa_result
)
2685 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2687 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2688 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2691 * 1. Rather than iterating, we use a simple numeric overestimate
2692 * to determine k = floor(log10(d)). We scale relevant
2693 * quantities using O(log2(k)) rather than O(k) multiplications.
2694 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2695 * try to generate digits strictly left to right. Instead, we
2696 * compute with fewer bits and propagate the carry if necessary
2697 * when rounding the final digit up. This is often faster.
2698 * 3. Under the assumption that input will be rounded nearest,
2699 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2700 * That is, we allow equality in stopping tests when the
2701 * round-nearest rule will give the same floating-point value
2702 * as would satisfaction of the stopping test with strict
2704 * 4. We remove common factors of powers of 2 from relevant
2706 * 5. When converting floating-point integers less than 1e16,
2707 * we use floating-point arithmetic rather than resorting
2708 * to multiple-precision integers.
2709 * 6. When asked to produce fewer than 15 digits, we first try
2710 * to get by with floating-point arithmetic; we resort to
2711 * multiple-precision integer arithmetic only if we cannot
2712 * guarantee that the floating-point calculation has given
2713 * the correctly rounded result. For k requested digits and
2714 * "uniformly" distributed input, the probability is
2715 * something like 10^(k-15) that we must resort to the Long
2722 (STATE_PARAM d
, mode
, ndigits
, decpt
, sign
, rve
)
2723 STATE_PARAM_DECL U d
; int mode
, ndigits
, *decpt
, *sign
; char **rve
;
2725 (STATE_PARAM U d
, int mode
, int ndigits
, int *decpt
, int *sign
, char **rve
)
2728 /* Arguments ndigits, decpt, sign are similar to those
2729 of ecvt and fcvt; trailing zeros are suppressed from
2730 the returned string. If not null, *rve is set to point
2731 to the end of the return value. If d is +-Infinity or NaN,
2732 then *decpt is set to 9999.
2735 0 ==> shortest string that yields d when read in
2736 and rounded to nearest.
2737 1 ==> like 0, but with Steele & White stopping rule;
2738 e.g. with IEEE P754 arithmetic , mode 0 gives
2739 1e23 whereas mode 1 gives 9.999999999999999e22.
2740 2 ==> max(1,ndigits) significant digits. This gives a
2741 return value similar to that of ecvt, except
2742 that trailing zeros are suppressed.
2743 3 ==> through ndigits past the decimal point. This
2744 gives a return value similar to that from fcvt,
2745 except that trailing zeros are suppressed, and
2746 ndigits can be negative.
2747 4,5 ==> similar to 2 and 3, respectively, but (in
2748 round-nearest mode) with the tests of mode 0 to
2749 possibly return a shorter string that rounds to d.
2750 With IEEE arithmetic and compilation with
2751 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2752 as modes 2 and 3 when FLT_ROUNDS != 1.
2753 6-9 ==> Debugging modes similar to mode - 4: don't try
2754 fast floating-point estimate (if applicable).
2756 Values of mode other than 0-9 are treated as mode 0.
2758 Sufficient space is allocated to the return value
2759 to hold the suppressed trailing zeros.
2762 int bbits
, b2
, b5
, be
, dig
, i
, ieps
, ilim
, ilim0
, ilim1
,
2763 j
, j1
, k
, k0
, k_check
, leftright
, m2
, m5
, s2
, s5
,
2764 spec_case
, try_quick
;
2766 #ifndef Sudden_Underflow
2770 Bigint
*b
, *b1
, *delta
, *mlo
, *mhi
, *S
;
2774 #ifdef Honor_FLT_ROUNDS
2778 int inexact
, oldinexact
;
2786 #ifdef USE_DTOA_RESULT
2788 freedtoa(PASS_STATE dtoa_result
);
2793 if (word0(d
) & Sign_bit
) {
2794 /* set sign for everything, including 0's and NaNs */
2796 word0(d
) &= ~Sign_bit
; /* clear sign bit */
2801 #if defined(IEEE_Arith) + defined(VAX)
2803 if ((word0(d
) & Exp_mask
) == Exp_mask
)
2805 if (word0(d
) == 0x8000)
2808 /* Infinity or NaN */
2811 if (!word1(d
) && !(word0(d
) & 0xfffff))
2812 return nrv_alloc(PASS_STATE
"Infinity", rve
, 8);
2814 return nrv_alloc(PASS_STATE
"NaN", rve
, 3);
2818 dval(d
) += 0; /* normalize */
2822 return nrv_alloc(PASS_STATE
"0", rve
, 1);
2826 try_quick
= oldinexact
= get_inexact();
2829 #ifdef Honor_FLT_ROUNDS
2830 if ((rounding
= Flt_Rounds
) >= 2) {
2832 rounding
= rounding
== 2 ? 0 : 2;
2839 b
= d2b(PASS_STATE d
, &be
, &bbits
);
2840 #ifdef Sudden_Underflow
2841 i
= (int)(word0(d
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
));
2843 if ((i
= (int)(word0(d
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
)))) {
2846 word0(d2
) &= Frac_mask1
;
2847 word0(d2
) |= Exp_11
;
2849 if (j
= 11 - hi0bits(word0(d2
) & Frac_mask
))
2853 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2854 * log10(x) = log(x) / log(10)
2855 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2856 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2858 * This suggests computing an approximation k to log10(d) by
2860 * k = (i - Bias)*0.301029995663981
2861 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2863 * We want k to be too large rather than too small.
2864 * The error in the first-order Taylor series approximation
2865 * is in our favor, so we just round up the constant enough
2866 * to compensate for any error in the multiplication of
2867 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2868 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2869 * adding 1e-13 to the constant term more than suffices.
2870 * Hence we adjust the constant term to 0.1760912590558.
2871 * (We could get a more accurate k by invoking log10,
2872 * but this is probably not worthwhile.)
2880 #ifndef Sudden_Underflow
2884 /* d is denormalized */
2886 i
= bbits
+ be
+ (Bias
+ (P
-1) - 1);
2887 x
= i
> 32 ? word0(d
) << (64 - i
) | word1(d
) >> (i
- 32)
2888 : word1(d
) << (32 - i
);
2890 word0(d2
) -= 31*Exp_msk1
; /* adjust exponent */
2891 i
-= (Bias
+ (P
-1) - 1) + 1;
2895 ds
= (dval(d2
)-1.5)*0.289529654602168 + 0.1760912590558 + i
*0.301029995663981;
2897 if (ds
< 0. && ds
!= k
)
2898 k
--; /* want k = floor(ds) */
2900 if (k
>= 0 && k
<= Ten_pmax
) {
2901 if (dval(d
) < tens
[k
])
2924 if (mode
< 0 || mode
> 9)
2928 #ifdef Check_FLT_ROUNDS
2929 try_quick
= Rounding
== 1;
2933 #endif /*SET_INEXACT*/
2953 ilim
= ilim1
= i
= ndigits
;
2959 i
= ndigits
+ k
+ 1;
2965 s
= s0
= rv_alloc(PASS_STATE i
);
2967 #ifdef Honor_FLT_ROUNDS
2968 if (mode
> 1 && rounding
!= 1)
2972 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
) {
2974 /* Try to get by with floating-point arithmetic. */
2980 ieps
= 2; /* conservative */
2985 /* prevent overflows */
2987 dval(d
) /= bigtens
[n_bigtens
-1];
2990 for(; j
; j
>>= 1, i
++)
2997 else if ((j1
= -k
)) {
2998 dval(d
) *= tens
[j1
& 0xf];
2999 for(j
= j1
>> 4; j
; j
>>= 1, i
++)
3002 dval(d
) *= bigtens
[i
];
3005 if (k_check
&& dval(d
) < 1. && ilim
> 0) {
3013 dval(eps
) = ieps
*dval(d
) + 7.;
3014 word0(eps
) -= (P
-1)*Exp_msk1
;
3018 if (dval(d
) > dval(eps
))
3020 if (dval(d
) < -dval(eps
))
3024 #ifndef No_leftright
3026 /* Use Steele & White method of only
3027 * generating digits needed.
3029 dval(eps
) = 0.5/tens
[ilim
-1] - dval(eps
);
3031 L
= (ULong
) dval(d
);
3033 *s
++ = '0' + (int)L
;
3034 if (dval(d
) < dval(eps
))
3036 if (1. - dval(d
) < dval(eps
))
3046 /* Generate ilim digits, then fix them up. */
3047 dval(eps
) *= tens
[ilim
-1];
3048 for(i
= 1;; i
++, dval(d
) *= 10.) {
3049 L
= (Long
)(dval(d
));
3050 if (!(dval(d
) -= L
))
3052 *s
++ = '0' + (int)L
;
3054 if (dval(d
) > 0.5 + dval(eps
))
3056 else if (dval(d
) < 0.5 - dval(eps
)) {
3064 #ifndef No_leftright
3074 /* Do we have a "small" integer? */
3076 if (be
>= 0 && k
<= Int_max
) {
3079 if (ndigits
< 0 && ilim
<= 0) {
3081 if (ilim
< 0 || dval(d
) < 5*ds
)
3085 for(i
= 1;; i
++, dval(d
) *= 10.) {
3086 L
= (Long
)(dval(d
) / ds
);
3088 #ifdef Check_FLT_ROUNDS
3089 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
3095 *s
++ = '0' + (int)L
;
3103 #ifdef Honor_FLT_ROUNDS
3107 case 2: goto bump_up
;
3111 if (dval(d
) > ds
|| (dval(d
) == ds
&& L
& 1)) {
3132 #ifndef Sudden_Underflow
3133 denorm
? be
+ (Bias
+ (P
-1) - 1 + 1) :
3136 1 + 4*P
- 3 - bbits
+ ((bbits
+ be
- 1) & 3);
3142 mhi
= i2b(PASS_STATE
1);
3144 if (m2
> 0 && s2
> 0) {
3145 i
= m2
< s2
? m2
: s2
;
3153 mhi
= pow5mult(PASS_STATE mhi
, m5
);
3154 b1
= mult(PASS_STATE mhi
, b
);
3155 Bfree(PASS_STATE b
);
3159 b
= pow5mult(PASS_STATE b
, j
);
3162 b
= pow5mult(PASS_STATE b
, b5
);
3164 S
= i2b(PASS_STATE
1);
3166 S
= pow5mult(PASS_STATE S
, s5
);
3168 /* Check for special case that d is a normalized power of 2. */
3171 if ((mode
< 2 || leftright
)
3172 #ifdef Honor_FLT_ROUNDS
3176 if (!word1(d
) && !(word0(d
) & Bndry_mask
)
3177 #ifndef Sudden_Underflow
3178 && word0(d
) & (Exp_mask
& ~Exp_msk1
)
3181 /* The special case */
3188 /* Arrange for convenient computation of quotients:
3189 * shift left if necessary so divisor has 4 leading 0 bits.
3191 * Perhaps we should just compute leading 28 bits of S once
3192 * and for all and pass them and a shift to quorem, so it
3193 * can do shifts and ors to compute the numerator for q.
3196 if ((i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0x1f))
3199 if (i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0xf)
3215 b
= lshift(PASS_STATE b
, b2
);
3217 S
= lshift(PASS_STATE S
, s2
);
3221 b
= multadd(PASS_STATE b
, 10, 0); /* we botched the k estimate */
3223 mhi
= multadd(PASS_STATE mhi
, 10, 0);
3227 if (ilim
<= 0 && (mode
== 3 || mode
== 5)) {
3228 if (ilim
< 0 || cmp(b
,S
= multadd(PASS_STATE S
,5,0)) < 0) {
3229 /* no digits, fcvt style */
3231 /* MOZILLA CHANGE: Always return a non-empty string. */
3243 mhi
= lshift(PASS_STATE mhi
, m2
);
3245 /* Compute mlo -- check for special case
3246 * that d is a normalized power of 2.
3251 mhi
= Balloc(PASS_STATE mhi
->k
);
3253 mhi
= lshift(PASS_STATE mhi
, Log2P
);
3257 dig
= quorem(b
,S
) + '0';
3258 /* Do we yet have the shortest decimal string
3259 * that will round to d?
3262 delta
= diff(PASS_STATE S
, mhi
);
3263 j1
= delta
->sign
? 1 : cmp(b
, delta
);
3264 Bfree(PASS_STATE delta
);
3265 #ifndef ROUND_BIASED
3266 if (j1
== 0 && mode
!= 1 && !(word1(d
) & 1)
3267 #ifdef Honor_FLT_ROUNDS
3276 else if (!b
->x
[0] && b
->wds
<= 1)
3283 if (j
< 0 || (j
== 0 && mode
!= 1
3284 #ifndef ROUND_BIASED
3288 if (!b
->x
[0] && b
->wds
<= 1) {
3294 #ifdef Honor_FLT_ROUNDS
3297 case 0: goto accept_dig
;
3298 case 2: goto keep_dig
;
3300 #endif /*Honor_FLT_ROUNDS*/
3302 b
= lshift(PASS_STATE b
, 1);
3304 if ((j1
> 0 || (j1
== 0 && dig
& 1))
3313 #ifdef Honor_FLT_ROUNDS
3317 if (dig
== '9') { /* possible if i == 1 */
3325 #ifdef Honor_FLT_ROUNDS
3331 b
= multadd(PASS_STATE b
, 10, 0);
3333 mlo
= mhi
= multadd(PASS_STATE mhi
, 10, 0);
3335 mlo
= multadd(PASS_STATE mlo
, 10, 0);
3336 mhi
= multadd(PASS_STATE mhi
, 10, 0);
3342 *s
++ = dig
= quorem(b
,S
) + '0';
3343 if (!b
->x
[0] && b
->wds
<= 1) {
3351 b
= multadd(PASS_STATE b
, 10, 0);
3354 /* Round off last digit */
3356 #ifdef Honor_FLT_ROUNDS
3358 case 0: goto trimzeros
;
3359 case 2: goto roundoff
;
3362 b
= lshift(PASS_STATE b
, 1);
3364 if (j
>= 0) { /* ECMA compatible rounding needed by Spidermonkey */
3375 #ifdef Honor_FLT_ROUNDS
3382 Bfree(PASS_STATE S
);
3384 if (mlo
&& mlo
!= mhi
)
3385 Bfree(PASS_STATE mlo
);
3386 Bfree(PASS_STATE mhi
);
3392 word0(d
) = Exp_1
+ (70 << Exp_shift
);
3397 else if (!oldinexact
)
3400 Bfree(PASS_STATE b
);