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