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