File indexing completed on 2024-05-12 15:43:19
0001 /**************************************************************** 0002 * 0003 * The author of this software is David M. Gay. 0004 * 0005 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. 0006 * 0007 * Permission to use, copy, modify, and distribute this software for any 0008 * purpose without fee is hereby granted, provided that this entire notice 0009 * is included in all copies of any software which is or includes a copy 0010 * or modification of this software and in all copies of the supporting 0011 * documentation for such software. 0012 * 0013 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 0014 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY 0015 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 0016 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 0017 * 0018 ***************************************************************/ 0019 0020 /* Please send bug reports to 0021 David M. Gay 0022 Bell Laboratories, Room 2C-463 0023 600 Mountain Avenue 0024 Murray Hill, NJ 07974-0636 0025 U.S.A. 0026 dmg@bell-labs.com 0027 */ 0028 0029 /* On a machine with IEEE extended-precision registers, it is 0030 * necessary to specify double-precision (53-bit) rounding precision 0031 * before invoking strtod or dtoa. If the machine uses (the equivalent 0032 * of) Intel 80x87 arithmetic, the call 0033 * _control87(PC_53, MCW_PC); 0034 * does this with many compilers. Whether this or another call is 0035 * appropriate depends on the compiler; for this to work, it may be 0036 * necessary to #include "float.h" or another system-dependent header 0037 * file. 0038 */ 0039 0040 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines. 0041 * 0042 * This strtod returns a nearest machine number to the input decimal 0043 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 0044 * broken by the IEEE round-even rule. Otherwise ties are broken by 0045 * biased rounding (add half and chop). 0046 * 0047 * Inspired loosely by William D. Clinger's paper "How to Read Floating 0048 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 0049 * 0050 * Modifications: 0051 * 0052 * 1. We only require IEEE, IBM, or VAX double-precision 0053 * arithmetic (not IEEE double-extended). 0054 * 2. We get by with floating-point arithmetic in a case that 0055 * Clinger missed -- when we're computing d * 10^n 0056 * for a small integer d and the integer n is not too 0057 * much larger than 22 (the maximum integer k for which 0058 * we can represent 10^k exactly), we may be able to 0059 * compute (d*10^k) * 10^(e-k) with just one roundoff. 0060 * 3. Rather than a bit-at-a-time adjustment of the binary 0061 * result in the hard case, we use floating-point 0062 * arithmetic to determine the adjustment to within 0063 * one bit; only in really hard cases do we need to 0064 * compute a second residual. 0065 * 4. Because of 3., we don't need a large table of powers of 10 0066 * for ten-to-e (just some small tables, e.g. of 10^k 0067 * for 0 <= k <= 22). 0068 */ 0069 0070 /* 0071 * #define IEEE_8087 for IEEE-arithmetic machines where the least 0072 * significant byte has the lowest address. 0073 * #define IEEE_MC68k for IEEE-arithmetic machines where the most 0074 * significant byte has the lowest address. 0075 * #define Long int on machines with 32-bit ints and 64-bit longs. 0076 * #define IBM for IBM mainframe-style floating-point arithmetic. 0077 * #define VAX for VAX-style floating-point arithmetic (D_floating). 0078 * #define No_leftright to omit left-right logic in fast floating-point 0079 * computation of dtoa. 0080 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 0081 * and strtod and dtoa should round accordingly. 0082 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 0083 * and Honor_FLT_ROUNDS is not #defined. 0084 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 0085 * that use extended-precision instructions to compute rounded 0086 * products and quotients) with IBM. 0087 * #define ROUND_BIASED for IEEE-format with biased rounding. 0088 * #define Inaccurate_Divide for IEEE-format with correctly rounded 0089 * products but inaccurate quotients, e.g., for Intel i860. 0090 * #define NO_LONG_LONG on machines that do not have a "long long" 0091 * integer type (of >= 64 bits). On such machines, you can 0092 * #define Just_16 to store 16 bits per 32-bit Long when doing 0093 * high-precision integer arithmetic. Whether this speeds things 0094 * up or slows things down depends on the machine and the number 0095 * being converted. If long long is available and the name is 0096 * something other than "long long", #define Llong to be the name, 0097 * and if "unsigned Llong" does not work as an unsigned version of 0098 * Llong, #define #ULLong to be the corresponding unsigned type. 0099 * #define KR_headers for old-style C function headers. 0100 * #define Bad_float_h if your system lacks a float.h or if it does not 0101 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, 0102 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. 0103 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) 0104 * if memory is available and otherwise does something you deem 0105 * appropriate. If MALLOC is undefined, malloc will be invoked 0106 * directly -- and assumed always to succeed. 0107 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making 0108 * memory allocations from a private pool of memory when possible. 0109 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes, 0110 * unless #defined to be a different length. This default length 0111 * suffices to get rid of MALLOC calls except for unusual cases, 0112 * such as decimal-to-binary conversion of a very long string of 0113 * digits. The longest string dtoa can return is about 751 bytes 0114 * long. For conversions by strtod of strings of 800 digits and 0115 * all dtoa conversions in single-threaded executions with 8-byte 0116 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte 0117 * pointers, PRIVATE_MEM >= 7112 appears adequate. 0118 * #define INFNAN_CHECK on IEEE systems to cause strtod to check for 0119 * Infinity and NaN (case insensitively). On some systems (e.g., 0120 * some HP systems), it may be necessary to #define NAN_WORD0 0121 * appropriately -- to the most significant word of a quiet NaN. 0122 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.) 0123 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined, 0124 * strtod also accepts (case insensitively) strings of the form 0125 * NaN(x), where x is a string of hexadecimal digits and spaces; 0126 * if there is only one string of hexadecimal digits, it is taken 0127 * for the 52 fraction bits of the resulting NaN; if there are two 0128 * or more strings of hex digits, the first is for the high 20 bits, 0129 * the second and subsequent for the low 32 bits, with intervening 0130 * white space ignored; but if this results in none of the 52 0131 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0 0132 * and NAN_WORD1 are used instead. 0133 * #define MULTIPLE_THREADS if the system offers preemptively scheduled 0134 * multiple threads. In this case, you must provide (or suitably 0135 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed 0136 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed 0137 * in pow5mult, ensures lazy evaluation of only one copy of high 0138 * powers of 5; omitting this lock would introduce a small 0139 * probability of wasting memory, but would otherwise be harmless.) 0140 * You must also invoke freedtoa(s) to free the value s returned by 0141 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. 0142 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that 0143 * avoids underflows on inputs whose result does not underflow. 0144 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format 0145 * floating-point numbers and flushes underflows to zero rather 0146 * than implementing gradual underflow, then you must also #define 0147 * Sudden_Underflow. 0148 * #define YES_ALIAS to permit aliasing certain double values with 0149 * arrays of ULongs. This leads to slightly better code with 0150 * some compilers and was always used prior to 19990916, but it 0151 * is not strictly legal and can cause trouble with aggressively 0152 * optimizing compilers (e.g., gcc 2.95.1 under -O2). 0153 * #define USE_LOCALE to use the current locale's decimal_point value. 0154 * #define SET_INEXACT if IEEE arithmetic is being used and extra 0155 * computation should be done to set the inexact flag when the 0156 * result is inexact and avoid setting inexact when the result 0157 * is exact. In this case, dtoa.c must be compiled in 0158 * an environment, perhaps provided by #include "dtoa.c" in a 0159 * suitable wrapper, that defines two functions, 0160 * int get_inexact(void); 0161 * void clear_inexact(void); 0162 * such that get_inexact() returns a nonzero value if the 0163 * inexact bit is already set, and clear_inexact() sets the 0164 * inexact bit to 0. When SET_INEXACT is #defined, strtod 0165 * also does extra computations to set the underflow and overflow 0166 * flags when appropriate (i.e., when the result is tiny and 0167 * inexact or when it is a numeric value rounded to +-infinity). 0168 * #define NO_ERRNO if strtod should not assign errno = ERANGE when 0169 * the result overflows to +-Infinity or underflows to 0. 0170 */ 0171 0172 #include "dtoa.h" 0173 0174 #include "global.h" 0175 0176 #if PLATFORM(BIG_ENDIAN) 0177 #define IEEE_MC68k 0178 #else 0179 #define IEEE_8087 0180 #endif 0181 #define INFNAN_CHECK 0182 0183 #ifndef Long 0184 #define Long int 0185 #endif 0186 #ifndef ULong 0187 typedef unsigned Long ULong; 0188 #endif 0189 0190 #ifdef DEBUG 0191 #include <stdio.h> 0192 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 0193 #endif 0194 0195 #include <stdlib.h> 0196 #include <string.h> 0197 0198 #ifdef USE_LOCALE 0199 #include <locale.h> 0200 #endif 0201 0202 #ifdef MALLOC 0203 extern void *MALLOC(size_t); 0204 #else 0205 #define MALLOC malloc 0206 #endif 0207 0208 #ifndef Omit_Private_Memory 0209 #ifndef PRIVATE_MEM 0210 #define PRIVATE_MEM 2304 0211 #endif 0212 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) 0213 static double private_mem[PRIVATE_mem], *pmem_next = private_mem; 0214 #endif 0215 0216 #undef IEEE_Arith 0217 #undef Avoid_Underflow 0218 #ifdef IEEE_MC68k 0219 #define IEEE_Arith 0220 #endif 0221 #ifdef IEEE_8087 0222 #define IEEE_Arith 0223 #endif 0224 0225 #if HAVE_ERRNO_H 0226 #include <errno.h> 0227 #endif 0228 0229 #ifdef Bad_float_h 0230 0231 #ifdef IEEE_Arith 0232 #define DBL_DIG 15 0233 #define DBL_MAX_10_EXP 308 0234 #define DBL_MAX_EXP 1024 0235 #define FLT_RADIX 2 0236 #endif /*IEEE_Arith*/ 0237 0238 #ifdef IBM 0239 #define DBL_DIG 16 0240 #define DBL_MAX_10_EXP 75 0241 #define DBL_MAX_EXP 63 0242 #define FLT_RADIX 16 0243 #define DBL_MAX 7.2370055773322621e+75 0244 #endif 0245 0246 #ifdef VAX 0247 #define DBL_DIG 16 0248 #define DBL_MAX_10_EXP 38 0249 #define DBL_MAX_EXP 127 0250 #define FLT_RADIX 2 0251 #define DBL_MAX 1.7014118346046923e+38 0252 #endif 0253 0254 #ifndef LONG_MAX 0255 #define LONG_MAX 2147483647 0256 #endif 0257 0258 #else /* ifndef Bad_float_h */ 0259 #include <float.h> 0260 #endif /* Bad_float_h */ 0261 0262 #ifndef __MATH_H__ 0263 #include <math.h> 0264 #endif 0265 0266 #define strtod kjs_strtod 0267 #define dtoa kjs_dtoa 0268 #define freedtoa kjs_freedtoa 0269 0270 #ifdef __cplusplus 0271 extern "C" { 0272 #endif 0273 0274 #ifndef CONST 0275 #define CONST const 0276 #endif 0277 0278 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 0279 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. 0280 #endif 0281 0282 typedef union { 0283 double d; 0284 ULong L[2]; 0285 } U; 0286 0287 #define dval(x) (x).d 0288 #ifdef IEEE_8087 0289 #define word0(x) (x).L[1] 0290 #define word1(x) (x).L[0] 0291 #else 0292 #define word0(x) (x).L[0] 0293 #define word1(x) (x).L[1] 0294 #endif 0295 0296 /* The following definition of Storeinc is appropriate for MIPS processors. 0297 * An alternative that might be better on some machines is 0298 */ 0299 #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 0300 0301 /* #define P DBL_MANT_DIG */ 0302 /* Ten_pmax = floor(P*log(2)/log(5)) */ 0303 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 0304 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 0305 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 0306 0307 #ifdef IEEE_Arith 0308 #define Exp_shift 20 0309 #define Exp_shift1 20 0310 #define Exp_msk1 0x100000 0311 #define Exp_msk11 0x100000 0312 #define Exp_mask 0x7ff00000 0313 #define P 53 0314 #define Bias 1023 0315 #define Emin (-1022) 0316 #define Exp_1 0x3ff00000 0317 #define Exp_11 0x3ff00000 0318 #define Ebits 11 0319 #define Frac_mask 0xfffff 0320 #define Frac_mask1 0xfffff 0321 #define Ten_pmax 22 0322 #define Bletch 0x10 0323 #define Bndry_mask 0xfffff 0324 #define Bndry_mask1 0xfffff 0325 #define LSB 1 0326 #define Sign_bit 0x80000000 0327 #define Log2P 1 0328 #define Tiny0 0 0329 #define Tiny1 1 0330 #define Quick_max 14 0331 #define Int_max 14 0332 #ifndef NO_IEEE_Scale 0333 #define Avoid_Underflow 0334 #ifdef Flush_Denorm /* debugging option */ 0335 #undef Sudden_Underflow 0336 #endif 0337 #endif 0338 0339 #ifndef Flt_Rounds 0340 #ifdef FLT_ROUNDS 0341 #define Flt_Rounds FLT_ROUNDS 0342 #else 0343 #define Flt_Rounds 1 0344 #endif 0345 #endif /*Flt_Rounds*/ 0346 0347 #ifdef Honor_FLT_ROUNDS 0348 #define Rounding rounding 0349 #undef Check_FLT_ROUNDS 0350 #define Check_FLT_ROUNDS 0351 #else 0352 #define Rounding Flt_Rounds 0353 #endif 0354 0355 #else /* ifndef IEEE_Arith */ 0356 #undef Check_FLT_ROUNDS 0357 #undef Honor_FLT_ROUNDS 0358 #undef SET_INEXACT 0359 #undef Sudden_Underflow 0360 #define Sudden_Underflow 0361 #ifdef IBM 0362 #undef Flt_Rounds 0363 #define Flt_Rounds 0 0364 #define Exp_shift 24 0365 #define Exp_shift1 24 0366 #define Exp_msk1 0x1000000 0367 #define Exp_msk11 0x1000000 0368 #define Exp_mask 0x7f000000 0369 #define P 14 0370 #define Bias 65 0371 #define Exp_1 0x41000000 0372 #define Exp_11 0x41000000 0373 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 0374 #define Frac_mask 0xffffff 0375 #define Frac_mask1 0xffffff 0376 #define Bletch 4 0377 #define Ten_pmax 22 0378 #define Bndry_mask 0xefffff 0379 #define Bndry_mask1 0xffffff 0380 #define LSB 1 0381 #define Sign_bit 0x80000000 0382 #define Log2P 4 0383 #define Tiny0 0x100000 0384 #define Tiny1 0 0385 #define Quick_max 14 0386 #define Int_max 15 0387 #else /* VAX */ 0388 #undef Flt_Rounds 0389 #define Flt_Rounds 1 0390 #define Exp_shift 23 0391 #define Exp_shift1 7 0392 #define Exp_msk1 0x80 0393 #define Exp_msk11 0x800000 0394 #define Exp_mask 0x7f80 0395 #define P 56 0396 #define Bias 129 0397 #define Exp_1 0x40800000 0398 #define Exp_11 0x4080 0399 #define Ebits 8 0400 #define Frac_mask 0x7fffff 0401 #define Frac_mask1 0xffff007f 0402 #define Ten_pmax 24 0403 #define Bletch 2 0404 #define Bndry_mask 0xffff007f 0405 #define Bndry_mask1 0xffff007f 0406 #define LSB 0x10000 0407 #define Sign_bit 0x8000 0408 #define Log2P 1 0409 #define Tiny0 0x80 0410 #define Tiny1 0 0411 #define Quick_max 15 0412 #define Int_max 15 0413 #endif /* IBM, VAX */ 0414 #endif /* IEEE_Arith */ 0415 0416 #ifndef IEEE_Arith 0417 #define ROUND_BIASED 0418 #endif 0419 0420 #ifdef RND_PRODQUOT 0421 #define rounded_product(a,b) a = rnd_prod(a, b) 0422 #define rounded_quotient(a,b) a = rnd_quot(a, b) 0423 extern double rnd_prod(double, double), rnd_quot(double, double); 0424 #else 0425 #define rounded_product(a,b) a *= b 0426 #define rounded_quotient(a,b) a /= b 0427 #endif 0428 0429 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 0430 #define Big1 0xffffffff 0431 0432 #ifndef Pack_32 0433 #define Pack_32 0434 #endif 0435 0436 #define FFFFFFFF 0xffffffffUL 0437 0438 #ifdef NO_LONG_LONG 0439 #undef ULLong 0440 #ifdef Just_16 0441 #undef Pack_32 0442 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long. 0443 * This makes some inner loops simpler and sometimes saves work 0444 * during multiplications, but it often seems to make things slightly 0445 * slower. Hence the default is now to store 32 bits per Long. 0446 */ 0447 #endif 0448 #else /* long long available */ 0449 #ifndef Llong 0450 #define Llong long long 0451 #endif 0452 #ifndef ULLong 0453 #define ULLong unsigned Llong 0454 #endif 0455 #endif /* NO_LONG_LONG */ 0456 0457 #ifndef MULTIPLE_THREADS 0458 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/ 0459 #define FREE_DTOA_LOCK(n) /*nothing*/ 0460 #endif 0461 0462 #define Kmax (sizeof(size_t) << 3) 0463 0464 struct 0465 Bigint { 0466 struct Bigint *next; 0467 int k, maxwds, sign, wds; 0468 ULong x[1]; 0469 }; 0470 0471 typedef struct Bigint Bigint; 0472 0473 static Bigint *freelist[Kmax + 1]; 0474 0475 static Bigint * 0476 Balloc 0477 (int k) 0478 { 0479 int x; 0480 Bigint *rv; 0481 #ifndef Omit_Private_Memory 0482 unsigned int len; 0483 #endif 0484 0485 ACQUIRE_DTOA_LOCK(0); 0486 if ((rv = freelist[k])) { 0487 freelist[k] = rv->next; 0488 } else { 0489 x = 1 << k; 0490 #ifdef Omit_Private_Memory 0491 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x - 1) * sizeof(ULong)); 0492 #else 0493 len = (sizeof(Bigint) + (x - 1) * sizeof(ULong) + sizeof(double) - 1) 0494 / sizeof(double); 0495 if (pmem_next - private_mem + len <= (unsigned)PRIVATE_mem) { 0496 rv = (Bigint *)pmem_next; 0497 pmem_next += len; 0498 } else { 0499 rv = (Bigint *)MALLOC(len * sizeof(double)); 0500 } 0501 #endif 0502 rv->k = k; 0503 rv->maxwds = x; 0504 } 0505 FREE_DTOA_LOCK(0); 0506 rv->sign = rv->wds = 0; 0507 return rv; 0508 } 0509 0510 static void 0511 Bfree 0512 (Bigint *v) 0513 { 0514 if (v) { 0515 ACQUIRE_DTOA_LOCK(0); 0516 v->next = freelist[v->k]; 0517 freelist[v->k] = v; 0518 FREE_DTOA_LOCK(0); 0519 } 0520 } 0521 0522 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 0523 y->wds*sizeof(Long) + 2*sizeof(int)) 0524 0525 static Bigint * 0526 multadd 0527 (Bigint *b, int m, int a) /* multiply by m and add a */ 0528 { 0529 int i, wds; 0530 #ifdef ULLong 0531 ULong *x; 0532 ULLong carry, y; 0533 #else 0534 ULong carry, *x, y; 0535 #ifdef Pack_32 0536 ULong xi, z; 0537 #endif 0538 #endif 0539 Bigint *b1; 0540 0541 wds = b->wds; 0542 x = b->x; 0543 i = 0; 0544 carry = a; 0545 do { 0546 #ifdef ULLong 0547 y = *x * (ULLong)m + carry; 0548 carry = y >> 32; 0549 *x++ = (ULong)y & FFFFFFFF; 0550 #else 0551 #ifdef Pack_32 0552 xi = *x; 0553 y = (xi & 0xffff) * m + carry; 0554 z = (xi >> 16) * m + (y >> 16); 0555 carry = z >> 16; 0556 *x++ = (z << 16) + (y & 0xffff); 0557 #else 0558 y = *x * m + carry; 0559 carry = y >> 16; 0560 *x++ = y & 0xffff; 0561 #endif 0562 #endif 0563 } while (++i < wds); 0564 if (carry) { 0565 if (wds >= b->maxwds) { 0566 b1 = Balloc(b->k + 1); 0567 Bcopy(b1, b); 0568 Bfree(b); 0569 b = b1; 0570 } 0571 b->x[wds++] = (ULong)carry; 0572 b->wds = wds; 0573 } 0574 return b; 0575 } 0576 0577 static Bigint * 0578 s2b 0579 (CONST char *s, int nd0, int nd, ULong y9) 0580 { 0581 Bigint *b; 0582 int i, k; 0583 Long x, y; 0584 0585 x = (nd + 8) / 9; 0586 for (k = 0, y = 1; x > y; y <<= 1, k++); 0587 #ifdef Pack_32 0588 b = Balloc(k); 0589 b->x[0] = y9; 0590 b->wds = 1; 0591 #else 0592 b = Balloc(k + 1); 0593 b->x[0] = y9 & 0xffff; 0594 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; 0595 #endif 0596 0597 i = 9; 0598 if (9 < nd0) { 0599 s += 9; 0600 do { 0601 b = multadd(b, 10, *s++ - '0'); 0602 } while (++i < nd0); 0603 s++; 0604 } else { 0605 s += 10; 0606 } 0607 for (; i < nd; i++) { 0608 b = multadd(b, 10, *s++ - '0'); 0609 } 0610 return b; 0611 } 0612 0613 static int 0614 hi0bits 0615 (ULong x) 0616 { 0617 int k = 0; 0618 0619 if (!(x & 0xffff0000)) { 0620 k = 16; 0621 x <<= 16; 0622 } 0623 if (!(x & 0xff000000)) { 0624 k += 8; 0625 x <<= 8; 0626 } 0627 if (!(x & 0xf0000000)) { 0628 k += 4; 0629 x <<= 4; 0630 } 0631 if (!(x & 0xc0000000)) { 0632 k += 2; 0633 x <<= 2; 0634 } 0635 if (!(x & 0x80000000)) { 0636 k++; 0637 if (!(x & 0x40000000)) { 0638 return 32; 0639 } 0640 } 0641 return k; 0642 } 0643 0644 static int 0645 lo0bits 0646 (ULong *y) 0647 { 0648 int k; 0649 ULong x = *y; 0650 0651 if (x & 7) { 0652 if (x & 1) { 0653 return 0; 0654 } 0655 if (x & 2) { 0656 *y = x >> 1; 0657 return 1; 0658 } 0659 *y = x >> 2; 0660 return 2; 0661 } 0662 k = 0; 0663 if (!(x & 0xffff)) { 0664 k = 16; 0665 x >>= 16; 0666 } 0667 if (!(x & 0xff)) { 0668 k += 8; 0669 x >>= 8; 0670 } 0671 if (!(x & 0xf)) { 0672 k += 4; 0673 x >>= 4; 0674 } 0675 if (!(x & 0x3)) { 0676 k += 2; 0677 x >>= 2; 0678 } 0679 if (!(x & 1)) { 0680 k++; 0681 x >>= 1; 0682 if (!x & 1) { 0683 return 32; 0684 } 0685 } 0686 *y = x; 0687 return k; 0688 } 0689 0690 static Bigint * 0691 i2b 0692 (int i) 0693 { 0694 Bigint *b; 0695 0696 b = Balloc(1); 0697 b->x[0] = i; 0698 b->wds = 1; 0699 return b; 0700 } 0701 0702 static Bigint * 0703 mult 0704 (Bigint *a, Bigint *b) 0705 { 0706 Bigint *c; 0707 int k, wa, wb, wc; 0708 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 0709 ULong y; 0710 #ifdef ULLong 0711 ULLong carry, z; 0712 #else 0713 ULong carry, z; 0714 #ifdef Pack_32 0715 ULong z2; 0716 #endif 0717 #endif 0718 0719 if (a->wds < b->wds) { 0720 c = a; 0721 a = b; 0722 b = c; 0723 } 0724 k = a->k; 0725 wa = a->wds; 0726 wb = b->wds; 0727 wc = wa + wb; 0728 if (wc > a->maxwds) { 0729 k++; 0730 } 0731 c = Balloc(k); 0732 for (x = c->x, xa = x + wc; x < xa; x++) { 0733 *x = 0; 0734 } 0735 xa = a->x; 0736 xae = xa + wa; 0737 xb = b->x; 0738 xbe = xb + wb; 0739 xc0 = c->x; 0740 #ifdef ULLong 0741 for (; xb < xbe; xc0++) { 0742 if ((y = *xb++)) { 0743 x = xa; 0744 xc = xc0; 0745 carry = 0; 0746 do { 0747 z = *x++ * (ULLong)y + *xc + carry; 0748 carry = z >> 32; 0749 *xc++ = (ULong)z & FFFFFFFF; 0750 } while (x < xae); 0751 *xc = (ULong)carry; 0752 } 0753 } 0754 #else 0755 #ifdef Pack_32 0756 for (; xb < xbe; xb++, xc0++) { 0757 if (y = *xb & 0xffff) { 0758 x = xa; 0759 xc = xc0; 0760 carry = 0; 0761 do { 0762 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 0763 carry = z >> 16; 0764 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 0765 carry = z2 >> 16; 0766 Storeinc(xc, z2, z); 0767 } while (x < xae); 0768 *xc = carry; 0769 } 0770 if (y = *xb >> 16) { 0771 x = xa; 0772 xc = xc0; 0773 carry = 0; 0774 z2 = *xc; 0775 do { 0776 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 0777 carry = z >> 16; 0778 Storeinc(xc, z, z2); 0779 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 0780 carry = z2 >> 16; 0781 } while (x < xae); 0782 *xc = z2; 0783 } 0784 } 0785 #else 0786 for (; xb < xbe; xc0++) { 0787 if (y = *xb++) { 0788 x = xa; 0789 xc = xc0; 0790 carry = 0; 0791 do { 0792 z = *x++ * y + *xc + carry; 0793 carry = z >> 16; 0794 *xc++ = z & 0xffff; 0795 } while (x < xae); 0796 *xc = carry; 0797 } 0798 } 0799 #endif 0800 #endif 0801 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc); 0802 c->wds = wc; 0803 return c; 0804 } 0805 0806 static Bigint *p5s; 0807 0808 static Bigint * 0809 pow5mult 0810 (Bigint *b, int k) 0811 { 0812 Bigint *b1, *p5, *p51; 0813 int i; 0814 static int p05[3] = { 5, 25, 125 }; 0815 0816 if ((i = k & 3)) { 0817 b = multadd(b, p05[i - 1], 0); 0818 } 0819 0820 if (!(k >>= 2)) { 0821 return b; 0822 } 0823 if (!(p5 = p5s)) { 0824 /* first time */ 0825 #ifdef MULTIPLE_THREADS 0826 ACQUIRE_DTOA_LOCK(1); 0827 if (!(p5 = p5s)) { 0828 p5 = p5s = i2b(625); 0829 p5->next = 0; 0830 } 0831 FREE_DTOA_LOCK(1); 0832 #else 0833 p5 = p5s = i2b(625); 0834 p5->next = nullptr; 0835 #endif 0836 } 0837 for (;;) { 0838 if (k & 1) { 0839 b1 = mult(b, p5); 0840 Bfree(b); 0841 b = b1; 0842 } 0843 if (!(k >>= 1)) { 0844 break; 0845 } 0846 if (!(p51 = p5->next)) { 0847 #ifdef MULTIPLE_THREADS 0848 ACQUIRE_DTOA_LOCK(1); 0849 if (!(p51 = p5->next)) { 0850 p51 = p5->next = mult(p5, p5); 0851 p51->next = 0; 0852 } 0853 FREE_DTOA_LOCK(1); 0854 #else 0855 p51 = p5->next = mult(p5, p5); 0856 p51->next = nullptr; 0857 #endif 0858 } 0859 p5 = p51; 0860 } 0861 return b; 0862 } 0863 0864 static Bigint * 0865 lshift 0866 (Bigint *b, int k) 0867 { 0868 int i, k1, n, n1; 0869 Bigint *b1; 0870 ULong *x, *x1, *xe, z; 0871 0872 #ifdef Pack_32 0873 n = k >> 5; 0874 #else 0875 n = k >> 4; 0876 #endif 0877 k1 = b->k; 0878 n1 = n + b->wds + 1; 0879 for (i = b->maxwds; n1 > i; i <<= 1) { 0880 k1++; 0881 } 0882 b1 = Balloc(k1); 0883 x1 = b1->x; 0884 for (i = 0; i < n; i++) { 0885 *x1++ = 0; 0886 } 0887 x = b->x; 0888 xe = x + b->wds; 0889 #ifdef Pack_32 0890 if (k &= 0x1f) { 0891 k1 = 32 - k; 0892 z = 0; 0893 do { 0894 *x1++ = *x << k | z; 0895 z = *x++ >> k1; 0896 } while (x < xe); 0897 if ((*x1 = z)) { 0898 ++n1; 0899 } 0900 } 0901 #else 0902 if (k &= 0xf) { 0903 k1 = 16 - k; 0904 z = 0; 0905 do { 0906 *x1++ = *x << k & 0xffff | z; 0907 z = *x++ >> k1; 0908 } while (x < xe); 0909 if (*x1 = z) { 0910 ++n1; 0911 } 0912 } 0913 #endif 0914 else do { 0915 *x1++ = *x++; 0916 } while (x < xe); 0917 b1->wds = n1 - 1; 0918 Bfree(b); 0919 return b1; 0920 } 0921 0922 static int 0923 cmp 0924 (Bigint *a, Bigint *b) 0925 { 0926 ULong *xa, *xa0, *xb, *xb0; 0927 int i, j; 0928 0929 i = a->wds; 0930 j = b->wds; 0931 #ifdef DEBUG 0932 if (i > 1 && !a->x[i - 1]) { 0933 Bug("cmp called with a->x[a->wds-1] == 0"); 0934 } 0935 if (j > 1 && !b->x[j - 1]) { 0936 Bug("cmp called with b->x[b->wds-1] == 0"); 0937 } 0938 #endif 0939 if (i -= j) { 0940 return i; 0941 } 0942 xa0 = a->x; 0943 xa = xa0 + j; 0944 xb0 = b->x; 0945 xb = xb0 + j; 0946 for (;;) { 0947 if (*--xa != *--xb) { 0948 return *xa < *xb ? -1 : 1; 0949 } 0950 if (xa <= xa0) { 0951 break; 0952 } 0953 } 0954 return 0; 0955 } 0956 0957 static Bigint * 0958 diff 0959 (Bigint *a, Bigint *b) 0960 { 0961 Bigint *c; 0962 int i, wa, wb; 0963 ULong *xa, *xae, *xb, *xbe, *xc; 0964 #ifdef ULLong 0965 ULLong borrow, y; 0966 #else 0967 ULong borrow, y; 0968 #ifdef Pack_32 0969 ULong z; 0970 #endif 0971 #endif 0972 0973 i = cmp(a, b); 0974 if (!i) { 0975 c = Balloc(0); 0976 c->wds = 1; 0977 c->x[0] = 0; 0978 return c; 0979 } 0980 if (i < 0) { 0981 c = a; 0982 a = b; 0983 b = c; 0984 i = 1; 0985 } else { 0986 i = 0; 0987 } 0988 c = Balloc(a->k); 0989 c->sign = i; 0990 wa = a->wds; 0991 xa = a->x; 0992 xae = xa + wa; 0993 wb = b->wds; 0994 xb = b->x; 0995 xbe = xb + wb; 0996 xc = c->x; 0997 borrow = 0; 0998 #ifdef ULLong 0999 do { 1000 y = (ULLong) * xa++ - *xb++ - borrow; 1001 borrow = y >> 32 & (ULong)1; 1002 *xc++ = (ULong)y & FFFFFFFF; 1003 } while (xb < xbe); 1004 while (xa < xae) { 1005 y = *xa++ - borrow; 1006 borrow = y >> 32 & (ULong)1; 1007 *xc++ = (ULong)y & FFFFFFFF; 1008 } 1009 #else 1010 #ifdef Pack_32 1011 do { 1012 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 1013 borrow = (y & 0x10000) >> 16; 1014 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 1015 borrow = (z & 0x10000) >> 16; 1016 Storeinc(xc, z, y); 1017 } while (xb < xbe); 1018 while (xa < xae) { 1019 y = (*xa & 0xffff) - borrow; 1020 borrow = (y & 0x10000) >> 16; 1021 z = (*xa++ >> 16) - borrow; 1022 borrow = (z & 0x10000) >> 16; 1023 Storeinc(xc, z, y); 1024 } 1025 #else 1026 do { 1027 y = *xa++ - *xb++ - borrow; 1028 borrow = (y & 0x10000) >> 16; 1029 *xc++ = y & 0xffff; 1030 } while (xb < xbe); 1031 while (xa < xae) { 1032 y = *xa++ - borrow; 1033 borrow = (y & 0x10000) >> 16; 1034 *xc++ = y & 0xffff; 1035 } 1036 #endif 1037 #endif 1038 while (!*--xc) { 1039 wa--; 1040 } 1041 c->wds = wa; 1042 return c; 1043 } 1044 1045 static double 1046 ulp 1047 (double dx) 1048 { 1049 Long L; 1050 U x, a; 1051 1052 dval(x) = dx; 1053 L = (word0(x) & Exp_mask) - (P - 1) * Exp_msk1; 1054 #ifndef Avoid_Underflow 1055 #ifndef Sudden_Underflow 1056 if (L > 0) { 1057 #endif 1058 #endif 1059 #ifdef IBM 1060 L |= Exp_msk1 >> 4; 1061 #endif 1062 word0(a) = L; 1063 word1(a) = 0; 1064 #ifndef Avoid_Underflow 1065 #ifndef Sudden_Underflow 1066 } else { 1067 L = -L >> Exp_shift; 1068 if (L < Exp_shift) { 1069 word0(a) = 0x80000 >> L; 1070 word1(a) = 0; 1071 } else { 1072 word0(a) = 0; 1073 L -= Exp_shift; 1074 word1(a) = L >= 31 ? 1 : 1 << 31 - L; 1075 } 1076 } 1077 #endif 1078 #endif 1079 return dval(a); 1080 } 1081 1082 static double 1083 b2d 1084 (Bigint *a, int *e) 1085 { 1086 ULong *xa, *xa0, w, y, z; 1087 int k; 1088 U d; 1089 #ifdef VAX 1090 ULong d0, d1; 1091 #else 1092 #define d0 word0(d) 1093 #define d1 word1(d) 1094 #endif 1095 1096 xa0 = a->x; 1097 xa = xa0 + a->wds; 1098 y = *--xa; 1099 #ifdef DEBUG 1100 if (!y) { 1101 Bug("zero y in b2d"); 1102 } 1103 #endif 1104 k = hi0bits(y); 1105 *e = 32 - k; 1106 #ifdef Pack_32 1107 if (k < Ebits) { 1108 d0 = Exp_1 | y >> (Ebits - k); 1109 w = xa > xa0 ? *--xa : 0; 1110 d1 = y << (32 - Ebits + k) | w >> (Ebits - k); 1111 goto ret_d; 1112 } 1113 z = xa > xa0 ? *--xa : 0; 1114 if (k -= Ebits) { 1115 d0 = Exp_1 | y << k | z >> (32 - k); 1116 y = xa > xa0 ? *--xa : 0; 1117 d1 = z << k | y >> (32 - k); 1118 } else { 1119 d0 = Exp_1 | y; 1120 d1 = z; 1121 } 1122 #else 1123 if (k < Ebits + 16) { 1124 z = xa > xa0 ? *--xa : 0; 1125 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 1126 w = xa > xa0 ? *--xa : 0; 1127 y = xa > xa0 ? *--xa : 0; 1128 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 1129 goto ret_d; 1130 } 1131 z = xa > xa0 ? *--xa : 0; 1132 w = xa > xa0 ? *--xa : 0; 1133 k -= Ebits + 16; 1134 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 1135 y = xa > xa0 ? *--xa : 0; 1136 d1 = w << k + 16 | y << k; 1137 #endif 1138 ret_d: 1139 #ifdef VAX 1140 word0(d) = d0 >> 16 | d0 << 16; 1141 word1(d) = d1 >> 16 | d1 << 16; 1142 #else 1143 #undef d0 1144 #undef d1 1145 #endif 1146 return dval(d); 1147 } 1148 1149 static Bigint * 1150 d2b 1151 (double dd, int *e, int *bits) 1152 { 1153 U d; 1154 Bigint *b; 1155 int de, k; 1156 ULong *x, y, z; 1157 #ifndef Sudden_Underflow 1158 int i; 1159 #endif 1160 #ifdef VAX 1161 ULong d0, d1; 1162 #endif 1163 dval(d) = dd; 1164 #ifdef VAX 1165 d0 = word0(d) >> 16 | word0(d) << 16; 1166 d1 = word1(d) >> 16 | word1(d) << 16; 1167 #else 1168 #define d0 word0(d) 1169 #define d1 word1(d) 1170 #endif 1171 1172 #ifdef Pack_32 1173 b = Balloc(1); 1174 #else 1175 b = Balloc(2); 1176 #endif 1177 x = b->x; 1178 1179 z = d0 & Frac_mask; 1180 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 1181 #ifdef Sudden_Underflow 1182 de = (int)(d0 >> Exp_shift); 1183 #ifndef IBM 1184 z |= Exp_msk11; 1185 #endif 1186 #else 1187 if ((de = (int)(d0 >> Exp_shift))) { 1188 z |= Exp_msk1; 1189 } 1190 #endif 1191 #ifdef Pack_32 1192 if ((y = d1)) { 1193 if ((k = lo0bits(&y))) { 1194 x[0] = y | z << (32 - k); 1195 z >>= k; 1196 } else { 1197 x[0] = y; 1198 } 1199 #ifndef Sudden_Underflow 1200 i = 1201 #endif 1202 b->wds = (x[1] = z) ? 2 : 1; 1203 } else { 1204 #ifdef DEBUG 1205 if (!z) { 1206 Bug("Zero passed to d2b"); 1207 } 1208 #endif 1209 k = lo0bits(&z); 1210 x[0] = z; 1211 #ifndef Sudden_Underflow 1212 i = 1213 #endif 1214 b->wds = 1; 1215 k += 32; 1216 } 1217 #else 1218 if (y = d1) { 1219 if (k = lo0bits(&y)) 1220 if (k >= 16) { 1221 x[0] = y | z << 32 - k & 0xffff; 1222 x[1] = z >> k - 16 & 0xffff; 1223 x[2] = z >> k; 1224 i = 2; 1225 } else { 1226 x[0] = y & 0xffff; 1227 x[1] = y >> 16 | z << 16 - k & 0xffff; 1228 x[2] = z >> k & 0xffff; 1229 x[3] = z >> k + 16; 1230 i = 3; 1231 } 1232 else { 1233 x[0] = y & 0xffff; 1234 x[1] = y >> 16; 1235 x[2] = z & 0xffff; 1236 x[3] = z >> 16; 1237 i = 3; 1238 } 1239 } else { 1240 #ifdef DEBUG 1241 if (!z) { 1242 Bug("Zero passed to d2b"); 1243 } 1244 #endif 1245 k = lo0bits(&z); 1246 if (k >= 16) { 1247 x[0] = z; 1248 i = 0; 1249 } else { 1250 x[0] = z & 0xffff; 1251 x[1] = z >> 16; 1252 i = 1; 1253 } 1254 k += 32; 1255 } 1256 while (!x[i]) { 1257 --i; 1258 } 1259 b->wds = i + 1; 1260 #endif 1261 #ifndef Sudden_Underflow 1262 if (de) { 1263 #endif 1264 #ifdef IBM 1265 *e = (de - Bias - (P - 1) << 2) + k; 1266 *bits = 4 * P + 8 - k - hi0bits(word0(d) & Frac_mask); 1267 #else 1268 *e = de - Bias - (P - 1) + k; 1269 *bits = P - k; 1270 #endif 1271 #ifndef Sudden_Underflow 1272 } else { 1273 *e = de - Bias - (P - 1) + 1 + k; 1274 #ifdef Pack_32 1275 *bits = 32 * i - hi0bits(x[i - 1]); 1276 #else 1277 *bits = (i + 2) * 16 - hi0bits(x[i]); 1278 #endif 1279 } 1280 #endif 1281 return b; 1282 } 1283 #undef d0 1284 #undef d1 1285 1286 static double 1287 ratio 1288 (Bigint *a, Bigint *b) 1289 { 1290 U da, db; 1291 int k, ka, kb; 1292 1293 dval(da) = b2d(a, &ka); 1294 dval(db) = b2d(b, &kb); 1295 #ifdef Pack_32 1296 k = ka - kb + 32 * (a->wds - b->wds); 1297 #else 1298 k = ka - kb + 16 * (a->wds - b->wds); 1299 #endif 1300 #ifdef IBM 1301 if (k > 0) { 1302 word0(da) += (k >> 2) * Exp_msk1; 1303 if (k &= 3) { 1304 dval(da) *= 1 << k; 1305 } 1306 } else { 1307 k = -k; 1308 word0(db) += (k >> 2) * Exp_msk1; 1309 if (k &= 3) { 1310 dval(db) *= 1 << k; 1311 } 1312 } 1313 #else 1314 if (k > 0) { 1315 word0(da) += k * Exp_msk1; 1316 } else { 1317 k = -k; 1318 word0(db) += k * Exp_msk1; 1319 } 1320 #endif 1321 return dval(da) / dval(db); 1322 } 1323 1324 static CONST double 1325 tens[] = { 1326 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1327 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1328 1e20, 1e21, 1e22 1329 #ifdef VAX 1330 , 1e23, 1e24 1331 #endif 1332 }; 1333 1334 static CONST double 1335 #ifdef IEEE_Arith 1336 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 1337 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1338 #ifdef Avoid_Underflow 1339 9007199254740992.*9007199254740992.e-256 1340 /* = 2^106 * 1e-53 */ 1341 #else 1342 1e-256 1343 #endif 1344 }; 1345 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 1346 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 1347 #define Scale_Bit 0x10 1348 #define n_bigtens 5 1349 #else 1350 #ifdef IBM 1351 bigtens[] = { 1e16, 1e32, 1e64 }; 1352 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 1353 #define n_bigtens 3 1354 #else 1355 bigtens[] = { 1e16, 1e32 }; 1356 static CONST double tinytens[] = { 1e-16, 1e-32 }; 1357 #define n_bigtens 2 1358 #endif 1359 #endif 1360 1361 #ifndef IEEE_Arith 1362 #undef INFNAN_CHECK 1363 #endif 1364 1365 #ifdef INFNAN_CHECK 1366 1367 #ifndef NAN_WORD0 1368 #define NAN_WORD0 0x7ff80000 1369 #endif 1370 1371 #ifndef NAN_WORD1 1372 #define NAN_WORD1 0 1373 #endif 1374 1375 static int 1376 match 1377 (CONST char **sp, CONST char *t) 1378 { 1379 int c, d; 1380 CONST char *s = *sp; 1381 1382 while ((d = *t++)) { 1383 if ((c = *++s) >= 'A' && c <= 'Z') { 1384 c += 'a' - 'A'; 1385 } 1386 if (c != d) { 1387 return 0; 1388 } 1389 } 1390 *sp = s + 1; 1391 return 1; 1392 } 1393 1394 #ifndef No_Hex_NaN 1395 static void 1396 hexnan 1397 (U *rvp, CONST char **sp) 1398 { 1399 ULong c, x[2]; 1400 CONST char *s; 1401 int havedig, udx0, xshift; 1402 1403 x[0] = x[1] = 0; 1404 havedig = xshift = 0; 1405 udx0 = 1; 1406 s = *sp; 1407 while ((c = *(CONST unsigned char *)++s)) { 1408 if (c >= '0' && c <= '9') { 1409 c -= '0'; 1410 } else if (c >= 'a' && c <= 'f') { 1411 c += 10 - 'a'; 1412 } else if (c >= 'A' && c <= 'F') { 1413 c += 10 - 'A'; 1414 } else if (c <= ' ') { 1415 if (udx0 && havedig) { 1416 udx0 = 0; 1417 xshift = 1; 1418 } 1419 continue; 1420 } else if (/*(*/ c == ')' && havedig) { 1421 *sp = s + 1; 1422 break; 1423 } else { 1424 return; /* invalid form: don't change *sp */ 1425 } 1426 havedig = 1; 1427 if (xshift) { 1428 xshift = 0; 1429 x[0] = x[1]; 1430 x[1] = 0; 1431 } 1432 if (udx0) { 1433 x[0] = (x[0] << 4) | (x[1] >> 28); 1434 } 1435 x[1] = (x[1] << 4) | c; 1436 } 1437 if ((x[0] &= 0xfffff) || x[1]) { 1438 word0(*rvp) = Exp_mask | x[0]; 1439 word1(*rvp) = x[1]; 1440 } 1441 } 1442 #endif /*No_Hex_NaN*/ 1443 #endif /* INFNAN_CHECK */ 1444 1445 double 1446 strtod 1447 (CONST char *s00, char **se) 1448 { 1449 #ifdef Avoid_Underflow 1450 int scale; 1451 #endif 1452 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, 1453 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 1454 CONST char *s, *s0, *s1; 1455 double aadj, aadj1, adj; 1456 U aadj2, rv, rv0; 1457 Long L; 1458 ULong y, z; 1459 Bigint *bb = nullptr, *bb1 = nullptr, *bd = nullptr, *bd0 = nullptr, *bs = nullptr, *delta = nullptr; 1460 #ifdef SET_INEXACT 1461 int inexact, oldinexact; 1462 #endif 1463 #ifdef Honor_FLT_ROUNDS 1464 int rounding; 1465 #endif 1466 #ifdef USE_LOCALE 1467 CONST char *s2; 1468 #endif 1469 1470 sign = nz0 = nz = 0; 1471 dval(rv) = 0.; 1472 for (s = s00;; s++) switch (*s) { 1473 case '-': 1474 sign = 1; 1475 /* no break */ 1476 case '+': 1477 if (*++s) { 1478 goto break2; 1479 } 1480 /* no break */ 1481 case 0: 1482 goto ret0; 1483 case '\t': 1484 case '\n': 1485 case '\v': 1486 case '\f': 1487 case '\r': 1488 case ' ': 1489 continue; 1490 default: 1491 goto break2; 1492 } 1493 break2: 1494 if (*s == '0') { 1495 nz0 = 1; 1496 while (*++s == '0'); 1497 if (!*s) { 1498 goto ret; 1499 } 1500 } 1501 s0 = s; 1502 y = z = 0; 1503 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 1504 if (nd < 9) { 1505 y = 10 * y + c - '0'; 1506 } else if (nd < 16) { 1507 z = 10 * z + c - '0'; 1508 } 1509 nd0 = nd; 1510 #ifdef USE_LOCALE 1511 s1 = localeconv()->decimal_point; 1512 if (c == *s1) { 1513 c = '.'; 1514 if (*++s1) { 1515 s2 = s; 1516 for (;;) { 1517 if (*++s2 != *s1) { 1518 c = 0; 1519 break; 1520 } 1521 if (!*++s1) { 1522 s = s2; 1523 break; 1524 } 1525 } 1526 } 1527 } 1528 #endif 1529 if (c == '.') { 1530 c = *++s; 1531 if (!nd) { 1532 for (; c == '0'; c = *++s) { 1533 nz++; 1534 } 1535 if (c > '0' && c <= '9') { 1536 s0 = s; 1537 nf += nz; 1538 nz = 0; 1539 goto have_dig; 1540 } 1541 goto dig_done; 1542 } 1543 for (; c >= '0' && c <= '9'; c = *++s) { 1544 have_dig: 1545 nz++; 1546 if (c -= '0') { 1547 nf += nz; 1548 for (i = 1; i < nz; i++) 1549 if (nd++ < 9) { 1550 y *= 10; 1551 } else if (nd <= DBL_DIG + 1) { 1552 z *= 10; 1553 } 1554 if (nd++ < 9) { 1555 y = 10 * y + c; 1556 } else if (nd <= DBL_DIG + 1) { 1557 z = 10 * z + c; 1558 } 1559 nz = 0; 1560 } 1561 } 1562 } 1563 dig_done: 1564 e = 0; 1565 if (c == 'e' || c == 'E') { 1566 if (!nd && !nz && !nz0) { 1567 goto ret0; 1568 } 1569 s00 = s; 1570 esign = 0; 1571 switch (c = *++s) { 1572 case '-': 1573 esign = 1; 1574 case '+': 1575 c = *++s; 1576 } 1577 if (c >= '0' && c <= '9') { 1578 while (c == '0') { 1579 c = *++s; 1580 } 1581 if (c > '0' && c <= '9') { 1582 L = c - '0'; 1583 s1 = s; 1584 while ((c = *++s) >= '0' && c <= '9') { 1585 L = 10 * L + c - '0'; 1586 } 1587 if (s - s1 > 8 || L > 19999) 1588 /* Avoid confusion from exponents 1589 * so large that e might overflow. 1590 */ 1591 { 1592 e = 19999; /* safe for 16 bit ints */ 1593 } else { 1594 e = (int)L; 1595 } 1596 if (esign) { 1597 e = -e; 1598 } 1599 } else { 1600 e = 0; 1601 } 1602 } else { 1603 s = s00; 1604 } 1605 } 1606 if (!nd) { 1607 if (!nz && !nz0) { 1608 #ifdef INFNAN_CHECK 1609 /* Check for Nan and Infinity */ 1610 switch (c) { 1611 case 'i': 1612 case 'I': 1613 if (match(&s, "nf")) { 1614 --s; 1615 if (!match(&s, "inity")) { 1616 ++s; 1617 } 1618 word0(rv) = 0x7ff00000; 1619 word1(rv) = 0; 1620 goto ret; 1621 } 1622 break; 1623 case 'n': 1624 case 'N': 1625 if (match(&s, "an")) { 1626 word0(rv) = NAN_WORD0; 1627 word1(rv) = NAN_WORD1; 1628 #ifndef No_Hex_NaN 1629 if (*s == '(') { /*)*/ 1630 hexnan(&rv, &s); 1631 } 1632 #endif 1633 goto ret; 1634 } 1635 } 1636 #endif /* INFNAN_CHECK */ 1637 ret0: 1638 s = s00; 1639 sign = 0; 1640 } 1641 goto ret; 1642 } 1643 e1 = e -= nf; 1644 1645 /* Now we have nd0 digits, starting at s0, followed by a 1646 * decimal point, followed by nd-nd0 digits. The number we're 1647 * after is the integer represented by those digits times 1648 * 10**e */ 1649 1650 if (!nd0) { 1651 nd0 = nd; 1652 } 1653 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 1654 dval(rv) = y; 1655 if (k > 9) { 1656 #ifdef SET_INEXACT 1657 if (k > DBL_DIG) { 1658 oldinexact = get_inexact(); 1659 } 1660 #endif 1661 dval(rv) = tens[k - 9] * dval(rv) + z; 1662 } 1663 bd0 = nullptr; 1664 if (nd <= DBL_DIG 1665 #ifndef RND_PRODQUOT 1666 #ifndef Honor_FLT_ROUNDS 1667 && Flt_Rounds == 1 1668 #endif 1669 #endif 1670 ) { 1671 if (!e) { 1672 goto ret; 1673 } 1674 if (e > 0) { 1675 if (e <= Ten_pmax) { 1676 #ifdef VAX 1677 goto vax_ovfl_check; 1678 #else 1679 #ifdef Honor_FLT_ROUNDS 1680 /* round correctly FLT_ROUNDS = 2 or 3 */ 1681 if (sign) { 1682 rv = -rv; 1683 sign = 0; 1684 } 1685 #endif 1686 /* rv = */ rounded_product(dval(rv), tens[e]); 1687 goto ret; 1688 #endif 1689 } 1690 i = DBL_DIG - nd; 1691 if (e <= Ten_pmax + i) { 1692 /* A fancier test would sometimes let us do 1693 * this for larger i values. 1694 */ 1695 #ifdef Honor_FLT_ROUNDS 1696 /* round correctly FLT_ROUNDS = 2 or 3 */ 1697 if (sign) { 1698 rv = -rv; 1699 sign = 0; 1700 } 1701 #endif 1702 e -= i; 1703 dval(rv) *= tens[i]; 1704 #ifdef VAX 1705 /* VAX exponent range is so narrow we must 1706 * worry about overflow here... 1707 */ 1708 vax_ovfl_check: 1709 word0(rv) -= P * Exp_msk1; 1710 /* rv = */ rounded_product(dval(rv), tens[e]); 1711 if ((word0(rv) & Exp_mask) 1712 > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) { 1713 goto ovfl; 1714 } 1715 word0(rv) += P * Exp_msk1; 1716 #else 1717 /* rv = */ rounded_product(dval(rv), tens[e]); 1718 #endif 1719 goto ret; 1720 } 1721 } 1722 #ifndef Inaccurate_Divide 1723 else if (e >= -Ten_pmax) { 1724 #ifdef Honor_FLT_ROUNDS 1725 /* round correctly FLT_ROUNDS = 2 or 3 */ 1726 if (sign) { 1727 rv = -rv; 1728 sign = 0; 1729 } 1730 #endif 1731 /* rv = */ rounded_quotient(dval(rv), tens[-e]); 1732 goto ret; 1733 } 1734 #endif 1735 } 1736 e1 += nd - k; 1737 1738 #ifdef IEEE_Arith 1739 #ifdef SET_INEXACT 1740 inexact = 1; 1741 if (k <= DBL_DIG) { 1742 oldinexact = get_inexact(); 1743 } 1744 #endif 1745 #ifdef Avoid_Underflow 1746 scale = 0; 1747 #endif 1748 #ifdef Honor_FLT_ROUNDS 1749 if ((rounding = Flt_Rounds) >= 2) { 1750 if (sign) { 1751 rounding = rounding == 2 ? 0 : 2; 1752 } else if (rounding != 2) { 1753 rounding = 0; 1754 } 1755 } 1756 #endif 1757 #endif /*IEEE_Arith*/ 1758 1759 /* Get starting approximation = rv * 10**e1 */ 1760 1761 if (e1 > 0) { 1762 if ((i = e1 & 15)) { 1763 dval(rv) *= tens[i]; 1764 } 1765 if (e1 &= ~15) { 1766 if (e1 > DBL_MAX_10_EXP) { 1767 ovfl: 1768 #if HAVE_ERRNO_H 1769 errno = ERANGE; 1770 #endif 1771 /* Can't trust HUGE_VAL */ 1772 #ifdef IEEE_Arith 1773 #ifdef Honor_FLT_ROUNDS 1774 switch (rounding) { 1775 case 0: /* toward 0 */ 1776 case 3: /* toward -infinity */ 1777 word0(rv) = Big0; 1778 word1(rv) = Big1; 1779 break; 1780 default: 1781 word0(rv) = Exp_mask; 1782 word1(rv) = 0; 1783 } 1784 #else /*Honor_FLT_ROUNDS*/ 1785 word0(rv) = Exp_mask; 1786 word1(rv) = 0; 1787 #endif /*Honor_FLT_ROUNDS*/ 1788 #ifdef SET_INEXACT 1789 /* set overflow bit */ 1790 dval(rv0) = 1e300; 1791 dval(rv0) *= dval(rv0); 1792 #endif 1793 #else /*IEEE_Arith*/ 1794 word0(rv) = Big0; 1795 word1(rv) = Big1; 1796 #endif /*IEEE_Arith*/ 1797 if (bd0) { 1798 goto retfree; 1799 } 1800 goto ret; 1801 } 1802 e1 >>= 4; 1803 for (j = 0; e1 > 1; j++, e1 >>= 1) 1804 if (e1 & 1) { 1805 dval(rv) *= bigtens[j]; 1806 } 1807 /* The last multiplication could overflow. */ 1808 word0(rv) -= P * Exp_msk1; 1809 dval(rv) *= bigtens[j]; 1810 if ((z = word0(rv) & Exp_mask) 1811 > Exp_msk1 * (DBL_MAX_EXP + Bias - P)) { 1812 goto ovfl; 1813 } 1814 if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) { 1815 /* set to largest number */ 1816 /* (Can't trust DBL_MAX) */ 1817 word0(rv) = Big0; 1818 word1(rv) = Big1; 1819 } else { 1820 word0(rv) += P * Exp_msk1; 1821 } 1822 } 1823 } else if (e1 < 0) { 1824 e1 = -e1; 1825 if ((i = e1 & 15)) { 1826 dval(rv) /= tens[i]; 1827 } 1828 if (e1 >>= 4) { 1829 if (e1 >= 1 << n_bigtens) { 1830 goto undfl; 1831 } 1832 #ifdef Avoid_Underflow 1833 if (e1 & Scale_Bit) { 1834 scale = 2 * P; 1835 } 1836 for (j = 0; e1 > 0; j++, e1 >>= 1) 1837 if (e1 & 1) { 1838 dval(rv) *= tinytens[j]; 1839 } 1840 if (scale && (j = 2 * P + 1 - ((word0(rv) & Exp_mask) 1841 >> Exp_shift)) > 0) { 1842 /* scaled rv is denormal; zap j low bits */ 1843 if (j >= 32) { 1844 word1(rv) = 0; 1845 if (j >= 53) { 1846 word0(rv) = (P + 2) * Exp_msk1; 1847 } else { 1848 word0(rv) &= 0xffffffff << (j - 32); 1849 } 1850 } else { 1851 word1(rv) &= 0xffffffff << j; 1852 } 1853 } 1854 #else 1855 for (j = 0; e1 > 1; j++, e1 >>= 1) 1856 if (e1 & 1) { 1857 dval(rv) *= tinytens[j]; 1858 } 1859 /* The last multiplication could underflow. */ 1860 dval(rv0) = dval(rv); 1861 dval(rv) *= tinytens[j]; 1862 if (!dval(rv)) { 1863 dval(rv) = 2.*dval(rv0); 1864 dval(rv) *= tinytens[j]; 1865 #endif 1866 if (!dval(rv)) { 1867 undfl: 1868 dval(rv) = 0.; 1869 #if HAVE_ERRNO_H 1870 errno = ERANGE; 1871 #endif 1872 if (bd0) { 1873 goto retfree; 1874 } 1875 goto ret; 1876 } 1877 #ifndef Avoid_Underflow 1878 word0(rv) = Tiny0; 1879 word1(rv) = Tiny1; 1880 /* The refinement below will clean 1881 * this approximation up. 1882 */ 1883 } 1884 #endif 1885 } 1886 } 1887 1888 /* Now the hard part -- adjusting rv to the correct value.*/ 1889 1890 /* Put digits into bd: true value = bd * 10^e */ 1891 1892 bd0 = s2b(s0, nd0, nd, y); 1893 1894 for (;;) { 1895 bd = Balloc(bd0->k); 1896 Bcopy(bd, bd0); 1897 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 1898 bs = i2b(1); 1899 1900 if (e >= 0) { 1901 bb2 = bb5 = 0; 1902 bd2 = bd5 = e; 1903 } else { 1904 bb2 = bb5 = -e; 1905 bd2 = bd5 = 0; 1906 } 1907 if (bbe >= 0) { 1908 bb2 += bbe; 1909 } else { 1910 bd2 -= bbe; 1911 } 1912 bs2 = bb2; 1913 #ifdef Honor_FLT_ROUNDS 1914 if (rounding != 1) { 1915 bs2++; 1916 } 1917 #endif 1918 #ifdef Avoid_Underflow 1919 j = bbe - scale; 1920 i = j + bbbits - 1; /* logb(rv) */ 1921 if (i < Emin) { /* denormal */ 1922 j += P - Emin; 1923 } else { 1924 j = P + 1 - bbbits; 1925 } 1926 #else /*Avoid_Underflow*/ 1927 #ifdef Sudden_Underflow 1928 #ifdef IBM 1929 j = 1 + 4 * P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 1930 #else 1931 j = P + 1 - bbbits; 1932 #endif 1933 #else /*Sudden_Underflow*/ 1934 j = bbe; 1935 i = j + bbbits - 1; /* logb(rv) */ 1936 if (i < Emin) { /* denormal */ 1937 j += P - Emin; 1938 } else { 1939 j = P + 1 - bbbits; 1940 } 1941 #endif /*Sudden_Underflow*/ 1942 #endif /*Avoid_Underflow*/ 1943 bb2 += j; 1944 bd2 += j; 1945 #ifdef Avoid_Underflow 1946 bd2 += scale; 1947 #endif 1948 i = bb2 < bd2 ? bb2 : bd2; 1949 if (i > bs2) { 1950 i = bs2; 1951 } 1952 if (i > 0) { 1953 bb2 -= i; 1954 bd2 -= i; 1955 bs2 -= i; 1956 } 1957 if (bb5 > 0) { 1958 bs = pow5mult(bs, bb5); 1959 bb1 = mult(bs, bb); 1960 Bfree(bb); 1961 bb = bb1; 1962 } 1963 if (bb2 > 0) { 1964 bb = lshift(bb, bb2); 1965 } 1966 if (bd5 > 0) { 1967 bd = pow5mult(bd, bd5); 1968 } 1969 if (bd2 > 0) { 1970 bd = lshift(bd, bd2); 1971 } 1972 if (bs2 > 0) { 1973 bs = lshift(bs, bs2); 1974 } 1975 delta = diff(bb, bd); 1976 dsign = delta->sign; 1977 delta->sign = 0; 1978 i = cmp(delta, bs); 1979 #ifdef Honor_FLT_ROUNDS 1980 if (rounding != 1) { 1981 if (i < 0) { 1982 /* Error is less than an ulp */ 1983 if (!delta->x[0] && delta->wds <= 1) { 1984 /* exact */ 1985 #ifdef SET_INEXACT 1986 inexact = 0; 1987 #endif 1988 break; 1989 } 1990 if (rounding) { 1991 if (dsign) { 1992 adj = 1.; 1993 goto apply_adj; 1994 } 1995 } else if (!dsign) { 1996 adj = -1.; 1997 if (!word1(rv) 1998 && !(word0(rv) & Frac_mask)) { 1999 y = word0(rv) & Exp_mask; 2000 #ifdef Avoid_Underflow 2001 if (!scale || y > 2 * P * Exp_msk1) 2002 #else 2003 if (y) 2004 #endif 2005 { 2006 delta = lshift(delta, Log2P); 2007 if (cmp(delta, bs) <= 0) { 2008 adj = -0.5; 2009 } 2010 } 2011 } 2012 apply_adj: 2013 #ifdef Avoid_Underflow 2014 if (scale && (y = word0(rv) & Exp_mask) 2015 <= 2 * P * Exp_msk1) { 2016 word0(adj) += (2 * P + 1) * Exp_msk1 - y; 2017 } 2018 #else 2019 #ifdef Sudden_Underflow 2020 if ((word0(rv) & Exp_mask) <= 2021 P * Exp_msk1) { 2022 word0(rv) += P * Exp_msk1; 2023 dval(rv) += adj * ulp(dval(rv)); 2024 word0(rv) -= P * Exp_msk1; 2025 } else 2026 #endif /*Sudden_Underflow*/ 2027 #endif /*Avoid_Underflow*/ 2028 dval(rv) += adj * ulp(dval(rv)); 2029 } 2030 break; 2031 } 2032 adj = ratio(delta, bs); 2033 if (adj < 1.) { 2034 adj = 1.; 2035 } 2036 if (adj <= 0x7ffffffe) { 2037 /* adj = rounding ? ceil(adj) : floor(adj); */ 2038 y = adj; 2039 if (y != adj) { 2040 if (!((rounding >> 1) ^ dsign)) { 2041 y++; 2042 } 2043 adj = y; 2044 } 2045 } 2046 #ifdef Avoid_Underflow 2047 if (scale && (y = word0(rv) & Exp_mask) <= 2 * P * Exp_msk1) { 2048 word0(adj) += (2 * P + 1) * Exp_msk1 - y; 2049 } 2050 #else 2051 #ifdef Sudden_Underflow 2052 if ((word0(rv) & Exp_mask) <= P * Exp_msk1) { 2053 word0(rv) += P * Exp_msk1; 2054 adj *= ulp(dval(rv)); 2055 if (dsign) { 2056 dval(rv) += adj; 2057 } else { 2058 dval(rv) -= adj; 2059 } 2060 word0(rv) -= P * Exp_msk1; 2061 goto cont; 2062 } 2063 #endif /*Sudden_Underflow*/ 2064 #endif /*Avoid_Underflow*/ 2065 adj *= ulp(dval(rv)); 2066 if (dsign) { 2067 dval(rv) += adj; 2068 } else { 2069 dval(rv) -= adj; 2070 } 2071 goto cont; 2072 } 2073 #endif /*Honor_FLT_ROUNDS*/ 2074 2075 if (i < 0) { 2076 /* Error is less than half an ulp -- check for 2077 * special case of mantissa a power of two. 2078 */ 2079 if (dsign || word1(rv) || word0(rv) & Bndry_mask 2080 #ifdef IEEE_Arith 2081 #ifdef Avoid_Underflow 2082 || (word0(rv) & Exp_mask) <= (2 * P + 1)*Exp_msk1 2083 #else 2084 || (word0(rv) & Exp_mask) <= Exp_msk1 2085 #endif 2086 #endif 2087 ) { 2088 #ifdef SET_INEXACT 2089 if (!delta->x[0] && delta->wds <= 1) { 2090 inexact = 0; 2091 } 2092 #endif 2093 break; 2094 } 2095 if (!delta->x[0] && delta->wds <= 1) { 2096 /* exact result */ 2097 #ifdef SET_INEXACT 2098 inexact = 0; 2099 #endif 2100 break; 2101 } 2102 delta = lshift(delta, Log2P); 2103 if (cmp(delta, bs) > 0) { 2104 goto drop_down; 2105 } 2106 break; 2107 } 2108 if (i == 0) { 2109 /* exactly half-way between */ 2110 if (dsign) { 2111 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 2112 && word1(rv) == ( 2113 #ifdef Avoid_Underflow 2114 (scale && (y = word0(rv) & Exp_mask) <= 2 * P * Exp_msk1) 2115 ? (0xffffffff & (0xffffffff << (2 * P + 1 - (y >> Exp_shift)))) : 2116 #endif 2117 0xffffffff)) { 2118 /*boundary case -- increment exponent*/ 2119 word0(rv) = (word0(rv) & Exp_mask) 2120 + Exp_msk1 2121 #ifdef IBM 2122 | Exp_msk1 >> 4 2123 #endif 2124 ; 2125 word1(rv) = 0; 2126 #ifdef Avoid_Underflow 2127 dsign = 0; 2128 #endif 2129 break; 2130 } 2131 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 2132 drop_down: 2133 /* boundary case -- decrement exponent */ 2134 #ifdef Sudden_Underflow /*{{*/ 2135 L = word0(rv) & Exp_mask; 2136 #ifdef IBM 2137 if (L < Exp_msk1) 2138 #else 2139 #ifdef Avoid_Underflow 2140 if (L <= (scale ? (2 * P + 1)*Exp_msk1 : Exp_msk1)) 2141 #else 2142 if (L <= Exp_msk1) 2143 #endif /*Avoid_Underflow*/ 2144 #endif /*IBM*/ 2145 goto undfl; 2146 L -= Exp_msk1; 2147 #else /*Sudden_Underflow}{*/ 2148 #ifdef Avoid_Underflow 2149 if (scale) { 2150 L = word0(rv) & Exp_mask; 2151 if (L <= (2 * P + 1)*Exp_msk1) { 2152 if (L > (P + 2)*Exp_msk1) 2153 /* round even ==> */ 2154 /* accept rv */ 2155 { 2156 break; 2157 } 2158 /* rv = smallest denormal */ 2159 goto undfl; 2160 } 2161 } 2162 #endif /*Avoid_Underflow*/ 2163 L = (word0(rv) & Exp_mask) - Exp_msk1; 2164 #endif /*Sudden_Underflow}}*/ 2165 word0(rv) = L | Bndry_mask1; 2166 word1(rv) = 0xffffffff; 2167 #ifdef IBM 2168 goto cont; 2169 #else 2170 break; 2171 #endif 2172 } 2173 #ifndef ROUND_BIASED 2174 if (!(word1(rv) & LSB)) { 2175 break; 2176 } 2177 #endif 2178 if (dsign) { 2179 dval(rv) += ulp(dval(rv)); 2180 } 2181 #ifndef ROUND_BIASED 2182 else { 2183 dval(rv) -= ulp(dval(rv)); 2184 #ifndef Sudden_Underflow 2185 if (!dval(rv)) { 2186 goto undfl; 2187 } 2188 #endif 2189 } 2190 #ifdef Avoid_Underflow 2191 dsign = 1 - dsign; 2192 #endif 2193 #endif 2194 break; 2195 } 2196 if ((aadj = ratio(delta, bs)) <= 2.) { 2197 if (dsign) { 2198 aadj = aadj1 = 1.; 2199 } else if (word1(rv) || word0(rv) & Bndry_mask) { 2200 #ifndef Sudden_Underflow 2201 if (word1(rv) == Tiny1 && !word0(rv)) { 2202 goto undfl; 2203 } 2204 #endif 2205 aadj = 1.; 2206 aadj1 = -1.; 2207 } else { 2208 /* special case -- power of FLT_RADIX to be */ 2209 /* rounded down... */ 2210 2211 if (aadj < 2. / FLT_RADIX) { 2212 aadj = 1. / FLT_RADIX; 2213 } else { 2214 aadj *= 0.5; 2215 } 2216 aadj1 = -aadj; 2217 } 2218 } else { 2219 aadj *= 0.5; 2220 aadj1 = dsign ? aadj : -aadj; 2221 #ifdef Check_FLT_ROUNDS 2222 switch (Rounding) { 2223 case 2: /* towards +infinity */ 2224 aadj1 -= 0.5; 2225 break; 2226 case 0: /* towards 0 */ 2227 case 3: /* towards -infinity */ 2228 aadj1 += 0.5; 2229 } 2230 #else 2231 if (Flt_Rounds == 0) { 2232 aadj1 += 0.5; 2233 } 2234 #endif /*Check_FLT_ROUNDS*/ 2235 } 2236 y = word0(rv) & Exp_mask; 2237 2238 /* Check for overflow */ 2239 2240 if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1)) { 2241 dval(rv0) = dval(rv); 2242 word0(rv) -= P * Exp_msk1; 2243 adj = aadj1 * ulp(dval(rv)); 2244 dval(rv) += adj; 2245 if ((word0(rv) & Exp_mask) >= 2246 Exp_msk1 * (DBL_MAX_EXP + Bias - P)) { 2247 if (word0(rv0) == Big0 && word1(rv0) == Big1) { 2248 goto ovfl; 2249 } 2250 word0(rv) = Big0; 2251 word1(rv) = Big1; 2252 goto cont; 2253 } else { 2254 word0(rv) += P * Exp_msk1; 2255 } 2256 } else { 2257 #ifdef Avoid_Underflow 2258 if (scale && y <= 2 * P * Exp_msk1) { 2259 if (aadj <= 0x7fffffff) { 2260 if ((z = (ULong)aadj) <= 0) { 2261 z = 1; 2262 } 2263 aadj = z; 2264 aadj1 = dsign ? aadj : -aadj; 2265 } 2266 dval(aadj2) = aadj1; 2267 word0(aadj2) += (2 * P + 1) * Exp_msk1 - y; 2268 aadj1 = dval(aadj2); 2269 } 2270 adj = aadj1 * ulp(dval(rv)); 2271 dval(rv) += adj; 2272 #else 2273 #ifdef Sudden_Underflow 2274 if ((word0(rv) & Exp_mask) <= P * Exp_msk1) { 2275 dval(rv0) = dval(rv); 2276 word0(rv) += P * Exp_msk1; 2277 adj = aadj1 * ulp(dval(rv)); 2278 dval(rv) += adj; 2279 #ifdef IBM 2280 if ((word0(rv) & Exp_mask) < P * Exp_msk1) 2281 #else 2282 if ((word0(rv) & Exp_mask) <= P * Exp_msk1) 2283 #endif 2284 { 2285 if (word0(rv0) == Tiny0 2286 && word1(rv0) == Tiny1) { 2287 goto undfl; 2288 } 2289 word0(rv) = Tiny0; 2290 word1(rv) = Tiny1; 2291 goto cont; 2292 } else { 2293 word0(rv) -= P * Exp_msk1; 2294 } 2295 } else { 2296 adj = aadj1 * ulp(dval(rv)); 2297 dval(rv) += adj; 2298 } 2299 #else /*Sudden_Underflow*/ 2300 /* Compute adj so that the IEEE rounding rules will 2301 * correctly round rv + adj in some half-way cases. 2302 * If rv * ulp(rv) is denormalized (i.e., 2303 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 2304 * trouble from bits lost to denormalization; 2305 * example: 1.2e-307 . 2306 */ 2307 if (y <= (P - 1)*Exp_msk1 && aadj > 1.) { 2308 aadj1 = (double)(int)(aadj + 0.5); 2309 if (!dsign) { 2310 aadj1 = -aadj1; 2311 } 2312 } 2313 adj = aadj1 * ulp(dval(rv)); 2314 dval(rv) += adj; 2315 #endif /*Sudden_Underflow*/ 2316 #endif /*Avoid_Underflow*/ 2317 } 2318 z = word0(rv) & Exp_mask; 2319 #ifndef SET_INEXACT 2320 #ifdef Avoid_Underflow 2321 if (!scale) 2322 #endif 2323 if (y == z) { 2324 /* Can we stop now? */ 2325 L = (Long)aadj; 2326 aadj -= L; 2327 /* The tolerances below are conservative. */ 2328 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 2329 if (aadj < .4999999 || aadj > .5000001) { 2330 break; 2331 } 2332 } else if (aadj < .4999999 / FLT_RADIX) { 2333 break; 2334 } 2335 } 2336 #endif 2337 cont: 2338 Bfree(bb); 2339 Bfree(bd); 2340 Bfree(bs); 2341 Bfree(delta); 2342 } 2343 #ifdef SET_INEXACT 2344 if (inexact) { 2345 if (!oldinexact) { 2346 word0(rv0) = Exp_1 + (70 << Exp_shift); 2347 word1(rv0) = 0; 2348 dval(rv0) += 1.; 2349 } 2350 } else if (!oldinexact) { 2351 clear_inexact(); 2352 } 2353 #endif 2354 #ifdef Avoid_Underflow 2355 if (scale) { 2356 word0(rv0) = Exp_1 - 2 * P * Exp_msk1; 2357 word1(rv0) = 0; 2358 dval(rv) *= dval(rv0); 2359 #if HAVE_ERRNO_H 2360 /* try to avoid the bug of testing an 8087 register value */ 2361 if (word0(rv) == 0 && word1(rv) == 0) { 2362 errno = ERANGE; 2363 } 2364 #endif 2365 } 2366 #endif /* Avoid_Underflow */ 2367 #ifdef SET_INEXACT 2368 if (inexact && !(word0(rv) & Exp_mask)) { 2369 /* set underflow bit */ 2370 dval(rv0) = 1e-300; 2371 dval(rv0) *= dval(rv0); 2372 } 2373 #endif 2374 retfree: 2375 Bfree(bb); 2376 Bfree(bd); 2377 Bfree(bs); 2378 Bfree(bd0); 2379 Bfree(delta); 2380 ret: 2381 if (se) { 2382 *se = (char *)s; 2383 } 2384 return sign ? -dval(rv) : dval(rv); 2385 } 2386 2387 static int 2388 quorem 2389 (Bigint *b, Bigint *S) 2390 { 2391 int n; 2392 ULong *bx, *bxe, q, *sx, *sxe; 2393 #ifdef ULLong 2394 ULLong borrow, carry, y, ys; 2395 #else 2396 ULong borrow, carry, y, ys; 2397 #ifdef Pack_32 2398 ULong si, z, zs; 2399 #endif 2400 #endif 2401 2402 n = S->wds; 2403 #ifdef DEBUG 2404 /*debug*/ if (b->wds > n) 2405 /*debug*/{ 2406 Bug("oversize b in quorem"); 2407 } 2408 #endif 2409 if (b->wds < n) { 2410 return 0; 2411 } 2412 sx = S->x; 2413 sxe = sx + --n; 2414 bx = b->x; 2415 bxe = bx + n; 2416 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 2417 #ifdef DEBUG 2418 /*debug*/ if (q > 9) 2419 /*debug*/{ 2420 Bug("oversized quotient in quorem"); 2421 } 2422 #endif 2423 if (q) { 2424 borrow = 0; 2425 carry = 0; 2426 do { 2427 #ifdef ULLong 2428 ys = *sx++ * (ULLong)q + carry; 2429 carry = ys >> 32; 2430 y = *bx - (ys & FFFFFFFF) - borrow; 2431 borrow = y >> 32 & (ULong)1; 2432 *bx++ = (ULong)y & FFFFFFFF; 2433 #else 2434 #ifdef Pack_32 2435 si = *sx++; 2436 ys = (si & 0xffff) * q + carry; 2437 zs = (si >> 16) * q + (ys >> 16); 2438 carry = zs >> 16; 2439 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 2440 borrow = (y & 0x10000) >> 16; 2441 z = (*bx >> 16) - (zs & 0xffff) - borrow; 2442 borrow = (z & 0x10000) >> 16; 2443 Storeinc(bx, z, y); 2444 #else 2445 ys = *sx++ * q + carry; 2446 carry = ys >> 16; 2447 y = *bx - (ys & 0xffff) - borrow; 2448 borrow = (y & 0x10000) >> 16; 2449 *bx++ = y & 0xffff; 2450 #endif 2451 #endif 2452 } while (sx <= sxe); 2453 if (!*bxe) { 2454 bx = b->x; 2455 while (--bxe > bx && !*bxe) { 2456 --n; 2457 } 2458 b->wds = n; 2459 } 2460 } 2461 if (cmp(b, S) >= 0) { 2462 q++; 2463 borrow = 0; 2464 carry = 0; 2465 bx = b->x; 2466 sx = S->x; 2467 do { 2468 #ifdef ULLong 2469 ys = *sx++ + carry; 2470 carry = ys >> 32; 2471 y = *bx - (ys & FFFFFFFF) - borrow; 2472 borrow = y >> 32 & (ULong)1; 2473 *bx++ = (ULong)y & FFFFFFFF; 2474 #else 2475 #ifdef Pack_32 2476 si = *sx++; 2477 ys = (si & 0xffff) + carry; 2478 zs = (si >> 16) + (ys >> 16); 2479 carry = zs >> 16; 2480 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 2481 borrow = (y & 0x10000) >> 16; 2482 z = (*bx >> 16) - (zs & 0xffff) - borrow; 2483 borrow = (z & 0x10000) >> 16; 2484 Storeinc(bx, z, y); 2485 #else 2486 ys = *sx++ + carry; 2487 carry = ys >> 16; 2488 y = *bx - (ys & 0xffff) - borrow; 2489 borrow = (y & 0x10000) >> 16; 2490 *bx++ = y & 0xffff; 2491 #endif 2492 #endif 2493 } while (sx <= sxe); 2494 bx = b->x; 2495 bxe = bx + n; 2496 if (!*bxe) { 2497 while (--bxe > bx && !*bxe) { 2498 --n; 2499 } 2500 b->wds = n; 2501 } 2502 } 2503 return q; 2504 } 2505 2506 #ifndef MULTIPLE_THREADS 2507 static char *dtoa_result; 2508 #endif 2509 2510 static char * 2511 rv_alloc(int i) 2512 { 2513 int j, k, *r; 2514 2515 j = sizeof(ULong); 2516 for (k = 0; 2517 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i; 2518 j <<= 1) { 2519 k++; 2520 } 2521 r = (int *)Balloc(k); 2522 *r = k; 2523 return 2524 #ifndef MULTIPLE_THREADS 2525 dtoa_result = 2526 #endif 2527 (char *)(r + 1); 2528 } 2529 2530 static char * 2531 nrv_alloc(CONST char *s, char **rve, int n) 2532 { 2533 char *rv, *t; 2534 2535 t = rv = rv_alloc(n); 2536 while ((*t = *s++)) { 2537 t++; 2538 } 2539 if (rve) { 2540 *rve = t; 2541 } 2542 return rv; 2543 } 2544 2545 /* freedtoa(s) must be used to free values s returned by dtoa 2546 * when MULTIPLE_THREADS is #defined. It should be used in all cases, 2547 * but for consistency with earlier versions of dtoa, it is optional 2548 * when MULTIPLE_THREADS is not defined. 2549 */ 2550 2551 void 2552 freedtoa(char *s) 2553 { 2554 Bigint *b = (Bigint *)((int *)s - 1); 2555 b->maxwds = 1 << (b->k = *(int *)b); 2556 Bfree(b); 2557 #ifndef MULTIPLE_THREADS 2558 if (s == dtoa_result) { 2559 dtoa_result = nullptr; 2560 } 2561 #endif 2562 } 2563 2564 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 2565 * 2566 * Inspired by "How to Print Floating-Point Numbers Accurately" by 2567 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101]. 2568 * 2569 * Modifications: 2570 * 1. Rather than iterating, we use a simple numeric overestimate 2571 * to determine k = floor(log10(d)). We scale relevant 2572 * quantities using O(log2(k)) rather than O(k) multiplications. 2573 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 2574 * try to generate digits strictly left to right. Instead, we 2575 * compute with fewer bits and propagate the carry if necessary 2576 * when rounding the final digit up. This is often faster. 2577 * 3. Under the assumption that input will be rounded nearest, 2578 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 2579 * That is, we allow equality in stopping tests when the 2580 * round-nearest rule will give the same floating-point value 2581 * as would satisfaction of the stopping test with strict 2582 * inequality. 2583 * 4. We remove common factors of powers of 2 from relevant 2584 * quantities. 2585 * 5. When converting floating-point integers less than 1e16, 2586 * we use floating-point arithmetic rather than resorting 2587 * to multiple-precision integers. 2588 * 6. When asked to produce fewer than 15 digits, we first try 2589 * to get by with floating-point arithmetic; we resort to 2590 * multiple-precision integer arithmetic only if we cannot 2591 * guarantee that the floating-point calculation has given 2592 * the correctly rounded result. For k requested digits and 2593 * "uniformly" distributed input, the probability is 2594 * something like 10^(k-15) that we must resort to the Long 2595 * calculation. 2596 */ 2597 2598 char * 2599 dtoa 2600 (double dd, int mode, int ndigits, int *decpt, int *sign, char **rve) 2601 { 2602 /* Arguments ndigits, decpt, sign are similar to those 2603 of ecvt and fcvt; trailing zeros are suppressed from 2604 the returned string. If not null, *rve is set to point 2605 to the end of the return value. If d is +-Infinity or NaN, 2606 then *decpt is set to 9999. 2607 2608 mode: 2609 0 ==> shortest string that yields d when read in 2610 and rounded to nearest. 2611 1 ==> like 0, but with Steele & White stopping rule; 2612 e.g. with IEEE P754 arithmetic , mode 0 gives 2613 1e23 whereas mode 1 gives 9.999999999999999e22. 2614 2 ==> max(1,ndigits) significant digits. This gives a 2615 return value similar to that of ecvt, except 2616 that trailing zeros are suppressed. 2617 3 ==> through ndigits past the decimal point. This 2618 gives a return value similar to that from fcvt, 2619 except that trailing zeros are suppressed, and 2620 ndigits can be negative. 2621 4,5 ==> similar to 2 and 3, respectively, but (in 2622 round-nearest mode) with the tests of mode 0 to 2623 possibly return a shorter string that rounds to d. 2624 With IEEE arithmetic and compilation with 2625 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same 2626 as modes 2 and 3 when FLT_ROUNDS != 1. 2627 6-9 ==> Debugging modes similar to mode - 4: don't try 2628 fast floating-point estimate (if applicable). 2629 2630 Values of mode other than 0-9 are treated as mode 0. 2631 2632 Sufficient space is allocated to the return value 2633 to hold the suppressed trailing zeros. 2634 */ 2635 2636 int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0, 2637 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, 2638 spec_case, try_quick; 2639 Long L; 2640 #ifndef Sudden_Underflow 2641 int denorm; 2642 ULong x; 2643 #endif 2644 Bigint *b, *b1, *delta, *mlo = nullptr, *mhi, *S; 2645 U d, d2, eps; 2646 double ds; 2647 char *s, *s0; 2648 #ifdef Honor_FLT_ROUNDS 2649 int rounding; 2650 #endif 2651 #ifdef SET_INEXACT 2652 int inexact, oldinexact; 2653 #endif 2654 2655 #ifndef MULTIPLE_THREADS 2656 if (dtoa_result) { 2657 freedtoa(dtoa_result); 2658 dtoa_result = nullptr; 2659 } 2660 #endif 2661 2662 dval(d) = dd; 2663 if (word0(d) & Sign_bit) { 2664 /* set sign for everything, including 0's and NaNs */ 2665 *sign = 1; 2666 word0(d) &= ~Sign_bit; /* clear sign bit */ 2667 } else { 2668 *sign = 0; 2669 } 2670 2671 #if defined(IEEE_Arith) + defined(VAX) 2672 #ifdef IEEE_Arith 2673 if ((word0(d) & Exp_mask) == Exp_mask) 2674 #else 2675 if (word0(d) == 0x8000) 2676 #endif 2677 { 2678 /* Infinity or NaN */ 2679 *decpt = 9999; 2680 #ifdef IEEE_Arith 2681 if (!word1(d) && !(word0(d) & 0xfffff)) { 2682 return nrv_alloc("Infinity", rve, 8); 2683 } 2684 #endif 2685 return nrv_alloc("NaN", rve, 3); 2686 } 2687 #endif 2688 #ifdef IBM 2689 dval(d) += 0; /* normalize */ 2690 #endif 2691 if (!dval(d)) { 2692 *decpt = 1; 2693 return nrv_alloc("0", rve, 1); 2694 } 2695 2696 #ifdef SET_INEXACT 2697 try_quick = oldinexact = get_inexact(); 2698 inexact = 1; 2699 #endif 2700 #ifdef Honor_FLT_ROUNDS 2701 if ((rounding = Flt_Rounds) >= 2) { 2702 if (*sign) { 2703 rounding = rounding == 2 ? 0 : 2; 2704 } else if (rounding != 2) { 2705 rounding = 0; 2706 } 2707 } 2708 #endif 2709 2710 b = d2b(dval(d), &be, &bbits); 2711 #ifdef Sudden_Underflow 2712 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask >> Exp_shift1)); 2713 #else 2714 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask >> Exp_shift1)))) { 2715 #endif 2716 dval(d2) = dval(d); 2717 word0(d2) &= Frac_mask1; 2718 word0(d2) |= Exp_11; 2719 #ifdef IBM 2720 if (j = 11 - hi0bits(word0(d2) & Frac_mask)) { 2721 dval(d2) /= 1 << j; 2722 } 2723 #endif 2724 2725 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 2726 * log10(x) = log(x) / log(10) 2727 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 2728 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 2729 * 2730 * This suggests computing an approximation k to log10(d) by 2731 * 2732 * k = (i - Bias)*0.301029995663981 2733 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 2734 * 2735 * We want k to be too large rather than too small. 2736 * The error in the first-order Taylor series approximation 2737 * is in our favor, so we just round up the constant enough 2738 * to compensate for any error in the multiplication of 2739 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 2740 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 2741 * adding 1e-13 to the constant term more than suffices. 2742 * Hence we adjust the constant term to 0.1760912590558. 2743 * (We could get a more accurate k by invoking log10, 2744 * but this is probably not worthwhile.) 2745 */ 2746 2747 i -= Bias; 2748 #ifdef IBM 2749 i <<= 2; 2750 i += j; 2751 #endif 2752 #ifndef Sudden_Underflow 2753 denorm = 0; 2754 } else { 2755 /* d is denormalized */ 2756 2757 i = bbits + be + (Bias + (P - 1) - 1); 2758 x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32) 2759 : word1(d) << (32 - i); 2760 dval(d2) = x; 2761 word0(d2) -= 31 * Exp_msk1; /* adjust exponent */ 2762 i -= (Bias + (P - 1) - 1) + 1; 2763 denorm = 1; 2764 } 2765 #endif 2766 ds = (dval(d2) - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981; 2767 k = (int)ds; 2768 if (ds < 0. && ds != k) { 2769 k--; /* want k = floor(ds) */ 2770 } 2771 k_check = 1; 2772 if (k >= 0 && k <= Ten_pmax) { 2773 if (dval(d) < tens[k]) { 2774 k--; 2775 } 2776 k_check = 0; 2777 } 2778 j = bbits - i - 1; 2779 if (j >= 0) { 2780 b2 = 0; 2781 s2 = j; 2782 } else { 2783 b2 = -j; 2784 s2 = 0; 2785 } 2786 if (k >= 0) { 2787 b5 = 0; 2788 s5 = k; 2789 s2 += k; 2790 } else { 2791 b2 -= k; 2792 b5 = -k; 2793 s5 = 0; 2794 } 2795 if (mode < 0 || mode > 9) { 2796 mode = 0; 2797 } 2798 2799 #ifndef SET_INEXACT 2800 #ifdef Check_FLT_ROUNDS 2801 try_quick = Rounding == 1; 2802 #else 2803 try_quick = 1; 2804 #endif 2805 #endif /*SET_INEXACT*/ 2806 2807 if (mode > 5) { 2808 mode -= 4; 2809 try_quick = 0; 2810 } 2811 leftright = 1; 2812 switch (mode) { 2813 case 0: 2814 case 1: 2815 ilim = ilim1 = -1; 2816 i = 18; 2817 ndigits = 0; 2818 break; 2819 case 2: 2820 leftright = 0; 2821 /* no break */ 2822 case 4: 2823 if (ndigits <= 0) { 2824 ndigits = 1; 2825 } 2826 ilim = ilim1 = i = ndigits; 2827 break; 2828 case 3: 2829 leftright = 0; 2830 /* no break */ 2831 case 5: 2832 i = ndigits + k + 1; 2833 ilim = i; 2834 ilim1 = i - 1; 2835 if (i <= 0) { 2836 i = 1; 2837 } 2838 } 2839 s = s0 = rv_alloc(i); 2840 2841 #ifdef Honor_FLT_ROUNDS 2842 if (mode > 1 && rounding != 1) { 2843 leftright = 0; 2844 } 2845 #endif 2846 2847 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 2848 2849 /* Try to get by with floating-point arithmetic. */ 2850 2851 i = 0; 2852 dval(d2) = dval(d); 2853 k0 = k; 2854 ilim0 = ilim; 2855 ieps = 2; /* conservative */ 2856 if (k > 0) { 2857 ds = tens[k & 0xf]; 2858 j = k >> 4; 2859 if (j & Bletch) { 2860 /* prevent overflows */ 2861 j &= Bletch - 1; 2862 dval(d) /= bigtens[n_bigtens - 1]; 2863 ieps++; 2864 } 2865 for (; j; j >>= 1, i++) 2866 if (j & 1) { 2867 ieps++; 2868 ds *= bigtens[i]; 2869 } 2870 dval(d) /= ds; 2871 } else if ((j1 = -k)) { 2872 dval(d) *= tens[j1 & 0xf]; 2873 for (j = j1 >> 4; j; j >>= 1, i++) 2874 if (j & 1) { 2875 ieps++; 2876 dval(d) *= bigtens[i]; 2877 } 2878 } 2879 if (k_check && dval(d) < 1. && ilim > 0) { 2880 if (ilim1 <= 0) { 2881 goto fast_failed; 2882 } 2883 ilim = ilim1; 2884 k--; 2885 dval(d) *= 10.; 2886 ieps++; 2887 } 2888 dval(eps) = ieps * dval(d) + 7.; 2889 word0(eps) -= (P - 1) * Exp_msk1; 2890 if (ilim == 0) { 2891 S = mhi = nullptr; 2892 dval(d) -= 5.; 2893 if (dval(d) > dval(eps)) { 2894 goto one_digit; 2895 } 2896 if (dval(d) < -dval(eps)) { 2897 goto no_digits; 2898 } 2899 goto fast_failed; 2900 } 2901 #ifndef No_leftright 2902 if (leftright) { 2903 /* Use Steele & White method of only 2904 * generating digits needed. 2905 */ 2906 dval(eps) = 0.5 / tens[ilim - 1] - dval(eps); 2907 for (i = 0;;) { 2908 L = (long int)dval(d); 2909 dval(d) -= L; 2910 *s++ = '0' + (int)L; 2911 if (dval(d) < dval(eps)) { 2912 goto ret1; 2913 } 2914 if (1. - dval(d) < dval(eps)) { 2915 goto bump_up; 2916 } 2917 if (++i >= ilim) { 2918 break; 2919 } 2920 dval(eps) *= 10.; 2921 dval(d) *= 10.; 2922 } 2923 } else { 2924 #endif 2925 /* Generate ilim digits, then fix them up. */ 2926 dval(eps) *= tens[ilim - 1]; 2927 for (i = 1;; i++, dval(d) *= 10.) { 2928 L = (Long)(dval(d)); 2929 if (!(dval(d) -= L)) { 2930 ilim = i; 2931 } 2932 *s++ = '0' + (int)L; 2933 if (i == ilim) { 2934 if (dval(d) > 0.5 + dval(eps)) { 2935 goto bump_up; 2936 } else if (dval(d) < 0.5 - dval(eps)) { 2937 while (*--s == '0') 2938 ; 2939 s++; 2940 goto ret1; 2941 } 2942 break; 2943 } 2944 } 2945 #ifndef No_leftright 2946 } 2947 #endif 2948 fast_failed: 2949 s = s0; 2950 dval(d) = dval(d2); 2951 k = k0; 2952 ilim = ilim0; 2953 } 2954 2955 /* Do we have a "small" integer? */ 2956 2957 if (be >= 0 && k <= Int_max) { 2958 /* Yes. */ 2959 ds = tens[k]; 2960 if (ndigits < 0 && ilim <= 0) { 2961 S = mhi = nullptr; 2962 if (ilim < 0 || dval(d) <= 5 * ds) { 2963 goto no_digits; 2964 } 2965 goto one_digit; 2966 } 2967 for (i = 1;; i++, dval(d) *= 10.) { 2968 L = (Long)(dval(d) / ds); 2969 dval(d) -= L * ds; 2970 #ifdef Check_FLT_ROUNDS 2971 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 2972 if (dval(d) < 0) { 2973 L--; 2974 dval(d) += ds; 2975 } 2976 #endif 2977 *s++ = '0' + (int)L; 2978 if (!dval(d)) { 2979 #ifdef SET_INEXACT 2980 inexact = 0; 2981 #endif 2982 break; 2983 } 2984 if (i == ilim) { 2985 #ifdef Honor_FLT_ROUNDS 2986 if (mode > 1) 2987 switch (rounding) { 2988 case 0: goto ret1; 2989 case 2: goto bump_up; 2990 } 2991 #endif 2992 dval(d) += dval(d); 2993 if (dval(d) > ds || (dval(d) == ds && L & 1)) { 2994 bump_up: 2995 while (*--s == '9') 2996 if (s == s0) { 2997 k++; 2998 *s = '0'; 2999 break; 3000 } 3001 ++*s++; 3002 } 3003 break; 3004 } 3005 } 3006 goto ret1; 3007 } 3008 3009 m2 = b2; 3010 m5 = b5; 3011 mhi = mlo = nullptr; 3012 if (leftright) { 3013 i = 3014 #ifndef Sudden_Underflow 3015 denorm ? be + (Bias + (P - 1) - 1 + 1) : 3016 #endif 3017 #ifdef IBM 3018 1 + 4 * P - 3 - bbits + ((bbits + be - 1) & 3); 3019 #else 3020 1 + P - bbits; 3021 #endif 3022 b2 += i; 3023 s2 += i; 3024 mhi = i2b(1); 3025 } 3026 if (m2 > 0 && s2 > 0) { 3027 i = m2 < s2 ? m2 : s2; 3028 b2 -= i; 3029 m2 -= i; 3030 s2 -= i; 3031 } 3032 if (b5 > 0) { 3033 if (leftright) { 3034 if (m5 > 0) { 3035 mhi = pow5mult(mhi, m5); 3036 b1 = mult(mhi, b); 3037 Bfree(b); 3038 b = b1; 3039 } 3040 if ((j = b5 - m5)) { 3041 b = pow5mult(b, j); 3042 } 3043 } else { 3044 b = pow5mult(b, b5); 3045 } 3046 } 3047 S = i2b(1); 3048 if (s5 > 0) { 3049 S = pow5mult(S, s5); 3050 } 3051 3052 /* Check for special case that d is a normalized power of 2. */ 3053 3054 spec_case = 0; 3055 if ((mode < 2 || leftright) 3056 #ifdef Honor_FLT_ROUNDS 3057 && rounding == 1 3058 #endif 3059 ) { 3060 if (!word1(d) && !(word0(d) & Bndry_mask) 3061 #ifndef Sudden_Underflow 3062 && word0(d) & (Exp_mask & ~Exp_msk1) 3063 #endif 3064 ) { 3065 /* The special case */ 3066 b2 += Log2P; 3067 s2 += Log2P; 3068 spec_case = 1; 3069 } 3070 } 3071 3072 /* Arrange for convenient computation of quotients: 3073 * shift left if necessary so divisor has 4 leading 0 bits. 3074 * 3075 * Perhaps we should just compute leading 28 bits of S once 3076 * and for all and pass them and a shift to quorem, so it 3077 * can do shifts and ors to compute the numerator for q. 3078 */ 3079 #ifdef Pack_32 3080 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds - 1]) : 1) + s2) & 0x1f)) { 3081 i = 32 - i; 3082 } 3083 #else 3084 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds - 1]) : 1) + s2) & 0xf) { 3085 i = 16 - i; 3086 } 3087 #endif 3088 if (i > 4) { 3089 i -= 4; 3090 b2 += i; 3091 m2 += i; 3092 s2 += i; 3093 } else if (i < 4) { 3094 i += 28; 3095 b2 += i; 3096 m2 += i; 3097 s2 += i; 3098 } 3099 if (b2 > 0) { 3100 b = lshift(b, b2); 3101 } 3102 if (s2 > 0) { 3103 S = lshift(S, s2); 3104 } 3105 if (k_check) { 3106 if (cmp(b, S) < 0) { 3107 k--; 3108 b = multadd(b, 10, 0); /* we botched the k estimate */ 3109 if (leftright) { 3110 mhi = multadd(mhi, 10, 0); 3111 } 3112 ilim = ilim1; 3113 } 3114 } 3115 if (ilim <= 0 && (mode == 3 || mode == 5)) { 3116 if (ilim < 0 || cmp(b, S = multadd(S, 5, 0)) <= 0) { 3117 /* no digits, fcvt style */ 3118 no_digits: 3119 k = -1 - ndigits; 3120 goto ret; 3121 } 3122 one_digit: 3123 *s++ = '1'; 3124 k++; 3125 goto ret; 3126 } 3127 if (leftright) { 3128 if (m2 > 0) { 3129 mhi = lshift(mhi, m2); 3130 } 3131 3132 /* Compute mlo -- check for special case 3133 * that d is a normalized power of 2. 3134 */ 3135 3136 mlo = mhi; 3137 if (spec_case) { 3138 mhi = Balloc(mhi->k); 3139 Bcopy(mhi, mlo); 3140 mhi = lshift(mhi, Log2P); 3141 } 3142 3143 for (i = 1;; i++) { 3144 dig = quorem(b, S) + '0'; 3145 /* Do we yet have the shortest decimal string 3146 * that will round to d? 3147 */ 3148 j = cmp(b, mlo); 3149 delta = diff(S, mhi); 3150 j1 = delta->sign ? 1 : cmp(b, delta); 3151 Bfree(delta); 3152 #ifndef ROUND_BIASED 3153 if (j1 == 0 && mode != 1 && !(word1(d) & 1) 3154 #ifdef Honor_FLT_ROUNDS 3155 && rounding >= 1 3156 #endif 3157 ) { 3158 if (dig == '9') { 3159 goto round_9_up; 3160 } 3161 if (j > 0) { 3162 dig++; 3163 } 3164 #ifdef SET_INEXACT 3165 else if (!b->x[0] && b->wds <= 1) { 3166 inexact = 0; 3167 } 3168 #endif 3169 *s++ = dig; 3170 goto ret; 3171 } 3172 #endif 3173 if (j < 0 || (j == 0 && mode != 1 3174 #ifndef ROUND_BIASED 3175 && !(word1(d) & 1) 3176 #endif 3177 )) { 3178 if (!b->x[0] && b->wds <= 1) { 3179 #ifdef SET_INEXACT 3180 inexact = 0; 3181 #endif 3182 goto accept_dig; 3183 } 3184 #ifdef Honor_FLT_ROUNDS 3185 if (mode > 1) 3186 switch (rounding) { 3187 case 0: goto accept_dig; 3188 case 2: goto keep_dig; 3189 } 3190 #endif /*Honor_FLT_ROUNDS*/ 3191 if (j1 > 0) { 3192 b = lshift(b, 1); 3193 j1 = cmp(b, S); 3194 if ((j1 > 0 || (j1 == 0 && dig & 1)) 3195 && dig++ == '9') { 3196 goto round_9_up; 3197 } 3198 } 3199 accept_dig: 3200 *s++ = dig; 3201 goto ret; 3202 } 3203 if (j1 > 0) { 3204 #ifdef Honor_FLT_ROUNDS 3205 if (!rounding) { 3206 goto accept_dig; 3207 } 3208 #endif 3209 if (dig == '9') { /* possible if i == 1 */ 3210 round_9_up: 3211 *s++ = '9'; 3212 goto roundoff; 3213 } 3214 *s++ = dig + 1; 3215 goto ret; 3216 } 3217 #ifdef Honor_FLT_ROUNDS 3218 keep_dig: 3219 #endif 3220 *s++ = dig; 3221 if (i == ilim) { 3222 break; 3223 } 3224 b = multadd(b, 10, 0); 3225 if (mlo == mhi) { 3226 mlo = mhi = multadd(mhi, 10, 0); 3227 } else { 3228 mlo = multadd(mlo, 10, 0); 3229 mhi = multadd(mhi, 10, 0); 3230 } 3231 } 3232 } else 3233 for (i = 1;; i++) { 3234 *s++ = dig = quorem(b, S) + '0'; 3235 if (!b->x[0] && b->wds <= 1) { 3236 #ifdef SET_INEXACT 3237 inexact = 0; 3238 #endif 3239 goto ret; 3240 } 3241 if (i >= ilim) { 3242 break; 3243 } 3244 b = multadd(b, 10, 0); 3245 } 3246 3247 /* Round off last digit */ 3248 3249 #ifdef Honor_FLT_ROUNDS 3250 switch (rounding) { 3251 case 0: goto trimzeros; 3252 case 2: goto roundoff; 3253 } 3254 #endif 3255 b = lshift(b, 1); 3256 j = cmp(b, S); 3257 if (j > 0 || (j == 0 && dig & 1)) { 3258 roundoff: 3259 while (*--s == '9') 3260 if (s == s0) { 3261 k++; 3262 *s++ = '1'; 3263 goto ret; 3264 } 3265 ++*s++; 3266 } else { 3267 #ifdef Honor_FLT_ROUNDS 3268 trimzeros: 3269 #endif 3270 while (*--s == '0') 3271 ; 3272 s++; 3273 } 3274 ret: 3275 Bfree(S); 3276 if (mhi) { 3277 if (mlo && mlo != mhi) { 3278 Bfree(mlo); 3279 } 3280 Bfree(mhi); 3281 } 3282 ret1: 3283 #ifdef SET_INEXACT 3284 if (inexact) { 3285 if (!oldinexact) { 3286 word0(d) = Exp_1 + (70 << Exp_shift); 3287 word1(d) = 0; 3288 dval(d) += 1.; 3289 } 3290 } else if (!oldinexact) { 3291 clear_inexact(); 3292 } 3293 #endif 3294 Bfree(b); 3295 *s = 0; 3296 *decpt = k + 1; 3297 if (rve) { 3298 *rve = s; 3299 } 3300 return s0; 3301 } 3302 #ifdef __cplusplus 3303 } 3304 #endif