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