Rev 4921 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
4349 | Serge | 1 | /**************************************************************** |
2 | * |
||
3 | * The author of this software is David M. Gay. |
||
4 | * |
||
5 | * Copyright (c) 1991 by AT&T. |
||
6 | * |
||
7 | * Permission to use, copy, modify, and distribute this software for any |
||
8 | * purpose without fee is hereby granted, provided that this entire notice |
||
9 | * is included in all copies of any software which is or includes a copy |
||
10 | * or modification of this software and in all copies of the supporting |
||
11 | * documentation for such software. |
||
12 | * |
||
13 | * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED |
||
14 | * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY |
||
15 | * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY |
||
16 | * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. |
||
17 | * |
||
18 | ***************************************************************/ |
||
19 | |||
20 | /* Please send bug reports to |
||
21 | David M. Gay |
||
22 | AT&T Bell Laboratories, Room 2C-463 |
||
23 | 600 Mountain Avenue |
||
24 | Murray Hill, NJ 07974-2070 |
||
25 | U.S.A. |
||
26 | dmg@research.att.com or research!dmg |
||
27 | */ |
||
28 | |||
29 | #include |
||
30 | #include |
||
31 | #include |
||
32 | #include |
||
33 | #include |
||
34 | #include |
||
35 | |||
36 | #ifdef __IEEE_LITTLE_ENDIAN |
||
37 | #define IEEE_8087 |
||
38 | #endif |
||
39 | |||
40 | #ifdef __IEEE_BIG_ENDIAN |
||
41 | #define IEEE_MC68k |
||
42 | #endif |
||
43 | |||
44 | #ifdef __Z8000__ |
||
45 | #define Just_16 |
||
46 | #endif |
||
47 | |||
48 | #ifdef DEBUG |
||
49 | #include "stdio.h" |
||
50 | #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} |
||
51 | #endif |
||
52 | |||
53 | #ifdef Unsigned_Shifts |
||
54 | #define Sign_Extend(a,b) if (b < 0) a |= (__uint32_t)0xffff0000; |
||
55 | #else |
||
56 | #define Sign_Extend(a,b) /*no-op*/ |
||
57 | #endif |
||
58 | |||
59 | #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 |
||
60 | Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. |
||
61 | #endif |
||
62 | |||
63 | /* If we are going to examine or modify specific bits in a double using |
||
64 | the word0 and/or word1 macros, then we must wrap the double inside |
||
65 | a union. This is necessary to avoid undefined behavior according to |
||
66 | the ANSI C spec. */ |
||
67 | union double_union |
||
68 | { |
||
69 | double d; |
||
70 | __uint32_t i[2]; |
||
71 | }; |
||
72 | |||
73 | #ifdef IEEE_8087 |
||
74 | #define word0(x) (x.i[1]) |
||
75 | #define word1(x) (x.i[0]) |
||
76 | #else |
||
77 | #define word0(x) (x.i[0]) |
||
78 | #define word1(x) (x.i[1]) |
||
79 | #endif |
||
80 | |||
81 | /* The following is taken from gdtoaimp.h for use with new strtod, but |
||
82 | adjusted to avoid invalid type-punning. */ |
||
83 | typedef __int32_t Long; |
||
84 | |||
85 | /* Unfortunately, because __ULong might be a different type than |
||
86 | __uint32_t, we can't re-use union double_union as-is without |
||
87 | further edits in strtod.c. */ |
||
88 | typedef union { double d; __ULong i[2]; } U; |
||
89 | |||
90 | #define dword0(x) word0(x) |
||
91 | #define dword1(x) word1(x) |
||
92 | #define dval(x) (x.d) |
||
93 | |||
94 | #undef SI |
||
95 | #ifdef Sudden_Underflow |
||
96 | #define SI 1 |
||
97 | #else |
||
98 | #define SI 0 |
||
99 | #endif |
||
100 | |||
4921 | Serge | 101 | #define Storeinc(a,b,c) (*(a)++ = ((b) << 16) | ((c) & 0xffff)) |
4349 | Serge | 102 | |
103 | /* #define P DBL_MANT_DIG */ |
||
104 | /* Ten_pmax = floor(P*log(2)/log(5)) */ |
||
105 | /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ |
||
106 | /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ |
||
107 | /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ |
||
108 | |||
109 | #if defined(IEEE_8087) + defined(IEEE_MC68k) |
||
110 | #if defined (_DOUBLE_IS_32BITS) |
||
111 | #define Exp_shift 23 |
||
112 | #define Exp_shift1 23 |
||
113 | #define Exp_msk1 ((__uint32_t)0x00800000L) |
||
114 | #define Exp_msk11 ((__uint32_t)0x00800000L) |
||
115 | #define Exp_mask ((__uint32_t)0x7f800000L) |
||
116 | #define P 24 |
||
117 | #define Bias 127 |
||
118 | #define NO_HEX_FP /* not supported in this case */ |
||
119 | #define IEEE_Arith |
||
120 | #define Emin (-126) |
||
121 | #define Exp_1 ((__uint32_t)0x3f800000L) |
||
122 | #define Exp_11 ((__uint32_t)0x3f800000L) |
||
123 | #define Ebits 8 |
||
124 | #define Frac_mask ((__uint32_t)0x007fffffL) |
||
125 | #define Frac_mask1 ((__uint32_t)0x007fffffL) |
||
126 | #define Ten_pmax 10 |
||
127 | #define Sign_bit ((__uint32_t)0x80000000L) |
||
128 | #define Ten_pmax 10 |
||
129 | #define Bletch 2 |
||
130 | #define Bndry_mask ((__uint32_t)0x007fffffL) |
||
131 | #define Bndry_mask1 ((__uint32_t)0x007fffffL) |
||
132 | #define LSB 1 |
||
133 | #define Sign_bit ((__uint32_t)0x80000000L) |
||
134 | #define Log2P 1 |
||
135 | #define Tiny0 0 |
||
136 | #define Tiny1 1 |
||
137 | #define Quick_max 5 |
||
138 | #define Int_max 6 |
||
139 | #define Infinite(x) (word0(x) == ((__uint32_t)0x7f800000L)) |
||
140 | #undef word0 |
||
141 | #undef word1 |
||
142 | #undef dword0 |
||
143 | #undef dword1 |
||
144 | |||
145 | #define word0(x) (x.i[0]) |
||
146 | #define word1(x) 0 |
||
147 | #define dword0(x) word0(x) |
||
148 | #define dword1(x) 0 |
||
149 | #else |
||
150 | |||
151 | #define Exp_shift 20 |
||
152 | #define Exp_shift1 20 |
||
153 | #define Exp_msk1 ((__uint32_t)0x100000L) |
||
154 | #define Exp_msk11 ((__uint32_t)0x100000L) |
||
155 | #define Exp_mask ((__uint32_t)0x7ff00000L) |
||
156 | #define P 53 |
||
157 | #define Bias 1023 |
||
158 | #define IEEE_Arith |
||
159 | #define Emin (-1022) |
||
160 | #define Exp_1 ((__uint32_t)0x3ff00000L) |
||
161 | #define Exp_11 ((__uint32_t)0x3ff00000L) |
||
162 | #define Ebits 11 |
||
163 | #define Frac_mask ((__uint32_t)0xfffffL) |
||
164 | #define Frac_mask1 ((__uint32_t)0xfffffL) |
||
165 | #define Ten_pmax 22 |
||
166 | #define Bletch 0x10 |
||
167 | #define Bndry_mask ((__uint32_t)0xfffffL) |
||
168 | #define Bndry_mask1 ((__uint32_t)0xfffffL) |
||
169 | #define LSB 1 |
||
170 | #define Sign_bit ((__uint32_t)0x80000000L) |
||
171 | #define Log2P 1 |
||
172 | #define Tiny0 0 |
||
173 | #define Tiny1 1 |
||
174 | #define Quick_max 14 |
||
175 | #define Int_max 14 |
||
176 | #define Infinite(x) (word0(x) == ((__uint32_t)0x7ff00000L)) /* sufficient test for here */ |
||
177 | |||
178 | #endif /* !_DOUBLE_IS_32BITS */ |
||
179 | |||
180 | #ifndef Flt_Rounds |
||
181 | #ifdef FLT_ROUNDS |
||
182 | #define Flt_Rounds FLT_ROUNDS |
||
183 | #else |
||
184 | #define Flt_Rounds 1 |
||
185 | #endif |
||
186 | #endif /*Flt_Rounds*/ |
||
187 | |||
188 | #else /* !IEEE_8087 && !IEEE_MC68k */ |
||
189 | #undef Sudden_Underflow |
||
190 | #define Sudden_Underflow |
||
191 | #ifdef IBM |
||
192 | #define Flt_Rounds 0 |
||
193 | #define Exp_shift 24 |
||
194 | #define Exp_shift1 24 |
||
195 | #define Exp_msk1 ((__uint32_t)0x1000000L) |
||
196 | #define Exp_msk11 ((__uint32_t)0x1000000L) |
||
197 | #define Exp_mask ((__uint32_t)0x7f000000L) |
||
198 | #define P 14 |
||
199 | #define Bias 65 |
||
200 | #define Exp_1 ((__uint32_t)0x41000000L) |
||
201 | #define Exp_11 ((__uint32_t)0x41000000L) |
||
202 | #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ |
||
203 | #define Frac_mask ((__uint32_t)0xffffffL) |
||
204 | #define Frac_mask1 ((__uint32_t)0xffffffL) |
||
205 | #define Bletch 4 |
||
206 | #define Ten_pmax 22 |
||
207 | #define Bndry_mask ((__uint32_t)0xefffffL) |
||
208 | #define Bndry_mask1 ((__uint32_t)0xffffffL) |
||
209 | #define LSB 1 |
||
210 | #define Sign_bit ((__uint32_t)0x80000000L) |
||
211 | #define Log2P 4 |
||
212 | #define Tiny0 ((__uint32_t)0x100000L) |
||
213 | #define Tiny1 0 |
||
214 | #define Quick_max 14 |
||
215 | #define Int_max 15 |
||
216 | #else /* VAX */ |
||
217 | #define Flt_Rounds 1 |
||
218 | #define Exp_shift 23 |
||
219 | #define Exp_shift1 7 |
||
220 | #define Exp_msk1 0x80 |
||
221 | #define Exp_msk11 ((__uint32_t)0x800000L) |
||
222 | #define Exp_mask ((__uint32_t)0x7f80L) |
||
223 | #define P 56 |
||
224 | #define Bias 129 |
||
225 | #define Exp_1 ((__uint32_t)0x40800000L) |
||
226 | #define Exp_11 ((__uint32_t)0x4080L) |
||
227 | #define Ebits 8 |
||
228 | #define Frac_mask ((__uint32_t)0x7fffffL) |
||
229 | #define Frac_mask1 ((__uint32_t)0xffff007fL) |
||
230 | #define Ten_pmax 24 |
||
231 | #define Bletch 2 |
||
232 | #define Bndry_mask ((__uint32_t)0xffff007fL) |
||
233 | #define Bndry_mask1 ((__uint32_t)0xffff007fL) |
||
234 | #define LSB ((__uint32_t)0x10000L) |
||
235 | #define Sign_bit ((__uint32_t)0x8000L) |
||
236 | #define Log2P 1 |
||
237 | #define Tiny0 0x80 |
||
238 | #define Tiny1 0 |
||
239 | #define Quick_max 15 |
||
240 | #define Int_max 15 |
||
241 | #endif |
||
242 | #endif |
||
243 | |||
244 | #ifndef IEEE_Arith |
||
245 | #define ROUND_BIASED |
||
246 | #else |
||
247 | #define Scale_Bit 0x10 |
||
248 | #if defined(_DOUBLE_IS_32BITS) && defined(__v800) |
||
249 | #define n_bigtens 2 |
||
250 | #else |
||
251 | #define n_bigtens 5 |
||
252 | #endif |
||
253 | #endif |
||
254 | |||
255 | #ifdef IBM |
||
256 | #define n_bigtens 3 |
||
257 | #endif |
||
258 | |||
259 | #ifdef VAX |
||
260 | #define n_bigtens 2 |
||
261 | #endif |
||
262 | |||
263 | #ifndef __NO_INFNAN_CHECK |
||
264 | #define INFNAN_CHECK |
||
265 | #endif |
||
266 | |||
267 | /* |
||
268 | * NAN_WORD0 and NAN_WORD1 are only referenced in strtod.c. Prior to |
||
269 | * 20050115, they used to be hard-wired here (to 0x7ff80000 and 0, |
||
270 | * respectively), but now are determined by compiling and running |
||
271 | * qnan.c to generate gd_qnan.h, which specifies d_QNAN0 and d_QNAN1. |
||
272 | * Formerly gdtoaimp.h recommended supplying suitable -DNAN_WORD0=... |
||
273 | * and -DNAN_WORD1=... values if necessary. This should still work. |
||
274 | * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.) |
||
275 | */ |
||
276 | #ifdef IEEE_Arith |
||
277 | #ifdef IEEE_MC68k |
||
278 | #define _0 0 |
||
279 | #define _1 1 |
||
280 | #ifndef NAN_WORD0 |
||
281 | #define NAN_WORD0 d_QNAN0 |
||
282 | #endif |
||
283 | #ifndef NAN_WORD1 |
||
284 | #define NAN_WORD1 d_QNAN1 |
||
285 | #endif |
||
286 | #else |
||
287 | #define _0 1 |
||
288 | #define _1 0 |
||
289 | #ifndef NAN_WORD0 |
||
290 | #define NAN_WORD0 d_QNAN1 |
||
291 | #endif |
||
292 | #ifndef NAN_WORD1 |
||
293 | #define NAN_WORD1 d_QNAN0 |
||
294 | #endif |
||
295 | #endif |
||
296 | #else |
||
297 | #undef INFNAN_CHECK |
||
298 | #endif |
||
299 | |||
300 | #ifdef RND_PRODQUOT |
||
301 | #define rounded_product(a,b) a = rnd_prod(a, b) |
||
302 | #define rounded_quotient(a,b) a = rnd_quot(a, b) |
||
303 | #ifdef KR_headers |
||
304 | extern double rnd_prod(), rnd_quot(); |
||
305 | #else |
||
306 | extern double rnd_prod(double, double), rnd_quot(double, double); |
||
307 | #endif |
||
308 | #else |
||
309 | #define rounded_product(a,b) a *= b |
||
310 | #define rounded_quotient(a,b) a /= b |
||
311 | #endif |
||
312 | |||
313 | #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) |
||
314 | #define Big1 ((__uint32_t)0xffffffffL) |
||
315 | |||
316 | #ifndef Just_16 |
||
317 | /* When Pack_32 is not defined, we store 16 bits per 32-bit long. |
||
318 | * This makes some inner loops simpler and sometimes saves work |
||
319 | * during multiplications, but it often seems to make things slightly |
||
320 | * slower. Hence the default is now to store 32 bits per long. |
||
321 | */ |
||
322 | |||
323 | #ifndef Pack_32 |
||
324 | #define Pack_32 |
||
325 | #endif |
||
326 | #else /* Just_16 */ |
||
327 | #ifndef Pack_16 |
||
328 | #define Pack_16 |
||
329 | #endif |
||
330 | #endif /* Just_16 */ |
||
331 | |||
332 | #ifdef Pack_32 |
||
333 | #define ULbits 32 |
||
334 | #define kshift 5 |
||
335 | #define kmask 31 |
||
336 | #define ALL_ON 0xffffffff |
||
337 | #else |
||
338 | #define ULbits 16 |
||
339 | #define kshift 4 |
||
340 | #define kmask 15 |
||
341 | #define ALL_ON 0xffff |
||
342 | #endif |
||
343 | |||
344 | #ifdef __cplusplus |
||
345 | extern "C" double strtod(const char *s00, char **se); |
||
346 | extern "C" char *dtoa(double d, int mode, int ndigits, |
||
347 | int *decpt, int *sign, char **rve); |
||
348 | #endif |
||
349 | |||
350 | |||
351 | typedef struct _Bigint _Bigint; |
||
352 | |||
353 | #define Balloc _Balloc |
||
354 | #define Bfree _Bfree |
||
355 | #define multadd __multadd |
||
356 | #define s2b __s2b |
||
357 | #define lo0bits __lo0bits |
||
358 | #define hi0bits __hi0bits |
||
359 | #define i2b __i2b |
||
360 | #define mult __multiply |
||
361 | #define pow5mult __pow5mult |
||
362 | #define lshift __lshift |
||
6099 | serge | 363 | #define match __match |
4349 | Serge | 364 | #define cmp __mcmp |
365 | #define diff __mdiff |
||
366 | #define ulp __ulp |
||
367 | #define b2d __b2d |
||
368 | #define d2b __d2b |
||
369 | #define ratio __ratio |
||
370 | #define any_on __any_on |
||
371 | #define gethex __gethex |
||
372 | #define copybits __copybits |
||
373 | #define hexnan __hexnan |
||
374 | |||
4921 | Serge | 375 | #if !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) |
376 | #define __get_hexdig(x) __hexdig[x] /* NOTE: must evaluate arg only once */ |
||
377 | #else /* !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) */ |
||
378 | #define __get_hexdig(x) __hexdig_fun(x) |
||
379 | #endif /* !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) */ |
||
4349 | Serge | 380 | |
381 | #define tens __mprec_tens |
||
382 | #define bigtens __mprec_bigtens |
||
383 | #define tinytens __mprec_tinytens |
||
384 | |||
385 | struct _reent ; |
||
386 | struct FPI; |
||
387 | double _EXFUN(ulp,(double x)); |
||
388 | double _EXFUN(b2d,(_Bigint *a , int *e)); |
||
389 | _Bigint * _EXFUN(Balloc,(struct _reent *p, int k)); |
||
390 | void _EXFUN(Bfree,(struct _reent *p, _Bigint *v)); |
||
391 | _Bigint * _EXFUN(multadd,(struct _reent *p, _Bigint *, int, int)); |
||
392 | _Bigint * _EXFUN(s2b,(struct _reent *, const char*, int, int, __ULong)); |
||
393 | _Bigint * _EXFUN(i2b,(struct _reent *,int)); |
||
394 | _Bigint * _EXFUN(mult, (struct _reent *, _Bigint *, _Bigint *)); |
||
395 | _Bigint * _EXFUN(pow5mult, (struct _reent *, _Bigint *, int k)); |
||
396 | int _EXFUN(hi0bits,(__ULong)); |
||
397 | int _EXFUN(lo0bits,(__ULong *)); |
||
398 | _Bigint * _EXFUN(d2b,(struct _reent *p, double d, int *e, int *bits)); |
||
399 | _Bigint * _EXFUN(lshift,(struct _reent *p, _Bigint *b, int k)); |
||
6099 | serge | 400 | int _EXFUN(match,(const char**, char*)); |
4349 | Serge | 401 | _Bigint * _EXFUN(diff,(struct _reent *p, _Bigint *a, _Bigint *b)); |
402 | int _EXFUN(cmp,(_Bigint *a, _Bigint *b)); |
||
4921 | Serge | 403 | int _EXFUN(gethex,(struct _reent *p, _CONST char **sp, _CONST struct FPI *fpi, Long *exp, _Bigint **bp, int sign)); |
4349 | Serge | 404 | double _EXFUN(ratio,(_Bigint *a, _Bigint *b)); |
405 | __ULong _EXFUN(any_on,(_Bigint *b, int k)); |
||
406 | void _EXFUN(copybits,(__ULong *c, int n, _Bigint *b)); |
||
6099 | serge | 407 | #if defined (_HAVE_LONG_DOUBLE) && !defined (_LDBL_EQ_DBL) |
408 | int _EXFUN(_strtorx_r,(struct _reent *, _CONST char *, char **, int, void *)); |
||
409 | int _EXFUN(_strtodg_r,(struct _reent *p, _CONST char *s00, char **se, struct FPI *fpi, Long *exp, __ULong *bits)); |
||
410 | #endif /* _HAVE_LONG_DOUBLE && !_LDBL_EQ_DBL */ |
||
411 | |||
4921 | Serge | 412 | #if defined(PREFER_SIZE_OVER_SPEED) || defined(__OPTIMIZE_SIZE__) || defined(_SMALL_HEXDIG) |
413 | unsigned char _EXFUN(__hexdig_fun,(unsigned char)); |
||
414 | #endif /* !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) */ |
||
4349 | Serge | 415 | #ifdef INFNAN_CHECK |
4921 | Serge | 416 | int _EXFUN(hexnan,(_CONST char **sp, _CONST struct FPI *fpi, __ULong *x0)); |
4349 | Serge | 417 | #endif |
418 | |||
419 | #define Bcopy(x,y) memcpy((char *)&x->_sign, (char *)&y->_sign, y->_wds*sizeof(__Long) + 2*sizeof(int)) |
||
420 | |||
421 | extern _CONST double tinytens[]; |
||
422 | extern _CONST double bigtens[]; |
||
423 | extern _CONST double tens[]; |
||
4921 | Serge | 424 | #if !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) |
425 | extern _CONST unsigned char __hexdig[]; |
||
426 | #endif /* !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) */ |
||
4349 | Serge | 427 | |
428 | |||
429 | double _EXFUN(_mprec_log10,(int));>><>> |