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