kjs Library API Documentation

dtoa.cpp

00001 /**************************************************************** 00002 * 00003 * The author of this software is David M. Gay. 00004 * 00005 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. 00006 * 00007 * Permission to use, copy, modify, and distribute this software for any 00008 * purpose without fee is hereby granted, provided that this entire notice 00009 * is included in all copies of any software which is or includes a copy 00010 * or modification of this software and in all copies of the supporting 00011 * documentation for such software. 00012 * 00013 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 00014 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY 00015 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 00016 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 00017 * 00018 ***************************************************************/ 00019 00020 /* Please send bug reports to 00021 David M. Gay 00022 Bell Laboratories, Room 2C-463 00023 600 Mountain Avenue 00024 Murray Hill, NJ 07974-0636 00025 U.S.A. 00026 dmg@bell-labs.com 00027 */ 00028 00029 /* On a machine with IEEE extended-precision registers, it is 00030 * necessary to specify double-precision (53-bit) rounding precision 00031 * before invoking strtod or dtoa. If the machine uses (the equivalent 00032 * of) Intel 80x87 arithmetic, the call 00033 * _control87(PC_53, MCW_PC); 00034 * does this with many compilers. Whether this or another call is 00035 * appropriate depends on the compiler; for this to work, it may be 00036 * necessary to #include "float.h" or another system-dependent header 00037 * file. 00038 */ 00039 00040 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines. 00041 * 00042 * This strtod returns a nearest machine number to the input decimal 00043 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 00044 * broken by the IEEE round-even rule. Otherwise ties are broken by 00045 * biased rounding (add half and chop). 00046 * 00047 * Inspired loosely by William D. Clinger's paper "How to Read Floating 00048 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 00049 * 00050 * Modifications: 00051 * 00052 * 1. We only require IEEE, IBM, or VAX double-precision 00053 * arithmetic (not IEEE double-extended). 00054 * 2. We get by with floating-point arithmetic in a case that 00055 * Clinger missed -- when we're computing d * 10^n 00056 * for a small integer d and the integer n is not too 00057 * much larger than 22 (the maximum integer k for which 00058 * we can represent 10^k exactly), we may be able to 00059 * compute (d*10^k) * 10^(e-k) with just one roundoff. 00060 * 3. Rather than a bit-at-a-time adjustment of the binary 00061 * result in the hard case, we use floating-point 00062 * arithmetic to determine the adjustment to within 00063 * one bit; only in really hard cases do we need to 00064 * compute a second residual. 00065 * 4. Because of 3., we don't need a large table of powers of 10 00066 * for ten-to-e (just some small tables, e.g. of 10^k 00067 * for 0 <= k <= 22). 00068 */ 00069 00070 /* 00071 * #define IEEE_8087 for IEEE-arithmetic machines where the least 00072 * significant byte has the lowest address. 00073 * #define IEEE_MC68k for IEEE-arithmetic machines where the most 00074 * significant byte has the lowest address. 00075 * #define Long int on machines with 32-bit ints and 64-bit longs. 00076 * #define IBM for IBM mainframe-style floating-point arithmetic. 00077 * #define VAX for VAX-style floating-point arithmetic (D_floating). 00078 * #define No_leftright to omit left-right logic in fast floating-point 00079 * computation of dtoa. 00080 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 00081 * and strtod and dtoa should round accordingly. 00082 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 00083 * and Honor_FLT_ROUNDS is not #defined. 00084 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 00085 * that use extended-precision instructions to compute rounded 00086 * products and quotients) with IBM. 00087 * #define ROUND_BIASED for IEEE-format with biased rounding. 00088 * #define Inaccurate_Divide for IEEE-format with correctly rounded 00089 * products but inaccurate quotients, e.g., for Intel i860. 00090 * #define NO_LONG_LONG on machines that do not have a "long long" 00091 * integer type (of >= 64 bits). On such machines, you can 00092 * #define Just_16 to store 16 bits per 32-bit Long when doing 00093 * high-precision integer arithmetic. Whether this speeds things 00094 * up or slows things down depends on the machine and the number 00095 * being converted. If long long is available and the name is 00096 * something other than "long long", #define Llong to be the name, 00097 * and if "unsigned Llong" does not work as an unsigned version of 00098 * Llong, #define #ULLong to be the corresponding unsigned type. 00099 * #define KR_headers for old-style C function headers. 00100 * #define Bad_float_h if your system lacks a float.h or if it does not 00101 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, 00102 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. 00103 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) 00104 * if memory is available and otherwise does something you deem 00105 * appropriate. If MALLOC is undefined, malloc will be invoked 00106 * directly -- and assumed always to succeed. 00107 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making 00108 * memory allocations from a private pool of memory when possible. 00109 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes, 00110 * unless #defined to be a different length. This default length 00111 * suffices to get rid of MALLOC calls except for unusual cases, 00112 * such as decimal-to-binary conversion of a very long string of 00113 * digits. The longest string dtoa can return is about 751 bytes 00114 * long. For conversions by strtod of strings of 800 digits and 00115 * all dtoa conversions in single-threaded executions with 8-byte 00116 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte 00117 * pointers, PRIVATE_MEM >= 7112 appears adequate. 00118 * #define INFNAN_CHECK on IEEE systems to cause strtod to check for 00119 * Infinity and NaN (case insensitively). On some systems (e.g., 00120 * some HP systems), it may be necessary to #define NAN_WORD0 00121 * appropriately -- to the most significant word of a quiet NaN. 00122 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.) 00123 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined, 00124 * strtod also accepts (case insensitively) strings of the form 00125 * NaN(x), where x is a string of hexadecimal digits and spaces; 00126 * if there is only one string of hexadecimal digits, it is taken 00127 * for the 52 fraction bits of the resulting NaN; if there are two 00128 * or more strings of hex digits, the first is for the high 20 bits, 00129 * the second and subsequent for the low 32 bits, with intervening 00130 * white space ignored; but if this results in none of the 52 00131 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0 00132 * and NAN_WORD1 are used instead. 00133 * #define MULTIPLE_THREADS if the system offers preemptively scheduled 00134 * multiple threads. In this case, you must provide (or suitably 00135 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed 00136 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed 00137 * in pow5mult, ensures lazy evaluation of only one copy of high 00138 * powers of 5; omitting this lock would introduce a small 00139 * probability of wasting memory, but would otherwise be harmless.) 00140 * You must also invoke freedtoa(s) to free the value s returned by 00141 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. 00142 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that 00143 * avoids underflows on inputs whose result does not underflow. 00144 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format 00145 * floating-point numbers and flushes underflows to zero rather 00146 * than implementing gradual underflow, then you must also #define 00147 * Sudden_Underflow. 00148 * #define YES_ALIAS to permit aliasing certain double values with 00149 * arrays of ULongs. This leads to slightly better code with 00150 * some compilers and was always used prior to 19990916, but it 00151 * is not strictly legal and can cause trouble with aggressively 00152 * optimizing compilers (e.g., gcc 2.95.1 under -O2). 00153 * #define USE_LOCALE to use the current locale's decimal_point value. 00154 * #define SET_INEXACT if IEEE arithmetic is being used and extra 00155 * computation should be done to set the inexact flag when the 00156 * result is inexact and avoid setting inexact when the result 00157 * is exact. In this case, dtoa.c must be compiled in 00158 * an environment, perhaps provided by #include "dtoa.c" in a 00159 * suitable wrapper, that defines two functions, 00160 * int get_inexact(void); 00161 * void clear_inexact(void); 00162 * such that get_inexact() returns a nonzero value if the 00163 * inexact bit is already set, and clear_inexact() sets the 00164 * inexact bit to 0. When SET_INEXACT is #defined, strtod 00165 * also does extra computations to set the underflow and overflow 00166 * flags when appropriate (i.e., when the result is tiny and 00167 * inexact or when it is a numeric value rounded to +-infinity). 00168 * #define NO_ERRNO if strtod should not assign errno = ERANGE when 00169 * the result overflows to +-Infinity or underflows to 0. 00170 */ 00171 00172 #include <config.h> 00173 00174 #include "stdlib.h" 00175 00176 #ifdef WORDS_BIGENDIAN 00177 #define IEEE_MC68k 00178 #else 00179 #define IEEE_8087 00180 #endif 00181 #define INFNAN_CHECK 00182 #include "dtoa.h" 00183 00184 00185 00186 #ifndef Long 00187 #define Long int 00188 #endif 00189 #ifndef ULong 00190 typedef unsigned Long ULong; 00191 #endif 00192 00193 #ifdef DEBUG 00194 #include "stdio.h" 00195 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 00196 #endif 00197 00198 #include "string.h" 00199 00200 #ifdef USE_LOCALE 00201 #include "locale.h" 00202 #endif 00203 00204 #ifdef MALLOC 00205 #ifdef KR_headers 00206 extern char *MALLOC(); 00207 #else 00208 extern void *MALLOC(size_t); 00209 #endif 00210 #else 00211 #define MALLOC malloc 00212 #endif 00213 00214 #ifndef Omit_Private_Memory 00215 #ifndef PRIVATE_MEM 00216 #define PRIVATE_MEM 2304 00217 #endif 00218 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) 00219 static double private_mem[PRIVATE_mem], *pmem_next = private_mem; 00220 #endif 00221 00222 #undef IEEE_Arith 00223 #undef Avoid_Underflow 00224 #ifdef IEEE_MC68k 00225 #define IEEE_Arith 00226 #endif 00227 #ifdef IEEE_8087 00228 #define IEEE_Arith 00229 #endif 00230 00231 #include "errno.h" 00232 00233 #ifdef Bad_float_h 00234 00235 #ifdef IEEE_Arith 00236 #define DBL_DIG 15 00237 #define DBL_MAX_10_EXP 308 00238 #define DBL_MAX_EXP 1024 00239 #define FLT_RADIX 2 00240 #endif /*IEEE_Arith*/ 00241 00242 #ifdef IBM 00243 #define DBL_DIG 16 00244 #define DBL_MAX_10_EXP 75 00245 #define DBL_MAX_EXP 63 00246 #define FLT_RADIX 16 00247 #define DBL_MAX 7.2370055773322621e+75 00248 #endif 00249 00250 #ifdef VAX 00251 #define DBL_DIG 16 00252 #define DBL_MAX_10_EXP 38 00253 #define DBL_MAX_EXP 127 00254 #define FLT_RADIX 2 00255 #define DBL_MAX 1.7014118346046923e+38 00256 #endif 00257 00258 #else /* ifndef Bad_float_h */ 00259 #include "float.h" 00260 #endif /* Bad_float_h */ 00261 00262 #ifndef __MATH_H__ 00263 #include "math.h" 00264 #endif 00265 00266 #ifdef __cplusplus 00267 extern "C" { 00268 #endif 00269 00270 #ifndef CONST 00271 #ifdef KR_headers 00272 #define CONST /* blank */ 00273 #else 00274 #define CONST const 00275 #endif 00276 #endif 00277 00278 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 00279 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. 00280 #endif 00281 00282 typedef union { double d; ULong L[2]; } U; 00283 00284 #ifdef YES_ALIAS 00285 #define dval(x) x 00286 #ifdef IEEE_8087 00287 #define word0(x) ((ULong *)&x)[1] 00288 #define word1(x) ((ULong *)&x)[0] 00289 #else 00290 #define word0(x) ((ULong *)&x)[0] 00291 #define word1(x) ((ULong *)&x)[1] 00292 #endif 00293 #else 00294 #ifdef IEEE_8087 00295 #define word0(x) ((U*)&x)->L[1] 00296 #define word1(x) ((U*)&x)->L[0] 00297 #else 00298 #define word0(x) ((U*)&x)->L[0] 00299 #define word1(x) ((U*)&x)->L[1] 00300 #endif 00301 #define dval(x) ((U*)&x)->d 00302 #endif 00303 00304 /* The following definition of Storeinc is appropriate for MIPS processors. 00305 * An alternative that might be better on some machines is 00306 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 00307 */ 00308 #if defined(IEEE_8087) + defined(VAX) 00309 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ 00310 ((unsigned short *)a)[0] = (unsigned short)c, a++) 00311 #else 00312 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ 00313 ((unsigned short *)a)[1] = (unsigned short)c, a++) 00314 #endif 00315 00316 /* #define P DBL_MANT_DIG */ 00317 /* Ten_pmax = floor(P*log(2)/log(5)) */ 00318 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 00319 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 00320 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 00321 00322 #ifdef IEEE_Arith 00323 #define Exp_shift 20 00324 #define Exp_shift1 20 00325 #define Exp_msk1 0x100000 00326 #define Exp_msk11 0x100000 00327 #define Exp_mask 0x7ff00000 00328 #define P 53 00329 #define Bias 1023 00330 #define Emin (-1022) 00331 #define Exp_1 0x3ff00000 00332 #define Exp_11 0x3ff00000 00333 #define Ebits 11 00334 #define Frac_mask 0xfffff 00335 #define Frac_mask1 0xfffff 00336 #define Ten_pmax 22 00337 #define Bletch 0x10 00338 #define Bndry_mask 0xfffff 00339 #define Bndry_mask1 0xfffff 00340 #define LSB 1 00341 #define Sign_bit 0x80000000 00342 #define Log2P 1 00343 #define Tiny0 0 00344 #define Tiny1 1 00345 #define Quick_max 14 00346 #define Int_max 14 00347 #ifndef NO_IEEE_Scale 00348 #define Avoid_Underflow 00349 #ifdef Flush_Denorm /* debugging option */ 00350 #undef Sudden_Underflow 00351 #endif 00352 #endif 00353 00354 #ifndef Flt_Rounds 00355 #ifdef FLT_ROUNDS 00356 #define Flt_Rounds FLT_ROUNDS 00357 #else 00358 #define Flt_Rounds 1 00359 #endif 00360 #endif /*Flt_Rounds*/ 00361 00362 #ifdef Honor_FLT_ROUNDS 00363 #define Rounding rounding 00364 #undef Check_FLT_ROUNDS 00365 #define Check_FLT_ROUNDS 00366 #else 00367 #define Rounding Flt_Rounds 00368 #endif 00369 00370 #else /* ifndef IEEE_Arith */ 00371 #undef Check_FLT_ROUNDS 00372 #undef Honor_FLT_ROUNDS 00373 #undef SET_INEXACT 00374 #undef Sudden_Underflow 00375 #define Sudden_Underflow 00376 #ifdef IBM 00377 #undef Flt_Rounds 00378 #define Flt_Rounds 0 00379 #define Exp_shift 24 00380 #define Exp_shift1 24 00381 #define Exp_msk1 0x1000000 00382 #define Exp_msk11 0x1000000 00383 #define Exp_mask 0x7f000000 00384 #define P 14 00385 #define Bias 65 00386 #define Exp_1 0x41000000 00387 #define Exp_11 0x41000000 00388 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 00389 #define Frac_mask 0xffffff 00390 #define Frac_mask1 0xffffff 00391 #define Bletch 4 00392 #define Ten_pmax 22 00393 #define Bndry_mask 0xefffff 00394 #define Bndry_mask1 0xffffff 00395 #define LSB 1 00396 #define Sign_bit 0x80000000 00397 #define Log2P 4 00398 #define Tiny0 0x100000 00399 #define Tiny1 0 00400 #define Quick_max 14 00401 #define Int_max 15 00402 #else /* VAX */ 00403 #undef Flt_Rounds 00404 #define Flt_Rounds 1 00405 #define Exp_shift 23 00406 #define Exp_shift1 7 00407 #define Exp_msk1 0x80 00408 #define Exp_msk11 0x800000 00409 #define Exp_mask 0x7f80 00410 #define P 56 00411 #define Bias 129 00412 #define Exp_1 0x40800000 00413 #define Exp_11 0x4080 00414 #define Ebits 8 00415 #define Frac_mask 0x7fffff 00416 #define Frac_mask1 0xffff007f 00417 #define Ten_pmax 24 00418 #define Bletch 2 00419 #define Bndry_mask 0xffff007f 00420 #define Bndry_mask1 0xffff007f 00421 #define LSB 0x10000 00422 #define Sign_bit 0x8000 00423 #define Log2P 1 00424 #define Tiny0 0x80 00425 #define Tiny1 0 00426 #define Quick_max 15 00427 #define Int_max 15 00428 #endif /* IBM, VAX */ 00429 #endif /* IEEE_Arith */ 00430 00431 #ifndef IEEE_Arith 00432 #define ROUND_BIASED 00433 #endif 00434 00435 #ifdef RND_PRODQUOT 00436 #define rounded_product(a,b) a = rnd_prod(a, b) 00437 #define rounded_quotient(a,b) a = rnd_quot(a, b) 00438 #ifdef KR_headers 00439 extern double rnd_prod(), rnd_quot(); 00440 #else 00441 extern double rnd_prod(double, double), rnd_quot(double, double); 00442 #endif 00443 #else 00444 #define rounded_product(a,b) a *= b 00445 #define rounded_quotient(a,b) a /= b 00446 #endif 00447 00448 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 00449 #define Big1 0xffffffff 00450 00451 #ifndef Pack_32 00452 #define Pack_32 00453 #endif 00454 00455 #ifdef KR_headers 00456 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff) 00457 #else 00458 #define FFFFFFFF 0xffffffffUL 00459 #endif 00460 00461 #ifdef NO_LONG_LONG 00462 #undef ULLong 00463 #ifdef Just_16 00464 #undef Pack_32 00465 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long. 00466 * This makes some inner loops simpler and sometimes saves work 00467 * during multiplications, but it often seems to make things slightly 00468 * slower. Hence the default is now to store 32 bits per Long. 00469 */ 00470 #endif 00471 #else /* long long available */ 00472 #ifndef Llong 00473 #define Llong long long 00474 #endif 00475 #ifndef ULLong 00476 #define ULLong unsigned Llong 00477 #endif 00478 #endif /* NO_LONG_LONG */ 00479 00480 #ifndef MULTIPLE_THREADS 00481 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/ 00482 #define FREE_DTOA_LOCK(n) /*nothing*/ 00483 #endif 00484 00485 #define Kmax 15 00486 00487 struct 00488 Bigint { 00489 struct Bigint *next; 00490 int k, maxwds, sign, wds; 00491 ULong x[1]; 00492 }; 00493 00494 typedef struct Bigint Bigint; 00495 00496 static Bigint *freelist[Kmax+1]; 00497 00498 static Bigint * 00499 Balloc 00500 #ifdef KR_headers 00501 (k) int k; 00502 #else 00503 (int k) 00504 #endif 00505 { 00506 int x; 00507 Bigint *rv; 00508 #ifndef Omit_Private_Memory 00509 unsigned int len; 00510 #endif 00511 00512 ACQUIRE_DTOA_LOCK(0); 00513 if ((rv = freelist[k])) { 00514 freelist[k] = rv->next; 00515 } 00516 else { 00517 x = 1 << k; 00518 #ifdef Omit_Private_Memory 00519 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); 00520 #else 00521 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) 00522 /sizeof(double); 00523 if (pmem_next - private_mem + len <= PRIVATE_mem) { 00524 rv = (Bigint*)pmem_next; 00525 pmem_next += len; 00526 } 00527 else 00528 rv = (Bigint*)MALLOC(len*sizeof(double)); 00529 #endif 00530 rv->k = k; 00531 rv->maxwds = x; 00532 } 00533 FREE_DTOA_LOCK(0); 00534 rv->sign = rv->wds = 0; 00535 return rv; 00536 } 00537 00538 static void 00539 Bfree 00540 #ifdef KR_headers 00541 (v) Bigint *v; 00542 #else 00543 (Bigint *v) 00544 #endif 00545 { 00546 if (v) { 00547 ACQUIRE_DTOA_LOCK(0); 00548 v->next = freelist[v->k]; 00549 freelist[v->k] = v; 00550 FREE_DTOA_LOCK(0); 00551 } 00552 } 00553 00554 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 00555 y->wds*sizeof(Long) + 2*sizeof(int)) 00556 00557 static Bigint * 00558 multadd 00559 #ifdef KR_headers 00560 (b, m, a) Bigint *b; int m, a; 00561 #else 00562 (Bigint *b, int m, int a) /* multiply by m and add a */ 00563 #endif 00564 { 00565 int i, wds; 00566 #ifdef ULLong 00567 ULong *x; 00568 ULLong carry, y; 00569 #else 00570 ULong carry, *x, y; 00571 #ifdef Pack_32 00572 ULong xi, z; 00573 #endif 00574 #endif 00575 Bigint *b1; 00576 00577 wds = b->wds; 00578 x = b->x; 00579 i = 0; 00580 carry = a; 00581 do { 00582 #ifdef ULLong 00583 y = *x * (ULLong)m + carry; 00584 carry = y >> 32; 00585 *x++ = y & FFFFFFFF; 00586 #else 00587 #ifdef Pack_32 00588 xi = *x; 00589 y = (xi & 0xffff) * m + carry; 00590 z = (xi >> 16) * m + (y >> 16); 00591 carry = z >> 16; 00592 *x++ = (z << 16) + (y & 0xffff); 00593 #else 00594 y = *x * m + carry; 00595 carry = y >> 16; 00596 *x++ = y & 0xffff; 00597 #endif 00598 #endif 00599 } 00600 while(++i < wds); 00601 if (carry) { 00602 if (wds >= b->maxwds) { 00603 b1 = Balloc(b->k+1); 00604 Bcopy(b1, b); 00605 Bfree(b); 00606 b = b1; 00607 } 00608 b->x[wds++] = carry; 00609 b->wds = wds; 00610 } 00611 return b; 00612 } 00613 00614 static Bigint * 00615 s2b 00616 #ifdef KR_headers 00617 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9; 00618 #else 00619 (CONST char *s, int nd0, int nd, ULong y9) 00620 #endif 00621 { 00622 Bigint *b; 00623 int i, k; 00624 Long x, y; 00625 00626 x = (nd + 8) / 9; 00627 for(k = 0, y = 1; x > y; y <<= 1, k++) ; 00628 #ifdef Pack_32 00629 b = Balloc(k); 00630 b->x[0] = y9; 00631 b->wds = 1; 00632 #else 00633 b = Balloc(k+1); 00634 b->x[0] = y9 & 0xffff; 00635 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; 00636 #endif 00637 00638 i = 9; 00639 if (9 < nd0) { 00640 s += 9; 00641 do b = multadd(b, 10, *s++ - '0'); 00642 while(++i < nd0); 00643 s++; 00644 } 00645 else 00646 s += 10; 00647 for(; i < nd; i++) 00648 b = multadd(b, 10, *s++ - '0'); 00649 return b; 00650 } 00651 00652 static int 00653 hi0bits 00654 #ifdef KR_headers 00655 (x) register ULong x; 00656 #else 00657 (register ULong x) 00658 #endif 00659 { 00660 register int k = 0; 00661 00662 if (!(x & 0xffff0000)) { 00663 k = 16; 00664 x <<= 16; 00665 } 00666 if (!(x & 0xff000000)) { 00667 k += 8; 00668 x <<= 8; 00669 } 00670 if (!(x & 0xf0000000)) { 00671 k += 4; 00672 x <<= 4; 00673 } 00674 if (!(x & 0xc0000000)) { 00675 k += 2; 00676 x <<= 2; 00677 } 00678 if (!(x & 0x80000000)) { 00679 k++; 00680 if (!(x & 0x40000000)) 00681 return 32; 00682 } 00683 return k; 00684 } 00685 00686 static int 00687 lo0bits 00688 #ifdef KR_headers 00689 (y) ULong *y; 00690 #else 00691 (ULong *y) 00692 #endif 00693 { 00694 register int k; 00695 register ULong x = *y; 00696 00697 if (x & 7) { 00698 if (x & 1) 00699 return 0; 00700 if (x & 2) { 00701 *y = x >> 1; 00702 return 1; 00703 } 00704 *y = x >> 2; 00705 return 2; 00706 } 00707 k = 0; 00708 if (!(x & 0xffff)) { 00709 k = 16; 00710 x >>= 16; 00711 } 00712 if (!(x & 0xff)) { 00713 k += 8; 00714 x >>= 8; 00715 } 00716 if (!(x & 0xf)) { 00717 k += 4; 00718 x >>= 4; 00719 } 00720 if (!(x & 0x3)) { 00721 k += 2; 00722 x >>= 2; 00723 } 00724 if (!(x & 1)) { 00725 k++; 00726 x >>= 1; 00727 if (!x & 1) 00728 return 32; 00729 } 00730 *y = x; 00731 return k; 00732 } 00733 00734 static Bigint * 00735 i2b 00736 #ifdef KR_headers 00737 (i) int i; 00738 #else 00739 (int i) 00740 #endif 00741 { 00742 Bigint *b; 00743 00744 b = Balloc(1); 00745 b->x[0] = i; 00746 b->wds = 1; 00747 return b; 00748 } 00749 00750 static Bigint * 00751 mult 00752 #ifdef KR_headers 00753 (a, b) Bigint *a, *b; 00754 #else 00755 (Bigint *a, Bigint *b) 00756 #endif 00757 { 00758 Bigint *c; 00759 int k, wa, wb, wc; 00760 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 00761 ULong y; 00762 #ifdef ULLong 00763 ULLong carry, z; 00764 #else 00765 ULong carry, z; 00766 #ifdef Pack_32 00767 ULong z2; 00768 #endif 00769 #endif 00770 00771 if (a->wds < b->wds) { 00772 c = a; 00773 a = b; 00774 b = c; 00775 } 00776 k = a->k; 00777 wa = a->wds; 00778 wb = b->wds; 00779 wc = wa + wb; 00780 if (wc > a->maxwds) 00781 k++; 00782 c = Balloc(k); 00783 for(x = c->x, xa = x + wc; x < xa; x++) 00784 *x = 0; 00785 xa = a->x; 00786 xae = xa + wa; 00787 xb = b->x; 00788 xbe = xb + wb; 00789 xc0 = c->x; 00790 #ifdef ULLong 00791 for(; xb < xbe; xc0++) { 00792 if ((y = *xb++)) { 00793 x = xa; 00794 xc = xc0; 00795 carry = 0; 00796 do { 00797 z = *x++ * (ULLong)y + *xc + carry; 00798 carry = z >> 32; 00799 *xc++ = z & FFFFFFFF; 00800 } 00801 while(x < xae); 00802 *xc = carry; 00803 } 00804 } 00805 #else 00806 #ifdef Pack_32 00807 for(; xb < xbe; xb++, xc0++) { 00808 if (y = *xb & 0xffff) { 00809 x = xa; 00810 xc = xc0; 00811 carry = 0; 00812 do { 00813 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 00814 carry = z >> 16; 00815 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 00816 carry = z2 >> 16; 00817 Storeinc(xc, z2, z); 00818 } 00819 while(x < xae); 00820 *xc = carry; 00821 } 00822 if (y = *xb >> 16) { 00823 x = xa; 00824 xc = xc0; 00825 carry = 0; 00826 z2 = *xc; 00827 do { 00828 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 00829 carry = z >> 16; 00830 Storeinc(xc, z, z2); 00831 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 00832 carry = z2 >> 16; 00833 } 00834 while(x < xae); 00835 *xc = z2; 00836 } 00837 } 00838 #else 00839 for(; xb < xbe; xc0++) { 00840 if (y = *xb++) { 00841 x = xa; 00842 xc = xc0; 00843 carry = 0; 00844 do { 00845 z = *x++ * y + *xc + carry; 00846 carry = z >> 16; 00847 *xc++ = z & 0xffff; 00848 } 00849 while(x < xae); 00850 *xc = carry; 00851 } 00852 } 00853 #endif 00854 #endif 00855 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 00856 c->wds = wc; 00857 return c; 00858 } 00859 00860 static Bigint *p5s; 00861 00862 static Bigint * 00863 pow5mult 00864 #ifdef KR_headers 00865 (b, k) Bigint *b; int k; 00866 #else 00867 (Bigint *b, int k) 00868 #endif 00869 { 00870 Bigint *b1, *p5, *p51; 00871 int i; 00872 static int p05[3] = { 5, 25, 125 }; 00873 00874 if ((i = k & 3)) 00875 b = multadd(b, p05[i-1], 0); 00876 00877 if (!(k >>= 2)) 00878 return b; 00879 if (!(p5 = p5s)) { 00880 /* first time */ 00881 #ifdef MULTIPLE_THREADS 00882 ACQUIRE_DTOA_LOCK(1); 00883 if (!(p5 = p5s)) { 00884 p5 = p5s = i2b(625); 00885 p5->next = 0; 00886 } 00887 FREE_DTOA_LOCK(1); 00888 #else 00889 p5 = p5s = i2b(625); 00890 p5->next = 0; 00891 #endif 00892 } 00893 for(;;) { 00894 if (k & 1) { 00895 b1 = mult(b, p5); 00896 Bfree(b); 00897 b = b1; 00898 } 00899 if (!(k >>= 1)) 00900 break; 00901 if (!(p51 = p5->next)) { 00902 #ifdef MULTIPLE_THREADS 00903 ACQUIRE_DTOA_LOCK(1); 00904 if (!(p51 = p5->next)) { 00905 p51 = p5->next = mult(p5,p5); 00906 p51->next = 0; 00907 } 00908 FREE_DTOA_LOCK(1); 00909 #else 00910 p51 = p5->next = mult(p5,p5); 00911 p51->next = 0; 00912 #endif 00913 } 00914 p5 = p51; 00915 } 00916 return b; 00917 } 00918 00919 static Bigint * 00920 lshift 00921 #ifdef KR_headers 00922 (b, k) Bigint *b; int k; 00923 #else 00924 (Bigint *b, int k) 00925 #endif 00926 { 00927 int i, k1, n, n1; 00928 Bigint *b1; 00929 ULong *x, *x1, *xe, z; 00930 00931 #ifdef Pack_32 00932 n = k >> 5; 00933 #else 00934 n = k >> 4; 00935 #endif 00936 k1 = b->k; 00937 n1 = n + b->wds + 1; 00938 for(i = b->maxwds; n1 > i; i <<= 1) 00939 k1++; 00940 b1 = Balloc(k1); 00941 x1 = b1->x; 00942 for(i = 0; i < n; i++) 00943 *x1++ = 0; 00944 x = b->x; 00945 xe = x + b->wds; 00946 #ifdef Pack_32 00947 if (k &= 0x1f) { 00948 k1 = 32 - k; 00949 z = 0; 00950 do { 00951 *x1++ = *x << k | z; 00952 z = *x++ >> k1; 00953 } 00954 while(x < xe); 00955 if ((*x1 = z)) 00956 ++n1; 00957 } 00958 #else 00959 if (k &= 0xf) { 00960 k1 = 16 - k; 00961 z = 0; 00962 do { 00963 *x1++ = *x << k & 0xffff | z; 00964 z = *x++ >> k1; 00965 } 00966 while(x < xe); 00967 if (*x1 = z) 00968 ++n1; 00969 } 00970 #endif 00971 else do 00972 *x1++ = *x++; 00973 while(x < xe); 00974 b1->wds = n1 - 1; 00975 Bfree(b); 00976 return b1; 00977 } 00978 00979 static int 00980 cmp 00981 #ifdef KR_headers 00982 (a, b) Bigint *a, *b; 00983 #else 00984 (Bigint *a, Bigint *b) 00985 #endif 00986 { 00987 ULong *xa, *xa0, *xb, *xb0; 00988 int i, j; 00989 00990 i = a->wds; 00991 j = b->wds; 00992 #ifdef DEBUG 00993 if (i > 1 && !a->x[i-1]) 00994 Bug("cmp called with a->x[a->wds-1] == 0"); 00995 if (j > 1 && !b->x[j-1]) 00996 Bug("cmp called with b->x[b->wds-1] == 0"); 00997 #endif 00998 if (i -= j) 00999 return i; 01000 xa0 = a->x; 01001 xa = xa0 + j; 01002 xb0 = b->x; 01003 xb = xb0 + j; 01004 for(;;) { 01005 if (*--xa != *--xb) 01006 return *xa < *xb ? -1 : 1; 01007 if (xa <= xa0) 01008 break; 01009 } 01010 return 0; 01011 } 01012 01013 static Bigint * 01014 diff 01015 #ifdef KR_headers 01016 (a, b) Bigint *a, *b; 01017 #else 01018 (Bigint *a, Bigint *b) 01019 #endif 01020 { 01021 Bigint *c; 01022 int i, wa, wb; 01023 ULong *xa, *xae, *xb, *xbe, *xc; 01024 #ifdef ULLong 01025 ULLong borrow, y; 01026 #else 01027 ULong borrow, y; 01028 #ifdef Pack_32 01029 ULong z; 01030 #endif 01031 #endif 01032 01033 i = cmp(a,b); 01034 if (!i) { 01035 c = Balloc(0); 01036 c->wds = 1; 01037 c->x[0] = 0; 01038 return c; 01039 } 01040 if (i < 0) { 01041 c = a; 01042 a = b; 01043 b = c; 01044 i = 1; 01045 } 01046 else 01047 i = 0; 01048 c = Balloc(a->k); 01049 c->sign = i; 01050 wa = a->wds; 01051 xa = a->x; 01052 xae = xa + wa; 01053 wb = b->wds; 01054 xb = b->x; 01055 xbe = xb + wb; 01056 xc = c->x; 01057 borrow = 0; 01058 #ifdef ULLong 01059 do { 01060 y = (ULLong)*xa++ - *xb++ - borrow; 01061 borrow = y >> 32 & (ULong)1; 01062 *xc++ = y & FFFFFFFF; 01063 } 01064 while(xb < xbe); 01065 while(xa < xae) { 01066 y = *xa++ - borrow; 01067 borrow = y >> 32 & (ULong)1; 01068 *xc++ = y & FFFFFFFF; 01069 } 01070 #else 01071 #ifdef Pack_32 01072 do { 01073 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 01074 borrow = (y & 0x10000) >> 16; 01075 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 01076 borrow = (z & 0x10000) >> 16; 01077 Storeinc(xc, z, y); 01078 } 01079 while(xb < xbe); 01080 while(xa < xae) { 01081 y = (*xa & 0xffff) - borrow; 01082 borrow = (y & 0x10000) >> 16; 01083 z = (*xa++ >> 16) - borrow; 01084 borrow = (z & 0x10000) >> 16; 01085 Storeinc(xc, z, y); 01086 } 01087 #else 01088 do { 01089 y = *xa++ - *xb++ - borrow; 01090 borrow = (y & 0x10000) >> 16; 01091 *xc++ = y & 0xffff; 01092 } 01093 while(xb < xbe); 01094 while(xa < xae) { 01095 y = *xa++ - borrow; 01096 borrow = (y & 0x10000) >> 16; 01097 *xc++ = y & 0xffff; 01098 } 01099 #endif 01100 #endif 01101 while(!*--xc) 01102 wa--; 01103 c->wds = wa; 01104 return c; 01105 } 01106 01107 static double 01108 ulp 01109 #ifdef KR_headers 01110 (x) double x; 01111 #else 01112 (double x) 01113 #endif 01114 { 01115 register Long L; 01116 double a; 01117 01118 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 01119 #ifndef Avoid_Underflow 01120 #ifndef Sudden_Underflow 01121 if (L > 0) { 01122 #endif 01123 #endif 01124 #ifdef IBM 01125 L |= Exp_msk1 >> 4; 01126 #endif 01127 word0(a) = L; 01128 word1(a) = 0; 01129 #ifndef Avoid_Underflow 01130 #ifndef Sudden_Underflow 01131 } 01132 else { 01133 L = -L >> Exp_shift; 01134 if (L < Exp_shift) { 01135 word0(a) = 0x80000 >> L; 01136 word1(a) = 0; 01137 } 01138 else { 01139 word0(a) = 0; 01140 L -= Exp_shift; 01141 word1(a) = L >= 31 ? 1 : 1 << 31 - L; 01142 } 01143 } 01144 #endif 01145 #endif 01146 return dval(a); 01147 } 01148 01149 static double 01150 b2d 01151 #ifdef KR_headers 01152 (a, e) Bigint *a; int *e; 01153 #else 01154 (Bigint *a, int *e) 01155 #endif 01156 { 01157 ULong *xa, *xa0, w, y, z; 01158 int k; 01159 double d; 01160 #ifdef VAX 01161 ULong d0, d1; 01162 #else 01163 #define d0 word0(d) 01164 #define d1 word1(d) 01165 #endif 01166 01167 xa0 = a->x; 01168 xa = xa0 + a->wds; 01169 y = *--xa; 01170 #ifdef DEBUG 01171 if (!y) Bug("zero y in b2d"); 01172 #endif 01173 k = hi0bits(y); 01174 *e = 32 - k; 01175 #ifdef Pack_32 01176 if (k < Ebits) { 01177 d0 = Exp_1 | y >> Ebits - k; 01178 w = xa > xa0 ? *--xa : 0; 01179 d1 = y << (32-Ebits) + k | w >> Ebits - k; 01180 goto ret_d; 01181 } 01182 z = xa > xa0 ? *--xa : 0; 01183 if (k -= Ebits) { 01184 d0 = Exp_1 | y << k | z >> 32 - k; 01185 y = xa > xa0 ? *--xa : 0; 01186 d1 = z << k | y >> 32 - k; 01187 } 01188 else { 01189 d0 = Exp_1 | y; 01190 d1 = z; 01191 } 01192 #else 01193 if (k < Ebits + 16) { 01194 z = xa > xa0 ? *--xa : 0; 01195 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 01196 w = xa > xa0 ? *--xa : 0; 01197 y = xa > xa0 ? *--xa : 0; 01198 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 01199 goto ret_d; 01200 } 01201 z = xa > xa0 ? *--xa : 0; 01202 w = xa > xa0 ? *--xa : 0; 01203 k -= Ebits + 16; 01204 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 01205 y = xa > xa0 ? *--xa : 0; 01206 d1 = w << k + 16 | y << k; 01207 #endif 01208 ret_d: 01209 #ifdef VAX 01210 word0(d) = d0 >> 16 | d0 << 16; 01211 word1(d) = d1 >> 16 | d1 << 16; 01212 #else 01213 #undef d0 01214 #undef d1 01215 #endif 01216 return dval(d); 01217 } 01218 01219 static Bigint * 01220 d2b 01221 #ifdef KR_headers 01222 (d, e, bits) double d; int *e, *bits; 01223 #else 01224 (double d, int *e, int *bits) 01225 #endif 01226 { 01227 Bigint *b; 01228 int de, k; 01229 ULong *x, y, z; 01230 #ifndef Sudden_Underflow 01231 int i; 01232 #endif 01233 #ifdef VAX 01234 ULong d0, d1; 01235 d0 = word0(d) >> 16 | word0(d) << 16; 01236 d1 = word1(d) >> 16 | word1(d) << 16; 01237 #else 01238 #define d0 word0(d) 01239 #define d1 word1(d) 01240 #endif 01241 01242 #ifdef Pack_32 01243 b = Balloc(1); 01244 #else 01245 b = Balloc(2); 01246 #endif 01247 x = b->x; 01248 01249 z = d0 & Frac_mask; 01250 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 01251 #ifdef Sudden_Underflow 01252 de = (int)(d0 >> Exp_shift); 01253 #ifndef IBM 01254 z |= Exp_msk11; 01255 #endif 01256 #else 01257 if ((de = (int)(d0 >> Exp_shift))) 01258 z |= Exp_msk1; 01259 #endif 01260 #ifdef Pack_32 01261 if ((y = d1)) { 01262 if ((k = lo0bits(&y))) { 01263 x[0] = y | z << 32 - k; 01264 z >>= k; 01265 } 01266 else 01267 x[0] = y; 01268 #ifndef Sudden_Underflow 01269 i = 01270 #endif 01271 b->wds = (x[1] = z) ? 2 : 1; 01272 } 01273 else { 01274 #ifdef DEBUG 01275 if (!z) 01276 Bug("Zero passed to d2b"); 01277 #endif 01278 k = lo0bits(&z); 01279 x[0] = z; 01280 #ifndef Sudden_Underflow 01281 i = 01282 #endif 01283 b->wds = 1; 01284 k += 32; 01285 } 01286 #else 01287 if (y = d1) { 01288 if (k = lo0bits(&y)) 01289 if (k >= 16) { 01290 x[0] = y | z << 32 - k & 0xffff; 01291 x[1] = z >> k - 16 & 0xffff; 01292 x[2] = z >> k; 01293 i = 2; 01294 } 01295 else { 01296 x[0] = y & 0xffff; 01297 x[1] = y >> 16 | z << 16 - k & 0xffff; 01298 x[2] = z >> k & 0xffff; 01299 x[3] = z >> k+16; 01300 i = 3; 01301 } 01302 else { 01303 x[0] = y & 0xffff; 01304 x[1] = y >> 16; 01305 x[2] = z & 0xffff; 01306 x[3] = z >> 16; 01307 i = 3; 01308 } 01309 } 01310 else { 01311 #ifdef DEBUG 01312 if (!z) 01313 Bug("Zero passed to d2b"); 01314 #endif 01315 k = lo0bits(&z); 01316 if (k >= 16) { 01317 x[0] = z; 01318 i = 0; 01319 } 01320 else { 01321 x[0] = z & 0xffff; 01322 x[1] = z >> 16; 01323 i = 1; 01324 } 01325 k += 32; 01326 } 01327 while(!x[i]) 01328 --i; 01329 b->wds = i + 1; 01330 #endif 01331 #ifndef Sudden_Underflow 01332 if (de) { 01333 #endif 01334 #ifdef IBM 01335 *e = (de - Bias - (P-1) << 2) + k; 01336 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 01337 #else 01338 *e = de - Bias - (P-1) + k; 01339 *bits = P - k; 01340 #endif 01341 #ifndef Sudden_Underflow 01342 } 01343 else { 01344 *e = de - Bias - (P-1) + 1 + k; 01345 #ifdef Pack_32 01346 *bits = 32*i - hi0bits(x[i-1]); 01347 #else 01348 *bits = (i+2)*16 - hi0bits(x[i]); 01349 #endif 01350 } 01351 #endif 01352 return b; 01353 } 01354 #undef d0 01355 #undef d1 01356 01357 static double 01358 ratio 01359 #ifdef KR_headers 01360 (a, b) Bigint *a, *b; 01361 #else 01362 (Bigint *a, Bigint *b) 01363 #endif 01364 { 01365 double da, db; 01366 int k, ka, kb; 01367 01368 dval(da) = b2d(a, &ka); 01369 dval(db) = b2d(b, &kb); 01370 #ifdef Pack_32 01371 k = ka - kb + 32*(a->wds - b->wds); 01372 #else 01373 k = ka - kb + 16*(a->wds - b->wds); 01374 #endif 01375 #ifdef IBM 01376 if (k > 0) { 01377 word0(da) += (k >> 2)*Exp_msk1; 01378 if (k &= 3) 01379 dval(da) *= 1 << k; 01380 } 01381 else { 01382 k = -k; 01383 word0(db) += (k >> 2)*Exp_msk1; 01384 if (k &= 3) 01385 dval(db) *= 1 << k; 01386 } 01387 #else 01388 if (k > 0) 01389 word0(da) += k*Exp_msk1; 01390 else { 01391 k = -k; 01392 word0(db) += k*Exp_msk1; 01393 } 01394 #endif 01395 return dval(da) / dval(db); 01396 } 01397 01398 static CONST double 01399 tens[] = { 01400 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 01401 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 01402 1e20, 1e21, 1e22 01403 #ifdef VAX 01404 , 1e23, 1e24 01405 #endif 01406 }; 01407 01408 static CONST double 01409 #ifdef IEEE_Arith 01410 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 01411 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 01412 #ifdef Avoid_Underflow 01413 9007199254740992.*9007199254740992.e-256 01414 /* = 2^106 * 1e-53 */ 01415 #else 01416 1e-256 01417 #endif 01418 }; 01419 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 01420 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 01421 #define Scale_Bit 0x10 01422 #define n_bigtens 5 01423 #else 01424 #ifdef IBM 01425 bigtens[] = { 1e16, 1e32, 1e64 }; 01426 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 01427 #define n_bigtens 3 01428 #else 01429 bigtens[] = { 1e16, 1e32 }; 01430 static CONST double tinytens[] = { 1e-16, 1e-32 }; 01431 #define n_bigtens 2 01432 #endif 01433 #endif 01434 01435 #ifndef IEEE_Arith 01436 #undef INFNAN_CHECK 01437 #endif 01438 01439 #ifdef INFNAN_CHECK 01440 01441 #ifndef NAN_WORD0 01442 #define NAN_WORD0 0x7ff80000 01443 #endif 01444 01445 #ifndef NAN_WORD1 01446 #define NAN_WORD1 0 01447 #endif 01448 01449 static int 01450 match 01451 #ifdef KR_headers 01452 (sp, t) char **sp, *t; 01453 #else 01454 (CONST char **sp, CONST char *t) 01455 #endif 01456 { 01457 int c, d; 01458 CONST char *s = *sp; 01459 01460 while((d = *t++)) { 01461 if ((c = *++s) >= 'A' && c <= 'Z') 01462 c += 'a' - 'A'; 01463 if (c != d) 01464 return 0; 01465 } 01466 *sp = s + 1; 01467 return 1; 01468 } 01469 01470 #ifndef No_Hex_NaN 01471 static void 01472 hexnan 01473 #ifdef KR_headers 01474 (rvp, sp) double *rvp; CONST char **sp; 01475 #else 01476 (double *rvp, CONST char **sp) 01477 #endif 01478 { 01479 ULong c, x[2]; 01480 CONST char *s; 01481 int havedig, udx0, xshift; 01482 01483 x[0] = x[1] = 0; 01484 havedig = xshift = 0; 01485 udx0 = 1; 01486 s = *sp; 01487 while((c = *(CONST unsigned char*)++s)) { 01488 if (c >= '0' && c <= '9') 01489 c -= '0'; 01490 else if (c >= 'a' && c <= 'f') 01491 c += 10 - 'a'; 01492 else if (c >= 'A' && c <= 'F') 01493 c += 10 - 'A'; 01494 else if (c <= ' ') { 01495 if (udx0 && havedig) { 01496 udx0 = 0; 01497 xshift = 1; 01498 } 01499 continue; 01500 } 01501 else if (/*(*/ c == ')' && havedig) { 01502 *sp = s + 1; 01503 break; 01504 } 01505 else 01506 return; /* invalid form: don't change *sp */ 01507 havedig = 1; 01508 if (xshift) { 01509 xshift = 0; 01510 x[0] = x[1]; 01511 x[1] = 0; 01512 } 01513 if (udx0) 01514 x[0] = (x[0] << 4) | (x[1] >> 28); 01515 x[1] = (x[1] << 4) | c; 01516 } 01517 if ((x[0] &= 0xfffff) || x[1]) { 01518 word0(*rvp) = Exp_mask | x[0]; 01519 word1(*rvp) = x[1]; 01520 } 01521 } 01522 #endif /*No_Hex_NaN*/ 01523 #endif /* INFNAN_CHECK */ 01524 01525 double 01526 kjs_strtod 01527 #ifdef KR_headers 01528 (s00, se) CONST char *s00; char **se; 01529 #else 01530 (CONST char *s00, char **se) 01531 #endif 01532 { 01533 #ifdef Avoid_Underflow 01534 int scale; 01535 #endif 01536 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, 01537 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 01538 CONST char *s, *s0, *s1; 01539 double aadj, aadj1, adj, rv, rv0; 01540 Long L; 01541 ULong y, z; 01542 Bigint *bb = NULL, *bb1 = NULL, *bd = NULL, *bd0 = NULL, *bs = NULL, *delta = NULL; 01543 #ifdef SET_INEXACT 01544 int inexact, oldinexact; 01545 #endif 01546 #ifdef Honor_FLT_ROUNDS 01547 int rounding; 01548 #endif 01549 #ifdef USE_LOCALE 01550 CONST char *s2; 01551 #endif 01552 01553 sign = nz0 = nz = 0; 01554 dval(rv) = 0.; 01555 for(s = s00;;s++) switch(*s) { 01556 case '-': 01557 sign = 1; 01558 /* no break */ 01559 case '+': 01560 if (*++s) 01561 goto break2; 01562 /* no break */ 01563 case 0: 01564 goto ret0; 01565 case '\t': 01566 case '\n': 01567 case '\v': 01568 case '\f': 01569 case '\r': 01570 case ' ': 01571 continue; 01572 default: 01573 goto break2; 01574 } 01575 break2: 01576 if (*s == '0') { 01577 nz0 = 1; 01578 while(*++s == '0') ; 01579 if (!*s) 01580 goto ret; 01581 } 01582 s0 = s; 01583 y = z = 0; 01584 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 01585 if (nd < 9) 01586 y = 10*y + c - '0'; 01587 else if (nd < 16) 01588 z = 10*z + c - '0'; 01589 nd0 = nd; 01590 #ifdef USE_LOCALE 01591 s1 = localeconv()->decimal_point; 01592 if (c == *s1) { 01593 c = '.'; 01594 if (*++s1) { 01595 s2 = s; 01596 for(;;) { 01597 if (*++s2 != *s1) { 01598 c = 0; 01599 break; 01600 } 01601 if (!*++s1) { 01602 s = s2; 01603 break; 01604 } 01605 } 01606 } 01607 } 01608 #endif 01609 if (c == '.') { 01610 c = *++s; 01611 if (!nd) { 01612 for(; c == '0'; c = *++s) 01613 nz++; 01614 if (c > '0' && c <= '9') { 01615 s0 = s; 01616 nf += nz; 01617 nz = 0; 01618 goto have_dig; 01619 } 01620 goto dig_done; 01621 } 01622 for(; c >= '0' && c <= '9'; c = *++s) { 01623 have_dig: 01624 nz++; 01625 if (c -= '0') { 01626 nf += nz; 01627 for(i = 1; i < nz; i++) 01628 if (nd++ < 9) 01629 y *= 10; 01630 else if (nd <= DBL_DIG + 1) 01631 z *= 10; 01632 if (nd++ < 9) 01633 y = 10*y + c; 01634 else if (nd <= DBL_DIG + 1) 01635 z = 10*z + c; 01636 nz = 0; 01637 } 01638 } 01639 } 01640 dig_done: 01641 e = 0; 01642 if (c == 'e' || c == 'E') { 01643 if (!nd && !nz && !nz0) { 01644 goto ret0; 01645 } 01646 s00 = s; 01647 esign = 0; 01648 switch(c = *++s) { 01649 case '-': 01650 esign = 1; 01651 case '+': 01652 c = *++s; 01653 } 01654 if (c >= '0' && c <= '9') { 01655 while(c == '0') 01656 c = *++s; 01657 if (c > '0' && c <= '9') { 01658 L = c - '0'; 01659 s1 = s; 01660 while((c = *++s) >= '0' && c <= '9') 01661 L = 10*L + c - '0'; 01662 if (s - s1 > 8 || L > 19999) 01663 /* Avoid confusion from exponents 01664 * so large that e might overflow. 01665 */ 01666 e = 19999; /* safe for 16 bit ints */ 01667 else 01668 e = (int)L; 01669 if (esign) 01670 e = -e; 01671 } 01672 else 01673 e = 0; 01674 } 01675 else 01676 s = s00; 01677 } 01678 if (!nd) { 01679 if (!nz && !nz0) { 01680 #ifdef INFNAN_CHECK 01681 /* Check for Nan and Infinity */ 01682 switch(c) { 01683 case 'i': 01684 case 'I': 01685 if (match(&s,"nf")) { 01686 --s; 01687 if (!match(&s,"inity")) 01688 ++s; 01689 word0(rv) = 0x7ff00000; 01690 word1(rv) = 0; 01691 goto ret; 01692 } 01693 break; 01694 case 'n': 01695 case 'N': 01696 if (match(&s, "an")) { 01697 word0(rv) = NAN_WORD0; 01698 word1(rv) = NAN_WORD1; 01699 #ifndef No_Hex_NaN 01700 if (*s == '(') /*)*/ 01701 hexnan(&rv, &s); 01702 #endif 01703 goto ret; 01704 } 01705 } 01706 #endif /* INFNAN_CHECK */ 01707 ret0: 01708 s = s00; 01709 sign = 0; 01710 } 01711 goto ret; 01712 } 01713 e1 = e -= nf; 01714 01715 /* Now we have nd0 digits, starting at s0, followed by a 01716 * decimal point, followed by nd-nd0 digits. The number we're 01717 * after is the integer represented by those digits times 01718 * 10**e */ 01719 01720 if (!nd0) 01721 nd0 = nd; 01722 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 01723 dval(rv) = y; 01724 if (k > 9) { 01725 #ifdef SET_INEXACT 01726 if (k > DBL_DIG) 01727 oldinexact = get_inexact(); 01728 #endif 01729 dval(rv) = tens[k - 9] * dval(rv) + z; 01730 } 01731 bd0 = 0; 01732 if (nd <= DBL_DIG 01733 #ifndef RND_PRODQUOT 01734 #ifndef Honor_FLT_ROUNDS 01735 && Flt_Rounds == 1 01736 #endif 01737 #endif 01738 ) { 01739 if (!e) 01740 goto ret; 01741 if (e > 0) { 01742 if (e <= Ten_pmax) { 01743 #ifdef VAX 01744 goto vax_ovfl_check; 01745 #else 01746 #ifdef Honor_FLT_ROUNDS 01747 /* round correctly FLT_ROUNDS = 2 or 3 */ 01748 if (sign) { 01749 rv = -rv; 01750 sign = 0; 01751 } 01752 #endif 01753 /* rv = */ rounded_product(dval(rv), tens[e]); 01754 goto ret; 01755 #endif 01756 } 01757 i = DBL_DIG - nd; 01758 if (e <= Ten_pmax + i) { 01759 /* A fancier test would sometimes let us do 01760 * this for larger i values. 01761 */ 01762 #ifdef Honor_FLT_ROUNDS 01763 /* round correctly FLT_ROUNDS = 2 or 3 */ 01764 if (sign) { 01765 rv = -rv; 01766 sign = 0; 01767 } 01768 #endif 01769 e -= i; 01770 dval(rv) *= tens[i]; 01771 #ifdef VAX 01772 /* VAX exponent range is so narrow we must 01773 * worry about overflow here... 01774 */ 01775 vax_ovfl_check: 01776 word0(rv) -= P*Exp_msk1; 01777 /* rv = */ rounded_product(dval(rv), tens[e]); 01778 if ((word0(rv) & Exp_mask) 01779 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 01780 goto ovfl; 01781 word0(rv) += P*Exp_msk1; 01782 #else 01783 /* rv = */ rounded_product(dval(rv), tens[e]); 01784 #endif 01785 goto ret; 01786 } 01787 } 01788 #ifndef Inaccurate_Divide 01789 else if (e >= -Ten_pmax) { 01790 #ifdef Honor_FLT_ROUNDS 01791 /* round correctly FLT_ROUNDS = 2 or 3 */ 01792 if (sign) { 01793 rv = -rv; 01794 sign = 0; 01795 } 01796 #endif 01797 /* rv = */ rounded_quotient(dval(rv), tens[-e]); 01798 goto ret; 01799 } 01800 #endif 01801 } 01802 e1 += nd - k; 01803 01804 #ifdef IEEE_Arith 01805 #ifdef SET_INEXACT 01806 inexact = 1; 01807 if (k <= DBL_DIG) 01808 oldinexact = get_inexact(); 01809 #endif 01810 #ifdef Avoid_Underflow 01811 scale = 0; 01812 #endif 01813 #ifdef Honor_FLT_ROUNDS 01814 if ((rounding = Flt_Rounds) >= 2) { 01815 if (sign) 01816 rounding = rounding == 2 ? 0 : 2; 01817 else 01818 if (rounding != 2) 01819 rounding = 0; 01820 } 01821 #endif 01822 #endif /*IEEE_Arith*/ 01823 01824 /* Get starting approximation = rv * 10**e1 */ 01825 01826 if (e1 > 0) { 01827 if ((i = e1 & 15)) 01828 dval(rv) *= tens[i]; 01829 if (e1 &= ~15) { 01830 if (e1 > DBL_MAX_10_EXP) { 01831 ovfl: 01832 #ifndef NO_ERRNO 01833 errno = ERANGE; 01834 #endif 01835 /* Can't trust HUGE_VAL */ 01836 #ifdef IEEE_Arith 01837 #ifdef Honor_FLT_ROUNDS 01838 switch(rounding) { 01839 case 0: /* toward 0 */ 01840 case 3: /* toward -infinity */ 01841 word0(rv) = Big0; 01842 word1(rv) = Big1; 01843 break; 01844 default: 01845 word0(rv) = Exp_mask; 01846 word1(rv) = 0; 01847 } 01848 #else /*Honor_FLT_ROUNDS*/ 01849 word0(rv) = Exp_mask; 01850 word1(rv) = 0; 01851 #endif /*Honor_FLT_ROUNDS*/ 01852 #ifdef SET_INEXACT 01853 /* set overflow bit */ 01854 dval(rv0) = 1e300; 01855 dval(rv0) *= dval(rv0); 01856 #endif 01857 #else /*IEEE_Arith*/ 01858 word0(rv) = Big0; 01859 word1(rv) = Big1; 01860 #endif /*IEEE_Arith*/ 01861 if (bd0) 01862 goto retfree; 01863 goto ret; 01864 } 01865 e1 >>= 4; 01866 for(j = 0; e1 > 1; j++, e1 >>= 1) 01867 if (e1 & 1) 01868 dval(rv) *= bigtens[j]; 01869 /* The last multiplication could overflow. */ 01870 word0(rv) -= P*Exp_msk1; 01871 dval(rv) *= bigtens[j]; 01872 if ((z = word0(rv) & Exp_mask) 01873 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 01874 goto ovfl; 01875 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 01876 /* set to largest number */ 01877 /* (Can't trust DBL_MAX) */ 01878 word0(rv) = Big0; 01879 word1(rv) = Big1; 01880 } 01881 else 01882 word0(rv) += P*Exp_msk1; 01883 } 01884 } 01885 else if (e1 < 0) { 01886 e1 = -e1; 01887 if ((i = e1 & 15)) 01888 dval(rv) /= tens[i]; 01889 if (e1 >>= 4) { 01890 if (e1 >= 1 << n_bigtens) 01891 goto undfl; 01892 #ifdef Avoid_Underflow 01893 if (e1 & Scale_Bit) 01894 scale = 2*P; 01895 for(j = 0; e1 > 0; j++, e1 >>= 1) 01896 if (e1 & 1) 01897 dval(rv) *= tinytens[j]; 01898 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask) 01899 >> Exp_shift)) > 0) { 01900 /* scaled rv is denormal; zap j low bits */ 01901 if (j >= 32) { 01902 word1(rv) = 0; 01903 if (j >= 53) 01904 word0(rv) = (P+2)*Exp_msk1; 01905 else 01906 word0(rv) &= 0xffffffff << j-32; 01907 } 01908 else 01909 word1(rv) &= 0xffffffff << j; 01910 } 01911 #else 01912 for(j = 0; e1 > 1; j++, e1 >>= 1) 01913 if (e1 & 1) 01914 dval(rv) *= tinytens[j]; 01915 /* The last multiplication could underflow. */ 01916 dval(rv0) = dval(rv); 01917 dval(rv) *= tinytens[j]; 01918 if (!dval(rv)) { 01919 dval(rv) = 2.*dval(rv0); 01920 dval(rv) *= tinytens[j]; 01921 #endif 01922 if (!dval(rv)) { 01923 undfl: 01924 dval(rv) = 0.; 01925 #ifndef NO_ERRNO 01926 errno = ERANGE; 01927 #endif 01928 if (bd0) 01929 goto retfree; 01930 goto ret; 01931 } 01932 #ifndef Avoid_Underflow 01933 word0(rv) = Tiny0; 01934 word1(rv) = Tiny1; 01935 /* The refinement below will clean 01936 * this approximation up. 01937 */ 01938 } 01939 #endif 01940 } 01941 } 01942 01943 /* Now the hard part -- adjusting rv to the correct value.*/ 01944 01945 /* Put digits into bd: true value = bd * 10^e */ 01946 01947 bd0 = s2b(s0, nd0, nd, y); 01948 01949 for(;;) { 01950 bd = Balloc(bd0->k); 01951 Bcopy(bd, bd0); 01952 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 01953 bs = i2b(1); 01954 01955 if (e >= 0) { 01956 bb2 = bb5 = 0; 01957 bd2 = bd5 = e; 01958 } 01959 else { 01960 bb2 = bb5 = -e; 01961 bd2 = bd5 = 0; 01962 } 01963 if (bbe >= 0) 01964 bb2 += bbe; 01965 else 01966 bd2 -= bbe; 01967 bs2 = bb2; 01968 #ifdef Honor_FLT_ROUNDS 01969 if (rounding != 1) 01970 bs2++; 01971 #endif 01972 #ifdef Avoid_Underflow 01973 j = bbe - scale; 01974 i = j + bbbits - 1; /* logb(rv) */ 01975 if (i < Emin) /* denormal */ 01976 j += P - Emin; 01977 else 01978 j = P + 1 - bbbits; 01979 #else /*Avoid_Underflow*/ 01980 #ifdef Sudden_Underflow 01981 #ifdef IBM 01982 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 01983 #else 01984 j = P + 1 - bbbits; 01985 #endif 01986 #else /*Sudden_Underflow*/ 01987 j = bbe; 01988 i = j + bbbits - 1; /* logb(rv) */ 01989 if (i < Emin) /* denormal */ 01990 j += P - Emin; 01991 else 01992 j = P + 1 - bbbits; 01993 #endif /*Sudden_Underflow*/ 01994 #endif /*Avoid_Underflow*/ 01995 bb2 += j; 01996 bd2 += j; 01997 #ifdef Avoid_Underflow 01998 bd2 += scale; 01999 #endif 02000 i = bb2 < bd2 ? bb2 : bd2; 02001 if (i > bs2) 02002 i = bs2; 02003 if (i > 0) { 02004 bb2 -= i; 02005 bd2 -= i; 02006 bs2 -= i; 02007 } 02008 if (bb5 > 0) { 02009 bs = pow5mult(bs, bb5); 02010 bb1 = mult(bs, bb); 02011 Bfree(bb); 02012 bb = bb1; 02013 } 02014 if (bb2 > 0) 02015 bb = lshift(bb, bb2); 02016 if (bd5 > 0) 02017 bd = pow5mult(bd, bd5); 02018 if (bd2 > 0) 02019 bd = lshift(bd, bd2); 02020 if (bs2 > 0) 02021 bs = lshift(bs, bs2); 02022 delta = diff(bb, bd); 02023 dsign = delta->sign; 02024 delta->sign = 0; 02025 i = cmp(delta, bs); 02026 #ifdef Honor_FLT_ROUNDS 02027 if (rounding != 1) { 02028 if (i < 0) { 02029 /* Error is less than an ulp */ 02030 if (!delta->x[0] && delta->wds <= 1) { 02031 /* exact */ 02032 #ifdef SET_INEXACT 02033 inexact = 0; 02034 #endif 02035 break; 02036 } 02037 if (rounding) { 02038 if (dsign) { 02039 adj = 1.; 02040 goto apply_adj; 02041 } 02042 } 02043 else if (!dsign) { 02044 adj = -1.; 02045 if (!word1(rv) 02046 && !(word0(rv) & Frac_mask)) { 02047 y = word0(rv) & Exp_mask; 02048 #ifdef Avoid_Underflow 02049 if (!scale || y > 2*P*Exp_msk1) 02050 #else 02051 if (y) 02052 #endif 02053 { 02054 delta = lshift(delta,Log2P); 02055 if (cmp(delta, bs) <= 0) 02056 adj = -0.5; 02057 } 02058 } 02059 apply_adj: 02060 #ifdef Avoid_Underflow 02061 if (scale && (y = word0(rv) & Exp_mask) 02062 <= 2*P*Exp_msk1) 02063 word0(adj) += (2*P+1)*Exp_msk1 - y; 02064 #else 02065 #ifdef Sudden_Underflow 02066 if ((word0(rv) & Exp_mask) <= 02067 P*Exp_msk1) { 02068 word0(rv) += P*Exp_msk1; 02069 dval(rv) += adj*ulp(dval(rv)); 02070 word0(rv) -= P*Exp_msk1; 02071 } 02072 else 02073 #endif /*Sudden_Underflow*/ 02074 #endif /*Avoid_Underflow*/ 02075 dval(rv) += adj*ulp(dval(rv)); 02076 } 02077 break; 02078 } 02079 adj = ratio(delta, bs); 02080 if (adj < 1.) 02081 adj = 1.; 02082 if (adj <= 0x7ffffffe) { 02083 /* adj = rounding ? ceil(adj) : floor(adj); */ 02084 y = adj; 02085 if (y != adj) { 02086 if (!((rounding>>1) ^ dsign)) 02087 y++; 02088 adj = y; 02089 } 02090 } 02091 #ifdef Avoid_Underflow 02092 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 02093 word0(adj) += (2*P+1)*Exp_msk1 - y; 02094 #else 02095 #ifdef Sudden_Underflow 02096 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 02097 word0(rv) += P*Exp_msk1; 02098 adj *= ulp(dval(rv)); 02099 if (dsign) 02100 dval(rv) += adj; 02101 else 02102 dval(rv) -= adj; 02103 word0(rv) -= P*Exp_msk1; 02104 goto cont; 02105 } 02106 #endif /*Sudden_Underflow*/ 02107 #endif /*Avoid_Underflow*/ 02108 adj *= ulp(dval(rv)); 02109 if (dsign) 02110 dval(rv) += adj; 02111 else 02112 dval(rv) -= adj; 02113 goto cont; 02114 } 02115 #endif /*Honor_FLT_ROUNDS*/ 02116 02117 if (i < 0) { 02118 /* Error is less than half an ulp -- check for 02119 * special case of mantissa a power of two. 02120 */ 02121 if (dsign || word1(rv) || word0(rv) & Bndry_mask 02122 #ifdef IEEE_Arith 02123 #ifdef Avoid_Underflow 02124 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 02125 #else 02126 || (word0(rv) & Exp_mask) <= Exp_msk1 02127 #endif 02128 #endif 02129 ) { 02130 #ifdef SET_INEXACT 02131 if (!delta->x[0] && delta->wds <= 1) 02132 inexact = 0; 02133 #endif 02134 break; 02135 } 02136 if (!delta->x[0] && delta->wds <= 1) { 02137 /* exact result */ 02138 #ifdef SET_INEXACT 02139 inexact = 0; 02140 #endif 02141 break; 02142 } 02143 delta = lshift(delta,Log2P); 02144 if (cmp(delta, bs) > 0) 02145 goto drop_down; 02146 break; 02147 } 02148 if (i == 0) { 02149 /* exactly half-way between */ 02150 if (dsign) { 02151 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 02152 && word1(rv) == ( 02153 #ifdef Avoid_Underflow 02154 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 02155 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 02156 #endif 02157 0xffffffff)) { 02158 /*boundary case -- increment exponent*/ 02159 word0(rv) = (word0(rv) & Exp_mask) 02160 + Exp_msk1 02161 #ifdef IBM 02162 | Exp_msk1 >> 4 02163 #endif 02164 ; 02165 word1(rv) = 0; 02166 #ifdef Avoid_Underflow 02167 dsign = 0; 02168 #endif 02169 break; 02170 } 02171 } 02172 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 02173 drop_down: 02174 /* boundary case -- decrement exponent */ 02175 #ifdef Sudden_Underflow /*{{*/ 02176 L = word0(rv) & Exp_mask; 02177 #ifdef IBM 02178 if (L < Exp_msk1) 02179 #else 02180 #ifdef Avoid_Underflow 02181 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 02182 #else 02183 if (L <= Exp_msk1) 02184 #endif /*Avoid_Underflow*/ 02185 #endif /*IBM*/ 02186 goto undfl; 02187 L -= Exp_msk1; 02188 #else /*Sudden_Underflow}{*/ 02189 #ifdef Avoid_Underflow 02190 if (scale) { 02191 L = word0(rv) & Exp_mask; 02192 if (L <= (2*P+1)*Exp_msk1) { 02193 if (L > (P+2)*Exp_msk1) 02194 /* round even ==> */ 02195 /* accept rv */ 02196 break; 02197 /* rv = smallest denormal */ 02198 goto undfl; 02199 } 02200 } 02201 #endif /*Avoid_Underflow*/ 02202 L = (word0(rv) & Exp_mask) - Exp_msk1; 02203 #endif /*Sudden_Underflow}}*/ 02204 word0(rv) = L | Bndry_mask1; 02205 word1(rv) = 0xffffffff; 02206 #ifdef IBM 02207 goto cont; 02208 #else 02209 break; 02210 #endif 02211 } 02212 #ifndef ROUND_BIASED 02213 if (!(word1(rv) & LSB)) 02214 break; 02215 #endif 02216 if (dsign) 02217 dval(rv) += ulp(dval(rv)); 02218 #ifndef ROUND_BIASED 02219 else { 02220 dval(rv) -= ulp(dval(rv)); 02221 #ifndef Sudden_Underflow 02222 if (!dval(rv)) 02223 goto undfl; 02224 #endif 02225 } 02226 #ifdef Avoid_Underflow 02227 dsign = 1 - dsign; 02228 #endif 02229 #endif 02230 break; 02231 } 02232 if ((aadj = ratio(delta, bs)) <= 2.) { 02233 if (dsign) 02234 aadj = aadj1 = 1.; 02235 else if (word1(rv) || word0(rv) & Bndry_mask) { 02236 #ifndef Sudden_Underflow 02237 if (word1(rv) == Tiny1 && !word0(rv)) 02238 goto undfl; 02239 #endif 02240 aadj = 1.; 02241 aadj1 = -1.; 02242 } 02243 else { 02244 /* special case -- power of FLT_RADIX to be */ 02245 /* rounded down... */ 02246 02247 if (aadj < 2./FLT_RADIX) 02248 aadj = 1./FLT_RADIX; 02249 else 02250 aadj *= 0.5; 02251 aadj1 = -aadj; 02252 } 02253 } 02254 else { 02255 aadj *= 0.5; 02256 aadj1 = dsign ? aadj : -aadj; 02257 #ifdef Check_FLT_ROUNDS 02258 switch(Rounding) { 02259 case 2: /* towards +infinity */ 02260 aadj1 -= 0.5; 02261 break; 02262 case 0: /* towards 0 */ 02263 case 3: /* towards -infinity */ 02264 aadj1 += 0.5; 02265 } 02266 #else 02267 if (Flt_Rounds == 0) 02268 aadj1 += 0.5; 02269 #endif /*Check_FLT_ROUNDS*/ 02270 } 02271 y = word0(rv) & Exp_mask; 02272 02273 /* Check for overflow */ 02274 02275 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 02276 dval(rv0) = dval(rv); 02277 word0(rv) -= P*Exp_msk1; 02278 adj = aadj1 * ulp(dval(rv)); 02279 dval(rv) += adj; 02280 if ((word0(rv) & Exp_mask) >= 02281 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 02282 if (word0(rv0) == Big0 && word1(rv0) == Big1) 02283 goto ovfl; 02284 word0(rv) = Big0; 02285 word1(rv) = Big1; 02286 goto cont; 02287 } 02288 else 02289 word0(rv) += P*Exp_msk1; 02290 } 02291 else { 02292 #ifdef Avoid_Underflow 02293 if (scale && y <= 2*P*Exp_msk1) { 02294 if (aadj <= 0x7fffffff) { 02295 if ((z = (ULong)aadj) <= 0) 02296 z = 1; 02297 aadj = z; 02298 aadj1 = dsign ? aadj : -aadj; 02299 } 02300 word0(aadj1) += (2*P+1)*Exp_msk1 - y; 02301 } 02302 adj = aadj1 * ulp(dval(rv)); 02303 dval(rv) += adj; 02304 #else 02305 #ifdef Sudden_Underflow 02306 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 02307 dval(rv0) = dval(rv); 02308 word0(rv) += P*Exp_msk1; 02309 adj = aadj1 * ulp(dval(rv)); 02310 dval(rv) += adj; 02311 #ifdef IBM 02312 if ((word0(rv) & Exp_mask) < P*Exp_msk1) 02313 #else 02314 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 02315 #endif 02316 { 02317 if (word0(rv0) == Tiny0 02318 && word1(rv0) == Tiny1) 02319 goto undfl; 02320 word0(rv) = Tiny0; 02321 word1(rv) = Tiny1; 02322 goto cont; 02323 } 02324 else 02325 word0(rv) -= P*Exp_msk1; 02326 } 02327 else { 02328 adj = aadj1 * ulp(dval(rv)); 02329 dval(rv) += adj; 02330 } 02331 #else /*Sudden_Underflow*/ 02332 /* Compute adj so that the IEEE rounding rules will 02333 * correctly round rv + adj in some half-way cases. 02334 * If rv * ulp(rv) is denormalized (i.e., 02335 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 02336 * trouble from bits lost to denormalization; 02337 * example: 1.2e-307 . 02338 */ 02339 if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 02340 aadj1 = (double)(int)(aadj + 0.5); 02341 if (!dsign) 02342 aadj1 = -aadj1; 02343 } 02344 adj = aadj1 * ulp(dval(rv)); 02345 dval(rv) += adj; 02346 #endif /*Sudden_Underflow*/ 02347 #endif /*Avoid_Underflow*/ 02348 } 02349 z = word0(rv) & Exp_mask; 02350 #ifndef SET_INEXACT 02351 #ifdef Avoid_Underflow 02352 if (!scale) 02353 #endif 02354 if (y == z) { 02355 /* Can we stop now? */ 02356 L = (Long)aadj; 02357 aadj -= L; 02358 /* The tolerances below are conservative. */ 02359 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 02360 if (aadj < .4999999 || aadj > .5000001) 02361 break; 02362 } 02363 else if (aadj < .4999999/FLT_RADIX) 02364 break; 02365 } 02366 #endif 02367 cont: 02368 Bfree(bb); 02369 Bfree(bd); 02370 Bfree(bs); 02371 Bfree(delta); 02372 } 02373 #ifdef SET_INEXACT 02374 if (inexact) { 02375 if (!oldinexact) { 02376 word0(rv0) = Exp_1 + (70 << Exp_shift); 02377 word1(rv0) = 0; 02378 dval(rv0) += 1.; 02379 } 02380 } 02381 else if (!oldinexact) 02382 clear_inexact(); 02383 #endif 02384 #ifdef Avoid_Underflow 02385 if (scale) { 02386 word0(rv0) = Exp_1 - 2*P*Exp_msk1; 02387 word1(rv0) = 0; 02388 dval(rv) *= dval(rv0); 02389 #ifndef NO_ERRNO 02390 /* try to avoid the bug of testing an 8087 register value */ 02391 if (word0(rv) == 0 && word1(rv) == 0) 02392 errno = ERANGE; 02393 #endif 02394 } 02395 #endif /* Avoid_Underflow */ 02396 #ifdef SET_INEXACT 02397 if (inexact && !(word0(rv) & Exp_mask)) { 02398 /* set underflow bit */ 02399 dval(rv0) = 1e-300; 02400 dval(rv0) *= dval(rv0); 02401 } 02402 #endif 02403 retfree: 02404 Bfree(bb); 02405 Bfree(bd); 02406 Bfree(bs); 02407 Bfree(bd0); 02408 Bfree(delta); 02409 ret: 02410 if (se) 02411 *se = (char *)s; 02412 return sign ? -dval(rv) : dval(rv); 02413 } 02414 02415 static int 02416 quorem 02417 #ifdef KR_headers 02418 (b, S) Bigint *b, *S; 02419 #else 02420 (Bigint *b, Bigint *S) 02421 #endif 02422 { 02423 int n; 02424 ULong *bx, *bxe, q, *sx, *sxe; 02425 #ifdef ULLong 02426 ULLong borrow, carry, y, ys; 02427 #else 02428 ULong borrow, carry, y, ys; 02429 #ifdef Pack_32 02430 ULong si, z, zs; 02431 #endif 02432 #endif 02433 02434 n = S->wds; 02435 #ifdef DEBUG 02436 /*debug*/ if (b->wds > n) 02437 /*debug*/ Bug("oversize b in quorem"); 02438 #endif 02439 if (b->wds < n) 02440 return 0; 02441 sx = S->x; 02442 sxe = sx + --n; 02443 bx = b->x; 02444 bxe = bx + n; 02445 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 02446 #ifdef DEBUG 02447 /*debug*/ if (q > 9) 02448 /*debug*/ Bug("oversized quotient in quorem"); 02449 #endif 02450 if (q) { 02451 borrow = 0; 02452 carry = 0; 02453 do { 02454 #ifdef ULLong 02455 ys = *sx++ * (ULLong)q + carry; 02456 carry = ys >> 32; 02457 y = *bx - (ys & FFFFFFFF) - borrow; 02458 borrow = y >> 32 & (ULong)1; 02459 *bx++ = y & FFFFFFFF; 02460 #else 02461 #ifdef Pack_32 02462 si = *sx++; 02463 ys = (si & 0xffff) * q + carry; 02464 zs = (si >> 16) * q + (ys >> 16); 02465 carry = zs >> 16; 02466 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 02467 borrow = (y & 0x10000) >> 16; 02468 z = (*bx >> 16) - (zs & 0xffff) - borrow; 02469 borrow = (z & 0x10000) >> 16; 02470 Storeinc(bx, z, y); 02471 #else 02472 ys = *sx++ * q + carry; 02473 carry = ys >> 16; 02474 y = *bx - (ys & 0xffff) - borrow; 02475 borrow = (y & 0x10000) >> 16; 02476 *bx++ = y & 0xffff; 02477 #endif 02478 #endif 02479 } 02480 while(sx <= sxe); 02481 if (!*bxe) { 02482 bx = b->x; 02483 while(--bxe > bx && !*bxe) 02484 --n; 02485 b->wds = n; 02486 } 02487 } 02488 if (cmp(b, S) >= 0) { 02489 q++; 02490 borrow = 0; 02491 carry = 0; 02492 bx = b->x; 02493 sx = S->x; 02494 do { 02495 #ifdef ULLong 02496 ys = *sx++ + carry; 02497 carry = ys >> 32; 02498 y = *bx - (ys & FFFFFFFF) - borrow; 02499 borrow = y >> 32 & (ULong)1; 02500 *bx++ = y & FFFFFFFF; 02501 #else 02502 #ifdef Pack_32 02503 si = *sx++; 02504 ys = (si & 0xffff) + carry; 02505 zs = (si >> 16) + (ys >> 16); 02506 carry = zs >> 16; 02507 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 02508 borrow = (y & 0x10000) >> 16; 02509 z = (*bx >> 16) - (zs & 0xffff) - borrow; 02510 borrow = (z & 0x10000) >> 16; 02511 Storeinc(bx, z, y); 02512 #else 02513 ys = *sx++ + carry; 02514 carry = ys >> 16; 02515 y = *bx - (ys & 0xffff) - borrow; 02516 borrow = (y & 0x10000) >> 16; 02517 *bx++ = y & 0xffff; 02518 #endif 02519 #endif 02520 } 02521 while(sx <= sxe); 02522 bx = b->x; 02523 bxe = bx + n; 02524 if (!*bxe) { 02525 while(--bxe > bx && !*bxe) 02526 --n; 02527 b->wds = n; 02528 } 02529 } 02530 return q; 02531 } 02532 02533 #ifndef MULTIPLE_THREADS 02534 static char *dtoa_result; 02535 #endif 02536 02537 static char * 02538 #ifdef KR_headers 02539 rv_alloc(i) int i; 02540 #else 02541 rv_alloc(int i) 02542 #endif 02543 { 02544 int j, k, *r; 02545 02546 j = sizeof(ULong); 02547 for(k = 0; 02548 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i; 02549 j <<= 1) 02550 k++; 02551 r = (int*)Balloc(k); 02552 *r = k; 02553 return 02554 #ifndef MULTIPLE_THREADS 02555 dtoa_result = 02556 #endif 02557 (char *)(r+1); 02558 } 02559 02560 static char * 02561 #ifdef KR_headers 02562 nrv_alloc(s, rve, n) char *s, **rve; int n; 02563 #else 02564 nrv_alloc(CONST char *s, char **rve, int n) 02565 #endif 02566 { 02567 char *rv, *t; 02568 02569 t = rv = rv_alloc(n); 02570 while((*t = *s++)) t++; 02571 if (rve) 02572 *rve = t; 02573 return rv; 02574 } 02575 02576 /* freedtoa(s) must be used to free values s returned by dtoa 02577 * when MULTIPLE_THREADS is #defined. It should be used in all cases, 02578 * but for consistency with earlier versions of dtoa, it is optional 02579 * when MULTIPLE_THREADS is not defined. 02580 */ 02581 02582 void 02583 #ifdef KR_headers 02584 kjs_freedtoa(s) char *s; 02585 #else 02586 kjs_freedtoa(char *s) 02587 #endif 02588 { 02589 Bigint *b = (Bigint *)((int *)s - 1); 02590 b->maxwds = 1 << (b->k = *(int*)b); 02591 Bfree(b); 02592 #ifndef MULTIPLE_THREADS 02593 if (s == dtoa_result) 02594 dtoa_result = 0; 02595 #endif 02596 } 02597 02598 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 02599 * 02600 * Inspired by "How to Print Floating-Point Numbers Accurately" by 02601 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101]. 02602 * 02603 * Modifications: 02604 * 1. Rather than iterating, we use a simple numeric overestimate 02605 * to determine k = floor(log10(d)). We scale relevant 02606 * quantities using O(log2(k)) rather than O(k) multiplications. 02607 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 02608 * try to generate digits strictly left to right. Instead, we 02609 * compute with fewer bits and propagate the carry if necessary 02610 * when rounding the final digit up. This is often faster. 02611 * 3. Under the assumption that input will be rounded nearest, 02612 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 02613 * That is, we allow equality in stopping tests when the 02614 * round-nearest rule will give the same floating-point value 02615 * as would satisfaction of the stopping test with strict 02616 * inequality. 02617 * 4. We remove common factors of powers of 2 from relevant 02618 * quantities. 02619 * 5. When converting floating-point integers less than 1e16, 02620 * we use floating-point arithmetic rather than resorting 02621 * to multiple-precision integers. 02622 * 6. When asked to produce fewer than 15 digits, we first try 02623 * to get by with floating-point arithmetic; we resort to 02624 * multiple-precision integer arithmetic only if we cannot 02625 * guarantee that the floating-point calculation has given 02626 * the correctly rounded result. For k requested digits and 02627 * "uniformly" distributed input, the probability is 02628 * something like 10^(k-15) that we must resort to the Long 02629 * calculation. 02630 */ 02631 02632 char * 02633 kjs_dtoa 02634 #ifdef KR_headers 02635 (d, mode, ndigits, decpt, sign, rve) 02636 double d; int mode, ndigits, *decpt, *sign; char **rve; 02637 #else 02638 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve) 02639 #endif 02640 { 02641 /* Arguments ndigits, decpt, sign are similar to those 02642 of ecvt and fcvt; trailing zeros are suppressed from 02643 the returned string. If not null, *rve is set to point 02644 to the end of the return value. If d is +-Infinity or NaN, 02645 then *decpt is set to 9999. 02646 02647 mode: 02648 0 ==> shortest string that yields d when read in 02649 and rounded to nearest. 02650 1 ==> like 0, but with Steele & White stopping rule; 02651 e.g. with IEEE P754 arithmetic , mode 0 gives 02652 1e23 whereas mode 1 gives 9.999999999999999e22. 02653 2 ==> max(1,ndigits) significant digits. This gives a 02654 return value similar to that of ecvt, except 02655 that trailing zeros are suppressed. 02656 3 ==> through ndigits past the decimal point. This 02657 gives a return value similar to that from fcvt, 02658 except that trailing zeros are suppressed, and 02659 ndigits can be negative. 02660 4,5 ==> similar to 2 and 3, respectively, but (in 02661 round-nearest mode) with the tests of mode 0 to 02662 possibly return a shorter string that rounds to d. 02663 With IEEE arithmetic and compilation with 02664 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same 02665 as modes 2 and 3 when FLT_ROUNDS != 1. 02666 6-9 ==> Debugging modes similar to mode - 4: don't try 02667 fast floating-point estimate (if applicable). 02668 02669 Values of mode other than 0-9 are treated as mode 0. 02670 02671 Sufficient space is allocated to the return value 02672 to hold the suppressed trailing zeros. 02673 */ 02674 02675 int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0, 02676 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, 02677 spec_case, try_quick; 02678 Long L; 02679 #ifndef Sudden_Underflow 02680 int denorm; 02681 ULong x; 02682 #endif 02683 Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S; 02684 double d2, ds, eps; 02685 char *s, *s0; 02686 #ifdef Honor_FLT_ROUNDS 02687 int rounding; 02688 #endif 02689 #ifdef SET_INEXACT 02690 int inexact, oldinexact; 02691 #endif 02692 02693 #ifndef MULTIPLE_THREADS 02694 if (dtoa_result) { 02695 kjs_freedtoa(dtoa_result); 02696 dtoa_result = 0; 02697 } 02698 #endif 02699 02700 if (word0(d) & Sign_bit) { 02701 /* set sign for everything, including 0's and NaNs */ 02702 *sign = 1; 02703 word0(d) &= ~Sign_bit; /* clear sign bit */ 02704 } 02705 else 02706 *sign = 0; 02707 02708 #if defined(IEEE_Arith) + defined(VAX) 02709 #ifdef IEEE_Arith 02710 if ((word0(d) & Exp_mask) == Exp_mask) 02711 #else 02712 if (word0(d) == 0x8000) 02713 #endif 02714 { 02715 /* Infinity or NaN */ 02716 *decpt = 9999; 02717 #ifdef IEEE_Arith 02718 if (!word1(d) && !(word0(d) & 0xfffff)) 02719 return nrv_alloc("Infinity", rve, 8); 02720 #endif 02721 return nrv_alloc("NaN", rve, 3); 02722 } 02723 #endif 02724 #ifdef IBM 02725 dval(d) += 0; /* normalize */ 02726 #endif 02727 if (!dval(d)) { 02728 *decpt = 1; 02729 return nrv_alloc("0", rve, 1); 02730 } 02731 02732 #ifdef SET_INEXACT 02733 try_quick = oldinexact = get_inexact(); 02734 inexact = 1; 02735 #endif 02736 #ifdef Honor_FLT_ROUNDS 02737 if ((rounding = Flt_Rounds) >= 2) { 02738 if (*sign) 02739 rounding = rounding == 2 ? 0 : 2; 02740 else 02741 if (rounding != 2) 02742 rounding = 0; 02743 } 02744 #endif 02745 02746 b = d2b(dval(d), &be, &bbits); 02747 #ifdef Sudden_Underflow 02748 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); 02749 #else 02750 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) { 02751 #endif 02752 dval(d2) = dval(d); 02753 word0(d2) &= Frac_mask1; 02754 word0(d2) |= Exp_11; 02755 #ifdef IBM 02756 if (j = 11 - hi0bits(word0(d2) & Frac_mask)) 02757 dval(d2) /= 1 << j; 02758 #endif 02759 02760 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 02761 * log10(x) = log(x) / log(10) 02762 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 02763 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 02764 * 02765 * This suggests computing an approximation k to log10(d) by 02766 * 02767 * k = (i - Bias)*0.301029995663981 02768 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 02769 * 02770 * We want k to be too large rather than too small. 02771 * The error in the first-order Taylor series approximation 02772 * is in our favor, so we just round up the constant enough 02773 * to compensate for any error in the multiplication of 02774 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 02775 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 02776 * adding 1e-13 to the constant term more than suffices. 02777 * Hence we adjust the constant term to 0.1760912590558. 02778 * (We could get a more accurate k by invoking log10, 02779 * but this is probably not worthwhile.) 02780 */ 02781 02782 i -= Bias; 02783 #ifdef IBM 02784 i <<= 2; 02785 i += j; 02786 #endif 02787 #ifndef Sudden_Underflow 02788 denorm = 0; 02789 } 02790 else { 02791 /* d is denormalized */ 02792 02793 i = bbits + be + (Bias + (P-1) - 1); 02794 x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32 02795 : word1(d) << 32 - i; 02796 dval(d2) = x; 02797 word0(d2) -= 31*Exp_msk1; /* adjust exponent */ 02798 i -= (Bias + (P-1) - 1) + 1; 02799 denorm = 1; 02800 } 02801 #endif 02802 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; 02803 k = (int)ds; 02804 if (ds < 0. && ds != k) 02805 k--; /* want k = floor(ds) */ 02806 k_check = 1; 02807 if (k >= 0 && k <= Ten_pmax) { 02808 if (dval(d) < tens[k]) 02809 k--; 02810 k_check = 0; 02811 } 02812 j = bbits - i - 1; 02813 if (j >= 0) { 02814 b2 = 0; 02815 s2 = j; 02816 } 02817 else { 02818 b2 = -j; 02819 s2 = 0; 02820 } 02821 if (k >= 0) { 02822 b5 = 0; 02823 s5 = k; 02824 s2 += k; 02825 } 02826 else { 02827 b2 -= k; 02828 b5 = -k; 02829 s5 = 0; 02830 } 02831 if (mode < 0 || mode > 9) 02832 mode = 0; 02833 02834 #ifndef SET_INEXACT 02835 #ifdef Check_FLT_ROUNDS 02836 try_quick = Rounding == 1; 02837 #else 02838 try_quick = 1; 02839 #endif 02840 #endif /*SET_INEXACT*/ 02841 02842 if (mode > 5) { 02843 mode -= 4; 02844 try_quick = 0; 02845 } 02846 leftright = 1; 02847 switch(mode) { 02848 case 0: 02849 case 1: 02850 ilim = ilim1 = -1; 02851 i = 18; 02852 ndigits = 0; 02853 break; 02854 case 2: 02855 leftright = 0; 02856 /* no break */ 02857 case 4: 02858 if (ndigits <= 0) 02859 ndigits = 1; 02860 ilim = ilim1 = i = ndigits; 02861 break; 02862 case 3: 02863 leftright = 0; 02864 /* no break */ 02865 case 5: 02866 i = ndigits + k + 1; 02867 ilim = i; 02868 ilim1 = i - 1; 02869 if (i <= 0) 02870 i = 1; 02871 } 02872 s = s0 = rv_alloc(i); 02873 02874 #ifdef Honor_FLT_ROUNDS 02875 if (mode > 1 && rounding != 1) 02876 leftright = 0; 02877 #endif 02878 02879 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 02880 02881 /* Try to get by with floating-point arithmetic. */ 02882 02883 i = 0; 02884 dval(d2) = dval(d); 02885 k0 = k; 02886 ilim0 = ilim; 02887 ieps = 2; /* conservative */ 02888 if (k > 0) { 02889 ds = tens[k&0xf]; 02890 j = k >> 4; 02891 if (j & Bletch) { 02892 /* prevent overflows */ 02893 j &= Bletch - 1; 02894 dval(d) /= bigtens[n_bigtens-1]; 02895 ieps++; 02896 } 02897 for(; j; j >>= 1, i++) 02898 if (j & 1) { 02899 ieps++; 02900 ds *= bigtens[i]; 02901 } 02902 dval(d) /= ds; 02903 } 02904 else if ((j1 = -k)) { 02905 dval(d) *= tens[j1 & 0xf]; 02906 for(j = j1 >> 4; j; j >>= 1, i++) 02907 if (j & 1) { 02908 ieps++; 02909 dval(d) *= bigtens[i]; 02910 } 02911 } 02912 if (k_check && dval(d) < 1. && ilim > 0) { 02913 if (ilim1 <= 0) 02914 goto fast_failed; 02915 ilim = ilim1; 02916 k--; 02917 dval(d) *= 10.; 02918 ieps++; 02919 } 02920 dval(eps) = ieps*dval(d) + 7.; 02921 word0(eps) -= (P-1)*Exp_msk1; 02922 if (ilim == 0) { 02923 S = mhi = 0; 02924 dval(d) -= 5.; 02925 if (dval(d) > dval(eps)) 02926 goto one_digit; 02927 if (dval(d) < -dval(eps)) 02928 goto no_digits; 02929 goto fast_failed; 02930 } 02931 #ifndef No_leftright 02932 if (leftright) { 02933 /* Use Steele & White method of only 02934 * generating digits needed. 02935 */ 02936 dval(eps) = 0.5/tens[ilim-1] - dval(eps); 02937 for(i = 0;;) { 02938 L = (long int)dval(d); 02939 dval(d) -= L; 02940 *s++ = '0' + (int)L; 02941 if (dval(d) < dval(eps)) 02942 goto ret1; 02943 if (1. - dval(d) < dval(eps)) 02944 goto bump_up; 02945 if (++i >= ilim) 02946 break; 02947 dval(eps) *= 10.; 02948 dval(d) *= 10.; 02949 } 02950 } 02951 else { 02952 #endif 02953 /* Generate ilim digits, then fix them up. */ 02954 dval(eps) *= tens[ilim-1]; 02955 for(i = 1;; i++, dval(d) *= 10.) { 02956 L = (Long)(dval(d)); 02957 if (!(dval(d) -= L)) 02958 ilim = i; 02959 *s++ = '0' + (int)L; 02960 if (i == ilim) { 02961 if (dval(d) > 0.5 + dval(eps)) 02962 goto bump_up; 02963 else if (dval(d) < 0.5 - dval(eps)) { 02964 while(*--s == '0'); 02965 s++; 02966 goto ret1; 02967 } 02968 break; 02969 } 02970 } 02971 #ifndef No_leftright 02972 } 02973 #endif 02974 fast_failed: 02975 s = s0; 02976 dval(d) = dval(d2); 02977 k = k0; 02978 ilim = ilim0; 02979 } 02980 02981 /* Do we have a "small" integer? */ 02982 02983 if (be >= 0 && k <= Int_max) { 02984 /* Yes. */ 02985 ds = tens[k]; 02986 if (ndigits < 0 && ilim <= 0) { 02987 S = mhi = 0; 02988 if (ilim < 0 || dval(d) <= 5*ds) 02989 goto no_digits; 02990 goto one_digit; 02991 } 02992 for(i = 1;; i++, dval(d) *= 10.) { 02993 L = (Long)(dval(d) / ds); 02994 dval(d) -= L*ds; 02995 #ifdef Check_FLT_ROUNDS 02996 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 02997 if (dval(d) < 0) { 02998 L--; 02999 dval(d) += ds; 03000 } 03001 #endif 03002 *s++ = '0' + (int)L; 03003 if (!dval(d)) { 03004 #ifdef SET_INEXACT 03005 inexact = 0; 03006 #endif 03007 break; 03008 } 03009 if (i == ilim) { 03010 #ifdef Honor_FLT_ROUNDS 03011 if (mode > 1) 03012 switch(rounding) { 03013 case 0: goto ret1; 03014 case 2: goto bump_up; 03015 } 03016 #endif 03017 dval(d) += dval(d); 03018 if (dval(d) > ds || dval(d) == ds && L & 1) { 03019 bump_up: 03020 while(*--s == '9') 03021 if (s == s0) { 03022 k++; 03023 *s = '0'; 03024 break; 03025 } 03026 ++*s++; 03027 } 03028 break; 03029 } 03030 } 03031 goto ret1; 03032 } 03033 03034 m2 = b2; 03035 m5 = b5; 03036 mhi = mlo = 0; 03037 if (leftright) { 03038 i = 03039 #ifndef Sudden_Underflow 03040 denorm ? be + (Bias + (P-1) - 1 + 1) : 03041 #endif 03042 #ifdef IBM 03043 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); 03044 #else 03045 1 + P - bbits; 03046 #endif 03047 b2 += i; 03048 s2 += i; 03049 mhi = i2b(1); 03050 } 03051 if (m2 > 0 && s2 > 0) { 03052 i = m2 < s2 ? m2 : s2; 03053 b2 -= i; 03054 m2 -= i; 03055 s2 -= i; 03056 } 03057 if (b5 > 0) { 03058 if (leftright) { 03059 if (m5 > 0) { 03060 mhi = pow5mult(mhi, m5); 03061 b1 = mult(mhi, b); 03062 Bfree(b); 03063 b = b1; 03064 } 03065 if ((j = b5 - m5)) 03066 b = pow5mult(b, j); 03067 } 03068 else 03069 b = pow5mult(b, b5); 03070 } 03071 S = i2b(1); 03072 if (s5 > 0) 03073 S = pow5mult(S, s5); 03074 03075 /* Check for special case that d is a normalized power of 2. */ 03076 03077 spec_case = 0; 03078 if ((mode < 2 || leftright) 03079 #ifdef Honor_FLT_ROUNDS 03080 && rounding == 1 03081 #endif 03082 ) { 03083 if (!word1(d) && !(word0(d) & Bndry_mask) 03084 #ifndef Sudden_Underflow 03085 && word0(d) & (Exp_mask & ~Exp_msk1) 03086 #endif 03087 ) { 03088 /* The special case */ 03089 b2 += Log2P; 03090 s2 += Log2P; 03091 spec_case = 1; 03092 } 03093 } 03094 03095 /* Arrange for convenient computation of quotients: 03096 * shift left if necessary so divisor has 4 leading 0 bits. 03097 * 03098 * Perhaps we should just compute leading 28 bits of S once 03099 * and for all and pass them and a shift to quorem, so it 03100 * can do shifts and ors to compute the numerator for q. 03101 */ 03102 #ifdef Pack_32 03103 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)) 03104 i = 32 - i; 03105 #else 03106 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) 03107 i = 16 - i; 03108 #endif 03109 if (i > 4) { 03110 i -= 4; 03111 b2 += i; 03112 m2 += i; 03113 s2 += i; 03114 } 03115 else if (i < 4) { 03116 i += 28; 03117 b2 += i; 03118 m2 += i; 03119 s2 += i; 03120 } 03121 if (b2 > 0) 03122 b = lshift(b, b2); 03123 if (s2 > 0) 03124 S = lshift(S, s2); 03125 if (k_check) { 03126 if (cmp(b,S) < 0) { 03127 k--; 03128 b = multadd(b, 10, 0); /* we botched the k estimate */ 03129 if (leftright) 03130 mhi = multadd(mhi, 10, 0); 03131 ilim = ilim1; 03132 } 03133 } 03134 if (ilim <= 0 && (mode == 3 || mode == 5)) { 03135 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 03136 /* no digits, fcvt style */ 03137 no_digits: 03138 k = -1 - ndigits; 03139 goto ret; 03140 } 03141 one_digit: 03142 *s++ = '1'; 03143 k++; 03144 goto ret; 03145 } 03146 if (leftright) { 03147 if (m2 > 0) 03148 mhi = lshift(mhi, m2); 03149 03150 /* Compute mlo -- check for special case 03151 * that d is a normalized power of 2. 03152 */ 03153 03154 mlo = mhi; 03155 if (spec_case) { 03156 mhi = Balloc(mhi->k); 03157 Bcopy(mhi, mlo); 03158 mhi = lshift(mhi, Log2P); 03159 } 03160 03161 for(i = 1;;i++) { 03162 dig = quorem(b,S) + '0'; 03163 /* Do we yet have the shortest decimal string 03164 * that will round to d? 03165 */ 03166 j = cmp(b, mlo); 03167 delta = diff(S, mhi); 03168 j1 = delta->sign ? 1 : cmp(b, delta); 03169 Bfree(delta); 03170 #ifndef ROUND_BIASED 03171 if (j1 == 0 && mode != 1 && !(word1(d) & 1) 03172 #ifdef Honor_FLT_ROUNDS 03173 && rounding >= 1 03174 #endif 03175 ) { 03176 if (dig == '9') 03177 goto round_9_up; 03178 if (j > 0) 03179 dig++; 03180 #ifdef SET_INEXACT 03181 else if (!b->x[0] && b->wds <= 1) 03182 inexact = 0; 03183 #endif 03184 *s++ = dig; 03185 goto ret; 03186 } 03187 #endif 03188 if (j < 0 || j == 0 && mode != 1 03189 #ifndef ROUND_BIASED 03190 && !(word1(d) & 1) 03191 #endif 03192 ) { 03193 if (!b->x[0] && b->wds <= 1) { 03194 #ifdef SET_INEXACT 03195 inexact = 0; 03196 #endif 03197 goto accept_dig; 03198 } 03199 #ifdef Honor_FLT_ROUNDS 03200 if (mode > 1) 03201 switch(rounding) { 03202 case 0: goto accept_dig; 03203 case 2: goto keep_dig; 03204 } 03205 #endif /*Honor_FLT_ROUNDS*/ 03206 if (j1 > 0) { 03207 b = lshift(b, 1); 03208 j1 = cmp(b, S); 03209 if ((j1 > 0 || j1 == 0 && dig & 1) 03210 && dig++ == '9') 03211 goto round_9_up; 03212 } 03213 accept_dig: 03214 *s++ = dig; 03215 goto ret; 03216 } 03217 if (j1 > 0) { 03218 #ifdef Honor_FLT_ROUNDS 03219 if (!rounding) 03220 goto accept_dig; 03221 #endif 03222 if (dig == '9') { /* possible if i == 1 */ 03223 round_9_up: 03224 *s++ = '9'; 03225 goto roundoff; 03226 } 03227 *s++ = dig + 1; 03228 goto ret; 03229 } 03230 #ifdef Honor_FLT_ROUNDS 03231 keep_dig: 03232 #endif 03233 *s++ = dig; 03234 if (i == ilim) 03235 break; 03236 b = multadd(b, 10, 0); 03237 if (mlo == mhi) 03238 mlo = mhi = multadd(mhi, 10, 0); 03239 else { 03240 mlo = multadd(mlo, 10, 0); 03241 mhi = multadd(mhi, 10, 0); 03242 } 03243 } 03244 } 03245 else 03246 for(i = 1;; i++) { 03247 *s++ = dig = quorem(b,S) + '0'; 03248 if (!b->x[0] && b->wds <= 1) { 03249 #ifdef SET_INEXACT 03250 inexact = 0; 03251 #endif 03252 goto ret; 03253 } 03254 if (i >= ilim) 03255 break; 03256 b = multadd(b, 10, 0); 03257 } 03258 03259 /* Round off last digit */ 03260 03261 #ifdef Honor_FLT_ROUNDS 03262 switch(rounding) { 03263 case 0: goto trimzeros; 03264 case 2: goto roundoff; 03265 } 03266 #endif 03267 b = lshift(b, 1); 03268 j = cmp(b, S); 03269 if (j > 0 || j == 0 && dig & 1) { 03270 roundoff: 03271 while(*--s == '9') 03272 if (s == s0) { 03273 k++; 03274 *s++ = '1'; 03275 goto ret; 03276 } 03277 ++*s++; 03278 } 03279 else { 03280 #ifdef Honor_FLT_ROUNDS 03281 trimzeros: 03282 #endif 03283 while(*--s == '0'); 03284 s++; 03285 } 03286 ret: 03287 Bfree(S); 03288 if (mhi) { 03289 if (mlo && mlo != mhi) 03290 Bfree(mlo); 03291 Bfree(mhi); 03292 } 03293 ret1: 03294 #ifdef SET_INEXACT 03295 if (inexact) { 03296 if (!oldinexact) { 03297 word0(d) = Exp_1 + (70 << Exp_shift); 03298 word1(d) = 0; 03299 dval(d) += 1.; 03300 } 03301 } 03302 else if (!oldinexact) 03303 clear_inexact(); 03304 #endif 03305 Bfree(b); 03306 *s = 0; 03307 *decpt = k + 1; 03308 if (rve) 03309 *rve = s; 03310 return s0; 03311 } 03312 #ifdef __cplusplus 03313 } 03314 #endif
KDE Logo
This file is part of the documentation for kjs Library Version 3.2.3.
Documentation copyright © 1996-2004 the KDE developers.
Generated on Fri Aug 20 09:48:58 2004 by doxygen 1.3.7 written by Dimitri van Heesch, © 1997-2003