2 |
- |
1 |
/****************************************************************
|
|
|
2 |
|
|
|
3 |
The author of this software (_dtoa, strtod) is David M. Gay.
|
|
|
4 |
Please send bug reports to
|
|
|
5 |
David M. Gay
|
|
|
6 |
Bell Laboratories, Room 2C-463
|
|
|
7 |
600 Mountain Avenue
|
|
|
8 |
Murray Hill, NJ 07974-2070
|
|
|
9 |
U.S.A.
|
|
|
10 |
dmg@research.bell-labs.com
|
|
|
11 |
*/
|
|
|
12 |
#include <stdlib.h>
|
|
|
13 |
#include <string.h>
|
|
|
14 |
#define _RESEARCH_SOURCE
|
|
|
15 |
#include <errno.h>
|
|
|
16 |
#include <float.h>
|
|
|
17 |
#include <math.h>
|
|
|
18 |
#define CONST const
|
|
|
19 |
|
|
|
20 |
/*
|
|
|
21 |
* #define IEEE_8087 for IEEE-arithmetic machines where the least
|
|
|
22 |
* significant byte has the lowest address.
|
|
|
23 |
* #define IEEE_MC68k for IEEE-arithmetic machines where the most
|
|
|
24 |
* significant byte has the lowest address.
|
|
|
25 |
* #define Sudden_Underflow for IEEE-format machines without gradual
|
|
|
26 |
* underflow (i.e., that flush to zero on underflow).
|
|
|
27 |
* #define IBM for IBM mainframe-style floating-point arithmetic.
|
|
|
28 |
* #define VAX for VAX-style floating-point arithmetic.
|
|
|
29 |
* #define Unsigned_Shifts if >> does treats its left operand as unsigned.
|
|
|
30 |
* #define No_leftright to omit left-right logic in fast floating-point
|
|
|
31 |
* computation of dtoa.
|
|
|
32 |
* #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
|
|
|
33 |
* #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
|
|
|
34 |
* that use extended-precision instructions to compute rounded
|
|
|
35 |
* products and quotients) with IBM.
|
|
|
36 |
* #define ROUND_BIASED for IEEE-format with biased rounding.
|
|
|
37 |
* #define Inaccurate_Divide for IEEE-format with correctly rounded
|
|
|
38 |
* products but inaccurate quotients, e.g., for Intel i860.
|
|
|
39 |
* #define Just_16 to store 16 bits per 32-bit long when doing high-precision
|
|
|
40 |
* integer arithmetic. Whether this speeds things up or slows things
|
|
|
41 |
* down depends on the machine and the number being converted.
|
|
|
42 |
*/
|
|
|
43 |
|
|
|
44 |
#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
|
|
|
45 |
Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
|
|
|
46 |
#endif
|
|
|
47 |
|
|
|
48 |
typedef union {
|
|
|
49 |
double d;
|
|
|
50 |
unsigned long ul[2];
|
|
|
51 |
} Dul;
|
|
|
52 |
|
|
|
53 |
#ifdef IEEE_8087
|
|
|
54 |
#define word0(x) ((x).ul[1])
|
|
|
55 |
#define word1(x) ((x).ul[0])
|
|
|
56 |
#else
|
|
|
57 |
#define word0(x) ((x).ul[0])
|
|
|
58 |
#define word1(x) ((x).ul[1])
|
|
|
59 |
#endif
|
|
|
60 |
|
|
|
61 |
#ifdef Unsigned_Shifts
|
|
|
62 |
#define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
|
|
|
63 |
#else
|
|
|
64 |
#define Sign_Extend(a,b) /*no-op*/
|
|
|
65 |
#endif
|
|
|
66 |
|
|
|
67 |
/* The following definition of Storeinc is appropriate for MIPS processors.
|
|
|
68 |
* An alternative that might be better on some machines is
|
|
|
69 |
* #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
|
|
|
70 |
*/
|
|
|
71 |
#if defined(IEEE_8087) + defined(VAX)
|
|
|
72 |
#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
|
|
|
73 |
((unsigned short *)a)[0] = (unsigned short)c, a++)
|
|
|
74 |
#else
|
|
|
75 |
#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
|
|
|
76 |
((unsigned short *)a)[1] = (unsigned short)c, a++)
|
|
|
77 |
#endif
|
|
|
78 |
|
|
|
79 |
/* #define P DBL_MANT_DIG */
|
|
|
80 |
/* Ten_pmax = floor(P*log(2)/log(5)) */
|
|
|
81 |
/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
|
|
|
82 |
/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
|
|
|
83 |
/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
|
|
|
84 |
|
|
|
85 |
#if defined(IEEE_8087) + defined(IEEE_MC68k)
|
|
|
86 |
#define Exp_shift 20
|
|
|
87 |
#define Exp_shift1 20
|
|
|
88 |
#define Exp_msk1 0x100000
|
|
|
89 |
#define Exp_msk11 0x100000
|
|
|
90 |
#define Exp_mask 0x7ff00000
|
|
|
91 |
#define P 53
|
|
|
92 |
#define Bias 1023
|
|
|
93 |
#define IEEE_Arith
|
|
|
94 |
#define Emin (-1022)
|
|
|
95 |
#define Exp_1 0x3ff00000
|
|
|
96 |
#define Exp_11 0x3ff00000
|
|
|
97 |
#define Ebits 11
|
|
|
98 |
#define Frac_mask 0xfffff
|
|
|
99 |
#define Frac_mask1 0xfffff
|
|
|
100 |
#define Ten_pmax 22
|
|
|
101 |
#define Bletch 0x10
|
|
|
102 |
#define Bndry_mask 0xfffff
|
|
|
103 |
#define Bndry_mask1 0xfffff
|
|
|
104 |
#define LSB 1
|
|
|
105 |
#define Sign_bit 0x80000000
|
|
|
106 |
#define Log2P 1
|
|
|
107 |
#define Tiny0 0
|
|
|
108 |
#define Tiny1 1
|
|
|
109 |
#define Quick_max 14
|
|
|
110 |
#define Int_max 14
|
|
|
111 |
#define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
|
|
|
112 |
#else
|
|
|
113 |
#undef Sudden_Underflow
|
|
|
114 |
#define Sudden_Underflow
|
|
|
115 |
#ifdef IBM
|
|
|
116 |
#define Exp_shift 24
|
|
|
117 |
#define Exp_shift1 24
|
|
|
118 |
#define Exp_msk1 0x1000000
|
|
|
119 |
#define Exp_msk11 0x1000000
|
|
|
120 |
#define Exp_mask 0x7f000000
|
|
|
121 |
#define P 14
|
|
|
122 |
#define Bias 65
|
|
|
123 |
#define Exp_1 0x41000000
|
|
|
124 |
#define Exp_11 0x41000000
|
|
|
125 |
#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
|
|
|
126 |
#define Frac_mask 0xffffff
|
|
|
127 |
#define Frac_mask1 0xffffff
|
|
|
128 |
#define Bletch 4
|
|
|
129 |
#define Ten_pmax 22
|
|
|
130 |
#define Bndry_mask 0xefffff
|
|
|
131 |
#define Bndry_mask1 0xffffff
|
|
|
132 |
#define LSB 1
|
|
|
133 |
#define Sign_bit 0x80000000
|
|
|
134 |
#define Log2P 4
|
|
|
135 |
#define Tiny0 0x100000
|
|
|
136 |
#define Tiny1 0
|
|
|
137 |
#define Quick_max 14
|
|
|
138 |
#define Int_max 15
|
|
|
139 |
#else /* VAX */
|
|
|
140 |
#define Exp_shift 23
|
|
|
141 |
#define Exp_shift1 7
|
|
|
142 |
#define Exp_msk1 0x80
|
|
|
143 |
#define Exp_msk11 0x800000
|
|
|
144 |
#define Exp_mask 0x7f80
|
|
|
145 |
#define P 56
|
|
|
146 |
#define Bias 129
|
|
|
147 |
#define Exp_1 0x40800000
|
|
|
148 |
#define Exp_11 0x4080
|
|
|
149 |
#define Ebits 8
|
|
|
150 |
#define Frac_mask 0x7fffff
|
|
|
151 |
#define Frac_mask1 0xffff007f
|
|
|
152 |
#define Ten_pmax 24
|
|
|
153 |
#define Bletch 2
|
|
|
154 |
#define Bndry_mask 0xffff007f
|
|
|
155 |
#define Bndry_mask1 0xffff007f
|
|
|
156 |
#define LSB 0x10000
|
|
|
157 |
#define Sign_bit 0x8000
|
|
|
158 |
#define Log2P 1
|
|
|
159 |
#define Tiny0 0x80
|
|
|
160 |
#define Tiny1 0
|
|
|
161 |
#define Quick_max 15
|
|
|
162 |
#define Int_max 15
|
|
|
163 |
#endif
|
|
|
164 |
#endif
|
|
|
165 |
|
|
|
166 |
#ifndef IEEE_Arith
|
|
|
167 |
#define ROUND_BIASED
|
|
|
168 |
#endif
|
|
|
169 |
|
|
|
170 |
#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
|
|
|
171 |
#define Big1 0xffffffff
|
|
|
172 |
|
|
|
173 |
#ifndef Just_16
|
|
|
174 |
/* When Pack_32 is not defined, we store 16 bits per 32-bit long.
|
|
|
175 |
* This makes some inner loops simpler and sometimes saves work
|
|
|
176 |
* during multiplications, but it often seems to make things slightly
|
|
|
177 |
* slower. Hence the default is now to store 32 bits per long.
|
|
|
178 |
*/
|
|
|
179 |
#ifndef Pack_32
|
|
|
180 |
#define Pack_32
|
|
|
181 |
#endif
|
|
|
182 |
#endif
|
|
|
183 |
|
|
|
184 |
#define Kmax 15
|
|
|
185 |
|
|
|
186 |
struct
|
|
|
187 |
Bigint {
|
|
|
188 |
struct Bigint *next;
|
|
|
189 |
int k, maxwds, sign, wds;
|
|
|
190 |
unsigned long x[1];
|
|
|
191 |
};
|
|
|
192 |
|
|
|
193 |
typedef struct Bigint Bigint;
|
|
|
194 |
|
|
|
195 |
/* This routines shouldn't be visible externally */
|
|
|
196 |
extern Bigint *_Balloc(int);
|
|
|
197 |
extern void _Bfree(Bigint *);
|
|
|
198 |
extern Bigint *_multadd(Bigint *, int, int);
|
|
|
199 |
extern int _hi0bits(unsigned long);
|
|
|
200 |
extern Bigint *_mult(Bigint *, Bigint *);
|
|
|
201 |
extern Bigint *_pow5mult(Bigint *, int);
|
|
|
202 |
extern Bigint *_lshift(Bigint *, int);
|
|
|
203 |
extern int _cmp(Bigint *, Bigint *);
|
|
|
204 |
extern Bigint *_diff(Bigint *, Bigint *);
|
|
|
205 |
extern Bigint *_d2b(double, int *, int *);
|
|
|
206 |
extern Bigint *_i2b(int);
|
|
|
207 |
|
|
|
208 |
extern double _tens[], _bigtens[], _tinytens[];
|
|
|
209 |
|
|
|
210 |
#ifdef IEEE_Arith
|
|
|
211 |
#define n_bigtens 5
|
|
|
212 |
#else
|
|
|
213 |
#ifdef IBM
|
|
|
214 |
#define n_bigtens 3
|
|
|
215 |
#else
|
|
|
216 |
#define n_bigtens 2
|
|
|
217 |
#endif
|
|
|
218 |
#endif
|
|
|
219 |
|
|
|
220 |
#define Balloc(x) _Balloc(x)
|
|
|
221 |
#define Bfree(x) _Bfree(x)
|
|
|
222 |
#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
|
|
|
223 |
y->wds*sizeof(long) + 2*sizeof(int))
|
|
|
224 |
#define multadd(x,y,z) _multadd(x,y,z)
|
|
|
225 |
#define hi0bits(x) _hi0bits(x)
|
|
|
226 |
#define i2b(x) _i2b(x)
|
|
|
227 |
#define mult(x,y) _mult(x,y)
|
|
|
228 |
#define pow5mult(x,y) _pow5mult(x,y)
|
|
|
229 |
#define lshift(x,y) _lshift(x,y)
|
|
|
230 |
#define cmp(x,y) _cmp(x,y)
|
|
|
231 |
#define diff(x,y) _diff(x,y)
|
|
|
232 |
#define d2b(x,y,z) _d2b(x,y,z)
|
|
|
233 |
|
|
|
234 |
#define tens _tens
|
|
|
235 |
#define bigtens _bigtens
|
|
|
236 |
#define tinytens _tinytens
|