OLD | NEW |
1 /**************************************************************** | 1 /**************************************************************** |
2 * | 2 * |
3 * The author of this software is David M. Gay. | 3 * The author of this software is David M. Gay. |
4 * | 4 * |
5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. | 5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. |
6 * | 6 * |
7 * Permission to use, copy, modify, and distribute this software for any | 7 * Permission to use, copy, modify, and distribute this software for any |
8 * purpose without fee is hereby granted, provided that this entire notice | 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 | 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 | 10 * or modification of this software and in all copies of the supporting |
(...skipping 85 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
96 * something other than "long long", #define Llong to be the name, | 96 * something other than "long long", #define Llong to be the name, |
97 * and if "unsigned Llong" does not work as an unsigned version of | 97 * and if "unsigned Llong" does not work as an unsigned version of |
98 * Llong, #define #ULLong to be the corresponding unsigned type. | 98 * Llong, #define #ULLong to be the corresponding unsigned type. |
99 * #define KR_headers for old-style C function headers. | 99 * #define KR_headers for old-style C function headers. |
100 * #define Bad_float_h if your system lacks a float.h or if it does not | 100 * #define Bad_float_h if your system lacks a float.h or if it does not |
101 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, | 101 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, |
102 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. | 102 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. |
103 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) | 103 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) |
104 * if memory is available and otherwise does something you deem | 104 * if memory is available and otherwise does something you deem |
105 * appropriate. If MALLOC is undefined, malloc will be invoked | 105 * appropriate. If MALLOC is undefined, malloc will be invoked |
106 *» directly -- and assumed always to succeed. | 106 *» directly -- and assumed always to succeed. Similarly, if you |
| 107 *» want something other than the system's free() to be called to |
| 108 *» recycle memory acquired from MALLOC, #define FREE to be the |
| 109 *» name of the alternate routine. (FREE or free is only called in |
| 110 *» pathological cases, e.g., in a dtoa call after a dtoa return in |
| 111 *» mode 3 with thousands of digits requested.) |
107 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making | 112 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making |
108 * memory allocations from a private pool of memory when possible. | 113 * memory allocations from a private pool of memory when possible. |
109 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes, | 114 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes, |
110 * unless #defined to be a different length. This default length | 115 * unless #defined to be a different length. This default length |
111 * suffices to get rid of MALLOC calls except for unusual cases, | 116 * suffices to get rid of MALLOC calls except for unusual cases, |
112 * such as decimal-to-binary conversion of a very long string of | 117 * such as decimal-to-binary conversion of a very long string of |
113 * digits. The longest string dtoa can return is about 751 bytes | 118 * digits. The longest string dtoa can return is about 751 bytes |
114 * long. For conversions by strtod of strings of 800 digits and | 119 * long. For conversions by strtod of strings of 800 digits and |
115 * all dtoa conversions in single-threaded executions with 8-byte | 120 * all dtoa conversions in single-threaded executions with 8-byte |
116 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte | 121 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte |
(...skipping 23 matching lines...) Expand all Loading... |
140 * powers of 5; omitting this lock would introduce a small | 145 * powers of 5; omitting this lock would introduce a small |
141 * probability of wasting memory, but would otherwise be harmless.) | 146 * probability of wasting memory, but would otherwise be harmless.) |
142 * You must also invoke freedtoa(s) to free the value s returned by | 147 * You must also invoke freedtoa(s) to free the value s returned by |
143 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. | 148 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. |
144 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that | 149 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that |
145 * avoids underflows on inputs whose result does not underflow. | 150 * avoids underflows on inputs whose result does not underflow. |
146 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format | 151 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format |
147 * floating-point numbers and flushes underflows to zero rather | 152 * floating-point numbers and flushes underflows to zero rather |
148 * than implementing gradual underflow, then you must also #define | 153 * than implementing gradual underflow, then you must also #define |
149 * Sudden_Underflow. | 154 * Sudden_Underflow. |
150 * #define YES_ALIAS to permit aliasing certain double values with | |
151 * arrays of ULongs. This leads to slightly better code with | |
152 * some compilers and was always used prior to 19990916, but it | |
153 * is not strictly legal and can cause trouble with aggressively | |
154 * optimizing compilers (e.g., gcc 2.95.1 under -O2). | |
155 * #define USE_LOCALE to use the current locale's decimal_point value. | 155 * #define USE_LOCALE to use the current locale's decimal_point value. |
156 * #define SET_INEXACT if IEEE arithmetic is being used and extra | 156 * #define SET_INEXACT if IEEE arithmetic is being used and extra |
157 * computation should be done to set the inexact flag when the | 157 * computation should be done to set the inexact flag when the |
158 * result is inexact and avoid setting inexact when the result | 158 * result is inexact and avoid setting inexact when the result |
159 * is exact. In this case, dtoa.c must be compiled in | 159 * is exact. In this case, dtoa.c must be compiled in |
160 * an environment, perhaps provided by #include "dtoa.c" in a | 160 * an environment, perhaps provided by #include "dtoa.c" in a |
161 * suitable wrapper, that defines two functions, | 161 * suitable wrapper, that defines two functions, |
162 * int get_inexact(void); | 162 * int get_inexact(void); |
163 * void clear_inexact(void); | 163 * void clear_inexact(void); |
164 * such that get_inexact() returns a nonzero value if the | 164 * such that get_inexact() returns a nonzero value if the |
165 * inexact bit is already set, and clear_inexact() sets the | 165 * inexact bit is already set, and clear_inexact() sets the |
166 * inexact bit to 0. When SET_INEXACT is #defined, strtod | 166 * inexact bit to 0. When SET_INEXACT is #defined, strtod |
167 * also does extra computations to set the underflow and overflow | 167 * also does extra computations to set the underflow and overflow |
168 * flags when appropriate (i.e., when the result is tiny and | 168 * flags when appropriate (i.e., when the result is tiny and |
169 * inexact or when it is a numeric value rounded to +-infinity). | 169 * inexact or when it is a numeric value rounded to +-infinity). |
170 * #define NO_ERRNO if strtod should not assign errno = ERANGE when | 170 * #define NO_ERRNO if strtod should not assign errno = ERANGE when |
171 * the result overflows to +-Infinity or underflows to 0. | 171 * the result overflows to +-Infinity or underflows to 0. |
| 172 * #define NO_HEX_FP to omit recognition of hexadecimal floating-point |
| 173 * values by strtod. |
| 174 * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now) |
| 175 * to disable logic for "fast" testing of very long input strings |
| 176 * to strtod. This testing proceeds by initially truncating the |
| 177 * input string, then if necessary comparing the whole string with |
| 178 * a decimal expansion to decide close cases. This logic is only |
| 179 * used for input more than STRTOD_DIGLIM digits long (default 40). |
172 */ | 180 */ |
173 | 181 |
174 #define IEEE_8087 | 182 #define IEEE_8087 |
175 #if defined(__GNUC__) | 183 #define NO_HEX_FP |
176 // Make gcc 4.3+ warnings about parentheses non-fatal warnings. | |
177 #pragma GCC diagnostic warning "-Wparentheses" | |
178 #endif | |
179 | 184 |
180 #ifndef Long | 185 #ifndef Long |
181 #define Long long | 186 #define Long long |
182 #endif | 187 #endif |
183 #ifndef ULong | 188 #ifndef ULong |
184 typedef unsigned Long ULong; | 189 typedef unsigned Long ULong; |
185 #endif | 190 #endif |
186 | 191 |
187 #ifdef DEBUG | 192 #ifdef DEBUG |
188 #include "stdio.h" | 193 #include "stdio.h" |
189 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} | 194 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} |
190 #endif | 195 #endif |
191 | 196 |
192 #include "stdlib.h" | 197 #include "stdlib.h" |
193 #include "string.h" | 198 #include "string.h" |
194 | 199 |
195 #ifdef USE_LOCALE | 200 #ifdef USE_LOCALE |
196 #include "locale.h" | 201 #include "locale.h" |
197 #endif | 202 #endif |
198 | 203 |
| 204 #ifdef Honor_FLT_ROUNDS |
| 205 #ifndef Trust_FLT_ROUNDS |
| 206 #include <fenv.h> |
| 207 #endif |
| 208 #endif |
| 209 |
199 #ifdef MALLOC | 210 #ifdef MALLOC |
200 #ifdef KR_headers | 211 #ifdef KR_headers |
201 extern char *MALLOC(); | 212 extern char *MALLOC(); |
202 #else | 213 #else |
203 extern void *MALLOC(size_t); | 214 extern void *MALLOC(size_t); |
204 #endif | 215 #endif |
205 #else | 216 #else |
206 #define MALLOC malloc | 217 #define MALLOC malloc |
207 #endif | 218 #endif |
208 | 219 |
(...skipping 14 matching lines...) Expand all Loading... |
223 #define IEEE_Arith | 234 #define IEEE_Arith |
224 #endif | 235 #endif |
225 | 236 |
226 #ifdef IEEE_Arith | 237 #ifdef IEEE_Arith |
227 #ifndef NO_INFNAN_CHECK | 238 #ifndef NO_INFNAN_CHECK |
228 #undef INFNAN_CHECK | 239 #undef INFNAN_CHECK |
229 #define INFNAN_CHECK | 240 #define INFNAN_CHECK |
230 #endif | 241 #endif |
231 #else | 242 #else |
232 #undef INFNAN_CHECK | 243 #undef INFNAN_CHECK |
| 244 #define NO_STRTOD_BIGCOMP |
233 #endif | 245 #endif |
234 | 246 |
235 #include "errno.h" | 247 #include "errno.h" |
236 | 248 |
237 #ifdef Bad_float_h | 249 #ifdef Bad_float_h |
238 | 250 |
239 #ifdef IEEE_Arith | 251 #ifdef IEEE_Arith |
240 #define DBL_DIG 15 | 252 #define DBL_DIG 15 |
241 #define DBL_MAX_10_EXP 308 | 253 #define DBL_MAX_10_EXP 308 |
242 #define DBL_MAX_EXP 1024 | 254 #define DBL_MAX_EXP 1024 |
(...skipping 37 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
280 #define CONST const | 292 #define CONST const |
281 #endif | 293 #endif |
282 #endif | 294 #endif |
283 | 295 |
284 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 | 296 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 |
285 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. | 297 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. |
286 #endif | 298 #endif |
287 | 299 |
288 typedef union { double d; ULong L[2]; } U; | 300 typedef union { double d; ULong L[2]; } U; |
289 | 301 |
290 #ifdef YES_ALIAS | |
291 #define dval(x) x | |
292 #ifdef IEEE_8087 | 302 #ifdef IEEE_8087 |
293 #define word0(x) ((ULong *)&x)[1] | 303 #define word0(x) (x)->L[1] |
294 #define word1(x) ((ULong *)&x)[0] | 304 #define word1(x) (x)->L[0] |
295 #else | 305 #else |
296 #define word0(x) ((ULong *)&x)[0] | 306 #define word0(x) (x)->L[0] |
297 #define word1(x) ((ULong *)&x)[1] | 307 #define word1(x) (x)->L[1] |
298 #endif | 308 #endif |
| 309 #define dval(x) (x)->d |
| 310 |
| 311 #ifndef STRTOD_DIGLIM |
| 312 #define STRTOD_DIGLIM 40 |
| 313 #endif |
| 314 |
| 315 #ifdef DIGLIM_DEBUG |
| 316 extern int strtod_diglim; |
299 #else | 317 #else |
300 #ifdef IEEE_8087 | 318 #define strtod_diglim STRTOD_DIGLIM |
301 #define word0(x) ((U*)&x)->L[1] | |
302 #define word1(x) ((U*)&x)->L[0] | |
303 #else | |
304 #define word0(x) ((U*)&x)->L[0] | |
305 #define word1(x) ((U*)&x)->L[1] | |
306 #endif | |
307 #define dval(x) ((U*)&x)->d | |
308 #endif | 319 #endif |
309 | 320 |
310 /* The following definition of Storeinc is appropriate for MIPS processors. | 321 /* The following definition of Storeinc is appropriate for MIPS processors. |
311 * An alternative that might be better on some machines is | 322 * An alternative that might be better on some machines is |
312 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) | 323 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) |
313 */ | 324 */ |
314 #if defined(IEEE_8087) + defined(VAX) | 325 #if defined(IEEE_8087) + defined(VAX) |
315 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ | 326 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ |
316 ((unsigned short *)a)[0] = (unsigned short)c, a++) | 327 ((unsigned short *)a)[0] = (unsigned short)c, a++) |
317 #else | 328 #else |
318 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ | 329 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ |
319 ((unsigned short *)a)[1] = (unsigned short)c, a++) | 330 ((unsigned short *)a)[1] = (unsigned short)c, a++) |
320 #endif | 331 #endif |
321 | 332 |
322 /* #define P DBL_MANT_DIG */ | 333 /* #define P DBL_MANT_DIG */ |
323 /* Ten_pmax = floor(P*log(2)/log(5)) */ | 334 /* Ten_pmax = floor(P*log(2)/log(5)) */ |
324 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ | 335 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ |
325 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ | 336 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ |
326 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ | 337 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ |
327 | 338 |
328 #ifdef IEEE_Arith | 339 #ifdef IEEE_Arith |
329 #define Exp_shift 20 | 340 #define Exp_shift 20 |
330 #define Exp_shift1 20 | 341 #define Exp_shift1 20 |
331 #define Exp_msk1 0x100000 | 342 #define Exp_msk1 0x100000 |
332 #define Exp_msk11 0x100000 | 343 #define Exp_msk11 0x100000 |
333 #define Exp_mask 0x7ff00000 | 344 #define Exp_mask 0x7ff00000 |
334 #define P 53 | 345 #define P 53 |
| 346 #define Nbits 53 |
335 #define Bias 1023 | 347 #define Bias 1023 |
| 348 #define Emax 1023 |
336 #define Emin (-1022) | 349 #define Emin (-1022) |
337 #define Exp_1 0x3ff00000 | 350 #define Exp_1 0x3ff00000 |
338 #define Exp_11 0x3ff00000 | 351 #define Exp_11 0x3ff00000 |
339 #define Ebits 11 | 352 #define Ebits 11 |
340 #define Frac_mask 0xfffff | 353 #define Frac_mask 0xfffff |
341 #define Frac_mask1 0xfffff | 354 #define Frac_mask1 0xfffff |
342 #define Ten_pmax 22 | 355 #define Ten_pmax 22 |
343 #define Bletch 0x10 | 356 #define Bletch 0x10 |
344 #define Bndry_mask 0xfffff | 357 #define Bndry_mask 0xfffff |
345 #define Bndry_mask1 0xfffff | 358 #define Bndry_mask1 0xfffff |
(...skipping 34 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
380 #define Sudden_Underflow | 393 #define Sudden_Underflow |
381 #ifdef IBM | 394 #ifdef IBM |
382 #undef Flt_Rounds | 395 #undef Flt_Rounds |
383 #define Flt_Rounds 0 | 396 #define Flt_Rounds 0 |
384 #define Exp_shift 24 | 397 #define Exp_shift 24 |
385 #define Exp_shift1 24 | 398 #define Exp_shift1 24 |
386 #define Exp_msk1 0x1000000 | 399 #define Exp_msk1 0x1000000 |
387 #define Exp_msk11 0x1000000 | 400 #define Exp_msk11 0x1000000 |
388 #define Exp_mask 0x7f000000 | 401 #define Exp_mask 0x7f000000 |
389 #define P 14 | 402 #define P 14 |
| 403 #define Nbits 56 |
390 #define Bias 65 | 404 #define Bias 65 |
| 405 #define Emax 248 |
| 406 #define Emin (-260) |
391 #define Exp_1 0x41000000 | 407 #define Exp_1 0x41000000 |
392 #define Exp_11 0x41000000 | 408 #define Exp_11 0x41000000 |
393 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ | 409 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ |
394 #define Frac_mask 0xffffff | 410 #define Frac_mask 0xffffff |
395 #define Frac_mask1 0xffffff | 411 #define Frac_mask1 0xffffff |
396 #define Bletch 4 | 412 #define Bletch 4 |
397 #define Ten_pmax 22 | 413 #define Ten_pmax 22 |
398 #define Bndry_mask 0xefffff | 414 #define Bndry_mask 0xefffff |
399 #define Bndry_mask1 0xffffff | 415 #define Bndry_mask1 0xffffff |
400 #define LSB 1 | 416 #define LSB 1 |
401 #define Sign_bit 0x80000000 | 417 #define Sign_bit 0x80000000 |
402 #define Log2P 4 | 418 #define Log2P 4 |
403 #define Tiny0 0x100000 | 419 #define Tiny0 0x100000 |
404 #define Tiny1 0 | 420 #define Tiny1 0 |
405 #define Quick_max 14 | 421 #define Quick_max 14 |
406 #define Int_max 15 | 422 #define Int_max 15 |
407 #else /* VAX */ | 423 #else /* VAX */ |
408 #undef Flt_Rounds | 424 #undef Flt_Rounds |
409 #define Flt_Rounds 1 | 425 #define Flt_Rounds 1 |
410 #define Exp_shift 23 | 426 #define Exp_shift 23 |
411 #define Exp_shift1 7 | 427 #define Exp_shift1 7 |
412 #define Exp_msk1 0x80 | 428 #define Exp_msk1 0x80 |
413 #define Exp_msk11 0x800000 | 429 #define Exp_msk11 0x800000 |
414 #define Exp_mask 0x7f80 | 430 #define Exp_mask 0x7f80 |
415 #define P 56 | 431 #define P 56 |
| 432 #define Nbits 56 |
416 #define Bias 129 | 433 #define Bias 129 |
| 434 #define Emax 126 |
| 435 #define Emin (-129) |
417 #define Exp_1 0x40800000 | 436 #define Exp_1 0x40800000 |
418 #define Exp_11 0x4080 | 437 #define Exp_11 0x4080 |
419 #define Ebits 8 | 438 #define Ebits 8 |
420 #define Frac_mask 0x7fffff | 439 #define Frac_mask 0x7fffff |
421 #define Frac_mask1 0xffff007f | 440 #define Frac_mask1 0xffff007f |
422 #define Ten_pmax 24 | 441 #define Ten_pmax 24 |
423 #define Bletch 2 | 442 #define Bletch 2 |
424 #define Bndry_mask 0xffff007f | 443 #define Bndry_mask 0xffff007f |
425 #define Bndry_mask1 0xffff007f | 444 #define Bndry_mask1 0xffff007f |
426 #define LSB 0x10000 | 445 #define LSB 0x10000 |
(...skipping 23 matching lines...) Expand all Loading... |
450 #define rounded_quotient(a,b) a /= b | 469 #define rounded_quotient(a,b) a /= b |
451 #endif | 470 #endif |
452 | 471 |
453 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) | 472 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) |
454 #define Big1 0xffffffff | 473 #define Big1 0xffffffff |
455 | 474 |
456 #ifndef Pack_32 | 475 #ifndef Pack_32 |
457 #define Pack_32 | 476 #define Pack_32 |
458 #endif | 477 #endif |
459 | 478 |
| 479 typedef struct BCinfo BCinfo; |
| 480 struct |
| 481 BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflc
hk; }; |
| 482 |
460 #ifdef KR_headers | 483 #ifdef KR_headers |
461 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff) | 484 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff) |
462 #else | 485 #else |
463 #define FFFFFFFF 0xffffffffUL | 486 #define FFFFFFFF 0xffffffffUL |
464 #endif | 487 #endif |
465 | 488 |
466 #ifdef NO_LONG_LONG | 489 #ifdef NO_LONG_LONG |
467 #undef ULLong | 490 #undef ULLong |
468 #ifdef Just_16 | 491 #ifdef Just_16 |
469 #undef Pack_32 | 492 #undef Pack_32 |
(...skipping 10 matching lines...) Expand all Loading... |
480 #ifndef ULLong | 503 #ifndef ULLong |
481 #define ULLong unsigned Llong | 504 #define ULLong unsigned Llong |
482 #endif | 505 #endif |
483 #endif /* NO_LONG_LONG */ | 506 #endif /* NO_LONG_LONG */ |
484 | 507 |
485 #ifndef MULTIPLE_THREADS | 508 #ifndef MULTIPLE_THREADS |
486 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/ | 509 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/ |
487 #define FREE_DTOA_LOCK(n) /*nothing*/ | 510 #define FREE_DTOA_LOCK(n) /*nothing*/ |
488 #endif | 511 #endif |
489 | 512 |
490 #define Kmax 15 | 513 #define Kmax 7 |
491 | 514 |
492 double strtod(const char *s00, char **se); | 515 double strtod(const char *s00, char **se); |
493 char *dtoa(double d, int mode, int ndigits, | 516 char *dtoa(double d, int mode, int ndigits, |
494 int *decpt, int *sign, char **rve); | 517 int *decpt, int *sign, char **rve); |
495 | 518 |
496 struct | 519 struct |
497 Bigint { | 520 Bigint { |
498 struct Bigint *next; | 521 struct Bigint *next; |
499 int k, maxwds, sign, wds; | 522 int k, maxwds, sign, wds; |
500 ULong x[1]; | 523 ULong x[1]; |
(...skipping 11 matching lines...) Expand all Loading... |
512 (int k) | 535 (int k) |
513 #endif | 536 #endif |
514 { | 537 { |
515 int x; | 538 int x; |
516 Bigint *rv; | 539 Bigint *rv; |
517 #ifndef Omit_Private_Memory | 540 #ifndef Omit_Private_Memory |
518 unsigned int len; | 541 unsigned int len; |
519 #endif | 542 #endif |
520 | 543 |
521 ACQUIRE_DTOA_LOCK(0); | 544 ACQUIRE_DTOA_LOCK(0); |
522 » if ((rv = freelist[k])) { | 545 » /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */ |
| 546 » /* but this case seems very unlikely. */ |
| 547 » if (k <= Kmax && (rv = freelist[k])) |
523 freelist[k] = rv->next; | 548 freelist[k] = rv->next; |
524 } | |
525 else { | 549 else { |
526 x = 1 << k; | 550 x = 1 << k; |
527 #ifdef Omit_Private_Memory | 551 #ifdef Omit_Private_Memory |
528 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); | 552 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); |
529 #else | 553 #else |
530 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1
) | 554 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1
) |
531 /sizeof(double); | 555 /sizeof(double); |
532 » » if (pmem_next - private_mem + len <= PRIVATE_mem) { | 556 » » if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) { |
533 rv = (Bigint*)pmem_next; | 557 rv = (Bigint*)pmem_next; |
534 pmem_next += len; | 558 pmem_next += len; |
535 } | 559 } |
536 else | 560 else |
537 rv = (Bigint*)MALLOC(len*sizeof(double)); | 561 rv = (Bigint*)MALLOC(len*sizeof(double)); |
538 #endif | 562 #endif |
539 rv->k = k; | 563 rv->k = k; |
540 rv->maxwds = x; | 564 rv->maxwds = x; |
541 } | 565 } |
542 FREE_DTOA_LOCK(0); | 566 FREE_DTOA_LOCK(0); |
543 rv->sign = rv->wds = 0; | 567 rv->sign = rv->wds = 0; |
544 return rv; | 568 return rv; |
545 } | 569 } |
546 | 570 |
547 static void | 571 static void |
548 Bfree | 572 Bfree |
549 #ifdef KR_headers | 573 #ifdef KR_headers |
550 (v) Bigint *v; | 574 (v) Bigint *v; |
551 #else | 575 #else |
552 (Bigint *v) | 576 (Bigint *v) |
553 #endif | 577 #endif |
554 { | 578 { |
555 if (v) { | 579 if (v) { |
556 » » ACQUIRE_DTOA_LOCK(0); | 580 » » if (v->k > Kmax) |
557 » » v->next = freelist[v->k]; | 581 #ifdef FREE |
558 » » freelist[v->k] = v; | 582 » » » FREE((void*)v); |
559 » » FREE_DTOA_LOCK(0); | 583 #else |
| 584 » » » free((void*)v); |
| 585 #endif |
| 586 » » else { |
| 587 » » » ACQUIRE_DTOA_LOCK(0); |
| 588 » » » v->next = freelist[v->k]; |
| 589 » » » freelist[v->k] = v; |
| 590 » » » FREE_DTOA_LOCK(0); |
| 591 » » » } |
560 } | 592 } |
561 } | 593 } |
562 | 594 |
563 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ | 595 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ |
564 y->wds*sizeof(Long) + 2*sizeof(int)) | 596 y->wds*sizeof(Long) + 2*sizeof(int)) |
565 | 597 |
566 static Bigint * | 598 static Bigint * |
567 multadd | 599 multadd |
568 #ifdef KR_headers | 600 #ifdef KR_headers |
569 (b, m, a) Bigint *b; int m, a; | 601 (b, m, a) Bigint *b; int m, a; |
(...skipping 46 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
616 } | 648 } |
617 b->x[wds++] = carry; | 649 b->x[wds++] = carry; |
618 b->wds = wds; | 650 b->wds = wds; |
619 } | 651 } |
620 return b; | 652 return b; |
621 } | 653 } |
622 | 654 |
623 static Bigint * | 655 static Bigint * |
624 s2b | 656 s2b |
625 #ifdef KR_headers | 657 #ifdef KR_headers |
626 » (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9; | 658 » (s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9; |
627 #else | 659 #else |
628 » (CONST char *s, int nd0, int nd, ULong y9) | 660 » (CONST char *s, int nd0, int nd, ULong y9, int dplen) |
629 #endif | 661 #endif |
630 { | 662 { |
631 Bigint *b; | 663 Bigint *b; |
632 int i, k; | 664 int i, k; |
633 Long x, y; | 665 Long x, y; |
634 | 666 |
635 x = (nd + 8) / 9; | 667 x = (nd + 8) / 9; |
636 for(k = 0, y = 1; x > y; y <<= 1, k++) ; | 668 for(k = 0, y = 1; x > y; y <<= 1, k++) ; |
637 #ifdef Pack_32 | 669 #ifdef Pack_32 |
638 b = Balloc(k); | 670 b = Balloc(k); |
639 b->x[0] = y9; | 671 b->x[0] = y9; |
640 b->wds = 1; | 672 b->wds = 1; |
641 #else | 673 #else |
642 b = Balloc(k+1); | 674 b = Balloc(k+1); |
643 b->x[0] = y9 & 0xffff; | 675 b->x[0] = y9 & 0xffff; |
644 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; | 676 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; |
645 #endif | 677 #endif |
646 | 678 |
647 i = 9; | 679 i = 9; |
648 if (9 < nd0) { | 680 if (9 < nd0) { |
649 s += 9; | 681 s += 9; |
650 do b = multadd(b, 10, *s++ - '0'); | 682 do b = multadd(b, 10, *s++ - '0'); |
651 while(++i < nd0); | 683 while(++i < nd0); |
652 » » s++; | 684 » » s += dplen; |
653 } | 685 } |
654 else | 686 else |
655 » » s += 10; | 687 » » s += dplen + 9; |
656 for(; i < nd; i++) | 688 for(; i < nd; i++) |
657 b = multadd(b, 10, *s++ - '0'); | 689 b = multadd(b, 10, *s++ - '0'); |
658 return b; | 690 return b; |
659 } | 691 } |
660 | 692 |
661 static int | 693 static int |
662 hi0bits | 694 hi0bits |
663 #ifdef KR_headers | 695 #ifdef KR_headers |
664 » (x) register ULong x; | 696 » (x) ULong x; |
665 #else | 697 #else |
666 » (register ULong x) | 698 » (ULong x) |
667 #endif | 699 #endif |
668 { | 700 { |
669 » register int k = 0; | 701 » int k = 0; |
670 | 702 |
671 if (!(x & 0xffff0000)) { | 703 if (!(x & 0xffff0000)) { |
672 k = 16; | 704 k = 16; |
673 x <<= 16; | 705 x <<= 16; |
674 } | 706 } |
675 if (!(x & 0xff000000)) { | 707 if (!(x & 0xff000000)) { |
676 k += 8; | 708 k += 8; |
677 x <<= 8; | 709 x <<= 8; |
678 } | 710 } |
679 if (!(x & 0xf0000000)) { | 711 if (!(x & 0xf0000000)) { |
(...skipping 13 matching lines...) Expand all Loading... |
693 } | 725 } |
694 | 726 |
695 static int | 727 static int |
696 lo0bits | 728 lo0bits |
697 #ifdef KR_headers | 729 #ifdef KR_headers |
698 (y) ULong *y; | 730 (y) ULong *y; |
699 #else | 731 #else |
700 (ULong *y) | 732 (ULong *y) |
701 #endif | 733 #endif |
702 { | 734 { |
703 » register int k; | 735 » int k; |
704 » register ULong x = *y; | 736 » ULong x = *y; |
705 | 737 |
706 if (x & 7) { | 738 if (x & 7) { |
707 if (x & 1) | 739 if (x & 1) |
708 return 0; | 740 return 0; |
709 if (x & 2) { | 741 if (x & 2) { |
710 *y = x >> 1; | 742 *y = x >> 1; |
711 return 1; | 743 return 1; |
712 } | 744 } |
713 *y = x >> 2; | 745 *y = x >> 2; |
714 return 2; | 746 return 2; |
(...skipping 394 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
1109 #endif | 1141 #endif |
1110 while(!*--xc) | 1142 while(!*--xc) |
1111 wa--; | 1143 wa--; |
1112 c->wds = wa; | 1144 c->wds = wa; |
1113 return c; | 1145 return c; |
1114 } | 1146 } |
1115 | 1147 |
1116 static double | 1148 static double |
1117 ulp | 1149 ulp |
1118 #ifdef KR_headers | 1150 #ifdef KR_headers |
1119 » (x) double x; | 1151 » (x) U *x; |
1120 #else | 1152 #else |
1121 » (double x) | 1153 » (U *x) |
1122 #endif | 1154 #endif |
1123 { | 1155 { |
1124 » register Long L; | 1156 » Long L; |
1125 » double a; | 1157 » U u; |
1126 | 1158 |
1127 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; | 1159 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; |
1128 #ifndef Avoid_Underflow | 1160 #ifndef Avoid_Underflow |
1129 #ifndef Sudden_Underflow | 1161 #ifndef Sudden_Underflow |
1130 if (L > 0) { | 1162 if (L > 0) { |
1131 #endif | 1163 #endif |
1132 #endif | 1164 #endif |
1133 #ifdef IBM | 1165 #ifdef IBM |
1134 L |= Exp_msk1 >> 4; | 1166 L |= Exp_msk1 >> 4; |
1135 #endif | 1167 #endif |
1136 » » word0(a) = L; | 1168 » » word0(&u) = L; |
1137 » » word1(a) = 0; | 1169 » » word1(&u) = 0; |
1138 #ifndef Avoid_Underflow | 1170 #ifndef Avoid_Underflow |
1139 #ifndef Sudden_Underflow | 1171 #ifndef Sudden_Underflow |
1140 } | 1172 } |
1141 else { | 1173 else { |
1142 L = -L >> Exp_shift; | 1174 L = -L >> Exp_shift; |
1143 if (L < Exp_shift) { | 1175 if (L < Exp_shift) { |
1144 » » » word0(a) = 0x80000 >> L; | 1176 » » » word0(&u) = 0x80000 >> L; |
1145 » » » word1(a) = 0; | 1177 » » » word1(&u) = 0; |
1146 } | 1178 } |
1147 else { | 1179 else { |
1148 » » » word0(a) = 0; | 1180 » » » word0(&u) = 0; |
1149 L -= Exp_shift; | 1181 L -= Exp_shift; |
1150 » » » word1(a) = L >= 31 ? 1 : 1 << 31 - L; | 1182 » » » word1(&u) = L >= 31 ? 1 : 1 << 31 - L; |
1151 } | 1183 } |
1152 } | 1184 } |
1153 #endif | 1185 #endif |
1154 #endif | 1186 #endif |
1155 » return dval(a); | 1187 » return dval(&u); |
1156 } | 1188 } |
1157 | 1189 |
1158 static double | 1190 static double |
1159 b2d | 1191 b2d |
1160 #ifdef KR_headers | 1192 #ifdef KR_headers |
1161 (a, e) Bigint *a; int *e; | 1193 (a, e) Bigint *a; int *e; |
1162 #else | 1194 #else |
1163 (Bigint *a, int *e) | 1195 (Bigint *a, int *e) |
1164 #endif | 1196 #endif |
1165 { | 1197 { |
1166 ULong *xa, *xa0, w, y, z; | 1198 ULong *xa, *xa0, w, y, z; |
1167 int k; | 1199 int k; |
1168 » double d; | 1200 » U d; |
1169 #ifdef VAX | 1201 #ifdef VAX |
1170 ULong d0, d1; | 1202 ULong d0, d1; |
1171 #else | 1203 #else |
1172 #define d0 word0(d) | 1204 #define d0 word0(&d) |
1173 #define d1 word1(d) | 1205 #define d1 word1(&d) |
1174 #endif | 1206 #endif |
1175 | 1207 |
1176 xa0 = a->x; | 1208 xa0 = a->x; |
1177 xa = xa0 + a->wds; | 1209 xa = xa0 + a->wds; |
1178 y = *--xa; | 1210 y = *--xa; |
1179 #ifdef DEBUG | 1211 #ifdef DEBUG |
1180 if (!y) Bug("zero y in b2d"); | 1212 if (!y) Bug("zero y in b2d"); |
1181 #endif | 1213 #endif |
1182 k = hi0bits(y); | 1214 k = hi0bits(y); |
1183 *e = 32 - k; | 1215 *e = 32 - k; |
1184 #ifdef Pack_32 | 1216 #ifdef Pack_32 |
1185 if (k < Ebits) { | 1217 if (k < Ebits) { |
1186 » » d0 = Exp_1 | y >> Ebits - k; | 1218 » » d0 = Exp_1 | y >> (Ebits - k); |
1187 w = xa > xa0 ? *--xa : 0; | 1219 w = xa > xa0 ? *--xa : 0; |
1188 » » d1 = y << (32-Ebits) + k | w >> Ebits - k; | 1220 » » d1 = y << ((32-Ebits) + k) | w >> (Ebits - k); |
1189 goto ret_d; | 1221 goto ret_d; |
1190 } | 1222 } |
1191 z = xa > xa0 ? *--xa : 0; | 1223 z = xa > xa0 ? *--xa : 0; |
1192 if (k -= Ebits) { | 1224 if (k -= Ebits) { |
1193 » » d0 = Exp_1 | y << k | z >> 32 - k; | 1225 » » d0 = Exp_1 | y << k | z >> (32 - k); |
1194 y = xa > xa0 ? *--xa : 0; | 1226 y = xa > xa0 ? *--xa : 0; |
1195 » » d1 = z << k | y >> 32 - k; | 1227 » » d1 = z << k | y >> (32 - k); |
1196 } | 1228 } |
1197 else { | 1229 else { |
1198 d0 = Exp_1 | y; | 1230 d0 = Exp_1 | y; |
1199 d1 = z; | 1231 d1 = z; |
1200 } | 1232 } |
1201 #else | 1233 #else |
1202 if (k < Ebits + 16) { | 1234 if (k < Ebits + 16) { |
1203 z = xa > xa0 ? *--xa : 0; | 1235 z = xa > xa0 ? *--xa : 0; |
1204 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; | 1236 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; |
1205 w = xa > xa0 ? *--xa : 0; | 1237 w = xa > xa0 ? *--xa : 0; |
1206 y = xa > xa0 ? *--xa : 0; | 1238 y = xa > xa0 ? *--xa : 0; |
1207 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; | 1239 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; |
1208 goto ret_d; | 1240 goto ret_d; |
1209 } | 1241 } |
1210 z = xa > xa0 ? *--xa : 0; | 1242 z = xa > xa0 ? *--xa : 0; |
1211 w = xa > xa0 ? *--xa : 0; | 1243 w = xa > xa0 ? *--xa : 0; |
1212 k -= Ebits + 16; | 1244 k -= Ebits + 16; |
1213 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; | 1245 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; |
1214 y = xa > xa0 ? *--xa : 0; | 1246 y = xa > xa0 ? *--xa : 0; |
1215 d1 = w << k + 16 | y << k; | 1247 d1 = w << k + 16 | y << k; |
1216 #endif | 1248 #endif |
1217 ret_d: | 1249 ret_d: |
1218 #ifdef VAX | 1250 #ifdef VAX |
1219 » word0(d) = d0 >> 16 | d0 << 16; | 1251 » word0(&d) = d0 >> 16 | d0 << 16; |
1220 » word1(d) = d1 >> 16 | d1 << 16; | 1252 » word1(&d) = d1 >> 16 | d1 << 16; |
1221 #else | 1253 #else |
1222 #undef d0 | 1254 #undef d0 |
1223 #undef d1 | 1255 #undef d1 |
1224 #endif | 1256 #endif |
1225 » return dval(d); | 1257 » return dval(&d); |
1226 } | 1258 } |
1227 | 1259 |
1228 static Bigint * | 1260 static Bigint * |
1229 d2b | 1261 d2b |
1230 #ifdef KR_headers | 1262 #ifdef KR_headers |
1231 » (d, e, bits) double d; int *e, *bits; | 1263 » (d, e, bits) U *d; int *e, *bits; |
1232 #else | 1264 #else |
1233 » (double d, int *e, int *bits) | 1265 » (U *d, int *e, int *bits) |
1234 #endif | 1266 #endif |
1235 { | 1267 { |
1236 Bigint *b; | 1268 Bigint *b; |
1237 int de, k; | 1269 int de, k; |
1238 ULong *x, y, z; | 1270 ULong *x, y, z; |
1239 #ifndef Sudden_Underflow | 1271 #ifndef Sudden_Underflow |
1240 int i; | 1272 int i; |
1241 #endif | 1273 #endif |
1242 #ifdef VAX | 1274 #ifdef VAX |
1243 ULong d0, d1; | 1275 ULong d0, d1; |
(...skipping 18 matching lines...) Expand all Loading... |
1262 #ifndef IBM | 1294 #ifndef IBM |
1263 z |= Exp_msk11; | 1295 z |= Exp_msk11; |
1264 #endif | 1296 #endif |
1265 #else | 1297 #else |
1266 if ((de = (int)(d0 >> Exp_shift))) | 1298 if ((de = (int)(d0 >> Exp_shift))) |
1267 z |= Exp_msk1; | 1299 z |= Exp_msk1; |
1268 #endif | 1300 #endif |
1269 #ifdef Pack_32 | 1301 #ifdef Pack_32 |
1270 if ((y = d1)) { | 1302 if ((y = d1)) { |
1271 if ((k = lo0bits(&y))) { | 1303 if ((k = lo0bits(&y))) { |
1272 » » » x[0] = y | z << 32 - k; | 1304 » » » x[0] = y | z << (32 - k); |
1273 z >>= k; | 1305 z >>= k; |
1274 } | 1306 } |
1275 else | 1307 else |
1276 x[0] = y; | 1308 x[0] = y; |
1277 #ifndef Sudden_Underflow | 1309 #ifndef Sudden_Underflow |
1278 i = | 1310 i = |
1279 #endif | 1311 #endif |
1280 b->wds = (x[1] = z) ? 2 : 1; | 1312 b->wds = (x[1] = z) ? 2 : 1; |
1281 } | 1313 } |
1282 else { | 1314 else { |
1283 #ifdef DEBUG | |
1284 if (!z) | |
1285 Bug("Zero passed to d2b"); | |
1286 #endif | |
1287 k = lo0bits(&z); | 1315 k = lo0bits(&z); |
1288 x[0] = z; | 1316 x[0] = z; |
1289 #ifndef Sudden_Underflow | 1317 #ifndef Sudden_Underflow |
1290 i = | 1318 i = |
1291 #endif | 1319 #endif |
1292 b->wds = 1; | 1320 b->wds = 1; |
1293 k += 32; | 1321 k += 32; |
1294 } | 1322 } |
1295 #else | 1323 #else |
1296 if (y = d1) { | 1324 if (y = d1) { |
(...skipping 67 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
1364 #undef d1 | 1392 #undef d1 |
1365 | 1393 |
1366 static double | 1394 static double |
1367 ratio | 1395 ratio |
1368 #ifdef KR_headers | 1396 #ifdef KR_headers |
1369 (a, b) Bigint *a, *b; | 1397 (a, b) Bigint *a, *b; |
1370 #else | 1398 #else |
1371 (Bigint *a, Bigint *b) | 1399 (Bigint *a, Bigint *b) |
1372 #endif | 1400 #endif |
1373 { | 1401 { |
1374 » double da, db; | 1402 » U da, db; |
1375 int k, ka, kb; | 1403 int k, ka, kb; |
1376 | 1404 |
1377 » dval(da) = b2d(a, &ka); | 1405 » dval(&da) = b2d(a, &ka); |
1378 » dval(db) = b2d(b, &kb); | 1406 » dval(&db) = b2d(b, &kb); |
1379 #ifdef Pack_32 | 1407 #ifdef Pack_32 |
1380 k = ka - kb + 32*(a->wds - b->wds); | 1408 k = ka - kb + 32*(a->wds - b->wds); |
1381 #else | 1409 #else |
1382 k = ka - kb + 16*(a->wds - b->wds); | 1410 k = ka - kb + 16*(a->wds - b->wds); |
1383 #endif | 1411 #endif |
1384 #ifdef IBM | 1412 #ifdef IBM |
1385 if (k > 0) { | 1413 if (k > 0) { |
1386 » » word0(da) += (k >> 2)*Exp_msk1; | 1414 » » word0(&da) += (k >> 2)*Exp_msk1; |
1387 if (k &= 3) | 1415 if (k &= 3) |
1388 » » » dval(da) *= 1 << k; | 1416 » » » dval(&da) *= 1 << k; |
1389 } | 1417 } |
1390 else { | 1418 else { |
1391 k = -k; | 1419 k = -k; |
1392 » » word0(db) += (k >> 2)*Exp_msk1; | 1420 » » word0(&db) += (k >> 2)*Exp_msk1; |
1393 if (k &= 3) | 1421 if (k &= 3) |
1394 » » » dval(db) *= 1 << k; | 1422 » » » dval(&db) *= 1 << k; |
1395 } | 1423 } |
1396 #else | 1424 #else |
1397 if (k > 0) | 1425 if (k > 0) |
1398 » » word0(da) += k*Exp_msk1; | 1426 » » word0(&da) += k*Exp_msk1; |
1399 else { | 1427 else { |
1400 k = -k; | 1428 k = -k; |
1401 » » word0(db) += k*Exp_msk1; | 1429 » » word0(&db) += k*Exp_msk1; |
1402 } | 1430 } |
1403 #endif | 1431 #endif |
1404 » return dval(da) / dval(db); | 1432 » return dval(&da) / dval(&db); |
1405 } | 1433 } |
1406 | 1434 |
1407 static CONST double | 1435 static CONST double |
1408 tens[] = { | 1436 tens[] = { |
1409 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, | 1437 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, |
1410 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, | 1438 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, |
1411 1e20, 1e21, 1e22 | 1439 1e20, 1e21, 1e22 |
1412 #ifdef VAX | 1440 #ifdef VAX |
1413 , 1e23, 1e24 | 1441 , 1e23, 1e24 |
1414 #endif | 1442 #endif |
(...skipping 19 matching lines...) Expand all Loading... |
1434 bigtens[] = { 1e16, 1e32, 1e64 }; | 1462 bigtens[] = { 1e16, 1e32, 1e64 }; |
1435 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; | 1463 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; |
1436 #define n_bigtens 3 | 1464 #define n_bigtens 3 |
1437 #else | 1465 #else |
1438 bigtens[] = { 1e16, 1e32 }; | 1466 bigtens[] = { 1e16, 1e32 }; |
1439 static CONST double tinytens[] = { 1e-16, 1e-32 }; | 1467 static CONST double tinytens[] = { 1e-16, 1e-32 }; |
1440 #define n_bigtens 2 | 1468 #define n_bigtens 2 |
1441 #endif | 1469 #endif |
1442 #endif | 1470 #endif |
1443 | 1471 |
| 1472 #undef Need_Hexdig |
| 1473 #ifdef INFNAN_CHECK |
| 1474 #ifndef No_Hex_NaN |
| 1475 #define Need_Hexdig |
| 1476 #endif |
| 1477 #endif |
| 1478 |
| 1479 #ifndef Need_Hexdig |
| 1480 #ifndef NO_HEX_FP |
| 1481 #define Need_Hexdig |
| 1482 #endif |
| 1483 #endif |
| 1484 |
| 1485 #ifdef Need_Hexdig /*{*/ |
| 1486 static unsigned char hexdig[256]; |
| 1487 |
| 1488 static void |
| 1489 #ifdef KR_headers |
| 1490 htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc; |
| 1491 #else |
| 1492 htinit(unsigned char *h, unsigned char *s, int inc) |
| 1493 #endif |
| 1494 { |
| 1495 int i, j; |
| 1496 for(i = 0; (j = s[i]) !=0; i++) |
| 1497 h[j] = i + inc; |
| 1498 } |
| 1499 |
| 1500 static void |
| 1501 #ifdef KR_headers |
| 1502 hexdig_init() |
| 1503 #else |
| 1504 hexdig_init(void) |
| 1505 #endif |
| 1506 { |
| 1507 #define USC (unsigned char *) |
| 1508 htinit(hexdig, USC "0123456789", 0x10); |
| 1509 htinit(hexdig, USC "abcdef", 0x10 + 10); |
| 1510 htinit(hexdig, USC "ABCDEF", 0x10 + 10); |
| 1511 } |
| 1512 #endif /* } Need_Hexdig */ |
| 1513 |
1444 #ifdef INFNAN_CHECK | 1514 #ifdef INFNAN_CHECK |
1445 | 1515 |
1446 #ifndef NAN_WORD0 | 1516 #ifndef NAN_WORD0 |
1447 #define NAN_WORD0 0x7ff80000 | 1517 #define NAN_WORD0 0x7ff80000 |
1448 #endif | 1518 #endif |
1449 | 1519 |
1450 #ifndef NAN_WORD1 | 1520 #ifndef NAN_WORD1 |
1451 #define NAN_WORD1 0 | 1521 #define NAN_WORD1 0 |
1452 #endif | 1522 #endif |
1453 | 1523 |
(...skipping 15 matching lines...) Expand all Loading... |
1469 return 0; | 1539 return 0; |
1470 } | 1540 } |
1471 *sp = s + 1; | 1541 *sp = s + 1; |
1472 return 1; | 1542 return 1; |
1473 } | 1543 } |
1474 | 1544 |
1475 #ifndef No_Hex_NaN | 1545 #ifndef No_Hex_NaN |
1476 static void | 1546 static void |
1477 hexnan | 1547 hexnan |
1478 #ifdef KR_headers | 1548 #ifdef KR_headers |
1479 » (rvp, sp) double *rvp; CONST char **sp; | 1549 » (rvp, sp) U *rvp; CONST char **sp; |
1480 #else | 1550 #else |
1481 » (double *rvp, CONST char **sp) | 1551 » (U *rvp, CONST char **sp) |
1482 #endif | 1552 #endif |
1483 { | 1553 { |
1484 ULong c, x[2]; | 1554 ULong c, x[2]; |
1485 CONST char *s; | 1555 CONST char *s; |
1486 » int havedig, udx0, xshift; | 1556 » int c1, havedig, udx0, xshift; |
1487 | 1557 |
| 1558 if (!hexdig['0']) |
| 1559 hexdig_init(); |
1488 x[0] = x[1] = 0; | 1560 x[0] = x[1] = 0; |
1489 havedig = xshift = 0; | 1561 havedig = xshift = 0; |
1490 udx0 = 1; | 1562 udx0 = 1; |
1491 s = *sp; | 1563 s = *sp; |
1492 /* allow optional initial 0x or 0X */ | 1564 /* allow optional initial 0x or 0X */ |
1493 while((c = *(CONST unsigned char*)(s+1)) && c <= ' ') | 1565 while((c = *(CONST unsigned char*)(s+1)) && c <= ' ') |
1494 ++s; | 1566 ++s; |
1495 if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')) | 1567 if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')) |
1496 s += 2; | 1568 s += 2; |
1497 while((c = *(CONST unsigned char*)++s)) { | 1569 while((c = *(CONST unsigned char*)++s)) { |
1498 » » if (c >= '0' && c <= '9') | 1570 » » if ((c1 = hexdig[c])) |
1499 » » » c -= '0'; | 1571 » » » c = c1 & 0xf; |
1500 » » else if (c >= 'a' && c <= 'f') | |
1501 » » » c += 10 - 'a'; | |
1502 » » else if (c >= 'A' && c <= 'F') | |
1503 » » » c += 10 - 'A'; | |
1504 else if (c <= ' ') { | 1572 else if (c <= ' ') { |
1505 if (udx0 && havedig) { | 1573 if (udx0 && havedig) { |
1506 udx0 = 0; | 1574 udx0 = 0; |
1507 xshift = 1; | 1575 xshift = 1; |
1508 } | 1576 } |
1509 continue; | 1577 continue; |
1510 } | 1578 } |
1511 #ifdef GDTOA_NON_PEDANTIC_NANCHECK | 1579 #ifdef GDTOA_NON_PEDANTIC_NANCHECK |
1512 else if (/*(*/ c == ')' && havedig) { | 1580 else if (/*(*/ c == ')' && havedig) { |
1513 *sp = s + 1; | 1581 *sp = s + 1; |
(...skipping 16 matching lines...) Expand all Loading... |
1530 if (xshift) { | 1598 if (xshift) { |
1531 xshift = 0; | 1599 xshift = 0; |
1532 x[0] = x[1]; | 1600 x[0] = x[1]; |
1533 x[1] = 0; | 1601 x[1] = 0; |
1534 } | 1602 } |
1535 if (udx0) | 1603 if (udx0) |
1536 x[0] = (x[0] << 4) | (x[1] >> 28); | 1604 x[0] = (x[0] << 4) | (x[1] >> 28); |
1537 x[1] = (x[1] << 4) | c; | 1605 x[1] = (x[1] << 4) | c; |
1538 } | 1606 } |
1539 if ((x[0] &= 0xfffff) || x[1]) { | 1607 if ((x[0] &= 0xfffff) || x[1]) { |
1540 » » word0(*rvp) = Exp_mask | x[0]; | 1608 » » word0(rvp) = Exp_mask | x[0]; |
1541 » » word1(*rvp) = x[1]; | 1609 » » word1(rvp) = x[1]; |
1542 } | 1610 } |
1543 } | 1611 } |
1544 #endif /*No_Hex_NaN*/ | 1612 #endif /*No_Hex_NaN*/ |
1545 #endif /* INFNAN_CHECK */ | 1613 #endif /* INFNAN_CHECK */ |
1546 | 1614 |
| 1615 #ifdef Pack_32 |
| 1616 #define ULbits 32 |
| 1617 #define kshift 5 |
| 1618 #define kmask 31 |
| 1619 #else |
| 1620 #define ULbits 16 |
| 1621 #define kshift 4 |
| 1622 #define kmask 15 |
| 1623 #endif |
| 1624 #ifndef NO_HEX_FP /*{*/ |
| 1625 |
| 1626 static void |
| 1627 #ifdef KR_headers |
| 1628 rshift(b, k) Bigint *b; int k; |
| 1629 #else |
| 1630 rshift(Bigint *b, int k) |
| 1631 #endif |
| 1632 { |
| 1633 ULong *x, *x1, *xe, y; |
| 1634 int n; |
| 1635 |
| 1636 x = x1 = b->x; |
| 1637 n = k >> kshift; |
| 1638 if (n < b->wds) { |
| 1639 xe = x + b->wds; |
| 1640 x += n; |
| 1641 if (k &= kmask) { |
| 1642 n = 32 - k; |
| 1643 y = *x++ >> k; |
| 1644 while(x < xe) { |
| 1645 *x1++ = (y | (*x << n)) & 0xffffffff; |
| 1646 y = *x++ >> k; |
| 1647 } |
| 1648 if ((*x1 = y) !=0) |
| 1649 x1++; |
| 1650 } |
| 1651 else |
| 1652 while(x < xe) |
| 1653 *x1++ = *x++; |
| 1654 } |
| 1655 if ((b->wds = x1 - b->x) == 0) |
| 1656 b->x[0] = 0; |
| 1657 } |
| 1658 |
| 1659 static ULong |
| 1660 #ifdef KR_headers |
| 1661 any_on(b, k) Bigint *b; int k; |
| 1662 #else |
| 1663 any_on(Bigint *b, int k) |
| 1664 #endif |
| 1665 { |
| 1666 int n, nwds; |
| 1667 ULong *x, *x0, x1, x2; |
| 1668 |
| 1669 x = b->x; |
| 1670 nwds = b->wds; |
| 1671 n = k >> kshift; |
| 1672 if (n > nwds) |
| 1673 n = nwds; |
| 1674 else if (n < nwds && (k &= kmask)) { |
| 1675 x1 = x2 = x[n]; |
| 1676 x1 >>= k; |
| 1677 x1 <<= k; |
| 1678 if (x1 != x2) |
| 1679 return 1; |
| 1680 } |
| 1681 x0 = x; |
| 1682 x += n; |
| 1683 while(x > x0) |
| 1684 if (*--x) |
| 1685 return 1; |
| 1686 return 0; |
| 1687 } |
| 1688 |
| 1689 enum { /* rounding values: same as FLT_ROUNDS */ |
| 1690 Round_zero = 0, |
| 1691 Round_near = 1, |
| 1692 Round_up = 2, |
| 1693 Round_down = 3 |
| 1694 }; |
| 1695 |
| 1696 static Bigint * |
| 1697 #ifdef KR_headers |
| 1698 increment(b) Bigint *b; |
| 1699 #else |
| 1700 increment(Bigint *b) |
| 1701 #endif |
| 1702 { |
| 1703 ULong *x, *xe; |
| 1704 Bigint *b1; |
| 1705 |
| 1706 x = b->x; |
| 1707 xe = x + b->wds; |
| 1708 do { |
| 1709 if (*x < (ULong)0xffffffffL) { |
| 1710 ++*x; |
| 1711 return b; |
| 1712 } |
| 1713 *x++ = 0; |
| 1714 } while(x < xe); |
| 1715 { |
| 1716 if (b->wds >= b->maxwds) { |
| 1717 b1 = Balloc(b->k+1); |
| 1718 Bcopy(b1,b); |
| 1719 Bfree(b); |
| 1720 b = b1; |
| 1721 } |
| 1722 b->x[b->wds++] = 1; |
| 1723 } |
| 1724 return b; |
| 1725 } |
| 1726 |
| 1727 void |
| 1728 #ifdef KR_headers |
| 1729 gethex(sp, rvp, rounding, sign) |
| 1730 CONST char **sp; U *rvp; int rounding, sign; |
| 1731 #else |
| 1732 gethex( CONST char **sp, U *rvp, int rounding, int sign) |
| 1733 #endif |
| 1734 { |
| 1735 Bigint *b; |
| 1736 CONST unsigned char *decpt, *s0, *s, *s1; |
| 1737 Long e, e1; |
| 1738 ULong L, lostbits, *x; |
| 1739 int big, denorm, esign, havedig, k, n, nbits, up, zret; |
| 1740 #ifdef IBM |
| 1741 int j; |
| 1742 #endif |
| 1743 enum { |
| 1744 #ifdef IEEE_Arith /*{{*/ |
| 1745 emax = 0x7fe - Bias - P + 1, |
| 1746 emin = Emin - P + 1 |
| 1747 #else /*}{*/ |
| 1748 emin = Emin - P, |
| 1749 #ifdef VAX |
| 1750 emax = 0x7ff - Bias - P + 1 |
| 1751 #endif |
| 1752 #ifdef IBM |
| 1753 emax = 0x7f - Bias - P |
| 1754 #endif |
| 1755 #endif /*}}*/ |
| 1756 }; |
| 1757 #ifdef USE_LOCALE |
| 1758 int i; |
| 1759 #ifdef NO_LOCALE_CACHE |
| 1760 const unsigned char *decimalpoint = (unsigned char*) |
| 1761 localeconv()->decimal_point; |
| 1762 #else |
| 1763 const unsigned char *decimalpoint; |
| 1764 static unsigned char *decimalpoint_cache; |
| 1765 if (!(s0 = decimalpoint_cache)) { |
| 1766 s0 = (unsigned char*)localeconv()->decimal_point; |
| 1767 if ((decimalpoint_cache = (unsigned char*) |
| 1768 MALLOC(strlen((CONST char*)s0) + 1))) { |
| 1769 strcpy((char*)decimalpoint_cache, (CONST char*)s0); |
| 1770 s0 = decimalpoint_cache; |
| 1771 } |
| 1772 } |
| 1773 decimalpoint = s0; |
| 1774 #endif |
| 1775 #endif |
| 1776 |
| 1777 if (!hexdig['0']) |
| 1778 hexdig_init(); |
| 1779 havedig = 0; |
| 1780 s0 = *(CONST unsigned char **)sp + 2; |
| 1781 while(s0[havedig] == '0') |
| 1782 havedig++; |
| 1783 s0 += havedig; |
| 1784 s = s0; |
| 1785 decpt = 0; |
| 1786 zret = 0; |
| 1787 e = 0; |
| 1788 if (hexdig[*s]) |
| 1789 havedig++; |
| 1790 else { |
| 1791 zret = 1; |
| 1792 #ifdef USE_LOCALE |
| 1793 for(i = 0; decimalpoint[i]; ++i) { |
| 1794 if (s[i] != decimalpoint[i]) |
| 1795 goto pcheck; |
| 1796 } |
| 1797 decpt = s += i; |
| 1798 #else |
| 1799 if (*s != '.') |
| 1800 goto pcheck; |
| 1801 decpt = ++s; |
| 1802 #endif |
| 1803 if (!hexdig[*s]) |
| 1804 goto pcheck; |
| 1805 while(*s == '0') |
| 1806 s++; |
| 1807 if (hexdig[*s]) |
| 1808 zret = 0; |
| 1809 havedig = 1; |
| 1810 s0 = s; |
| 1811 } |
| 1812 while(hexdig[*s]) |
| 1813 s++; |
| 1814 #ifdef USE_LOCALE |
| 1815 if (*s == *decimalpoint && !decpt) { |
| 1816 for(i = 1; decimalpoint[i]; ++i) { |
| 1817 if (s[i] != decimalpoint[i]) |
| 1818 goto pcheck; |
| 1819 } |
| 1820 decpt = s += i; |
| 1821 #else |
| 1822 if (*s == '.' && !decpt) { |
| 1823 decpt = ++s; |
| 1824 #endif |
| 1825 while(hexdig[*s]) |
| 1826 s++; |
| 1827 }/*}*/ |
| 1828 if (decpt) |
| 1829 e = -(((Long)(s-decpt)) << 2); |
| 1830 pcheck: |
| 1831 s1 = s; |
| 1832 big = esign = 0; |
| 1833 switch(*s) { |
| 1834 case 'p': |
| 1835 case 'P': |
| 1836 switch(*++s) { |
| 1837 case '-': |
| 1838 esign = 1; |
| 1839 /* no break */ |
| 1840 case '+': |
| 1841 s++; |
| 1842 } |
| 1843 if ((n = hexdig[*s]) == 0 || n > 0x19) { |
| 1844 s = s1; |
| 1845 break; |
| 1846 } |
| 1847 e1 = n - 0x10; |
| 1848 while((n = hexdig[*++s]) !=0 && n <= 0x19) { |
| 1849 if (e1 & 0xf8000000) |
| 1850 big = 1; |
| 1851 e1 = 10*e1 + n - 0x10; |
| 1852 } |
| 1853 if (esign) |
| 1854 e1 = -e1; |
| 1855 e += e1; |
| 1856 } |
| 1857 *sp = (char*)s; |
| 1858 if (!havedig) |
| 1859 *sp = (char*)s0 - 1; |
| 1860 if (zret) |
| 1861 goto retz1; |
| 1862 if (big) { |
| 1863 if (esign) { |
| 1864 #ifdef IEEE_Arith |
| 1865 switch(rounding) { |
| 1866 case Round_up: |
| 1867 if (sign) |
| 1868 break; |
| 1869 goto ret_tiny; |
| 1870 case Round_down: |
| 1871 if (!sign) |
| 1872 break; |
| 1873 goto ret_tiny; |
| 1874 } |
| 1875 #endif |
| 1876 goto retz; |
| 1877 #ifdef IEEE_Arith |
| 1878 ret_tiny: |
| 1879 #ifndef NO_ERRNO |
| 1880 errno = ERANGE; |
| 1881 #endif |
| 1882 word0(rvp) = 0; |
| 1883 word1(rvp) = 1; |
| 1884 return; |
| 1885 #endif /* IEEE_Arith */ |
| 1886 } |
| 1887 switch(rounding) { |
| 1888 case Round_near: |
| 1889 goto ovfl1; |
| 1890 case Round_up: |
| 1891 if (!sign) |
| 1892 goto ovfl1; |
| 1893 goto ret_big; |
| 1894 case Round_down: |
| 1895 if (sign) |
| 1896 goto ovfl1; |
| 1897 goto ret_big; |
| 1898 } |
| 1899 ret_big: |
| 1900 word0(rvp) = Big0; |
| 1901 word1(rvp) = Big1; |
| 1902 return; |
| 1903 } |
| 1904 n = s1 - s0 - 1; |
| 1905 for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1) |
| 1906 k++; |
| 1907 b = Balloc(k); |
| 1908 x = b->x; |
| 1909 n = 0; |
| 1910 L = 0; |
| 1911 #ifdef USE_LOCALE |
| 1912 for(i = 0; decimalpoint[i+1]; ++i); |
| 1913 #endif |
| 1914 while(s1 > s0) { |
| 1915 #ifdef USE_LOCALE |
| 1916 if (*--s1 == decimalpoint[i]) { |
| 1917 s1 -= i; |
| 1918 continue; |
| 1919 } |
| 1920 #else |
| 1921 if (*--s1 == '.') |
| 1922 continue; |
| 1923 #endif |
| 1924 if (n == ULbits) { |
| 1925 *x++ = L; |
| 1926 L = 0; |
| 1927 n = 0; |
| 1928 } |
| 1929 L |= (hexdig[*s1] & 0x0f) << n; |
| 1930 n += 4; |
| 1931 } |
| 1932 *x++ = L; |
| 1933 b->wds = n = x - b->x; |
| 1934 n = ULbits*n - hi0bits(L); |
| 1935 nbits = Nbits; |
| 1936 lostbits = 0; |
| 1937 x = b->x; |
| 1938 if (n > nbits) { |
| 1939 n -= nbits; |
| 1940 if (any_on(b,n)) { |
| 1941 lostbits = 1; |
| 1942 k = n - 1; |
| 1943 if (x[k>>kshift] & 1 << (k & kmask)) { |
| 1944 lostbits = 2; |
| 1945 if (k > 0 && any_on(b,k)) |
| 1946 lostbits = 3; |
| 1947 } |
| 1948 } |
| 1949 rshift(b, n); |
| 1950 e += n; |
| 1951 } |
| 1952 else if (n < nbits) { |
| 1953 n = nbits - n; |
| 1954 b = lshift(b, n); |
| 1955 e -= n; |
| 1956 x = b->x; |
| 1957 } |
| 1958 if (e > Emax) { |
| 1959 ovfl: |
| 1960 Bfree(b); |
| 1961 ovfl1: |
| 1962 #ifndef NO_ERRNO |
| 1963 errno = ERANGE; |
| 1964 #endif |
| 1965 word0(rvp) = Exp_mask; |
| 1966 word1(rvp) = 0; |
| 1967 return; |
| 1968 } |
| 1969 denorm = 0; |
| 1970 if (e < emin) { |
| 1971 denorm = 1; |
| 1972 n = emin - e; |
| 1973 if (n >= nbits) { |
| 1974 #ifdef IEEE_Arith /*{*/ |
| 1975 switch (rounding) { |
| 1976 case Round_near: |
| 1977 if (n == nbits && (n < 2 || any_on(b,n-1))) |
| 1978 goto ret_tiny; |
| 1979 break; |
| 1980 case Round_up: |
| 1981 if (!sign) |
| 1982 goto ret_tiny; |
| 1983 break; |
| 1984 case Round_down: |
| 1985 if (sign) |
| 1986 goto ret_tiny; |
| 1987 } |
| 1988 #endif /* } IEEE_Arith */ |
| 1989 Bfree(b); |
| 1990 retz: |
| 1991 #ifndef NO_ERRNO |
| 1992 errno = ERANGE; |
| 1993 #endif |
| 1994 retz1: |
| 1995 rvp->d = 0.; |
| 1996 return; |
| 1997 } |
| 1998 k = n - 1; |
| 1999 if (lostbits) |
| 2000 lostbits = 1; |
| 2001 else if (k > 0) |
| 2002 lostbits = any_on(b,k); |
| 2003 if (x[k>>kshift] & 1 << (k & kmask)) |
| 2004 lostbits |= 2; |
| 2005 nbits -= n; |
| 2006 rshift(b,n); |
| 2007 e = emin; |
| 2008 } |
| 2009 if (lostbits) { |
| 2010 up = 0; |
| 2011 switch(rounding) { |
| 2012 case Round_zero: |
| 2013 break; |
| 2014 case Round_near: |
| 2015 if (lostbits & 2 |
| 2016 && (lostbits & 1) | (x[0] & 1)) |
| 2017 up = 1; |
| 2018 break; |
| 2019 case Round_up: |
| 2020 up = 1 - sign; |
| 2021 break; |
| 2022 case Round_down: |
| 2023 up = sign; |
| 2024 } |
| 2025 if (up) { |
| 2026 k = b->wds; |
| 2027 b = increment(b); |
| 2028 x = b->x; |
| 2029 if (denorm) { |
| 2030 #if 0 |
| 2031 if (nbits == Nbits - 1 |
| 2032 && x[nbits >> kshift] & 1 << (nbits & kmask)) |
| 2033 denorm = 0; /* not currently used */ |
| 2034 #endif |
| 2035 } |
| 2036 else if (b->wds > k |
| 2037 || ((n = nbits & kmask) !=0 |
| 2038 && hi0bits(x[k-1]) < 32-n)) { |
| 2039 rshift(b,1); |
| 2040 if (++e > Emax) |
| 2041 goto ovfl; |
| 2042 } |
| 2043 } |
| 2044 } |
| 2045 #ifdef IEEE_Arith |
| 2046 if (denorm) |
| 2047 word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0; |
| 2048 else |
| 2049 word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20); |
| 2050 word1(rvp) = b->x[0]; |
| 2051 #endif |
| 2052 #ifdef IBM |
| 2053 if ((j = e & 3)) { |
| 2054 k = b->x[0] & ((1 << j) - 1); |
| 2055 rshift(b,j); |
| 2056 if (k) { |
| 2057 switch(rounding) { |
| 2058 case Round_up: |
| 2059 if (!sign) |
| 2060 increment(b); |
| 2061 break; |
| 2062 case Round_down: |
| 2063 if (sign) |
| 2064 increment(b); |
| 2065 break; |
| 2066 case Round_near: |
| 2067 j = 1 << (j-1); |
| 2068 if (k & j && ((k & (j-1)) | lostbits)) |
| 2069 increment(b); |
| 2070 } |
| 2071 } |
| 2072 } |
| 2073 e >>= 2; |
| 2074 word0(rvp) = b->x[1] | ((e + 65 + 13) << 24); |
| 2075 word1(rvp) = b->x[0]; |
| 2076 #endif |
| 2077 #ifdef VAX |
| 2078 /* The next two lines ignore swap of low- and high-order 2 bytes. */ |
| 2079 /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */ |
| 2080 /* word1(rvp) = b->x[0]; */ |
| 2081 word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b-
>x[1] << 16); |
| 2082 word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16); |
| 2083 #endif |
| 2084 Bfree(b); |
| 2085 } |
| 2086 #endif /*}!NO_HEX_FP*/ |
| 2087 |
| 2088 static int |
| 2089 #ifdef KR_headers |
| 2090 dshift(b, p2) Bigint *b; int p2; |
| 2091 #else |
| 2092 dshift(Bigint *b, int p2) |
| 2093 #endif |
| 2094 { |
| 2095 int rv = hi0bits(b->x[b->wds-1]) - 4; |
| 2096 if (p2 > 0) |
| 2097 rv -= p2; |
| 2098 return rv & kmask; |
| 2099 } |
| 2100 |
| 2101 static int |
| 2102 quorem |
| 2103 #ifdef KR_headers |
| 2104 (b, S) Bigint *b, *S; |
| 2105 #else |
| 2106 (Bigint *b, Bigint *S) |
| 2107 #endif |
| 2108 { |
| 2109 int n; |
| 2110 ULong *bx, *bxe, q, *sx, *sxe; |
| 2111 #ifdef ULLong |
| 2112 ULLong borrow, carry, y, ys; |
| 2113 #else |
| 2114 ULong borrow, carry, y, ys; |
| 2115 #ifdef Pack_32 |
| 2116 ULong si, z, zs; |
| 2117 #endif |
| 2118 #endif |
| 2119 |
| 2120 n = S->wds; |
| 2121 #ifdef DEBUG |
| 2122 /*debug*/ if (b->wds > n) |
| 2123 /*debug*/ Bug("oversize b in quorem"); |
| 2124 #endif |
| 2125 if (b->wds < n) |
| 2126 return 0; |
| 2127 sx = S->x; |
| 2128 sxe = sx + --n; |
| 2129 bx = b->x; |
| 2130 bxe = bx + n; |
| 2131 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ |
| 2132 #ifdef DEBUG |
| 2133 /*debug*/ if (q > 9) |
| 2134 /*debug*/ Bug("oversized quotient in quorem"); |
| 2135 #endif |
| 2136 if (q) { |
| 2137 borrow = 0; |
| 2138 carry = 0; |
| 2139 do { |
| 2140 #ifdef ULLong |
| 2141 ys = *sx++ * (ULLong)q + carry; |
| 2142 carry = ys >> 32; |
| 2143 y = *bx - (ys & FFFFFFFF) - borrow; |
| 2144 borrow = y >> 32 & (ULong)1; |
| 2145 *bx++ = y & FFFFFFFF; |
| 2146 #else |
| 2147 #ifdef Pack_32 |
| 2148 si = *sx++; |
| 2149 ys = (si & 0xffff) * q + carry; |
| 2150 zs = (si >> 16) * q + (ys >> 16); |
| 2151 carry = zs >> 16; |
| 2152 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; |
| 2153 borrow = (y & 0x10000) >> 16; |
| 2154 z = (*bx >> 16) - (zs & 0xffff) - borrow; |
| 2155 borrow = (z & 0x10000) >> 16; |
| 2156 Storeinc(bx, z, y); |
| 2157 #else |
| 2158 ys = *sx++ * q + carry; |
| 2159 carry = ys >> 16; |
| 2160 y = *bx - (ys & 0xffff) - borrow; |
| 2161 borrow = (y & 0x10000) >> 16; |
| 2162 *bx++ = y & 0xffff; |
| 2163 #endif |
| 2164 #endif |
| 2165 } |
| 2166 while(sx <= sxe); |
| 2167 if (!*bxe) { |
| 2168 bx = b->x; |
| 2169 while(--bxe > bx && !*bxe) |
| 2170 --n; |
| 2171 b->wds = n; |
| 2172 } |
| 2173 } |
| 2174 if (cmp(b, S) >= 0) { |
| 2175 q++; |
| 2176 borrow = 0; |
| 2177 carry = 0; |
| 2178 bx = b->x; |
| 2179 sx = S->x; |
| 2180 do { |
| 2181 #ifdef ULLong |
| 2182 ys = *sx++ + carry; |
| 2183 carry = ys >> 32; |
| 2184 y = *bx - (ys & FFFFFFFF) - borrow; |
| 2185 borrow = y >> 32 & (ULong)1; |
| 2186 *bx++ = y & FFFFFFFF; |
| 2187 #else |
| 2188 #ifdef Pack_32 |
| 2189 si = *sx++; |
| 2190 ys = (si & 0xffff) + carry; |
| 2191 zs = (si >> 16) + (ys >> 16); |
| 2192 carry = zs >> 16; |
| 2193 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; |
| 2194 borrow = (y & 0x10000) >> 16; |
| 2195 z = (*bx >> 16) - (zs & 0xffff) - borrow; |
| 2196 borrow = (z & 0x10000) >> 16; |
| 2197 Storeinc(bx, z, y); |
| 2198 #else |
| 2199 ys = *sx++ + carry; |
| 2200 carry = ys >> 16; |
| 2201 y = *bx - (ys & 0xffff) - borrow; |
| 2202 borrow = (y & 0x10000) >> 16; |
| 2203 *bx++ = y & 0xffff; |
| 2204 #endif |
| 2205 #endif |
| 2206 } |
| 2207 while(sx <= sxe); |
| 2208 bx = b->x; |
| 2209 bxe = bx + n; |
| 2210 if (!*bxe) { |
| 2211 while(--bxe > bx && !*bxe) |
| 2212 --n; |
| 2213 b->wds = n; |
| 2214 } |
| 2215 } |
| 2216 return q; |
| 2217 } |
| 2218 |
| 2219 #ifndef NO_STRTOD_BIGCOMP |
| 2220 |
| 2221 static void |
| 2222 bigcomp |
| 2223 #ifdef KR_headers |
| 2224 (rv, s0, bc) |
| 2225 U *rv; CONST char *s0; BCinfo *bc; |
| 2226 #else |
| 2227 (U *rv, CONST char *s0, BCinfo *bc) |
| 2228 #endif |
| 2229 { |
| 2230 Bigint *b, *d; |
| 2231 int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase; |
| 2232 |
| 2233 dsign = bc->dsign; |
| 2234 nd = bc->nd; |
| 2235 nd0 = bc->nd0; |
| 2236 p5 = nd + bc->e0 - 1; |
| 2237 speccase = 0; |
| 2238 #ifndef Sudden_Underflow |
| 2239 if (rv->d == 0.) { /* special case: value near underflow-to-zero */ |
| 2240 /* threshold was rounded to zero */ |
| 2241 b = i2b(1); |
| 2242 p2 = Emin - P + 1; |
| 2243 bbits = 1; |
| 2244 #ifdef Avoid_Underflow |
| 2245 word0(rv) = (P+2) << Exp_shift; |
| 2246 #else |
| 2247 word1(rv) = 1; |
| 2248 #endif |
| 2249 i = 0; |
| 2250 #ifdef Honor_FLT_ROUNDS |
| 2251 if (bc->rounding == 1) |
| 2252 #endif |
| 2253 { |
| 2254 speccase = 1; |
| 2255 --p2; |
| 2256 dsign = 0; |
| 2257 goto have_i; |
| 2258 } |
| 2259 } |
| 2260 else |
| 2261 #endif |
| 2262 b = d2b(rv, &p2, &bbits); |
| 2263 #ifdef Avoid_Underflow |
| 2264 p2 -= bc->scale; |
| 2265 #endif |
| 2266 /* floor(log2(rv)) == bbits - 1 + p2 */ |
| 2267 /* Check for denormal case. */ |
| 2268 i = P - bbits; |
| 2269 if (i > (j = P - Emin - 1 + p2)) { |
| 2270 #ifdef Sudden_Underflow |
| 2271 Bfree(b); |
| 2272 b = i2b(1); |
| 2273 p2 = Emin; |
| 2274 i = P - 1; |
| 2275 #ifdef Avoid_Underflow |
| 2276 word0(rv) = (1 + bc->scale) << Exp_shift; |
| 2277 #else |
| 2278 word0(rv) = Exp_msk1; |
| 2279 #endif |
| 2280 word1(rv) = 0; |
| 2281 #else |
| 2282 i = j; |
| 2283 #endif |
| 2284 } |
| 2285 #ifdef Honor_FLT_ROUNDS |
| 2286 if (bc->rounding != 1) { |
| 2287 if (i > 0) |
| 2288 b = lshift(b, i); |
| 2289 if (dsign) |
| 2290 b = increment(b); |
| 2291 } |
| 2292 else |
| 2293 #endif |
| 2294 { |
| 2295 b = lshift(b, ++i); |
| 2296 b->x[0] |= 1; |
| 2297 } |
| 2298 #ifndef Sudden_Underflow |
| 2299 have_i: |
| 2300 #endif |
| 2301 p2 -= p5 + i; |
| 2302 d = i2b(1); |
| 2303 /* Arrange for convenient computation of quotients: |
| 2304 * shift left if necessary so divisor has 4 leading 0 bits. |
| 2305 */ |
| 2306 if (p5 > 0) |
| 2307 d = pow5mult(d, p5); |
| 2308 else if (p5 < 0) |
| 2309 b = pow5mult(b, -p5); |
| 2310 if (p2 > 0) { |
| 2311 b2 = p2; |
| 2312 d2 = 0; |
| 2313 } |
| 2314 else { |
| 2315 b2 = 0; |
| 2316 d2 = -p2; |
| 2317 } |
| 2318 i = dshift(d, d2); |
| 2319 if ((b2 += i) > 0) |
| 2320 b = lshift(b, b2); |
| 2321 if ((d2 += i) > 0) |
| 2322 d = lshift(d, d2); |
| 2323 |
| 2324 /* Now b/d = exactly half-way between the two floating-point values */ |
| 2325 /* on either side of the input string. Compute first digit of b/d. */ |
| 2326 |
| 2327 if (!(dig = quorem(b,d))) { |
| 2328 b = multadd(b, 10, 0); /* very unlikely */ |
| 2329 dig = quorem(b,d); |
| 2330 } |
| 2331 |
| 2332 /* Compare b/d with s0 */ |
| 2333 |
| 2334 for(i = 0; i < nd0; ) { |
| 2335 if ((dd = s0[i++] - '0' - dig)) |
| 2336 goto ret; |
| 2337 if (!b->x[0] && b->wds == 1) { |
| 2338 if (i < nd) |
| 2339 dd = 1; |
| 2340 goto ret; |
| 2341 } |
| 2342 b = multadd(b, 10, 0); |
| 2343 dig = quorem(b,d); |
| 2344 } |
| 2345 for(j = bc->dp1; i++ < nd;) { |
| 2346 if ((dd = s0[j++] - '0' - dig)) |
| 2347 goto ret; |
| 2348 if (!b->x[0] && b->wds == 1) { |
| 2349 if (i < nd) |
| 2350 dd = 1; |
| 2351 goto ret; |
| 2352 } |
| 2353 b = multadd(b, 10, 0); |
| 2354 dig = quorem(b,d); |
| 2355 } |
| 2356 if (b->x[0] || b->wds > 1) |
| 2357 dd = -1; |
| 2358 ret: |
| 2359 Bfree(b); |
| 2360 Bfree(d); |
| 2361 #ifdef Honor_FLT_ROUNDS |
| 2362 if (bc->rounding != 1) { |
| 2363 if (dd < 0) { |
| 2364 if (bc->rounding == 0) { |
| 2365 if (!dsign) |
| 2366 goto retlow1; |
| 2367 } |
| 2368 else if (dsign) |
| 2369 goto rethi1; |
| 2370 } |
| 2371 else if (dd > 0) { |
| 2372 if (bc->rounding == 0) { |
| 2373 if (dsign) |
| 2374 goto rethi1; |
| 2375 goto ret1; |
| 2376 } |
| 2377 if (!dsign) |
| 2378 goto rethi1; |
| 2379 dval(rv) += 2.*ulp(rv); |
| 2380 } |
| 2381 else { |
| 2382 bc->inexact = 0; |
| 2383 if (dsign) |
| 2384 goto rethi1; |
| 2385 } |
| 2386 } |
| 2387 else |
| 2388 #endif |
| 2389 if (speccase) { |
| 2390 if (dd <= 0) |
| 2391 rv->d = 0.; |
| 2392 } |
| 2393 else if (dd < 0) { |
| 2394 if (!dsign) /* does not happen for round-near */ |
| 2395 retlow1: |
| 2396 dval(rv) -= ulp(rv); |
| 2397 } |
| 2398 else if (dd > 0) { |
| 2399 if (dsign) { |
| 2400 rethi1: |
| 2401 dval(rv) += ulp(rv); |
| 2402 } |
| 2403 } |
| 2404 else { |
| 2405 /* Exact half-way case: apply round-even rule. */ |
| 2406 if (word1(rv) & 1) { |
| 2407 if (dsign) |
| 2408 goto rethi1; |
| 2409 goto retlow1; |
| 2410 } |
| 2411 } |
| 2412 |
| 2413 #ifdef Honor_FLT_ROUNDS |
| 2414 ret1: |
| 2415 #endif |
| 2416 return; |
| 2417 } |
| 2418 #endif /* NO_STRTOD_BIGCOMP */ |
| 2419 |
1547 double | 2420 double |
1548 strtod | 2421 strtod |
1549 #ifdef KR_headers | 2422 #ifdef KR_headers |
1550 (s00, se) CONST char *s00; char **se; | 2423 (s00, se) CONST char *s00; char **se; |
1551 #else | 2424 #else |
1552 (CONST char *s00, char **se) | 2425 (CONST char *s00, char **se) |
1553 #endif | 2426 #endif |
1554 { | 2427 { |
1555 #ifdef Avoid_Underflow | 2428 » int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1; |
1556 » int scale; | 2429 » int esign, i, j, k, nd, nd0, nf, nz, nz0, sign; |
1557 #endif | |
1558 » int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, | |
1559 » » e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; | |
1560 CONST char *s, *s0, *s1; | 2430 CONST char *s, *s0, *s1; |
1561 » double aadj, aadj1, adj, rv, rv0; | 2431 » double aadj, aadj1; |
1562 Long L; | 2432 Long L; |
| 2433 U aadj2, adj, rv, rv0; |
1563 ULong y, z; | 2434 ULong y, z; |
| 2435 BCinfo bc; |
1564 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; | 2436 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; |
1565 #ifdef SET_INEXACT | 2437 #ifdef SET_INEXACT |
1566 » int inexact, oldinexact; | 2438 » int oldinexact; |
1567 #endif | 2439 #endif |
1568 #ifdef Honor_FLT_ROUNDS /*{*/ | 2440 #ifdef Honor_FLT_ROUNDS /*{*/ |
1569 int Rounding; | |
1570 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ | 2441 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ |
1571 » Rounding = Flt_Rounds; | 2442 » bc.rounding = Flt_Rounds; |
1572 #else /*}{*/ | 2443 #else /*}{*/ |
1573 » Rounding = 1; | 2444 » bc.rounding = 1; |
1574 switch(fegetround()) { | 2445 switch(fegetround()) { |
1575 » case FE_TOWARDZERO:» Rounding = 0; break; | 2446 » case FE_TOWARDZERO:» bc.rounding = 0; break; |
1576 » case FE_UPWARD:» Rounding = 2; break; | 2447 » case FE_UPWARD:» bc.rounding = 2; break; |
1577 » case FE_DOWNWARD:» Rounding = 3; | 2448 » case FE_DOWNWARD:» bc.rounding = 3; |
1578 } | 2449 } |
1579 #endif /*}}*/ | 2450 #endif /*}}*/ |
1580 #endif /*}*/ | 2451 #endif /*}*/ |
1581 #ifdef USE_LOCALE | 2452 #ifdef USE_LOCALE |
1582 CONST char *s2; | 2453 CONST char *s2; |
1583 #endif | 2454 #endif |
1584 | 2455 |
1585 » sign = nz0 = nz = 0; | 2456 » sign = nz0 = nz = bc.dplen = bc.uflchk = 0; |
1586 » dval(rv) = 0.; | 2457 » dval(&rv) = 0.; |
1587 for(s = s00;;s++) switch(*s) { | 2458 for(s = s00;;s++) switch(*s) { |
1588 case '-': | 2459 case '-': |
1589 sign = 1; | 2460 sign = 1; |
1590 /* no break */ | 2461 /* no break */ |
1591 case '+': | 2462 case '+': |
1592 if (*++s) | 2463 if (*++s) |
1593 goto break2; | 2464 goto break2; |
1594 /* no break */ | 2465 /* no break */ |
1595 case 0: | 2466 case 0: |
1596 goto ret0; | 2467 goto ret0; |
1597 case '\t': | 2468 case '\t': |
1598 case '\n': | 2469 case '\n': |
1599 case '\v': | 2470 case '\v': |
1600 case '\f': | 2471 case '\f': |
1601 case '\r': | 2472 case '\r': |
1602 case ' ': | 2473 case ' ': |
1603 continue; | 2474 continue; |
1604 default: | 2475 default: |
1605 goto break2; | 2476 goto break2; |
1606 } | 2477 } |
1607 break2: | 2478 break2: |
1608 if (*s == '0') { | 2479 if (*s == '0') { |
| 2480 #ifndef NO_HEX_FP /*{*/ |
| 2481 switch(s[1]) { |
| 2482 case 'x': |
| 2483 case 'X': |
| 2484 #ifdef Honor_FLT_ROUNDS |
| 2485 gethex(&s, &rv, bc.rounding, sign); |
| 2486 #else |
| 2487 gethex(&s, &rv, 1, sign); |
| 2488 #endif |
| 2489 goto ret; |
| 2490 } |
| 2491 #endif /*}*/ |
1609 nz0 = 1; | 2492 nz0 = 1; |
1610 while(*++s == '0') ; | 2493 while(*++s == '0') ; |
1611 if (!*s) | 2494 if (!*s) |
1612 goto ret; | 2495 goto ret; |
1613 } | 2496 } |
1614 s0 = s; | 2497 s0 = s; |
1615 y = z = 0; | 2498 y = z = 0; |
1616 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) | 2499 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) |
1617 if (nd < 9) | 2500 if (nd < 9) |
1618 y = 10*y + c - '0'; | 2501 y = 10*y + c - '0'; |
1619 else if (nd < 16) | 2502 else if (nd < 16) |
1620 z = 10*z + c - '0'; | 2503 z = 10*z + c - '0'; |
1621 nd0 = nd; | 2504 nd0 = nd; |
| 2505 bc.dp0 = bc.dp1 = s - s0; |
1622 #ifdef USE_LOCALE | 2506 #ifdef USE_LOCALE |
1623 s1 = localeconv()->decimal_point; | 2507 s1 = localeconv()->decimal_point; |
1624 if (c == *s1) { | 2508 if (c == *s1) { |
1625 c = '.'; | 2509 c = '.'; |
1626 if (*++s1) { | 2510 if (*++s1) { |
1627 s2 = s; | 2511 s2 = s; |
1628 for(;;) { | 2512 for(;;) { |
1629 if (*++s2 != *s1) { | 2513 if (*++s2 != *s1) { |
1630 c = 0; | 2514 c = 0; |
1631 break; | 2515 break; |
1632 } | 2516 } |
1633 if (!*++s1) { | 2517 if (!*++s1) { |
1634 s = s2; | 2518 s = s2; |
1635 break; | 2519 break; |
1636 } | 2520 } |
1637 } | 2521 } |
1638 } | 2522 } |
1639 } | 2523 } |
1640 #endif | 2524 #endif |
1641 if (c == '.') { | 2525 if (c == '.') { |
1642 c = *++s; | 2526 c = *++s; |
| 2527 bc.dp1 = s - s0; |
| 2528 bc.dplen = bc.dp1 - bc.dp0; |
1643 if (!nd) { | 2529 if (!nd) { |
1644 for(; c == '0'; c = *++s) | 2530 for(; c == '0'; c = *++s) |
1645 nz++; | 2531 nz++; |
1646 if (c > '0' && c <= '9') { | 2532 if (c > '0' && c <= '9') { |
1647 s0 = s; | 2533 s0 = s; |
1648 nf += nz; | 2534 nf += nz; |
1649 nz = 0; | 2535 nz = 0; |
1650 goto have_dig; | 2536 goto have_dig; |
1651 } | 2537 } |
1652 goto dig_done; | 2538 goto dig_done; |
(...skipping 51 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
1704 else | 2590 else |
1705 e = 0; | 2591 e = 0; |
1706 } | 2592 } |
1707 else | 2593 else |
1708 s = s00; | 2594 s = s00; |
1709 } | 2595 } |
1710 if (!nd) { | 2596 if (!nd) { |
1711 if (!nz && !nz0) { | 2597 if (!nz && !nz0) { |
1712 #ifdef INFNAN_CHECK | 2598 #ifdef INFNAN_CHECK |
1713 /* Check for Nan and Infinity */ | 2599 /* Check for Nan and Infinity */ |
1714 » » » switch(c) { | 2600 » » » if (!bc.dplen) |
| 2601 » » » switch(c) { |
1715 case 'i': | 2602 case 'i': |
1716 case 'I': | 2603 case 'I': |
1717 if (match(&s,"nf")) { | 2604 if (match(&s,"nf")) { |
1718 --s; | 2605 --s; |
1719 if (!match(&s,"inity")) | 2606 if (!match(&s,"inity")) |
1720 ++s; | 2607 ++s; |
1721 » » » » » word0(rv) = 0x7ff00000; | 2608 » » » » » word0(&rv) = 0x7ff00000; |
1722 » » » » » word1(rv) = 0; | 2609 » » » » » word1(&rv) = 0; |
1723 goto ret; | 2610 goto ret; |
1724 } | 2611 } |
1725 break; | 2612 break; |
1726 case 'n': | 2613 case 'n': |
1727 case 'N': | 2614 case 'N': |
1728 if (match(&s, "an")) { | 2615 if (match(&s, "an")) { |
1729 » » » » » word0(rv) = NAN_WORD0; | 2616 » » » » » word0(&rv) = NAN_WORD0; |
1730 » » » » » word1(rv) = NAN_WORD1; | 2617 » » » » » word1(&rv) = NAN_WORD1; |
1731 #ifndef No_Hex_NaN | 2618 #ifndef No_Hex_NaN |
1732 if (*s == '(') /*)*/ | 2619 if (*s == '(') /*)*/ |
1733 hexnan(&rv, &s); | 2620 hexnan(&rv, &s); |
1734 #endif | 2621 #endif |
1735 goto ret; | 2622 goto ret; |
1736 } | 2623 } |
1737 } | 2624 } |
1738 #endif /* INFNAN_CHECK */ | 2625 #endif /* INFNAN_CHECK */ |
1739 ret0: | 2626 ret0: |
1740 s = s00; | 2627 s = s00; |
1741 sign = 0; | 2628 sign = 0; |
1742 } | 2629 } |
1743 goto ret; | 2630 goto ret; |
1744 } | 2631 } |
1745 » e1 = e -= nf; | 2632 » bc.e0 = e1 = e -= nf; |
1746 | 2633 |
1747 /* Now we have nd0 digits, starting at s0, followed by a | 2634 /* Now we have nd0 digits, starting at s0, followed by a |
1748 * decimal point, followed by nd-nd0 digits. The number we're | 2635 * decimal point, followed by nd-nd0 digits. The number we're |
1749 * after is the integer represented by those digits times | 2636 * after is the integer represented by those digits times |
1750 * 10**e */ | 2637 * 10**e */ |
1751 | 2638 |
1752 if (!nd0) | 2639 if (!nd0) |
1753 nd0 = nd; | 2640 nd0 = nd; |
1754 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; | 2641 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; |
1755 » dval(rv) = y; | 2642 » dval(&rv) = y; |
1756 if (k > 9) { | 2643 if (k > 9) { |
1757 #ifdef SET_INEXACT | 2644 #ifdef SET_INEXACT |
1758 if (k > DBL_DIG) | 2645 if (k > DBL_DIG) |
1759 oldinexact = get_inexact(); | 2646 oldinexact = get_inexact(); |
1760 #endif | 2647 #endif |
1761 » » dval(rv) = tens[k - 9] * dval(rv) + z; | 2648 » » dval(&rv) = tens[k - 9] * dval(&rv) + z; |
1762 } | 2649 } |
1763 bd0 = 0; | 2650 bd0 = 0; |
1764 if (nd <= DBL_DIG | 2651 if (nd <= DBL_DIG |
1765 #ifndef RND_PRODQUOT | 2652 #ifndef RND_PRODQUOT |
1766 #ifndef Honor_FLT_ROUNDS | 2653 #ifndef Honor_FLT_ROUNDS |
1767 && Flt_Rounds == 1 | 2654 && Flt_Rounds == 1 |
1768 #endif | 2655 #endif |
1769 #endif | 2656 #endif |
1770 ) { | 2657 ) { |
1771 if (!e) | 2658 if (!e) |
1772 goto ret; | 2659 goto ret; |
1773 if (e > 0) { | 2660 if (e > 0) { |
1774 if (e <= Ten_pmax) { | 2661 if (e <= Ten_pmax) { |
1775 #ifdef VAX | 2662 #ifdef VAX |
1776 goto vax_ovfl_check; | 2663 goto vax_ovfl_check; |
1777 #else | 2664 #else |
1778 #ifdef Honor_FLT_ROUNDS | 2665 #ifdef Honor_FLT_ROUNDS |
1779 /* round correctly FLT_ROUNDS = 2 or 3 */ | 2666 /* round correctly FLT_ROUNDS = 2 or 3 */ |
1780 if (sign) { | 2667 if (sign) { |
1781 » » » » » rv = -rv; | 2668 » » » » » rv.d = -rv.d; |
1782 sign = 0; | 2669 sign = 0; |
1783 } | 2670 } |
1784 #endif | 2671 #endif |
1785 » » » » /* rv = */ rounded_product(dval(rv), tens[e]); | 2672 » » » » /* rv = */ rounded_product(dval(&rv), tens[e]); |
1786 goto ret; | 2673 goto ret; |
1787 #endif | 2674 #endif |
1788 } | 2675 } |
1789 i = DBL_DIG - nd; | 2676 i = DBL_DIG - nd; |
1790 if (e <= Ten_pmax + i) { | 2677 if (e <= Ten_pmax + i) { |
1791 /* A fancier test would sometimes let us do | 2678 /* A fancier test would sometimes let us do |
1792 * this for larger i values. | 2679 * this for larger i values. |
1793 */ | 2680 */ |
1794 #ifdef Honor_FLT_ROUNDS | 2681 #ifdef Honor_FLT_ROUNDS |
1795 /* round correctly FLT_ROUNDS = 2 or 3 */ | 2682 /* round correctly FLT_ROUNDS = 2 or 3 */ |
1796 if (sign) { | 2683 if (sign) { |
1797 » » » » » rv = -rv; | 2684 » » » » » rv.d = -rv.d; |
1798 sign = 0; | 2685 sign = 0; |
1799 } | 2686 } |
1800 #endif | 2687 #endif |
1801 e -= i; | 2688 e -= i; |
1802 » » » » dval(rv) *= tens[i]; | 2689 » » » » dval(&rv) *= tens[i]; |
1803 #ifdef VAX | 2690 #ifdef VAX |
1804 /* VAX exponent range is so narrow we must | 2691 /* VAX exponent range is so narrow we must |
1805 * worry about overflow here... | 2692 * worry about overflow here... |
1806 */ | 2693 */ |
1807 vax_ovfl_check: | 2694 vax_ovfl_check: |
1808 » » » » word0(rv) -= P*Exp_msk1; | 2695 » » » » word0(&rv) -= P*Exp_msk1; |
1809 » » » » /* rv = */ rounded_product(dval(rv), tens[e]); | 2696 » » » » /* rv = */ rounded_product(dval(&rv), tens[e]); |
1810 » » » » if ((word0(rv) & Exp_mask) | 2697 » » » » if ((word0(&rv) & Exp_mask) |
1811 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) | 2698 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) |
1812 goto ovfl; | 2699 goto ovfl; |
1813 » » » » word0(rv) += P*Exp_msk1; | 2700 » » » » word0(&rv) += P*Exp_msk1; |
1814 #else | 2701 #else |
1815 » » » » /* rv = */ rounded_product(dval(rv), tens[e]); | 2702 » » » » /* rv = */ rounded_product(dval(&rv), tens[e]); |
1816 #endif | 2703 #endif |
1817 goto ret; | 2704 goto ret; |
1818 } | 2705 } |
1819 } | 2706 } |
1820 #ifndef Inaccurate_Divide | 2707 #ifndef Inaccurate_Divide |
1821 else if (e >= -Ten_pmax) { | 2708 else if (e >= -Ten_pmax) { |
1822 #ifdef Honor_FLT_ROUNDS | 2709 #ifdef Honor_FLT_ROUNDS |
1823 /* round correctly FLT_ROUNDS = 2 or 3 */ | 2710 /* round correctly FLT_ROUNDS = 2 or 3 */ |
1824 if (sign) { | 2711 if (sign) { |
1825 » » » » rv = -rv; | 2712 » » » » rv.d = -rv.d; |
1826 sign = 0; | 2713 sign = 0; |
1827 } | 2714 } |
1828 #endif | 2715 #endif |
1829 » » » /* rv = */ rounded_quotient(dval(rv), tens[-e]); | 2716 » » » /* rv = */ rounded_quotient(dval(&rv), tens[-e]); |
1830 goto ret; | 2717 goto ret; |
1831 } | 2718 } |
1832 #endif | 2719 #endif |
1833 } | 2720 } |
1834 e1 += nd - k; | 2721 e1 += nd - k; |
1835 | 2722 |
1836 #ifdef IEEE_Arith | 2723 #ifdef IEEE_Arith |
1837 #ifdef SET_INEXACT | 2724 #ifdef SET_INEXACT |
1838 » inexact = 1; | 2725 » bc.inexact = 1; |
1839 if (k <= DBL_DIG) | 2726 if (k <= DBL_DIG) |
1840 oldinexact = get_inexact(); | 2727 oldinexact = get_inexact(); |
1841 #endif | 2728 #endif |
1842 #ifdef Avoid_Underflow | 2729 #ifdef Avoid_Underflow |
1843 » scale = 0; | 2730 » bc.scale = 0; |
1844 #endif | 2731 #endif |
1845 #ifdef Honor_FLT_ROUNDS | 2732 #ifdef Honor_FLT_ROUNDS |
1846 » if (Rounding >= 2) { | 2733 » if (bc.rounding >= 2) { |
1847 if (sign) | 2734 if (sign) |
1848 » » » Rounding = Rounding == 2 ? 0 : 2; | 2735 » » » bc.rounding = bc.rounding == 2 ? 0 : 2; |
1849 else | 2736 else |
1850 » » » if (Rounding != 2) | 2737 » » » if (bc.rounding != 2) |
1851 » » » » Rounding = 0; | 2738 » » » » bc.rounding = 0; |
1852 } | 2739 } |
1853 #endif | 2740 #endif |
1854 #endif /*IEEE_Arith*/ | 2741 #endif /*IEEE_Arith*/ |
1855 | 2742 |
1856 /* Get starting approximation = rv * 10**e1 */ | 2743 /* Get starting approximation = rv * 10**e1 */ |
1857 | 2744 |
1858 if (e1 > 0) { | 2745 if (e1 > 0) { |
1859 if ((i = e1 & 15)) | 2746 if ((i = e1 & 15)) |
1860 » » » dval(rv) *= tens[i]; | 2747 » » » dval(&rv) *= tens[i]; |
1861 if (e1 &= ~15) { | 2748 if (e1 &= ~15) { |
1862 if (e1 > DBL_MAX_10_EXP) { | 2749 if (e1 > DBL_MAX_10_EXP) { |
1863 ovfl: | 2750 ovfl: |
1864 #ifndef NO_ERRNO | 2751 #ifndef NO_ERRNO |
1865 errno = ERANGE; | 2752 errno = ERANGE; |
1866 #endif | 2753 #endif |
1867 /* Can't trust HUGE_VAL */ | 2754 /* Can't trust HUGE_VAL */ |
1868 #ifdef IEEE_Arith | 2755 #ifdef IEEE_Arith |
1869 #ifdef Honor_FLT_ROUNDS | 2756 #ifdef Honor_FLT_ROUNDS |
1870 » » » » switch(Rounding) { | 2757 » » » » switch(bc.rounding) { |
1871 case 0: /* toward 0 */ | 2758 case 0: /* toward 0 */ |
1872 case 3: /* toward -infinity */ | 2759 case 3: /* toward -infinity */ |
1873 » » » » » word0(rv) = Big0; | 2760 » » » » » word0(&rv) = Big0; |
1874 » » » » » word1(rv) = Big1; | 2761 » » » » » word1(&rv) = Big1; |
1875 break; | 2762 break; |
1876 default: | 2763 default: |
1877 » » » » » word0(rv) = Exp_mask; | 2764 » » » » » word0(&rv) = Exp_mask; |
1878 » » » » » word1(rv) = 0; | 2765 » » » » » word1(&rv) = 0; |
1879 } | 2766 } |
1880 #else /*Honor_FLT_ROUNDS*/ | 2767 #else /*Honor_FLT_ROUNDS*/ |
1881 » » » » word0(rv) = Exp_mask; | 2768 » » » » word0(&rv) = Exp_mask; |
1882 » » » » word1(rv) = 0; | 2769 » » » » word1(&rv) = 0; |
1883 #endif /*Honor_FLT_ROUNDS*/ | 2770 #endif /*Honor_FLT_ROUNDS*/ |
1884 #ifdef SET_INEXACT | 2771 #ifdef SET_INEXACT |
1885 /* set overflow bit */ | 2772 /* set overflow bit */ |
1886 » » » » dval(rv0) = 1e300; | 2773 » » » » dval(&rv0) = 1e300; |
1887 » » » » dval(rv0) *= dval(rv0); | 2774 » » » » dval(&rv0) *= dval(&rv0); |
1888 #endif | 2775 #endif |
1889 #else /*IEEE_Arith*/ | 2776 #else /*IEEE_Arith*/ |
1890 » » » » word0(rv) = Big0; | 2777 » » » » word0(&rv) = Big0; |
1891 » » » » word1(rv) = Big1; | 2778 » » » » word1(&rv) = Big1; |
1892 #endif /*IEEE_Arith*/ | 2779 #endif /*IEEE_Arith*/ |
1893 if (bd0) | |
1894 goto retfree; | |
1895 goto ret; | 2780 goto ret; |
1896 } | 2781 } |
1897 e1 >>= 4; | 2782 e1 >>= 4; |
1898 for(j = 0; e1 > 1; j++, e1 >>= 1) | 2783 for(j = 0; e1 > 1; j++, e1 >>= 1) |
1899 if (e1 & 1) | 2784 if (e1 & 1) |
1900 » » » » » dval(rv) *= bigtens[j]; | 2785 » » » » » dval(&rv) *= bigtens[j]; |
1901 /* The last multiplication could overflow. */ | 2786 /* The last multiplication could overflow. */ |
1902 » » » word0(rv) -= P*Exp_msk1; | 2787 » » » word0(&rv) -= P*Exp_msk1; |
1903 » » » dval(rv) *= bigtens[j]; | 2788 » » » dval(&rv) *= bigtens[j]; |
1904 » » » if ((z = word0(rv) & Exp_mask) | 2789 » » » if ((z = word0(&rv) & Exp_mask) |
1905 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) | 2790 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) |
1906 goto ovfl; | 2791 goto ovfl; |
1907 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { | 2792 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { |
1908 /* set to largest number */ | 2793 /* set to largest number */ |
1909 /* (Can't trust DBL_MAX) */ | 2794 /* (Can't trust DBL_MAX) */ |
1910 » » » » word0(rv) = Big0; | 2795 » » » » word0(&rv) = Big0; |
1911 » » » » word1(rv) = Big1; | 2796 » » » » word1(&rv) = Big1; |
1912 } | 2797 } |
1913 else | 2798 else |
1914 » » » » word0(rv) += P*Exp_msk1; | 2799 » » » » word0(&rv) += P*Exp_msk1; |
1915 } | 2800 } |
1916 } | 2801 } |
1917 else if (e1 < 0) { | 2802 else if (e1 < 0) { |
1918 e1 = -e1; | 2803 e1 = -e1; |
1919 if ((i = e1 & 15)) | 2804 if ((i = e1 & 15)) |
1920 » » » dval(rv) /= tens[i]; | 2805 » » » dval(&rv) /= tens[i]; |
1921 if (e1 >>= 4) { | 2806 if (e1 >>= 4) { |
1922 if (e1 >= 1 << n_bigtens) | 2807 if (e1 >= 1 << n_bigtens) |
1923 goto undfl; | 2808 goto undfl; |
1924 #ifdef Avoid_Underflow | 2809 #ifdef Avoid_Underflow |
1925 if (e1 & Scale_Bit) | 2810 if (e1 & Scale_Bit) |
1926 » » » » scale = 2*P; | 2811 » » » » bc.scale = 2*P; |
1927 for(j = 0; e1 > 0; j++, e1 >>= 1) | 2812 for(j = 0; e1 > 0; j++, e1 >>= 1) |
1928 if (e1 & 1) | 2813 if (e1 & 1) |
1929 » » » » » dval(rv) *= tinytens[j]; | 2814 » » » » » dval(&rv) *= tinytens[j]; |
1930 » » » if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask) | 2815 » » » if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask) |
1931 >> Exp_shift)) > 0) { | 2816 >> Exp_shift)) > 0) { |
1932 /* scaled rv is denormal; clear j low bits */ | 2817 /* scaled rv is denormal; clear j low bits */ |
1933 if (j >= 32) { | 2818 if (j >= 32) { |
1934 » » » » » word1(rv) = 0; | 2819 » » » » » word1(&rv) = 0; |
1935 if (j >= 53) | 2820 if (j >= 53) |
1936 » » » » » word0(rv) = (P+2)*Exp_msk1; | 2821 » » » » » word0(&rv) = (P+2)*Exp_msk1; |
1937 else | 2822 else |
1938 » » » » » word0(rv) &= 0xffffffff << j-32; | 2823 » » » » » word0(&rv) &= 0xffffffff << (j-32); |
1939 } | 2824 } |
1940 else | 2825 else |
1941 » » » » » word1(rv) &= 0xffffffff << j; | 2826 » » » » » word1(&rv) &= 0xffffffff << j; |
1942 } | 2827 } |
1943 #else | 2828 #else |
1944 for(j = 0; e1 > 1; j++, e1 >>= 1) | 2829 for(j = 0; e1 > 1; j++, e1 >>= 1) |
1945 if (e1 & 1) | 2830 if (e1 & 1) |
1946 » » » » » dval(rv) *= tinytens[j]; | 2831 » » » » » dval(&rv) *= tinytens[j]; |
1947 /* The last multiplication could underflow. */ | 2832 /* The last multiplication could underflow. */ |
1948 » » » dval(rv0) = dval(rv); | 2833 » » » dval(&rv0) = dval(&rv); |
1949 » » » dval(rv) *= tinytens[j]; | 2834 » » » dval(&rv) *= tinytens[j]; |
1950 » » » if (!dval(rv)) { | 2835 » » » if (!dval(&rv)) { |
1951 » » » » dval(rv) = 2.*dval(rv0); | 2836 » » » » dval(&rv) = 2.*dval(&rv0); |
1952 » » » » dval(rv) *= tinytens[j]; | 2837 » » » » dval(&rv) *= tinytens[j]; |
1953 #endif | 2838 #endif |
1954 » » » » if (!dval(rv)) { | 2839 » » » » if (!dval(&rv)) { |
1955 undfl: | 2840 undfl: |
1956 » » » » » dval(rv) = 0.; | 2841 » » » » » dval(&rv) = 0.; |
1957 #ifndef NO_ERRNO | 2842 #ifndef NO_ERRNO |
1958 errno = ERANGE; | 2843 errno = ERANGE; |
1959 #endif | 2844 #endif |
1960 if (bd0) | |
1961 goto retfree; | |
1962 goto ret; | 2845 goto ret; |
1963 } | 2846 } |
1964 #ifndef Avoid_Underflow | 2847 #ifndef Avoid_Underflow |
1965 » » » » word0(rv) = Tiny0; | 2848 » » » » word0(&rv) = Tiny0; |
1966 » » » » word1(rv) = Tiny1; | 2849 » » » » word1(&rv) = Tiny1; |
1967 /* The refinement below will clean | 2850 /* The refinement below will clean |
1968 * this approximation up. | 2851 * this approximation up. |
1969 */ | 2852 */ |
1970 } | 2853 } |
1971 #endif | 2854 #endif |
1972 } | 2855 } |
1973 } | 2856 } |
1974 | 2857 |
1975 /* Now the hard part -- adjusting rv to the correct value.*/ | 2858 /* Now the hard part -- adjusting rv to the correct value.*/ |
1976 | 2859 |
1977 /* Put digits into bd: true value = bd * 10^e */ | 2860 /* Put digits into bd: true value = bd * 10^e */ |
1978 | 2861 |
1979 » bd0 = s2b(s0, nd0, nd, y); | 2862 » bc.nd = nd; |
| 2863 #ifndef NO_STRTOD_BIGCOMP |
| 2864 » bc.nd0 = nd0;» /* Only needed if nd > strtod_diglim, but done here */ |
| 2865 » » » /* to silence an erroneous warning about bc.nd0 */ |
| 2866 » » » /* possibly not being initialized. */ |
| 2867 » if (nd > strtod_diglim) { |
| 2868 » » /* ASSERT(strtod_diglim >= 18); 18 == one more than the */ |
| 2869 » » /* minimum number of decimal digits to distinguish double values
*/ |
| 2870 » » /* in IEEE arithmetic. */ |
| 2871 » » i = j = 18; |
| 2872 » » if (i > nd0) |
| 2873 » » » j += bc.dplen; |
| 2874 » » for(;;) { |
| 2875 » » » if (--j <= bc.dp1 && j >= bc.dp0) |
| 2876 » » » » j = bc.dp0 - 1; |
| 2877 » » » if (s0[j] != '0') |
| 2878 » » » » break; |
| 2879 » » » --i; |
| 2880 » » » } |
| 2881 » » e += nd - i; |
| 2882 » » nd = i; |
| 2883 » » if (nd0 > nd) |
| 2884 » » » nd0 = nd; |
| 2885 » » if (nd < 9) { /* must recompute y */ |
| 2886 » » » y = 0; |
| 2887 » » » for(i = 0; i < nd0; ++i) |
| 2888 » » » » y = 10*y + s0[i] - '0'; |
| 2889 » » » for(j = bc.dp1; i < nd; ++i) |
| 2890 » » » » y = 10*y + s0[j++] - '0'; |
| 2891 » » » } |
| 2892 » » } |
| 2893 #endif |
| 2894 » bd0 = s2b(s0, nd0, nd, y, bc.dplen); |
1980 | 2895 |
1981 for(;;) { | 2896 for(;;) { |
1982 bd = Balloc(bd0->k); | 2897 bd = Balloc(bd0->k); |
1983 Bcopy(bd, bd0); | 2898 Bcopy(bd, bd0); |
1984 » » bb = d2b(dval(rv), &bbe, &bbbits);» /* rv = bb * 2^bbe */ | 2899 » » bb = d2b(&rv, &bbe, &bbbits);» /* rv = bb * 2^bbe */ |
1985 bs = i2b(1); | 2900 bs = i2b(1); |
1986 | 2901 |
1987 if (e >= 0) { | 2902 if (e >= 0) { |
1988 bb2 = bb5 = 0; | 2903 bb2 = bb5 = 0; |
1989 bd2 = bd5 = e; | 2904 bd2 = bd5 = e; |
1990 } | 2905 } |
1991 else { | 2906 else { |
1992 bb2 = bb5 = -e; | 2907 bb2 = bb5 = -e; |
1993 bd2 = bd5 = 0; | 2908 bd2 = bd5 = 0; |
1994 } | 2909 } |
1995 if (bbe >= 0) | 2910 if (bbe >= 0) |
1996 bb2 += bbe; | 2911 bb2 += bbe; |
1997 else | 2912 else |
1998 bd2 -= bbe; | 2913 bd2 -= bbe; |
1999 bs2 = bb2; | 2914 bs2 = bb2; |
2000 #ifdef Honor_FLT_ROUNDS | 2915 #ifdef Honor_FLT_ROUNDS |
2001 » » if (Rounding != 1) | 2916 » » if (bc.rounding != 1) |
2002 bs2++; | 2917 bs2++; |
2003 #endif | 2918 #endif |
2004 #ifdef Avoid_Underflow | 2919 #ifdef Avoid_Underflow |
2005 » » j = bbe - scale; | 2920 » » j = bbe - bc.scale; |
2006 i = j + bbbits - 1; /* logb(rv) */ | 2921 i = j + bbbits - 1; /* logb(rv) */ |
2007 if (i < Emin) /* denormal */ | 2922 if (i < Emin) /* denormal */ |
2008 j += P - Emin; | 2923 j += P - Emin; |
2009 else | 2924 else |
2010 j = P + 1 - bbbits; | 2925 j = P + 1 - bbbits; |
2011 #else /*Avoid_Underflow*/ | 2926 #else /*Avoid_Underflow*/ |
2012 #ifdef Sudden_Underflow | 2927 #ifdef Sudden_Underflow |
2013 #ifdef IBM | 2928 #ifdef IBM |
2014 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); | 2929 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); |
2015 #else | 2930 #else |
2016 j = P + 1 - bbbits; | 2931 j = P + 1 - bbbits; |
2017 #endif | 2932 #endif |
2018 #else /*Sudden_Underflow*/ | 2933 #else /*Sudden_Underflow*/ |
2019 j = bbe; | 2934 j = bbe; |
2020 i = j + bbbits - 1; /* logb(rv) */ | 2935 i = j + bbbits - 1; /* logb(rv) */ |
2021 if (i < Emin) /* denormal */ | 2936 if (i < Emin) /* denormal */ |
2022 j += P - Emin; | 2937 j += P - Emin; |
2023 else | 2938 else |
2024 j = P + 1 - bbbits; | 2939 j = P + 1 - bbbits; |
2025 #endif /*Sudden_Underflow*/ | 2940 #endif /*Sudden_Underflow*/ |
2026 #endif /*Avoid_Underflow*/ | 2941 #endif /*Avoid_Underflow*/ |
2027 bb2 += j; | 2942 bb2 += j; |
2028 bd2 += j; | 2943 bd2 += j; |
2029 #ifdef Avoid_Underflow | 2944 #ifdef Avoid_Underflow |
2030 » » bd2 += scale; | 2945 » » bd2 += bc.scale; |
2031 #endif | 2946 #endif |
2032 i = bb2 < bd2 ? bb2 : bd2; | 2947 i = bb2 < bd2 ? bb2 : bd2; |
2033 if (i > bs2) | 2948 if (i > bs2) |
2034 i = bs2; | 2949 i = bs2; |
2035 if (i > 0) { | 2950 if (i > 0) { |
2036 bb2 -= i; | 2951 bb2 -= i; |
2037 bd2 -= i; | 2952 bd2 -= i; |
2038 bs2 -= i; | 2953 bs2 -= i; |
2039 } | 2954 } |
2040 if (bb5 > 0) { | 2955 if (bb5 > 0) { |
2041 bs = pow5mult(bs, bb5); | 2956 bs = pow5mult(bs, bb5); |
2042 bb1 = mult(bs, bb); | 2957 bb1 = mult(bs, bb); |
2043 Bfree(bb); | 2958 Bfree(bb); |
2044 bb = bb1; | 2959 bb = bb1; |
2045 } | 2960 } |
2046 if (bb2 > 0) | 2961 if (bb2 > 0) |
2047 bb = lshift(bb, bb2); | 2962 bb = lshift(bb, bb2); |
2048 if (bd5 > 0) | 2963 if (bd5 > 0) |
2049 bd = pow5mult(bd, bd5); | 2964 bd = pow5mult(bd, bd5); |
2050 if (bd2 > 0) | 2965 if (bd2 > 0) |
2051 bd = lshift(bd, bd2); | 2966 bd = lshift(bd, bd2); |
2052 if (bs2 > 0) | 2967 if (bs2 > 0) |
2053 bs = lshift(bs, bs2); | 2968 bs = lshift(bs, bs2); |
2054 delta = diff(bb, bd); | 2969 delta = diff(bb, bd); |
2055 » » dsign = delta->sign; | 2970 » » bc.dsign = delta->sign; |
2056 delta->sign = 0; | 2971 delta->sign = 0; |
2057 i = cmp(delta, bs); | 2972 i = cmp(delta, bs); |
| 2973 #ifndef NO_STRTOD_BIGCOMP |
| 2974 if (bc.nd > nd && i <= 0) { |
| 2975 if (bc.dsign) |
| 2976 break; /* Must use bigcomp(). */ |
2058 #ifdef Honor_FLT_ROUNDS | 2977 #ifdef Honor_FLT_ROUNDS |
2059 » » if (Rounding != 1) { | 2978 » » » if (bc.rounding != 1) { |
| 2979 » » » » if (i < 0) |
| 2980 » » » » » break; |
| 2981 » » » » } |
| 2982 » » » else |
| 2983 #endif |
| 2984 » » » » { |
| 2985 » » » » bc.nd = nd; |
| 2986 » » » » i = -1;»/* Discarded digits make delta smaller.
*/ |
| 2987 » » » » } |
| 2988 » » » } |
| 2989 #endif |
| 2990 #ifdef Honor_FLT_ROUNDS |
| 2991 » » if (bc.rounding != 1) { |
2060 if (i < 0) { | 2992 if (i < 0) { |
2061 /* Error is less than an ulp */ | 2993 /* Error is less than an ulp */ |
2062 if (!delta->x[0] && delta->wds <= 1) { | 2994 if (!delta->x[0] && delta->wds <= 1) { |
2063 /* exact */ | 2995 /* exact */ |
2064 #ifdef SET_INEXACT | 2996 #ifdef SET_INEXACT |
2065 » » » » » inexact = 0; | 2997 » » » » » bc.inexact = 0; |
2066 #endif | 2998 #endif |
2067 break; | 2999 break; |
2068 } | 3000 } |
2069 » » » » if (Rounding) { | 3001 » » » » if (bc.rounding) { |
2070 » » » » » if (dsign) { | 3002 » » » » » if (bc.dsign) { |
2071 » » » » » » adj = 1.; | 3003 » » » » » » adj.d = 1.; |
2072 goto apply_adj; | 3004 goto apply_adj; |
2073 } | 3005 } |
2074 } | 3006 } |
2075 » » » » else if (!dsign) { | 3007 » » » » else if (!bc.dsign) { |
2076 » » » » » adj = -1.; | 3008 » » » » » adj.d = -1.; |
2077 » » » » » if (!word1(rv) | 3009 » » » » » if (!word1(&rv) |
2078 » » » » » && !(word0(rv) & Frac_mask)) { | 3010 » » » » » && !(word0(&rv) & Frac_mask)) { |
2079 » » » » » » y = word0(rv) & Exp_mask; | 3011 » » » » » » y = word0(&rv) & Exp_mask; |
2080 #ifdef Avoid_Underflow | 3012 #ifdef Avoid_Underflow |
2081 » » » » » » if (!scale || y > 2*P*Exp_msk1) | 3013 » » » » » » if (!bc.scale || y > 2*P*Exp_msk
1) |
2082 #else | 3014 #else |
2083 if (y) | 3015 if (y) |
2084 #endif | 3016 #endif |
2085 { | 3017 { |
2086 delta = lshift(delta,Log2P); | 3018 delta = lshift(delta,Log2P); |
2087 if (cmp(delta, bs) <= 0) | 3019 if (cmp(delta, bs) <= 0) |
2088 » » » » » » » adj = -0.5; | 3020 » » » » » » » adj.d = -0.5; |
2089 } | 3021 } |
2090 } | 3022 } |
2091 apply_adj: | 3023 apply_adj: |
2092 #ifdef Avoid_Underflow | 3024 #ifdef Avoid_Underflow |
2093 » » » » » if (scale && (y = word0(rv) & Exp_mask) | 3025 » » » » » if (bc.scale && (y = word0(&rv) & Exp_ma
sk) |
2094 <= 2*P*Exp_msk1) | 3026 <= 2*P*Exp_msk1) |
2095 » » » » » word0(adj) += (2*P+1)*Exp_msk1 - y; | 3027 » » » » » word0(&adj) += (2*P+1)*Exp_msk1 - y; |
2096 #else | 3028 #else |
2097 #ifdef Sudden_Underflow | 3029 #ifdef Sudden_Underflow |
2098 » » » » » if ((word0(rv) & Exp_mask) <= | 3030 » » » » » if ((word0(&rv) & Exp_mask) <= |
2099 P*Exp_msk1) { | 3031 P*Exp_msk1) { |
2100 » » » » » » word0(rv) += P*Exp_msk1; | 3032 » » » » » » word0(&rv) += P*Exp_msk1; |
2101 » » » » » » dval(rv) += adj*ulp(dval(rv)); | 3033 » » » » » » dval(&rv) += adj.d*ulp(dval(&rv)
); |
2102 » » » » » » word0(rv) -= P*Exp_msk1; | 3034 » » » » » » word0(&rv) -= P*Exp_msk1; |
2103 } | 3035 } |
2104 else | 3036 else |
2105 #endif /*Sudden_Underflow*/ | 3037 #endif /*Sudden_Underflow*/ |
2106 #endif /*Avoid_Underflow*/ | 3038 #endif /*Avoid_Underflow*/ |
2107 » » » » » dval(rv) += adj*ulp(dval(rv)); | 3039 » » » » » dval(&rv) += adj.d*ulp(&rv); |
2108 } | 3040 } |
2109 break; | 3041 break; |
2110 } | 3042 } |
2111 » » » adj = ratio(delta, bs); | 3043 » » » adj.d = ratio(delta, bs); |
2112 » » » if (adj < 1.) | 3044 » » » if (adj.d < 1.) |
2113 » » » » adj = 1.; | 3045 » » » » adj.d = 1.; |
2114 » » » if (adj <= 0x7ffffffe) { | 3046 » » » if (adj.d <= 0x7ffffffe) { |
2115 /* adj = rounding ? ceil(adj) : floor(adj); */ | 3047 /* adj = rounding ? ceil(adj) : floor(adj); */ |
2116 » » » » y = adj; | 3048 » » » » y = adj.d; |
2117 » » » » if (y != adj) { | 3049 » » » » if (y != adj.d) { |
2118 » » » » » if (!((Rounding>>1) ^ dsign)) | 3050 » » » » » if (!((bc.rounding>>1) ^ bc.dsign)) |
2119 y++; | 3051 y++; |
2120 » » » » » adj = y; | 3052 » » » » » adj.d = y; |
2121 } | 3053 } |
2122 } | 3054 } |
2123 #ifdef Avoid_Underflow | 3055 #ifdef Avoid_Underflow |
2124 » » » if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) | 3056 » » » if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_m
sk1) |
2125 » » » » word0(adj) += (2*P+1)*Exp_msk1 - y; | 3057 » » » » word0(&adj) += (2*P+1)*Exp_msk1 - y; |
2126 #else | 3058 #else |
2127 #ifdef Sudden_Underflow | 3059 #ifdef Sudden_Underflow |
2128 » » » if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { | 3060 » » » if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { |
2129 » » » » word0(rv) += P*Exp_msk1; | 3061 » » » » word0(&rv) += P*Exp_msk1; |
2130 » » » » adj *= ulp(dval(rv)); | 3062 » » » » adj.d *= ulp(dval(&rv)); |
2131 » » » » if (dsign) | 3063 » » » » if (bc.dsign) |
2132 » » » » » dval(rv) += adj; | 3064 » » » » » dval(&rv) += adj.d; |
2133 else | 3065 else |
2134 » » » » » dval(rv) -= adj; | 3066 » » » » » dval(&rv) -= adj.d; |
2135 » » » » word0(rv) -= P*Exp_msk1; | 3067 » » » » word0(&rv) -= P*Exp_msk1; |
2136 goto cont; | 3068 goto cont; |
2137 } | 3069 } |
2138 #endif /*Sudden_Underflow*/ | 3070 #endif /*Sudden_Underflow*/ |
2139 #endif /*Avoid_Underflow*/ | 3071 #endif /*Avoid_Underflow*/ |
2140 » » » adj *= ulp(dval(rv)); | 3072 » » » adj.d *= ulp(&rv); |
2141 » » » if (dsign) { | 3073 » » » if (bc.dsign) { |
2142 » » » » if (word0(rv) == Big0 && word1(rv) == Big1) | 3074 » » » » if (word0(&rv) == Big0 && word1(&rv) == Big1) |
2143 goto ovfl; | 3075 goto ovfl; |
2144 » » » » dval(rv) += adj; | 3076 » » » » dval(&rv) += adj.d; |
2145 } | 3077 } |
2146 else | 3078 else |
2147 » » » » dval(rv) -= adj; | 3079 » » » » dval(&rv) -= adj.d; |
2148 goto cont; | 3080 goto cont; |
2149 } | 3081 } |
2150 #endif /*Honor_FLT_ROUNDS*/ | 3082 #endif /*Honor_FLT_ROUNDS*/ |
2151 | 3083 |
2152 if (i < 0) { | 3084 if (i < 0) { |
2153 /* Error is less than half an ulp -- check for | 3085 /* Error is less than half an ulp -- check for |
2154 * special case of mantissa a power of two. | 3086 * special case of mantissa a power of two. |
2155 */ | 3087 */ |
2156 » » » if (dsign || word1(rv) || word0(rv) & Bndry_mask | 3088 » » » if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask |
2157 #ifdef IEEE_Arith | 3089 #ifdef IEEE_Arith |
2158 #ifdef Avoid_Underflow | 3090 #ifdef Avoid_Underflow |
2159 » » » || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 | 3091 » » » || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 |
2160 #else | 3092 #else |
2161 » » » || (word0(rv) & Exp_mask) <= Exp_msk1 | 3093 » » » || (word0(&rv) & Exp_mask) <= Exp_msk1 |
2162 #endif | 3094 #endif |
2163 #endif | 3095 #endif |
2164 ) { | 3096 ) { |
2165 #ifdef SET_INEXACT | 3097 #ifdef SET_INEXACT |
2166 if (!delta->x[0] && delta->wds <= 1) | 3098 if (!delta->x[0] && delta->wds <= 1) |
2167 » » » » » inexact = 0; | 3099 » » » » » bc.inexact = 0; |
2168 #endif | 3100 #endif |
2169 break; | 3101 break; |
2170 } | 3102 } |
2171 if (!delta->x[0] && delta->wds <= 1) { | 3103 if (!delta->x[0] && delta->wds <= 1) { |
2172 /* exact result */ | 3104 /* exact result */ |
2173 #ifdef SET_INEXACT | 3105 #ifdef SET_INEXACT |
2174 » » » » inexact = 0; | 3106 » » » » bc.inexact = 0; |
2175 #endif | 3107 #endif |
2176 break; | 3108 break; |
2177 } | 3109 } |
2178 delta = lshift(delta,Log2P); | 3110 delta = lshift(delta,Log2P); |
2179 if (cmp(delta, bs) > 0) | 3111 if (cmp(delta, bs) > 0) |
2180 goto drop_down; | 3112 goto drop_down; |
2181 break; | 3113 break; |
2182 } | 3114 } |
2183 if (i == 0) { | 3115 if (i == 0) { |
2184 /* exactly half-way between */ | 3116 /* exactly half-way between */ |
2185 » » » if (dsign) { | 3117 » » » if (bc.dsign) { |
2186 » » » » if ((word0(rv) & Bndry_mask1) == Bndry_mask1 | 3118 » » » » if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 |
2187 » » » » && word1(rv) == ( | 3119 » » » » && word1(&rv) == ( |
2188 #ifdef Avoid_Underflow | 3120 #ifdef Avoid_Underflow |
2189 » » » (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) | 3121 » » » (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1
) |
2190 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : | 3122 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : |
2191 #endif | 3123 #endif |
2192 0xffffffff)) { | 3124 0xffffffff)) { |
2193 /*boundary case -- increment exponent*/ | 3125 /*boundary case -- increment exponent*/ |
2194 » » » » » word0(rv) = (word0(rv) & Exp_mask) | 3126 » » » » » word0(&rv) = (word0(&rv) & Exp_mask) |
2195 + Exp_msk1 | 3127 + Exp_msk1 |
2196 #ifdef IBM | 3128 #ifdef IBM |
2197 | Exp_msk1 >> 4 | 3129 | Exp_msk1 >> 4 |
2198 #endif | 3130 #endif |
2199 ; | 3131 ; |
2200 » » » » » word1(rv) = 0; | 3132 » » » » » word1(&rv) = 0; |
2201 #ifdef Avoid_Underflow | 3133 #ifdef Avoid_Underflow |
2202 » » » » » dsign = 0; | 3134 » » » » » bc.dsign = 0; |
2203 #endif | 3135 #endif |
2204 break; | 3136 break; |
2205 } | 3137 } |
2206 } | 3138 } |
2207 » » » else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { | 3139 » » » else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) { |
2208 drop_down: | 3140 drop_down: |
2209 /* boundary case -- decrement exponent */ | 3141 /* boundary case -- decrement exponent */ |
2210 #ifdef Sudden_Underflow /*{{*/ | 3142 #ifdef Sudden_Underflow /*{{*/ |
2211 » » » » L = word0(rv) & Exp_mask; | 3143 » » » » L = word0(&rv) & Exp_mask; |
2212 #ifdef IBM | 3144 #ifdef IBM |
2213 if (L < Exp_msk1) | 3145 if (L < Exp_msk1) |
2214 #else | 3146 #else |
2215 #ifdef Avoid_Underflow | 3147 #ifdef Avoid_Underflow |
2216 » » » » if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) | 3148 » » » » if (L <= (bc.scale ? (2*P+1)*Exp_msk1 : Exp_msk1
)) |
2217 #else | 3149 #else |
2218 if (L <= Exp_msk1) | 3150 if (L <= Exp_msk1) |
2219 #endif /*Avoid_Underflow*/ | 3151 #endif /*Avoid_Underflow*/ |
2220 #endif /*IBM*/ | 3152 #endif /*IBM*/ |
| 3153 { |
| 3154 if (bc.nd >nd) { |
| 3155 bc.uflchk = 1; |
| 3156 break; |
| 3157 } |
2221 goto undfl; | 3158 goto undfl; |
| 3159 } |
2222 L -= Exp_msk1; | 3160 L -= Exp_msk1; |
2223 #else /*Sudden_Underflow}{*/ | 3161 #else /*Sudden_Underflow}{*/ |
2224 #ifdef Avoid_Underflow | 3162 #ifdef Avoid_Underflow |
2225 » » » » if (scale) { | 3163 » » » » if (bc.scale) { |
2226 » » » » » L = word0(rv) & Exp_mask; | 3164 » » » » » L = word0(&rv) & Exp_mask; |
2227 if (L <= (2*P+1)*Exp_msk1) { | 3165 if (L <= (2*P+1)*Exp_msk1) { |
2228 if (L > (P+2)*Exp_msk1) | 3166 if (L > (P+2)*Exp_msk1) |
2229 /* round even ==> */ | 3167 /* round even ==> */ |
2230 /* accept rv */ | 3168 /* accept rv */ |
2231 break; | 3169 break; |
2232 /* rv = smallest denormal */ | 3170 /* rv = smallest denormal */ |
| 3171 if (bc.nd >nd) { |
| 3172 bc.uflchk = 1; |
| 3173 break; |
| 3174 } |
2233 goto undfl; | 3175 goto undfl; |
2234 } | 3176 } |
2235 } | 3177 } |
2236 #endif /*Avoid_Underflow*/ | 3178 #endif /*Avoid_Underflow*/ |
2237 » » » » L = (word0(rv) & Exp_mask) - Exp_msk1; | 3179 » » » » L = (word0(&rv) & Exp_mask) - Exp_msk1; |
2238 #endif /*Sudden_Underflow}}*/ | 3180 #endif /*Sudden_Underflow}}*/ |
2239 » » » » word0(rv) = L | Bndry_mask1; | 3181 » » » » word0(&rv) = L | Bndry_mask1; |
2240 » » » » word1(rv) = 0xffffffff; | 3182 » » » » word1(&rv) = 0xffffffff; |
2241 #ifdef IBM | 3183 #ifdef IBM |
2242 goto cont; | 3184 goto cont; |
2243 #else | 3185 #else |
2244 break; | 3186 break; |
2245 #endif | 3187 #endif |
2246 } | 3188 } |
2247 #ifndef ROUND_BIASED | 3189 #ifndef ROUND_BIASED |
2248 » » » if (!(word1(rv) & LSB)) | 3190 » » » if (!(word1(&rv) & LSB)) |
2249 break; | 3191 break; |
2250 #endif | 3192 #endif |
2251 » » » if (dsign) | 3193 » » » if (bc.dsign) |
2252 » » » » dval(rv) += ulp(dval(rv)); | 3194 » » » » dval(&rv) += ulp(&rv); |
2253 #ifndef ROUND_BIASED | 3195 #ifndef ROUND_BIASED |
2254 else { | 3196 else { |
2255 » » » » dval(rv) -= ulp(dval(rv)); | 3197 » » » » dval(&rv) -= ulp(&rv); |
2256 #ifndef Sudden_Underflow | 3198 #ifndef Sudden_Underflow |
2257 » » » » if (!dval(rv)) | 3199 » » » » if (!dval(&rv)) { |
| 3200 » » » » » if (bc.nd >nd) { |
| 3201 » » » » » » bc.uflchk = 1; |
| 3202 » » » » » » break; |
| 3203 » » » » » » } |
2258 goto undfl; | 3204 goto undfl; |
| 3205 } |
2259 #endif | 3206 #endif |
2260 } | 3207 } |
2261 #ifdef Avoid_Underflow | 3208 #ifdef Avoid_Underflow |
2262 » » » dsign = 1 - dsign; | 3209 » » » bc.dsign = 1 - bc.dsign; |
2263 #endif | 3210 #endif |
2264 #endif | 3211 #endif |
2265 break; | 3212 break; |
2266 } | 3213 } |
2267 if ((aadj = ratio(delta, bs)) <= 2.) { | 3214 if ((aadj = ratio(delta, bs)) <= 2.) { |
2268 » » » if (dsign) | 3215 » » » if (bc.dsign) |
2269 aadj = aadj1 = 1.; | 3216 aadj = aadj1 = 1.; |
2270 » » » else if (word1(rv) || word0(rv) & Bndry_mask) { | 3217 » » » else if (word1(&rv) || word0(&rv) & Bndry_mask) { |
2271 #ifndef Sudden_Underflow | 3218 #ifndef Sudden_Underflow |
2272 » » » » if (word1(rv) == Tiny1 && !word0(rv)) | 3219 » » » » if (word1(&rv) == Tiny1 && !word0(&rv)) { |
| 3220 » » » » » if (bc.nd >nd) { |
| 3221 » » » » » » bc.uflchk = 1; |
| 3222 » » » » » » break; |
| 3223 » » » » » » } |
2273 goto undfl; | 3224 goto undfl; |
| 3225 } |
2274 #endif | 3226 #endif |
2275 aadj = 1.; | 3227 aadj = 1.; |
2276 aadj1 = -1.; | 3228 aadj1 = -1.; |
2277 } | 3229 } |
2278 else { | 3230 else { |
2279 /* special case -- power of FLT_RADIX to be */ | 3231 /* special case -- power of FLT_RADIX to be */ |
2280 /* rounded down... */ | 3232 /* rounded down... */ |
2281 | 3233 |
2282 if (aadj < 2./FLT_RADIX) | 3234 if (aadj < 2./FLT_RADIX) |
2283 aadj = 1./FLT_RADIX; | 3235 aadj = 1./FLT_RADIX; |
2284 else | 3236 else |
2285 aadj *= 0.5; | 3237 aadj *= 0.5; |
2286 aadj1 = -aadj; | 3238 aadj1 = -aadj; |
2287 } | 3239 } |
2288 } | 3240 } |
2289 else { | 3241 else { |
2290 aadj *= 0.5; | 3242 aadj *= 0.5; |
2291 » » » aadj1 = dsign ? aadj : -aadj; | 3243 » » » aadj1 = bc.dsign ? aadj : -aadj; |
2292 #ifdef Check_FLT_ROUNDS | 3244 #ifdef Check_FLT_ROUNDS |
2293 » » » switch(Rounding) { | 3245 » » » switch(bc.rounding) { |
2294 case 2: /* towards +infinity */ | 3246 case 2: /* towards +infinity */ |
2295 aadj1 -= 0.5; | 3247 aadj1 -= 0.5; |
2296 break; | 3248 break; |
2297 case 0: /* towards 0 */ | 3249 case 0: /* towards 0 */ |
2298 case 3: /* towards -infinity */ | 3250 case 3: /* towards -infinity */ |
2299 aadj1 += 0.5; | 3251 aadj1 += 0.5; |
2300 } | 3252 } |
2301 #else | 3253 #else |
2302 if (Flt_Rounds == 0) | 3254 if (Flt_Rounds == 0) |
2303 aadj1 += 0.5; | 3255 aadj1 += 0.5; |
2304 #endif /*Check_FLT_ROUNDS*/ | 3256 #endif /*Check_FLT_ROUNDS*/ |
2305 } | 3257 } |
2306 » » y = word0(rv) & Exp_mask; | 3258 » » y = word0(&rv) & Exp_mask; |
2307 | 3259 |
2308 /* Check for overflow */ | 3260 /* Check for overflow */ |
2309 | 3261 |
2310 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { | 3262 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { |
2311 » » » dval(rv0) = dval(rv); | 3263 » » » dval(&rv0) = dval(&rv); |
2312 » » » word0(rv) -= P*Exp_msk1; | 3264 » » » word0(&rv) -= P*Exp_msk1; |
2313 » » » adj = aadj1 * ulp(dval(rv)); | 3265 » » » adj.d = aadj1 * ulp(&rv); |
2314 » » » dval(rv) += adj; | 3266 » » » dval(&rv) += adj.d; |
2315 » » » if ((word0(rv) & Exp_mask) >= | 3267 » » » if ((word0(&rv) & Exp_mask) >= |
2316 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { | 3268 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { |
2317 » » » » if (word0(rv0) == Big0 && word1(rv0) == Big1) | 3269 » » » » if (word0(&rv0) == Big0 && word1(&rv0) == Big1) |
2318 goto ovfl; | 3270 goto ovfl; |
2319 » » » » word0(rv) = Big0; | 3271 » » » » word0(&rv) = Big0; |
2320 » » » » word1(rv) = Big1; | 3272 » » » » word1(&rv) = Big1; |
2321 goto cont; | 3273 goto cont; |
2322 } | 3274 } |
2323 else | 3275 else |
2324 » » » » word0(rv) += P*Exp_msk1; | 3276 » » » » word0(&rv) += P*Exp_msk1; |
2325 } | 3277 } |
2326 else { | 3278 else { |
2327 #ifdef Avoid_Underflow | 3279 #ifdef Avoid_Underflow |
2328 » » » if (scale && y <= 2*P*Exp_msk1) { | 3280 » » » if (bc.scale && y <= 2*P*Exp_msk1) { |
2329 if (aadj <= 0x7fffffff) { | 3281 if (aadj <= 0x7fffffff) { |
2330 if ((z = aadj) <= 0) | 3282 if ((z = aadj) <= 0) |
2331 z = 1; | 3283 z = 1; |
2332 aadj = z; | 3284 aadj = z; |
2333 » » » » » aadj1 = dsign ? aadj : -aadj; | 3285 » » » » » aadj1 = bc.dsign ? aadj : -aadj; |
2334 } | 3286 } |
2335 » » » » word0(aadj1) += (2*P+1)*Exp_msk1 - y; | 3287 » » » » dval(&aadj2) = aadj1; |
| 3288 » » » » word0(&aadj2) += (2*P+1)*Exp_msk1 - y; |
| 3289 » » » » aadj1 = dval(&aadj2); |
2336 } | 3290 } |
2337 » » » adj = aadj1 * ulp(dval(rv)); | 3291 » » » adj.d = aadj1 * ulp(&rv); |
2338 » » » dval(rv) += adj; | 3292 » » » dval(&rv) += adj.d; |
2339 #else | 3293 #else |
2340 #ifdef Sudden_Underflow | 3294 #ifdef Sudden_Underflow |
2341 » » » if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { | 3295 » » » if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { |
2342 » » » » dval(rv0) = dval(rv); | 3296 » » » » dval(&rv0) = dval(&rv); |
2343 » » » » word0(rv) += P*Exp_msk1; | 3297 » » » » word0(&rv) += P*Exp_msk1; |
2344 » » » » adj = aadj1 * ulp(dval(rv)); | 3298 » » » » adj.d = aadj1 * ulp(&rv); |
2345 » » » » dval(rv) += adj; | 3299 » » » » dval(&rv) += adj.d; |
2346 #ifdef IBM | 3300 #ifdef IBM |
2347 » » » » if ((word0(rv) & Exp_mask) < P*Exp_msk1) | 3301 » » » » if ((word0(&rv) & Exp_mask) < P*Exp_msk1) |
2348 #else | 3302 #else |
2349 » » » » if ((word0(rv) & Exp_mask) <= P*Exp_msk1) | 3303 » » » » if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) |
2350 #endif | 3304 #endif |
2351 { | 3305 { |
2352 » » » » » if (word0(rv0) == Tiny0 | 3306 » » » » » if (word0(&rv0) == Tiny0 |
2353 » » » » » && word1(rv0) == Tiny1) | 3307 » » » » » && word1(&rv0) == Tiny1) { |
| 3308 » » » » » » if (bc.nd >nd) { |
| 3309 » » » » » » » bc.uflchk = 1; |
| 3310 » » » » » » » break; |
| 3311 » » » » » » » } |
2354 goto undfl; | 3312 goto undfl; |
2355 » » » » » word0(rv) = Tiny0; | 3313 » » » » » » } |
2356 » » » » » word1(rv) = Tiny1; | 3314 » » » » » word0(&rv) = Tiny0; |
| 3315 » » » » » word1(&rv) = Tiny1; |
2357 goto cont; | 3316 goto cont; |
2358 } | 3317 } |
2359 else | 3318 else |
2360 » » » » » word0(rv) -= P*Exp_msk1; | 3319 » » » » » word0(&rv) -= P*Exp_msk1; |
2361 } | 3320 } |
2362 else { | 3321 else { |
2363 » » » » adj = aadj1 * ulp(dval(rv)); | 3322 » » » » adj.d = aadj1 * ulp(&rv); |
2364 » » » » dval(rv) += adj; | 3323 » » » » dval(&rv) += adj.d; |
2365 } | 3324 } |
2366 #else /*Sudden_Underflow*/ | 3325 #else /*Sudden_Underflow*/ |
2367 /* Compute adj so that the IEEE rounding rules will | 3326 /* Compute adj so that the IEEE rounding rules will |
2368 * correctly round rv + adj in some half-way cases. | 3327 * correctly round rv + adj in some half-way cases. |
2369 * If rv * ulp(rv) is denormalized (i.e., | 3328 * If rv * ulp(rv) is denormalized (i.e., |
2370 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid | 3329 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid |
2371 * trouble from bits lost to denormalization; | 3330 * trouble from bits lost to denormalization; |
2372 * example: 1.2e-307 . | 3331 * example: 1.2e-307 . |
2373 */ | 3332 */ |
2374 if (y <= (P-1)*Exp_msk1 && aadj > 1.) { | 3333 if (y <= (P-1)*Exp_msk1 && aadj > 1.) { |
2375 aadj1 = (double)(int)(aadj + 0.5); | 3334 aadj1 = (double)(int)(aadj + 0.5); |
2376 » » » » if (!dsign) | 3335 » » » » if (!bc.dsign) |
2377 aadj1 = -aadj1; | 3336 aadj1 = -aadj1; |
2378 } | 3337 } |
2379 » » » adj = aadj1 * ulp(dval(rv)); | 3338 » » » adj.d = aadj1 * ulp(&rv); |
2380 » » » dval(rv) += adj; | 3339 » » » dval(&rv) += adj.d; |
2381 #endif /*Sudden_Underflow*/ | 3340 #endif /*Sudden_Underflow*/ |
2382 #endif /*Avoid_Underflow*/ | 3341 #endif /*Avoid_Underflow*/ |
2383 } | 3342 } |
2384 » » z = word0(rv) & Exp_mask; | 3343 » » z = word0(&rv) & Exp_mask; |
2385 #ifndef SET_INEXACT | 3344 #ifndef SET_INEXACT |
| 3345 if (bc.nd == nd) { |
2386 #ifdef Avoid_Underflow | 3346 #ifdef Avoid_Underflow |
2387 » » if (!scale) | 3347 » » if (!bc.scale) |
2388 #endif | 3348 #endif |
2389 if (y == z) { | 3349 if (y == z) { |
2390 /* Can we stop now? */ | 3350 /* Can we stop now? */ |
2391 L = (Long)aadj; | 3351 L = (Long)aadj; |
2392 aadj -= L; | 3352 aadj -= L; |
2393 /* The tolerances below are conservative. */ | 3353 /* The tolerances below are conservative. */ |
2394 » » » if (dsign || word1(rv) || word0(rv) & Bndry_mask) { | 3354 » » » if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) { |
2395 if (aadj < .4999999 || aadj > .5000001) | 3355 if (aadj < .4999999 || aadj > .5000001) |
2396 break; | 3356 break; |
2397 } | 3357 } |
2398 else if (aadj < .4999999/FLT_RADIX) | 3358 else if (aadj < .4999999/FLT_RADIX) |
2399 break; | 3359 break; |
2400 } | 3360 } |
| 3361 } |
2401 #endif | 3362 #endif |
2402 cont: | 3363 cont: |
2403 Bfree(bb); | 3364 Bfree(bb); |
2404 Bfree(bd); | 3365 Bfree(bd); |
2405 Bfree(bs); | 3366 Bfree(bs); |
2406 Bfree(delta); | 3367 Bfree(delta); |
2407 } | 3368 } |
| 3369 Bfree(bb); |
| 3370 Bfree(bd); |
| 3371 Bfree(bs); |
| 3372 Bfree(bd0); |
| 3373 Bfree(delta); |
| 3374 #ifndef NO_STRTOD_BIGCOMP |
| 3375 if (bc.nd > nd) |
| 3376 bigcomp(&rv, s0, &bc); |
| 3377 #endif |
2408 #ifdef SET_INEXACT | 3378 #ifdef SET_INEXACT |
2409 » if (inexact) { | 3379 » if (bc.inexact) { |
2410 if (!oldinexact) { | 3380 if (!oldinexact) { |
2411 » » » word0(rv0) = Exp_1 + (70 << Exp_shift); | 3381 » » » word0(&rv0) = Exp_1 + (70 << Exp_shift); |
2412 » » » word1(rv0) = 0; | 3382 » » » word1(&rv0) = 0; |
2413 » » » dval(rv0) += 1.; | 3383 » » » dval(&rv0) += 1.; |
2414 } | 3384 } |
2415 } | 3385 } |
2416 else if (!oldinexact) | 3386 else if (!oldinexact) |
2417 clear_inexact(); | 3387 clear_inexact(); |
2418 #endif | 3388 #endif |
2419 #ifdef Avoid_Underflow | 3389 #ifdef Avoid_Underflow |
2420 » if (scale) { | 3390 » if (bc.scale) { |
2421 » » word0(rv0) = Exp_1 - 2*P*Exp_msk1; | 3391 » » word0(&rv0) = Exp_1 - 2*P*Exp_msk1; |
2422 » » word1(rv0) = 0; | 3392 » » word1(&rv0) = 0; |
2423 » » dval(rv) *= dval(rv0); | 3393 » » dval(&rv) *= dval(&rv0); |
2424 #ifndef NO_ERRNO | 3394 #ifndef NO_ERRNO |
2425 /* try to avoid the bug of testing an 8087 register value */ | 3395 /* try to avoid the bug of testing an 8087 register value */ |
2426 #ifdef IEEE_Arith | 3396 #ifdef IEEE_Arith |
2427 » » if (!(word0(rv) & Exp_mask)) | 3397 » » if (!(word0(&rv) & Exp_mask)) |
2428 #else | 3398 #else |
2429 » » if (word0(rv) == 0 && word1(rv) == 0) | 3399 » » if (word0(&rv) == 0 && word1(&rv) == 0) |
2430 #endif | 3400 #endif |
2431 errno = ERANGE; | 3401 errno = ERANGE; |
2432 #endif | 3402 #endif |
2433 } | 3403 } |
2434 #endif /* Avoid_Underflow */ | 3404 #endif /* Avoid_Underflow */ |
2435 #ifdef SET_INEXACT | 3405 #ifdef SET_INEXACT |
2436 » if (inexact && !(word0(rv) & Exp_mask)) { | 3406 » if (bc.inexact && !(word0(&rv) & Exp_mask)) { |
2437 /* set underflow bit */ | 3407 /* set underflow bit */ |
2438 » » dval(rv0) = 1e-300; | 3408 » » dval(&rv0) = 1e-300; |
2439 » » dval(rv0) *= dval(rv0); | 3409 » » dval(&rv0) *= dval(&rv0); |
2440 } | 3410 } |
2441 #endif | 3411 #endif |
2442 retfree: | |
2443 Bfree(bb); | |
2444 Bfree(bd); | |
2445 Bfree(bs); | |
2446 Bfree(bd0); | |
2447 Bfree(delta); | |
2448 ret: | 3412 ret: |
2449 if (se) | 3413 if (se) |
2450 *se = (char *)s; | 3414 *se = (char *)s; |
2451 » return sign ? -dval(rv) : dval(rv); | 3415 » return sign ? -dval(&rv) : dval(&rv); |
2452 » } | |
2453 | |
2454 static int | |
2455 quorem | |
2456 #ifdef KR_headers | |
2457 » (b, S) Bigint *b, *S; | |
2458 #else | |
2459 » (Bigint *b, Bigint *S) | |
2460 #endif | |
2461 { | |
2462 » int n; | |
2463 » ULong *bx, *bxe, q, *sx, *sxe; | |
2464 #ifdef ULLong | |
2465 » ULLong borrow, carry, y, ys; | |
2466 #else | |
2467 » ULong borrow, carry, y, ys; | |
2468 #ifdef Pack_32 | |
2469 » ULong si, z, zs; | |
2470 #endif | |
2471 #endif | |
2472 | |
2473 » n = S->wds; | |
2474 #ifdef DEBUG | |
2475 » /*debug*/ if (b->wds > n) | |
2476 » /*debug*/» Bug("oversize b in quorem"); | |
2477 #endif | |
2478 » if (b->wds < n) | |
2479 » » return 0; | |
2480 » sx = S->x; | |
2481 » sxe = sx + --n; | |
2482 » bx = b->x; | |
2483 » bxe = bx + n; | |
2484 » q = *bxe / (*sxe + 1);» /* ensure q <= true quotient */ | |
2485 #ifdef DEBUG | |
2486 » /*debug*/ if (q > 9) | |
2487 » /*debug*/» Bug("oversized quotient in quorem"); | |
2488 #endif | |
2489 » if (q) { | |
2490 » » borrow = 0; | |
2491 » » carry = 0; | |
2492 » » do { | |
2493 #ifdef ULLong | |
2494 » » » ys = *sx++ * (ULLong)q + carry; | |
2495 » » » carry = ys >> 32; | |
2496 » » » y = *bx - (ys & FFFFFFFF) - borrow; | |
2497 » » » borrow = y >> 32 & (ULong)1; | |
2498 » » » *bx++ = y & FFFFFFFF; | |
2499 #else | |
2500 #ifdef Pack_32 | |
2501 » » » si = *sx++; | |
2502 » » » ys = (si & 0xffff) * q + carry; | |
2503 » » » zs = (si >> 16) * q + (ys >> 16); | |
2504 » » » carry = zs >> 16; | |
2505 » » » y = (*bx & 0xffff) - (ys & 0xffff) - borrow; | |
2506 » » » borrow = (y & 0x10000) >> 16; | |
2507 » » » z = (*bx >> 16) - (zs & 0xffff) - borrow; | |
2508 » » » borrow = (z & 0x10000) >> 16; | |
2509 » » » Storeinc(bx, z, y); | |
2510 #else | |
2511 » » » ys = *sx++ * q + carry; | |
2512 » » » carry = ys >> 16; | |
2513 » » » y = *bx - (ys & 0xffff) - borrow; | |
2514 » » » borrow = (y & 0x10000) >> 16; | |
2515 » » » *bx++ = y & 0xffff; | |
2516 #endif | |
2517 #endif | |
2518 » » » } | |
2519 » » » while(sx <= sxe); | |
2520 » » if (!*bxe) { | |
2521 » » » bx = b->x; | |
2522 » » » while(--bxe > bx && !*bxe) | |
2523 » » » » --n; | |
2524 » » » b->wds = n; | |
2525 » » » } | |
2526 » » } | |
2527 » if (cmp(b, S) >= 0) { | |
2528 » » q++; | |
2529 » » borrow = 0; | |
2530 » » carry = 0; | |
2531 » » bx = b->x; | |
2532 » » sx = S->x; | |
2533 » » do { | |
2534 #ifdef ULLong | |
2535 » » » ys = *sx++ + carry; | |
2536 » » » carry = ys >> 32; | |
2537 » » » y = *bx - (ys & FFFFFFFF) - borrow; | |
2538 » » » borrow = y >> 32 & (ULong)1; | |
2539 » » » *bx++ = y & FFFFFFFF; | |
2540 #else | |
2541 #ifdef Pack_32 | |
2542 » » » si = *sx++; | |
2543 » » » ys = (si & 0xffff) + carry; | |
2544 » » » zs = (si >> 16) + (ys >> 16); | |
2545 » » » carry = zs >> 16; | |
2546 » » » y = (*bx & 0xffff) - (ys & 0xffff) - borrow; | |
2547 » » » borrow = (y & 0x10000) >> 16; | |
2548 » » » z = (*bx >> 16) - (zs & 0xffff) - borrow; | |
2549 » » » borrow = (z & 0x10000) >> 16; | |
2550 » » » Storeinc(bx, z, y); | |
2551 #else | |
2552 » » » ys = *sx++ + carry; | |
2553 » » » carry = ys >> 16; | |
2554 » » » y = *bx - (ys & 0xffff) - borrow; | |
2555 » » » borrow = (y & 0x10000) >> 16; | |
2556 » » » *bx++ = y & 0xffff; | |
2557 #endif | |
2558 #endif | |
2559 » » » } | |
2560 » » » while(sx <= sxe); | |
2561 » » bx = b->x; | |
2562 » » bxe = bx + n; | |
2563 » » if (!*bxe) { | |
2564 » » » while(--bxe > bx && !*bxe) | |
2565 » » » » --n; | |
2566 » » » b->wds = n; | |
2567 » » » } | |
2568 » » } | |
2569 » return q; | |
2570 } | 3416 } |
2571 | 3417 |
2572 #ifndef MULTIPLE_THREADS | 3418 #ifndef MULTIPLE_THREADS |
2573 static char *dtoa_result; | 3419 static char *dtoa_result; |
2574 #endif | 3420 #endif |
2575 | 3421 |
2576 static char * | 3422 static char * |
2577 #ifdef KR_headers | 3423 #ifdef KR_headers |
2578 rv_alloc(i) int i; | 3424 rv_alloc(i) int i; |
2579 #else | 3425 #else |
(...skipping 84 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2664 * guarantee that the floating-point calculation has given | 3510 * guarantee that the floating-point calculation has given |
2665 * the correctly rounded result. For k requested digits and | 3511 * the correctly rounded result. For k requested digits and |
2666 * "uniformly" distributed input, the probability is | 3512 * "uniformly" distributed input, the probability is |
2667 * something like 10^(k-15) that we must resort to the Long | 3513 * something like 10^(k-15) that we must resort to the Long |
2668 * calculation. | 3514 * calculation. |
2669 */ | 3515 */ |
2670 | 3516 |
2671 char * | 3517 char * |
2672 dtoa | 3518 dtoa |
2673 #ifdef KR_headers | 3519 #ifdef KR_headers |
2674 » (d, mode, ndigits, decpt, sign, rve) | 3520 » (dd, mode, ndigits, decpt, sign, rve) |
2675 » double d; int mode, ndigits, *decpt, *sign; char **rve; | 3521 » double dd; int mode, ndigits, *decpt, *sign; char **rve; |
2676 #else | 3522 #else |
2677 » (double d, int mode, int ndigits, int *decpt, int *sign, char **rve) | 3523 » (double dd, int mode, int ndigits, int *decpt, int *sign, char **rve) |
2678 #endif | 3524 #endif |
2679 { | 3525 { |
2680 /* Arguments ndigits, decpt, sign are similar to those | 3526 /* Arguments ndigits, decpt, sign are similar to those |
2681 of ecvt and fcvt; trailing zeros are suppressed from | 3527 of ecvt and fcvt; trailing zeros are suppressed from |
2682 the returned string. If not null, *rve is set to point | 3528 the returned string. If not null, *rve is set to point |
2683 to the end of the return value. If d is +-Infinity or NaN, | 3529 to the end of the return value. If d is +-Infinity or NaN, |
2684 then *decpt is set to 9999. | 3530 then *decpt is set to 9999. |
2685 | 3531 |
2686 mode: | 3532 mode: |
2687 0 ==> shortest string that yields d when read in | 3533 0 ==> shortest string that yields d when read in |
(...skipping 16 matching lines...) Expand all Loading... |
2704 as modes 2 and 3 when FLT_ROUNDS != 1. | 3550 as modes 2 and 3 when FLT_ROUNDS != 1. |
2705 6-9 ==> Debugging modes similar to mode - 4: don't try | 3551 6-9 ==> Debugging modes similar to mode - 4: don't try |
2706 fast floating-point estimate (if applicable). | 3552 fast floating-point estimate (if applicable). |
2707 | 3553 |
2708 Values of mode other than 0-9 are treated as mode 0. | 3554 Values of mode other than 0-9 are treated as mode 0. |
2709 | 3555 |
2710 Sufficient space is allocated to the return value | 3556 Sufficient space is allocated to the return value |
2711 to hold the suppressed trailing zeros. | 3557 to hold the suppressed trailing zeros. |
2712 */ | 3558 */ |
2713 | 3559 |
2714 » int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0, | 3560 » int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, |
2715 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, | 3561 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, |
2716 spec_case, try_quick; | 3562 spec_case, try_quick; |
2717 Long L; | 3563 Long L; |
2718 #ifndef Sudden_Underflow | 3564 #ifndef Sudden_Underflow |
2719 int denorm; | 3565 int denorm; |
2720 ULong x; | 3566 ULong x; |
2721 #endif | 3567 #endif |
2722 Bigint *b, *b1, *delta, *mlo, *mhi, *S; | 3568 Bigint *b, *b1, *delta, *mlo, *mhi, *S; |
2723 » double d2, ds, eps; | 3569 » U d2, eps, u; |
| 3570 » double ds; |
2724 char *s, *s0; | 3571 char *s, *s0; |
2725 #ifdef SET_INEXACT | 3572 #ifdef SET_INEXACT |
2726 int inexact, oldinexact; | 3573 int inexact, oldinexact; |
2727 #endif | 3574 #endif |
2728 #ifdef Honor_FLT_ROUNDS /*{*/ | 3575 #ifdef Honor_FLT_ROUNDS /*{*/ |
2729 int Rounding; | 3576 int Rounding; |
2730 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ | 3577 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ |
2731 Rounding = Flt_Rounds; | 3578 Rounding = Flt_Rounds; |
2732 #else /*}{*/ | 3579 #else /*}{*/ |
2733 Rounding = 1; | 3580 Rounding = 1; |
2734 switch(fegetround()) { | 3581 switch(fegetround()) { |
2735 case FE_TOWARDZERO: Rounding = 0; break; | 3582 case FE_TOWARDZERO: Rounding = 0; break; |
2736 case FE_UPWARD: Rounding = 2; break; | 3583 case FE_UPWARD: Rounding = 2; break; |
2737 case FE_DOWNWARD: Rounding = 3; | 3584 case FE_DOWNWARD: Rounding = 3; |
2738 } | 3585 } |
2739 #endif /*}}*/ | 3586 #endif /*}}*/ |
2740 #endif /*}*/ | 3587 #endif /*}*/ |
2741 | 3588 |
2742 #ifndef MULTIPLE_THREADS | 3589 #ifndef MULTIPLE_THREADS |
2743 if (dtoa_result) { | 3590 if (dtoa_result) { |
2744 freedtoa(dtoa_result); | 3591 freedtoa(dtoa_result); |
2745 dtoa_result = 0; | 3592 dtoa_result = 0; |
2746 } | 3593 } |
2747 #endif | 3594 #endif |
2748 | 3595 |
2749 » if (word0(d) & Sign_bit) { | 3596 » u.d = dd; |
| 3597 » if (word0(&u) & Sign_bit) { |
2750 /* set sign for everything, including 0's and NaNs */ | 3598 /* set sign for everything, including 0's and NaNs */ |
2751 *sign = 1; | 3599 *sign = 1; |
2752 » » word0(d) &= ~Sign_bit;» /* clear sign bit */ | 3600 » » word0(&u) &= ~Sign_bit;»/* clear sign bit */ |
2753 } | 3601 } |
2754 else | 3602 else |
2755 *sign = 0; | 3603 *sign = 0; |
2756 | 3604 |
2757 #if defined(IEEE_Arith) + defined(VAX) | 3605 #if defined(IEEE_Arith) + defined(VAX) |
2758 #ifdef IEEE_Arith | 3606 #ifdef IEEE_Arith |
2759 » if ((word0(d) & Exp_mask) == Exp_mask) | 3607 » if ((word0(&u) & Exp_mask) == Exp_mask) |
2760 #else | 3608 #else |
2761 » if (word0(d) == 0x8000) | 3609 » if (word0(&u) == 0x8000) |
2762 #endif | 3610 #endif |
2763 { | 3611 { |
2764 /* Infinity or NaN */ | 3612 /* Infinity or NaN */ |
2765 *decpt = 9999; | 3613 *decpt = 9999; |
2766 #ifdef IEEE_Arith | 3614 #ifdef IEEE_Arith |
2767 » » if (!word1(d) && !(word0(d) & 0xfffff)) | 3615 » » if (!word1(&u) && !(word0(&u) & 0xfffff)) |
2768 return nrv_alloc("Infinity", rve, 8); | 3616 return nrv_alloc("Infinity", rve, 8); |
2769 #endif | 3617 #endif |
2770 return nrv_alloc("NaN", rve, 3); | 3618 return nrv_alloc("NaN", rve, 3); |
2771 } | 3619 } |
2772 #endif | 3620 #endif |
2773 #ifdef IBM | 3621 #ifdef IBM |
2774 » dval(d) += 0; /* normalize */ | 3622 » dval(&u) += 0; /* normalize */ |
2775 #endif | 3623 #endif |
2776 » if (!dval(d)) { | 3624 » if (!dval(&u)) { |
2777 *decpt = 1; | 3625 *decpt = 1; |
2778 return nrv_alloc("0", rve, 1); | 3626 return nrv_alloc("0", rve, 1); |
2779 } | 3627 } |
2780 | 3628 |
2781 #ifdef SET_INEXACT | 3629 #ifdef SET_INEXACT |
2782 try_quick = oldinexact = get_inexact(); | 3630 try_quick = oldinexact = get_inexact(); |
2783 inexact = 1; | 3631 inexact = 1; |
2784 #endif | 3632 #endif |
2785 #ifdef Honor_FLT_ROUNDS | 3633 #ifdef Honor_FLT_ROUNDS |
2786 if (Rounding >= 2) { | 3634 if (Rounding >= 2) { |
2787 if (*sign) | 3635 if (*sign) |
2788 Rounding = Rounding == 2 ? 0 : 2; | 3636 Rounding = Rounding == 2 ? 0 : 2; |
2789 else | 3637 else |
2790 if (Rounding != 2) | 3638 if (Rounding != 2) |
2791 Rounding = 0; | 3639 Rounding = 0; |
2792 } | 3640 } |
2793 #endif | 3641 #endif |
2794 | 3642 |
2795 » b = d2b(dval(d), &be, &bbits); | 3643 » b = d2b(&u, &be, &bbits); |
2796 #ifdef Sudden_Underflow | 3644 #ifdef Sudden_Underflow |
2797 » i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); | 3645 » i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); |
2798 #else | 3646 #else |
2799 » if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) { | 3647 » if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) { |
2800 #endif | 3648 #endif |
2801 » » dval(d2) = dval(d); | 3649 » » dval(&d2) = dval(&u); |
2802 » » word0(d2) &= Frac_mask1; | 3650 » » word0(&d2) &= Frac_mask1; |
2803 » » word0(d2) |= Exp_11; | 3651 » » word0(&d2) |= Exp_11; |
2804 #ifdef IBM | 3652 #ifdef IBM |
2805 » » if (j = 11 - hi0bits(word0(d2) & Frac_mask)) | 3653 » » if (j = 11 - hi0bits(word0(&d2) & Frac_mask)) |
2806 » » » dval(d2) /= 1 << j; | 3654 » » » dval(&d2) /= 1 << j; |
2807 #endif | 3655 #endif |
2808 | 3656 |
2809 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 | 3657 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 |
2810 * log10(x) = log(x) / log(10) | 3658 * log10(x) = log(x) / log(10) |
2811 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) | 3659 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) |
2812 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) | 3660 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) |
2813 * | 3661 * |
2814 * This suggests computing an approximation k to log10(d) by | 3662 * This suggests computing an approximation k to log10(d) by |
2815 * | 3663 * |
2816 * k = (i - Bias)*0.301029995663981 | 3664 * k = (i - Bias)*0.301029995663981 |
(...skipping 16 matching lines...) Expand all Loading... |
2833 i <<= 2; | 3681 i <<= 2; |
2834 i += j; | 3682 i += j; |
2835 #endif | 3683 #endif |
2836 #ifndef Sudden_Underflow | 3684 #ifndef Sudden_Underflow |
2837 denorm = 0; | 3685 denorm = 0; |
2838 } | 3686 } |
2839 else { | 3687 else { |
2840 /* d is denormalized */ | 3688 /* d is denormalized */ |
2841 | 3689 |
2842 i = bbits + be + (Bias + (P-1) - 1); | 3690 i = bbits + be + (Bias + (P-1) - 1); |
2843 » » x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32 | 3691 » » x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32) |
2844 » » » : word1(d) << 32 - i; | 3692 » » » : word1(&u) << (32 - i); |
2845 » » dval(d2) = x; | 3693 » » dval(&d2) = x; |
2846 » » word0(d2) -= 31*Exp_msk1; /* adjust exponent */ | 3694 » » word0(&d2) -= 31*Exp_msk1; /* adjust exponent */ |
2847 i -= (Bias + (P-1) - 1) + 1; | 3695 i -= (Bias + (P-1) - 1) + 1; |
2848 denorm = 1; | 3696 denorm = 1; |
2849 } | 3697 } |
2850 #endif | 3698 #endif |
2851 » ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.3010299956
63981; | 3699 » ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995
663981; |
2852 k = (int)ds; | 3700 k = (int)ds; |
2853 if (ds < 0. && ds != k) | 3701 if (ds < 0. && ds != k) |
2854 k--; /* want k = floor(ds) */ | 3702 k--; /* want k = floor(ds) */ |
2855 k_check = 1; | 3703 k_check = 1; |
2856 if (k >= 0 && k <= Ten_pmax) { | 3704 if (k >= 0 && k <= Ten_pmax) { |
2857 » » if (dval(d) < tens[k]) | 3705 » » if (dval(&u) < tens[k]) |
2858 k--; | 3706 k--; |
2859 k_check = 0; | 3707 k_check = 0; |
2860 } | 3708 } |
2861 j = bbits - i - 1; | 3709 j = bbits - i - 1; |
2862 if (j >= 0) { | 3710 if (j >= 0) { |
2863 b2 = 0; | 3711 b2 = 0; |
2864 s2 = j; | 3712 s2 = j; |
2865 } | 3713 } |
2866 else { | 3714 else { |
2867 b2 = -j; | 3715 b2 = -j; |
(...skipping 18 matching lines...) Expand all Loading... |
2886 #else | 3734 #else |
2887 try_quick = 1; | 3735 try_quick = 1; |
2888 #endif | 3736 #endif |
2889 #endif /*SET_INEXACT*/ | 3737 #endif /*SET_INEXACT*/ |
2890 | 3738 |
2891 if (mode > 5) { | 3739 if (mode > 5) { |
2892 mode -= 4; | 3740 mode -= 4; |
2893 try_quick = 0; | 3741 try_quick = 0; |
2894 } | 3742 } |
2895 leftright = 1; | 3743 leftright = 1; |
| 3744 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */ |
| 3745 /* silence erroneous "gcc -Wall" warning. */ |
2896 switch(mode) { | 3746 switch(mode) { |
2897 case 0: | 3747 case 0: |
2898 case 1: | 3748 case 1: |
2899 ilim = ilim1 = -1; | |
2900 i = 18; | 3749 i = 18; |
2901 ndigits = 0; | 3750 ndigits = 0; |
2902 break; | 3751 break; |
2903 case 2: | 3752 case 2: |
2904 leftright = 0; | 3753 leftright = 0; |
2905 /* no break */ | 3754 /* no break */ |
2906 case 4: | 3755 case 4: |
2907 if (ndigits <= 0) | 3756 if (ndigits <= 0) |
2908 ndigits = 1; | 3757 ndigits = 1; |
2909 ilim = ilim1 = i = ndigits; | 3758 ilim = ilim1 = i = ndigits; |
(...skipping 13 matching lines...) Expand all Loading... |
2923 #ifdef Honor_FLT_ROUNDS | 3772 #ifdef Honor_FLT_ROUNDS |
2924 if (mode > 1 && Rounding != 1) | 3773 if (mode > 1 && Rounding != 1) |
2925 leftright = 0; | 3774 leftright = 0; |
2926 #endif | 3775 #endif |
2927 | 3776 |
2928 if (ilim >= 0 && ilim <= Quick_max && try_quick) { | 3777 if (ilim >= 0 && ilim <= Quick_max && try_quick) { |
2929 | 3778 |
2930 /* Try to get by with floating-point arithmetic. */ | 3779 /* Try to get by with floating-point arithmetic. */ |
2931 | 3780 |
2932 i = 0; | 3781 i = 0; |
2933 » » dval(d2) = dval(d); | 3782 » » dval(&d2) = dval(&u); |
2934 k0 = k; | 3783 k0 = k; |
2935 ilim0 = ilim; | 3784 ilim0 = ilim; |
2936 ieps = 2; /* conservative */ | 3785 ieps = 2; /* conservative */ |
2937 if (k > 0) { | 3786 if (k > 0) { |
2938 ds = tens[k&0xf]; | 3787 ds = tens[k&0xf]; |
2939 j = k >> 4; | 3788 j = k >> 4; |
2940 if (j & Bletch) { | 3789 if (j & Bletch) { |
2941 /* prevent overflows */ | 3790 /* prevent overflows */ |
2942 j &= Bletch - 1; | 3791 j &= Bletch - 1; |
2943 » » » » dval(d) /= bigtens[n_bigtens-1]; | 3792 » » » » dval(&u) /= bigtens[n_bigtens-1]; |
2944 ieps++; | 3793 ieps++; |
2945 } | 3794 } |
2946 for(; j; j >>= 1, i++) | 3795 for(; j; j >>= 1, i++) |
2947 if (j & 1) { | 3796 if (j & 1) { |
2948 ieps++; | 3797 ieps++; |
2949 ds *= bigtens[i]; | 3798 ds *= bigtens[i]; |
2950 } | 3799 } |
2951 » » » dval(d) /= ds; | 3800 » » » dval(&u) /= ds; |
2952 } | 3801 } |
2953 else if ((j1 = -k)) { | 3802 else if ((j1 = -k)) { |
2954 » » » dval(d) *= tens[j1 & 0xf]; | 3803 » » » dval(&u) *= tens[j1 & 0xf]; |
2955 for(j = j1 >> 4; j; j >>= 1, i++) | 3804 for(j = j1 >> 4; j; j >>= 1, i++) |
2956 if (j & 1) { | 3805 if (j & 1) { |
2957 ieps++; | 3806 ieps++; |
2958 » » » » » dval(d) *= bigtens[i]; | 3807 » » » » » dval(&u) *= bigtens[i]; |
2959 } | 3808 } |
2960 } | 3809 } |
2961 » » if (k_check && dval(d) < 1. && ilim > 0) { | 3810 » » if (k_check && dval(&u) < 1. && ilim > 0) { |
2962 if (ilim1 <= 0) | 3811 if (ilim1 <= 0) |
2963 goto fast_failed; | 3812 goto fast_failed; |
2964 ilim = ilim1; | 3813 ilim = ilim1; |
2965 k--; | 3814 k--; |
2966 » » » dval(d) *= 10.; | 3815 » » » dval(&u) *= 10.; |
2967 ieps++; | 3816 ieps++; |
2968 } | 3817 } |
2969 » » dval(eps) = ieps*dval(d) + 7.; | 3818 » » dval(&eps) = ieps*dval(&u) + 7.; |
2970 » » word0(eps) -= (P-1)*Exp_msk1; | 3819 » » word0(&eps) -= (P-1)*Exp_msk1; |
2971 if (ilim == 0) { | 3820 if (ilim == 0) { |
2972 S = mhi = 0; | 3821 S = mhi = 0; |
2973 » » » dval(d) -= 5.; | 3822 » » » dval(&u) -= 5.; |
2974 » » » if (dval(d) > dval(eps)) | 3823 » » » if (dval(&u) > dval(&eps)) |
2975 goto one_digit; | 3824 goto one_digit; |
2976 » » » if (dval(d) < -dval(eps)) | 3825 » » » if (dval(&u) < -dval(&eps)) |
2977 goto no_digits; | 3826 goto no_digits; |
2978 goto fast_failed; | 3827 goto fast_failed; |
2979 } | 3828 } |
2980 #ifndef No_leftright | 3829 #ifndef No_leftright |
2981 if (leftright) { | 3830 if (leftright) { |
2982 /* Use Steele & White method of only | 3831 /* Use Steele & White method of only |
2983 * generating digits needed. | 3832 * generating digits needed. |
2984 */ | 3833 */ |
2985 » » » dval(eps) = 0.5/tens[ilim-1] - dval(eps); | 3834 » » » dval(&eps) = 0.5/tens[ilim-1] - dval(&eps); |
2986 for(i = 0;;) { | 3835 for(i = 0;;) { |
2987 » » » » L = dval(d); | 3836 » » » » L = dval(&u); |
2988 » » » » dval(d) -= L; | 3837 » » » » dval(&u) -= L; |
2989 *s++ = '0' + (int)L; | 3838 *s++ = '0' + (int)L; |
2990 » » » » if (dval(d) < dval(eps)) | 3839 » » » » if (dval(&u) < dval(&eps)) |
2991 goto ret1; | 3840 goto ret1; |
2992 » » » » if (1. - dval(d) < dval(eps)) | 3841 » » » » if (1. - dval(&u) < dval(&eps)) |
2993 goto bump_up; | 3842 goto bump_up; |
2994 if (++i >= ilim) | 3843 if (++i >= ilim) |
2995 break; | 3844 break; |
2996 » » » » dval(eps) *= 10.; | 3845 » » » » dval(&eps) *= 10.; |
2997 » » » » dval(d) *= 10.; | 3846 » » » » dval(&u) *= 10.; |
2998 } | 3847 } |
2999 } | 3848 } |
3000 else { | 3849 else { |
3001 #endif | 3850 #endif |
3002 /* Generate ilim digits, then fix them up. */ | 3851 /* Generate ilim digits, then fix them up. */ |
3003 » » » dval(eps) *= tens[ilim-1]; | 3852 » » » dval(&eps) *= tens[ilim-1]; |
3004 » » » for(i = 1;; i++, dval(d) *= 10.) { | 3853 » » » for(i = 1;; i++, dval(&u) *= 10.) { |
3005 » » » » L = (Long)(dval(d)); | 3854 » » » » L = (Long)(dval(&u)); |
3006 » » » » if (!(dval(d) -= L)) | 3855 » » » » if (!(dval(&u) -= L)) |
3007 ilim = i; | 3856 ilim = i; |
3008 *s++ = '0' + (int)L; | 3857 *s++ = '0' + (int)L; |
3009 if (i == ilim) { | 3858 if (i == ilim) { |
3010 » » » » » if (dval(d) > 0.5 + dval(eps)) | 3859 » » » » » if (dval(&u) > 0.5 + dval(&eps)) |
3011 goto bump_up; | 3860 goto bump_up; |
3012 » » » » » else if (dval(d) < 0.5 - dval(eps)) { | 3861 » » » » » else if (dval(&u) < 0.5 - dval(&eps)) { |
3013 while(*--s == '0'); | 3862 while(*--s == '0'); |
3014 s++; | 3863 s++; |
3015 goto ret1; | 3864 goto ret1; |
3016 } | 3865 } |
3017 break; | 3866 break; |
3018 } | 3867 } |
3019 } | 3868 } |
3020 #ifndef No_leftright | 3869 #ifndef No_leftright |
3021 } | 3870 } |
3022 #endif | 3871 #endif |
3023 fast_failed: | 3872 fast_failed: |
3024 s = s0; | 3873 s = s0; |
3025 » » dval(d) = dval(d2); | 3874 » » dval(&u) = dval(&d2); |
3026 k = k0; | 3875 k = k0; |
3027 ilim = ilim0; | 3876 ilim = ilim0; |
3028 } | 3877 } |
3029 | 3878 |
3030 /* Do we have a "small" integer? */ | 3879 /* Do we have a "small" integer? */ |
3031 | 3880 |
3032 if (be >= 0 && k <= Int_max) { | 3881 if (be >= 0 && k <= Int_max) { |
3033 /* Yes. */ | 3882 /* Yes. */ |
3034 ds = tens[k]; | 3883 ds = tens[k]; |
3035 if (ndigits < 0 && ilim <= 0) { | 3884 if (ndigits < 0 && ilim <= 0) { |
3036 S = mhi = 0; | 3885 S = mhi = 0; |
3037 » » » if (ilim < 0 || dval(d) <= 5*ds) | 3886 » » » if (ilim < 0 || dval(&u) <= 5*ds) |
3038 goto no_digits; | 3887 goto no_digits; |
3039 goto one_digit; | 3888 goto one_digit; |
3040 } | 3889 } |
3041 » » for(i = 1;; i++, dval(d) *= 10.) { | 3890 » » for(i = 1;; i++, dval(&u) *= 10.) { |
3042 » » » L = (Long)(dval(d) / ds); | 3891 » » » L = (Long)(dval(&u) / ds); |
3043 » » » dval(d) -= L*ds; | 3892 » » » dval(&u) -= L*ds; |
3044 #ifdef Check_FLT_ROUNDS | 3893 #ifdef Check_FLT_ROUNDS |
3045 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ | 3894 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ |
3046 » » » if (dval(d) < 0) { | 3895 » » » if (dval(&u) < 0) { |
3047 L--; | 3896 L--; |
3048 » » » » dval(d) += ds; | 3897 » » » » dval(&u) += ds; |
3049 } | 3898 } |
3050 #endif | 3899 #endif |
3051 *s++ = '0' + (int)L; | 3900 *s++ = '0' + (int)L; |
3052 » » » if (!dval(d)) { | 3901 » » » if (!dval(&u)) { |
3053 #ifdef SET_INEXACT | 3902 #ifdef SET_INEXACT |
3054 inexact = 0; | 3903 inexact = 0; |
3055 #endif | 3904 #endif |
3056 break; | 3905 break; |
3057 } | 3906 } |
3058 if (i == ilim) { | 3907 if (i == ilim) { |
3059 #ifdef Honor_FLT_ROUNDS | 3908 #ifdef Honor_FLT_ROUNDS |
3060 if (mode > 1) | 3909 if (mode > 1) |
3061 switch(Rounding) { | 3910 switch(Rounding) { |
3062 case 0: goto ret1; | 3911 case 0: goto ret1; |
3063 case 2: goto bump_up; | 3912 case 2: goto bump_up; |
3064 } | 3913 } |
3065 #endif | 3914 #endif |
3066 » » » » dval(d) += dval(d); | 3915 » » » » dval(&u) += dval(&u); |
3067 » » » » if (dval(d) > ds || dval(d) == ds && L & 1) { | 3916 » » » » if (dval(&u) > ds || (dval(&u) == ds && L & 1))
{ |
3068 bump_up: | 3917 bump_up: |
3069 while(*--s == '9') | 3918 while(*--s == '9') |
3070 if (s == s0) { | 3919 if (s == s0) { |
3071 k++; | 3920 k++; |
3072 *s = '0'; | 3921 *s = '0'; |
3073 break; | 3922 break; |
3074 } | 3923 } |
3075 ++*s++; | 3924 ++*s++; |
3076 } | 3925 } |
3077 break; | 3926 break; |
(...skipping 44 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
3122 S = pow5mult(S, s5); | 3971 S = pow5mult(S, s5); |
3123 | 3972 |
3124 /* Check for special case that d is a normalized power of 2. */ | 3973 /* Check for special case that d is a normalized power of 2. */ |
3125 | 3974 |
3126 spec_case = 0; | 3975 spec_case = 0; |
3127 if ((mode < 2 || leftright) | 3976 if ((mode < 2 || leftright) |
3128 #ifdef Honor_FLT_ROUNDS | 3977 #ifdef Honor_FLT_ROUNDS |
3129 && Rounding == 1 | 3978 && Rounding == 1 |
3130 #endif | 3979 #endif |
3131 ) { | 3980 ) { |
3132 » » if (!word1(d) && !(word0(d) & Bndry_mask) | 3981 » » if (!word1(&u) && !(word0(&u) & Bndry_mask) |
3133 #ifndef Sudden_Underflow | 3982 #ifndef Sudden_Underflow |
3134 » » && word0(d) & (Exp_mask & ~Exp_msk1) | 3983 » » && word0(&u) & (Exp_mask & ~Exp_msk1) |
3135 #endif | 3984 #endif |
3136 ) { | 3985 ) { |
3137 /* The special case */ | 3986 /* The special case */ |
3138 b2 += Log2P; | 3987 b2 += Log2P; |
3139 s2 += Log2P; | 3988 s2 += Log2P; |
3140 spec_case = 1; | 3989 spec_case = 1; |
3141 } | 3990 } |
3142 } | 3991 } |
3143 | 3992 |
3144 /* Arrange for convenient computation of quotients: | 3993 /* Arrange for convenient computation of quotients: |
3145 * shift left if necessary so divisor has 4 leading 0 bits. | 3994 * shift left if necessary so divisor has 4 leading 0 bits. |
3146 * | 3995 * |
3147 * Perhaps we should just compute leading 28 bits of S once | 3996 * Perhaps we should just compute leading 28 bits of S once |
3148 * and for all and pass them and a shift to quorem, so it | 3997 * and for all and pass them and a shift to quorem, so it |
3149 * can do shifts and ors to compute the numerator for q. | 3998 * can do shifts and ors to compute the numerator for q. |
3150 */ | 3999 */ |
3151 #ifdef Pack_32 | 4000 #ifdef Pack_32 |
3152 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)) | 4001 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)) |
3153 i = 32 - i; | 4002 i = 32 - i; |
| 4003 #define iInc 28 |
3154 #else | 4004 #else |
3155 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) | 4005 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) |
3156 i = 16 - i; | 4006 i = 16 - i; |
| 4007 #define iInc 12 |
3157 #endif | 4008 #endif |
3158 » if (i > 4) { | 4009 » i = dshift(S, s2); |
3159 » » i -= 4; | 4010 » b2 += i; |
3160 » » b2 += i; | 4011 » m2 += i; |
3161 » » m2 += i; | 4012 » s2 += i; |
3162 » » s2 += i; | |
3163 » » } | |
3164 » else if (i < 4) { | |
3165 » » i += 28; | |
3166 » » b2 += i; | |
3167 » » m2 += i; | |
3168 » » s2 += i; | |
3169 » » } | |
3170 if (b2 > 0) | 4013 if (b2 > 0) |
3171 b = lshift(b, b2); | 4014 b = lshift(b, b2); |
3172 if (s2 > 0) | 4015 if (s2 > 0) |
3173 S = lshift(S, s2); | 4016 S = lshift(S, s2); |
3174 if (k_check) { | 4017 if (k_check) { |
3175 if (cmp(b,S) < 0) { | 4018 if (cmp(b,S) < 0) { |
3176 k--; | 4019 k--; |
3177 b = multadd(b, 10, 0); /* we botched the k estimate */ | 4020 b = multadd(b, 10, 0); /* we botched the k estimate */ |
3178 if (leftright) | 4021 if (leftright) |
3179 mhi = multadd(mhi, 10, 0); | 4022 mhi = multadd(mhi, 10, 0); |
(...skipping 30 matching lines...) Expand all Loading... |
3210 for(i = 1;;i++) { | 4053 for(i = 1;;i++) { |
3211 dig = quorem(b,S) + '0'; | 4054 dig = quorem(b,S) + '0'; |
3212 /* Do we yet have the shortest decimal string | 4055 /* Do we yet have the shortest decimal string |
3213 * that will round to d? | 4056 * that will round to d? |
3214 */ | 4057 */ |
3215 j = cmp(b, mlo); | 4058 j = cmp(b, mlo); |
3216 delta = diff(S, mhi); | 4059 delta = diff(S, mhi); |
3217 j1 = delta->sign ? 1 : cmp(b, delta); | 4060 j1 = delta->sign ? 1 : cmp(b, delta); |
3218 Bfree(delta); | 4061 Bfree(delta); |
3219 #ifndef ROUND_BIASED | 4062 #ifndef ROUND_BIASED |
3220 » » » if (j1 == 0 && mode != 1 && !(word1(d) & 1) | 4063 » » » if (j1 == 0 && mode != 1 && !(word1(&u) & 1) |
3221 #ifdef Honor_FLT_ROUNDS | 4064 #ifdef Honor_FLT_ROUNDS |
3222 && Rounding >= 1 | 4065 && Rounding >= 1 |
3223 #endif | 4066 #endif |
3224 ) { | 4067 ) { |
3225 if (dig == '9') | 4068 if (dig == '9') |
3226 goto round_9_up; | 4069 goto round_9_up; |
3227 if (j > 0) | 4070 if (j > 0) |
3228 dig++; | 4071 dig++; |
3229 #ifdef SET_INEXACT | 4072 #ifdef SET_INEXACT |
3230 else if (!b->x[0] && b->wds <= 1) | 4073 else if (!b->x[0] && b->wds <= 1) |
3231 inexact = 0; | 4074 inexact = 0; |
3232 #endif | 4075 #endif |
3233 *s++ = dig; | 4076 *s++ = dig; |
3234 goto ret; | 4077 goto ret; |
3235 } | 4078 } |
3236 #endif | 4079 #endif |
3237 » » » if (j < 0 || j == 0 && mode != 1 | 4080 » » » if (j < 0 || (j == 0 && mode != 1 |
3238 #ifndef ROUND_BIASED | 4081 #ifndef ROUND_BIASED |
3239 » » » » » » » && !(word1(d) & 1) | 4082 » » » » » » » && !(word1(&u) & 1) |
3240 #endif | 4083 #endif |
3241 » » » » » ) { | 4084 » » » » » )) { |
3242 if (!b->x[0] && b->wds <= 1) { | 4085 if (!b->x[0] && b->wds <= 1) { |
3243 #ifdef SET_INEXACT | 4086 #ifdef SET_INEXACT |
3244 inexact = 0; | 4087 inexact = 0; |
3245 #endif | 4088 #endif |
3246 goto accept_dig; | 4089 goto accept_dig; |
3247 } | 4090 } |
3248 #ifdef Honor_FLT_ROUNDS | 4091 #ifdef Honor_FLT_ROUNDS |
3249 if (mode > 1) | 4092 if (mode > 1) |
3250 switch(Rounding) { | 4093 switch(Rounding) { |
3251 case 0: goto accept_dig; | 4094 case 0: goto accept_dig; |
3252 case 2: goto keep_dig; | 4095 case 2: goto keep_dig; |
3253 } | 4096 } |
3254 #endif /*Honor_FLT_ROUNDS*/ | 4097 #endif /*Honor_FLT_ROUNDS*/ |
3255 if (j1 > 0) { | 4098 if (j1 > 0) { |
3256 b = lshift(b, 1); | 4099 b = lshift(b, 1); |
3257 j1 = cmp(b, S); | 4100 j1 = cmp(b, S); |
3258 » » » » » if ((j1 > 0 || j1 == 0 && dig & 1) | 4101 » » » » » if ((j1 > 0 || (j1 == 0 && dig & 1)) |
3259 && dig++ == '9') | 4102 && dig++ == '9') |
3260 goto round_9_up; | 4103 goto round_9_up; |
3261 } | 4104 } |
3262 accept_dig: | 4105 accept_dig: |
3263 *s++ = dig; | 4106 *s++ = dig; |
3264 goto ret; | 4107 goto ret; |
3265 } | 4108 } |
3266 if (j1 > 0) { | 4109 if (j1 > 0) { |
3267 #ifdef Honor_FLT_ROUNDS | 4110 #ifdef Honor_FLT_ROUNDS |
3268 if (!Rounding) | 4111 if (!Rounding) |
(...skipping 39 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
3308 /* Round off last digit */ | 4151 /* Round off last digit */ |
3309 | 4152 |
3310 #ifdef Honor_FLT_ROUNDS | 4153 #ifdef Honor_FLT_ROUNDS |
3311 switch(Rounding) { | 4154 switch(Rounding) { |
3312 case 0: goto trimzeros; | 4155 case 0: goto trimzeros; |
3313 case 2: goto roundoff; | 4156 case 2: goto roundoff; |
3314 } | 4157 } |
3315 #endif | 4158 #endif |
3316 b = lshift(b, 1); | 4159 b = lshift(b, 1); |
3317 j = cmp(b, S); | 4160 j = cmp(b, S); |
3318 » if (j > 0 || j == 0 && dig & 1) { | 4161 » if (j > 0 || (j == 0 && dig & 1)) { |
3319 roundoff: | 4162 roundoff: |
3320 while(*--s == '9') | 4163 while(*--s == '9') |
3321 if (s == s0) { | 4164 if (s == s0) { |
3322 k++; | 4165 k++; |
3323 *s++ = '1'; | 4166 *s++ = '1'; |
3324 goto ret; | 4167 goto ret; |
3325 } | 4168 } |
3326 ++*s++; | 4169 ++*s++; |
3327 } | 4170 } |
3328 else { | 4171 else { |
3329 #ifdef Honor_FLT_ROUNDS | 4172 #ifdef Honor_FLT_ROUNDS |
3330 trimzeros: | 4173 trimzeros: |
3331 #endif | 4174 #endif |
3332 while(*--s == '0'); | 4175 while(*--s == '0'); |
3333 s++; | 4176 s++; |
3334 } | 4177 } |
3335 ret: | 4178 ret: |
3336 Bfree(S); | 4179 Bfree(S); |
3337 if (mhi) { | 4180 if (mhi) { |
3338 if (mlo && mlo != mhi) | 4181 if (mlo && mlo != mhi) |
3339 Bfree(mlo); | 4182 Bfree(mlo); |
3340 Bfree(mhi); | 4183 Bfree(mhi); |
3341 } | 4184 } |
3342 ret1: | 4185 ret1: |
3343 #ifdef SET_INEXACT | 4186 #ifdef SET_INEXACT |
3344 if (inexact) { | 4187 if (inexact) { |
3345 if (!oldinexact) { | 4188 if (!oldinexact) { |
3346 » » » word0(d) = Exp_1 + (70 << Exp_shift); | 4189 » » » word0(&u) = Exp_1 + (70 << Exp_shift); |
3347 » » » word1(d) = 0; | 4190 » » » word1(&u) = 0; |
3348 » » » dval(d) += 1.; | 4191 » » » dval(&u) += 1.; |
3349 } | 4192 } |
3350 } | 4193 } |
3351 else if (!oldinexact) | 4194 else if (!oldinexact) |
3352 clear_inexact(); | 4195 clear_inexact(); |
3353 #endif | 4196 #endif |
3354 Bfree(b); | 4197 Bfree(b); |
3355 *s = 0; | 4198 *s = 0; |
3356 *decpt = k + 1; | 4199 *decpt = k + 1; |
3357 if (rve) | 4200 if (rve) |
3358 *rve = s; | 4201 *rve = s; |
3359 return s0; | 4202 return s0; |
3360 } | 4203 } |
3361 | 4204 |
3362 } // namespace dmg_fp | 4205 } // namespace dmg_fp |
OLD | NEW |