1/********************************************************************** 2 3 util.c - 4 5 $Author: nagachika $ 6 created at: Fri Mar 10 17:22:34 JST 1995 7 8 Copyright (C) 1993-2008 Yukihiro Matsumoto 9 10**********************************************************************/ 11 12#include "ruby/ruby.h" 13#include "internal.h" 14 15#include <ctype.h> 16#include <stdio.h> 17#include <errno.h> 18#include <math.h> 19#include <float.h> 20 21#ifdef _WIN32 22#include "missing/file.h" 23#endif 24 25#include "ruby/util.h" 26 27unsigned long 28ruby_scan_oct(const char *start, size_t len, size_t *retlen) 29{ 30 register const char *s = start; 31 register unsigned long retval = 0; 32 33 while (len-- && *s >= '0' && *s <= '7') { 34 retval <<= 3; 35 retval |= *s++ - '0'; 36 } 37 *retlen = (int)(s - start); /* less than len */ 38 return retval; 39} 40 41unsigned long 42ruby_scan_hex(const char *start, size_t len, size_t *retlen) 43{ 44 static const char hexdigit[] = "0123456789abcdef0123456789ABCDEF"; 45 register const char *s = start; 46 register unsigned long retval = 0; 47 const char *tmp; 48 49 while (len-- && *s && (tmp = strchr(hexdigit, *s))) { 50 retval <<= 4; 51 retval |= (tmp - hexdigit) & 15; 52 s++; 53 } 54 *retlen = (int)(s - start); /* less than len */ 55 return retval; 56} 57 58static unsigned long 59scan_digits(const char *str, int base, size_t *retlen, int *overflow) 60{ 61 static signed char table[] = { 62 /* 0 1 2 3 4 5 6 7 8 9 a b c d e f */ 63 /*0*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 64 /*1*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 65 /*2*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 66 /*3*/ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,-1,-1,-1,-1,-1,-1, 67 /*4*/ -1,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24, 68 /*5*/ 25,26,27,28,29,30,31,32,33,34,35,-1,-1,-1,-1,-1, 69 /*6*/ -1,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24, 70 /*7*/ 25,26,27,28,29,30,31,32,33,34,35,-1,-1,-1,-1,-1, 71 /*8*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 72 /*9*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 73 /*a*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 74 /*b*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 75 /*c*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 76 /*d*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 77 /*e*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 78 /*f*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 79 }; 80 81 const char *start = str; 82 unsigned long ret = 0, x; 83 unsigned long mul_overflow = (~(unsigned long)0) / base; 84 int c; 85 *overflow = 0; 86 87 while ((c = (unsigned char)*str++) != '\0') { 88 int d = table[c]; 89 if (d == -1 || base <= d) { 90 *retlen = (str-1) - start; 91 return ret; 92 } 93 if (mul_overflow < ret) 94 *overflow = 1; 95 ret *= base; 96 x = ret; 97 ret += d; 98 if (ret < x) 99 *overflow = 1; 100 } 101 *retlen = (str-1) - start; 102 return ret; 103} 104 105unsigned long 106ruby_strtoul(const char *str, char **endptr, int base) 107{ 108 int c, b, overflow; 109 int sign = 0; 110 size_t len; 111 unsigned long ret; 112 const char *subject_found = str; 113 114 if (base == 1 || 36 < base) { 115 errno = EINVAL; 116 return 0; 117 } 118 119 while ((c = *str) && ISSPACE(c)) 120 str++; 121 122 if (c == '+') { 123 sign = 1; 124 str++; 125 } 126 else if (c == '-') { 127 sign = -1; 128 str++; 129 } 130 131 if (str[0] == '0') { 132 subject_found = str+1; 133 if (base == 0 || base == 16) { 134 if (str[1] == 'x' || str[1] == 'X') { 135 b = 16; 136 str += 2; 137 } 138 else { 139 b = base == 0 ? 8 : 16; 140 str++; 141 } 142 } 143 else { 144 b = base; 145 str++; 146 } 147 } 148 else { 149 b = base == 0 ? 10 : base; 150 } 151 152 ret = scan_digits(str, b, &len, &overflow); 153 154 if (0 < len) 155 subject_found = str+len; 156 157 if (endptr) 158 *endptr = (char*)subject_found; 159 160 if (overflow) { 161 errno = ERANGE; 162 return ULONG_MAX; 163 } 164 165 if (sign < 0) { 166 ret = (unsigned long)(-(long)ret); 167 return ret; 168 } 169 else { 170 return ret; 171 } 172} 173 174#include <sys/types.h> 175#include <sys/stat.h> 176#ifdef HAVE_UNISTD_H 177#include <unistd.h> 178#endif 179#if defined(HAVE_FCNTL_H) 180#include <fcntl.h> 181#endif 182 183#ifndef S_ISDIR 184# define S_ISDIR(m) (((m) & S_IFMT) == S_IFDIR) 185#endif 186 187 188/* mm.c */ 189 190#define mmtype long 191#define mmcount (16 / SIZEOF_LONG) 192#define A ((mmtype*)a) 193#define B ((mmtype*)b) 194#define C ((mmtype*)c) 195#define D ((mmtype*)d) 196 197#define mmstep (sizeof(mmtype) * mmcount) 198#define mmprepare(base, size) do {\ 199 if (((VALUE)(base) % sizeof(mmtype)) == 0 && ((size) % sizeof(mmtype)) == 0) \ 200 if ((size) >= mmstep) mmkind = 1;\ 201 else mmkind = 0;\ 202 else mmkind = -1;\ 203 high = ((size) / mmstep) * mmstep;\ 204 low = ((size) % mmstep);\ 205} while (0)\ 206 207#define mmarg mmkind, size, high, low 208#define mmargdecl int mmkind, size_t size, size_t high, size_t low 209 210static void mmswap_(register char *a, register char *b, mmargdecl) 211{ 212 if (a == b) return; 213 if (mmkind >= 0) { 214 register mmtype s; 215#if mmcount > 1 216 if (mmkind > 0) { 217 register char *t = a + high; 218 do { 219 s = A[0]; A[0] = B[0]; B[0] = s; 220 s = A[1]; A[1] = B[1]; B[1] = s; 221#if mmcount > 2 222 s = A[2]; A[2] = B[2]; B[2] = s; 223#if mmcount > 3 224 s = A[3]; A[3] = B[3]; B[3] = s; 225#endif 226#endif 227 a += mmstep; b += mmstep; 228 } while (a < t); 229 } 230#endif 231 if (low != 0) { s = A[0]; A[0] = B[0]; B[0] = s; 232#if mmcount > 2 233 if (low >= 2 * sizeof(mmtype)) { s = A[1]; A[1] = B[1]; B[1] = s; 234#if mmcount > 3 235 if (low >= 3 * sizeof(mmtype)) {s = A[2]; A[2] = B[2]; B[2] = s;} 236#endif 237 } 238#endif 239 } 240 } 241 else { 242 register char *t = a + size, s; 243 do {s = *a; *a++ = *b; *b++ = s;} while (a < t); 244 } 245} 246#define mmswap(a,b) mmswap_((a),(b),mmarg) 247 248/* a, b, c = b, c, a */ 249static void mmrot3_(register char *a, register char *b, register char *c, mmargdecl) 250{ 251 if (mmkind >= 0) { 252 register mmtype s; 253#if mmcount > 1 254 if (mmkind > 0) { 255 register char *t = a + high; 256 do { 257 s = A[0]; A[0] = B[0]; B[0] = C[0]; C[0] = s; 258 s = A[1]; A[1] = B[1]; B[1] = C[1]; C[1] = s; 259#if mmcount > 2 260 s = A[2]; A[2] = B[2]; B[2] = C[2]; C[2] = s; 261#if mmcount > 3 262 s = A[3]; A[3] = B[3]; B[3] = C[3]; C[3] = s; 263#endif 264#endif 265 a += mmstep; b += mmstep; c += mmstep; 266 } while (a < t); 267 } 268#endif 269 if (low != 0) { s = A[0]; A[0] = B[0]; B[0] = C[0]; C[0] = s; 270#if mmcount > 2 271 if (low >= 2 * sizeof(mmtype)) { s = A[1]; A[1] = B[1]; B[1] = C[1]; C[1] = s; 272#if mmcount > 3 273 if (low == 3 * sizeof(mmtype)) {s = A[2]; A[2] = B[2]; B[2] = C[2]; C[2] = s;} 274#endif 275 } 276#endif 277 } 278 } 279 else { 280 register char *t = a + size, s; 281 do {s = *a; *a++ = *b; *b++ = *c; *c++ = s;} while (a < t); 282 } 283} 284#define mmrot3(a,b,c) mmrot3_((a),(b),(c),mmarg) 285 286/* qs6.c */ 287/*****************************************************/ 288/* */ 289/* qs6 (Quick sort function) */ 290/* */ 291/* by Tomoyuki Kawamura 1995.4.21 */ 292/* kawamura@tokuyama.ac.jp */ 293/*****************************************************/ 294 295typedef struct { char *LL, *RR; } stack_node; /* Stack structure for L,l,R,r */ 296#define PUSH(ll,rr) do { top->LL = (ll); top->RR = (rr); ++top; } while (0) /* Push L,l,R,r */ 297#define POP(ll,rr) do { --top; (ll) = top->LL; (rr) = top->RR; } while (0) /* Pop L,l,R,r */ 298 299#define med3(a,b,c) ((*cmp)((a),(b),d)<0 ? \ 300 ((*cmp)((b),(c),d)<0 ? (b) : ((*cmp)((a),(c),d)<0 ? (c) : (a))) : \ 301 ((*cmp)((b),(c),d)>0 ? (b) : ((*cmp)((a),(c),d)<0 ? (a) : (c)))) 302 303typedef int (cmpfunc_t)(const void*, const void*, void*); 304void 305ruby_qsort(void* base, const size_t nel, const size_t size, cmpfunc_t *cmp, void *d) 306{ 307 register char *l, *r, *m; /* l,r:left,right group m:median point */ 308 register int t, eq_l, eq_r; /* eq_l: all items in left group are equal to S */ 309 char *L = base; /* left end of current region */ 310 char *R = (char*)base + size*(nel-1); /* right end of current region */ 311 size_t chklim = 63; /* threshold of ordering element check */ 312 enum {size_bits = sizeof(size) * CHAR_BIT}; 313 stack_node stack[size_bits]; /* enough for size_t size */ 314 stack_node *top = stack; 315 int mmkind; 316 size_t high, low, n; 317 318 if (nel <= 1) return; /* need not to sort */ 319 mmprepare(base, size); 320 goto start; 321 322 nxt: 323 if (stack == top) return; /* return if stack is empty */ 324 POP(L,R); 325 326 for (;;) { 327 start: 328 if (L + size == R) { /* 2 elements */ 329 if ((*cmp)(L,R,d) > 0) mmswap(L,R); goto nxt; 330 } 331 332 l = L; r = R; 333 n = (r - l + size) / size; /* number of elements */ 334 m = l + size * (n >> 1); /* calculate median value */ 335 336 if (n >= 60) { 337 register char *m1; 338 register char *m3; 339 if (n >= 200) { 340 n = size*(n>>3); /* number of bytes in splitting 8 */ 341 { 342 register char *p1 = l + n; 343 register char *p2 = p1 + n; 344 register char *p3 = p2 + n; 345 m1 = med3(p1, p2, p3); 346 p1 = m + n; 347 p2 = p1 + n; 348 p3 = p2 + n; 349 m3 = med3(p1, p2, p3); 350 } 351 } 352 else { 353 n = size*(n>>2); /* number of bytes in splitting 4 */ 354 m1 = l + n; 355 m3 = m + n; 356 } 357 m = med3(m1, m, m3); 358 } 359 360 if ((t = (*cmp)(l,m,d)) < 0) { /*3-5-?*/ 361 if ((t = (*cmp)(m,r,d)) < 0) { /*3-5-7*/ 362 if (chklim && nel >= chklim) { /* check if already ascending order */ 363 char *p; 364 chklim = 0; 365 for (p=l; p<r; p+=size) if ((*cmp)(p,p+size,d) > 0) goto fail; 366 goto nxt; 367 } 368 fail: goto loopA; /*3-5-7*/ 369 } 370 if (t > 0) { 371 if ((*cmp)(l,r,d) <= 0) {mmswap(m,r); goto loopA;} /*3-5-4*/ 372 mmrot3(r,m,l); goto loopA; /*3-5-2*/ 373 } 374 goto loopB; /*3-5-5*/ 375 } 376 377 if (t > 0) { /*7-5-?*/ 378 if ((t = (*cmp)(m,r,d)) > 0) { /*7-5-3*/ 379 if (chklim && nel >= chklim) { /* check if already ascending order */ 380 char *p; 381 chklim = 0; 382 for (p=l; p<r; p+=size) if ((*cmp)(p,p+size,d) < 0) goto fail2; 383 while (l<r) {mmswap(l,r); l+=size; r-=size;} /* reverse region */ 384 goto nxt; 385 } 386 fail2: mmswap(l,r); goto loopA; /*7-5-3*/ 387 } 388 if (t < 0) { 389 if ((*cmp)(l,r,d) <= 0) {mmswap(l,m); goto loopB;} /*7-5-8*/ 390 mmrot3(l,m,r); goto loopA; /*7-5-6*/ 391 } 392 mmswap(l,r); goto loopA; /*7-5-5*/ 393 } 394 395 if ((t = (*cmp)(m,r,d)) < 0) {goto loopA;} /*5-5-7*/ 396 if (t > 0) {mmswap(l,r); goto loopB;} /*5-5-3*/ 397 398 /* determining splitting type in case 5-5-5 */ /*5-5-5*/ 399 for (;;) { 400 if ((l += size) == r) goto nxt; /*5-5-5*/ 401 if (l == m) continue; 402 if ((t = (*cmp)(l,m,d)) > 0) {mmswap(l,r); l = L; goto loopA;}/*575-5*/ 403 if (t < 0) {mmswap(L,l); l = L; goto loopB;} /*535-5*/ 404 } 405 406 loopA: eq_l = 1; eq_r = 1; /* splitting type A */ /* left <= median < right */ 407 for (;;) { 408 for (;;) { 409 if ((l += size) == r) 410 {l -= size; if (l != m) mmswap(m,l); l -= size; goto fin;} 411 if (l == m) continue; 412 if ((t = (*cmp)(l,m,d)) > 0) {eq_r = 0; break;} 413 if (t < 0) eq_l = 0; 414 } 415 for (;;) { 416 if (l == (r -= size)) 417 {l -= size; if (l != m) mmswap(m,l); l -= size; goto fin;} 418 if (r == m) {m = l; break;} 419 if ((t = (*cmp)(r,m,d)) < 0) {eq_l = 0; break;} 420 if (t == 0) break; 421 } 422 mmswap(l,r); /* swap left and right */ 423 } 424 425 loopB: eq_l = 1; eq_r = 1; /* splitting type B */ /* left < median <= right */ 426 for (;;) { 427 for (;;) { 428 if (l == (r -= size)) 429 {r += size; if (r != m) mmswap(r,m); r += size; goto fin;} 430 if (r == m) continue; 431 if ((t = (*cmp)(r,m,d)) < 0) {eq_l = 0; break;} 432 if (t > 0) eq_r = 0; 433 } 434 for (;;) { 435 if ((l += size) == r) 436 {r += size; if (r != m) mmswap(r,m); r += size; goto fin;} 437 if (l == m) {m = r; break;} 438 if ((t = (*cmp)(l,m,d)) > 0) {eq_r = 0; break;} 439 if (t == 0) break; 440 } 441 mmswap(l,r); /* swap left and right */ 442 } 443 444 fin: 445 if (eq_l == 0) /* need to sort left side */ 446 if (eq_r == 0) /* need to sort right side */ 447 if (l-L < R-r) {PUSH(r,R); R = l;} /* sort left side first */ 448 else {PUSH(L,l); L = r;} /* sort right side first */ 449 else R = l; /* need to sort left side only */ 450 else if (eq_r == 0) L = r; /* need to sort right side only */ 451 else goto nxt; /* need not to sort both sides */ 452 } 453} 454 455char * 456ruby_strdup(const char *str) 457{ 458 char *tmp; 459 size_t len = strlen(str) + 1; 460 461 tmp = xmalloc(len); 462 memcpy(tmp, str, len); 463 464 return tmp; 465} 466 467#ifdef __native_client__ 468char * 469ruby_getcwd(void) 470{ 471 char *buf = xmalloc(2); 472 strcpy(buf, "."); 473 return buf; 474} 475#else 476char * 477ruby_getcwd(void) 478{ 479#ifdef HAVE_GETCWD 480 int size = 200; 481 char *buf = xmalloc(size); 482 483 while (!getcwd(buf, size)) { 484 if (errno != ERANGE) { 485 xfree(buf); 486 rb_sys_fail("getcwd"); 487 } 488 size *= 2; 489 buf = xrealloc(buf, size); 490 } 491#else 492# ifndef PATH_MAX 493# define PATH_MAX 8192 494# endif 495 char *buf = xmalloc(PATH_MAX+1); 496 497 if (!getwd(buf)) { 498 xfree(buf); 499 rb_sys_fail("getwd"); 500 } 501#endif 502 return buf; 503} 504#endif 505 506/**************************************************************** 507 * 508 * The author of this software is David M. Gay. 509 * 510 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. 511 * 512 * Permission to use, copy, modify, and distribute this software for any 513 * purpose without fee is hereby granted, provided that this entire notice 514 * is included in all copies of any software which is or includes a copy 515 * or modification of this software and in all copies of the supporting 516 * documentation for such software. 517 * 518 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 519 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY 520 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 521 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 522 * 523 ***************************************************************/ 524 525/* Please send bug reports to David M. Gay (dmg at acm dot org, 526 * with " at " changed at "@" and " dot " changed to "."). */ 527 528/* On a machine with IEEE extended-precision registers, it is 529 * necessary to specify double-precision (53-bit) rounding precision 530 * before invoking strtod or dtoa. If the machine uses (the equivalent 531 * of) Intel 80x87 arithmetic, the call 532 * _control87(PC_53, MCW_PC); 533 * does this with many compilers. Whether this or another call is 534 * appropriate depends on the compiler; for this to work, it may be 535 * necessary to #include "float.h" or another system-dependent header 536 * file. 537 */ 538 539/* strtod for IEEE-, VAX-, and IBM-arithmetic machines. 540 * 541 * This strtod returns a nearest machine number to the input decimal 542 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 543 * broken by the IEEE round-even rule. Otherwise ties are broken by 544 * biased rounding (add half and chop). 545 * 546 * Inspired loosely by William D. Clinger's paper "How to Read Floating 547 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 548 * 549 * Modifications: 550 * 551 * 1. We only require IEEE, IBM, or VAX double-precision 552 * arithmetic (not IEEE double-extended). 553 * 2. We get by with floating-point arithmetic in a case that 554 * Clinger missed -- when we're computing d * 10^n 555 * for a small integer d and the integer n is not too 556 * much larger than 22 (the maximum integer k for which 557 * we can represent 10^k exactly), we may be able to 558 * compute (d*10^k) * 10^(e-k) with just one roundoff. 559 * 3. Rather than a bit-at-a-time adjustment of the binary 560 * result in the hard case, we use floating-point 561 * arithmetic to determine the adjustment to within 562 * one bit; only in really hard cases do we need to 563 * compute a second residual. 564 * 4. Because of 3., we don't need a large table of powers of 10 565 * for ten-to-e (just some small tables, e.g. of 10^k 566 * for 0 <= k <= 22). 567 */ 568 569/* 570 * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least 571 * significant byte has the lowest address. 572 * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most 573 * significant byte has the lowest address. 574 * #define Long int on machines with 32-bit ints and 64-bit longs. 575 * #define IBM for IBM mainframe-style floating-point arithmetic. 576 * #define VAX for VAX-style floating-point arithmetic (D_floating). 577 * #define No_leftright to omit left-right logic in fast floating-point 578 * computation of dtoa. 579 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 580 * and strtod and dtoa should round accordingly. 581 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 582 * and Honor_FLT_ROUNDS is not #defined. 583 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 584 * that use extended-precision instructions to compute rounded 585 * products and quotients) with IBM. 586 * #define ROUND_BIASED for IEEE-format with biased rounding. 587 * #define Inaccurate_Divide for IEEE-format with correctly rounded 588 * products but inaccurate quotients, e.g., for Intel i860. 589 * #define NO_LONG_LONG on machines that do not have a "long long" 590 * integer type (of >= 64 bits). On such machines, you can 591 * #define Just_16 to store 16 bits per 32-bit Long when doing 592 * high-precision integer arithmetic. Whether this speeds things 593 * up or slows things down depends on the machine and the number 594 * being converted. If long long is available and the name is 595 * something other than "long long", #define Llong to be the name, 596 * and if "unsigned Llong" does not work as an unsigned version of 597 * Llong, #define #ULLong to be the corresponding unsigned type. 598 * #define KR_headers for old-style C function headers. 599 * #define Bad_float_h if your system lacks a float.h or if it does not 600 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, 601 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. 602 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) 603 * if memory is available and otherwise does something you deem 604 * appropriate. If MALLOC is undefined, malloc will be invoked 605 * directly -- and assumed always to succeed. 606 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making 607 * memory allocations from a private pool of memory when possible. 608 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes, 609 * unless #defined to be a different length. This default length 610 * suffices to get rid of MALLOC calls except for unusual cases, 611 * such as decimal-to-binary conversion of a very long string of 612 * digits. The longest string dtoa can return is about 751 bytes 613 * long. For conversions by strtod of strings of 800 digits and 614 * all dtoa conversions in single-threaded executions with 8-byte 615 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte 616 * pointers, PRIVATE_MEM >= 7112 appears adequate. 617 * #define INFNAN_CHECK on IEEE systems to cause strtod to check for 618 * Infinity and NaN (case insensitively). On some systems (e.g., 619 * some HP systems), it may be necessary to #define NAN_WORD0 620 * appropriately -- to the most significant word of a quiet NaN. 621 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.) 622 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined, 623 * strtod also accepts (case insensitively) strings of the form 624 * NaN(x), where x is a string of hexadecimal digits and spaces; 625 * if there is only one string of hexadecimal digits, it is taken 626 * for the 52 fraction bits of the resulting NaN; if there are two 627 * or more strings of hex digits, the first is for the high 20 bits, 628 * the second and subsequent for the low 32 bits, with intervening 629 * white space ignored; but if this results in none of the 52 630 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0 631 * and NAN_WORD1 are used instead. 632 * #define MULTIPLE_THREADS if the system offers preemptively scheduled 633 * multiple threads. In this case, you must provide (or suitably 634 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed 635 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed 636 * in pow5mult, ensures lazy evaluation of only one copy of high 637 * powers of 5; omitting this lock would introduce a small 638 * probability of wasting memory, but would otherwise be harmless.) 639 * You must also invoke freedtoa(s) to free the value s returned by 640 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. 641 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that 642 * avoids underflows on inputs whose result does not underflow. 643 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format 644 * floating-point numbers and flushes underflows to zero rather 645 * than implementing gradual underflow, then you must also #define 646 * Sudden_Underflow. 647 * #define YES_ALIAS to permit aliasing certain double values with 648 * arrays of ULongs. This leads to slightly better code with 649 * some compilers and was always used prior to 19990916, but it 650 * is not strictly legal and can cause trouble with aggressively 651 * optimizing compilers (e.g., gcc 2.95.1 under -O2). 652 * #define USE_LOCALE to use the current locale's decimal_point value. 653 * #define SET_INEXACT if IEEE arithmetic is being used and extra 654 * computation should be done to set the inexact flag when the 655 * result is inexact and avoid setting inexact when the result 656 * is exact. In this case, dtoa.c must be compiled in 657 * an environment, perhaps provided by #include "dtoa.c" in a 658 * suitable wrapper, that defines two functions, 659 * int get_inexact(void); 660 * void clear_inexact(void); 661 * such that get_inexact() returns a nonzero value if the 662 * inexact bit is already set, and clear_inexact() sets the 663 * inexact bit to 0. When SET_INEXACT is #defined, strtod 664 * also does extra computations to set the underflow and overflow 665 * flags when appropriate (i.e., when the result is tiny and 666 * inexact or when it is a numeric value rounded to +-infinity). 667 * #define NO_ERRNO if strtod should not assign errno = ERANGE when 668 * the result overflows to +-Infinity or underflows to 0. 669 */ 670 671#ifdef WORDS_BIGENDIAN 672#define IEEE_BIG_ENDIAN 673#else 674#define IEEE_LITTLE_ENDIAN 675#endif 676 677#ifdef __vax__ 678#define VAX 679#undef IEEE_BIG_ENDIAN 680#undef IEEE_LITTLE_ENDIAN 681#endif 682 683#if defined(__arm__) && !defined(__VFP_FP__) 684#define IEEE_BIG_ENDIAN 685#undef IEEE_LITTLE_ENDIAN 686#endif 687 688#undef Long 689#undef ULong 690 691#if SIZEOF_INT == 4 692#define Long int 693#define ULong unsigned int 694#elif SIZEOF_LONG == 4 695#define Long long int 696#define ULong unsigned long int 697#endif 698 699#if HAVE_LONG_LONG 700#define Llong LONG_LONG 701#endif 702 703#ifdef DEBUG 704#include "stdio.h" 705#define Bug(x) {fprintf(stderr, "%s\n", (x)); exit(EXIT_FAILURE);} 706#endif 707 708#include "stdlib.h" 709#include "string.h" 710 711#ifdef USE_LOCALE 712#include "locale.h" 713#endif 714 715#ifdef MALLOC 716extern void *MALLOC(size_t); 717#else 718#define MALLOC malloc 719#endif 720#ifdef FREE 721extern void FREE(void*); 722#else 723#define FREE free 724#endif 725 726#ifndef Omit_Private_Memory 727#ifndef PRIVATE_MEM 728#define PRIVATE_MEM 2304 729#endif 730#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) 731static double private_mem[PRIVATE_mem], *pmem_next = private_mem; 732#endif 733 734#undef IEEE_Arith 735#undef Avoid_Underflow 736#ifdef IEEE_BIG_ENDIAN 737#define IEEE_Arith 738#endif 739#ifdef IEEE_LITTLE_ENDIAN 740#define IEEE_Arith 741#endif 742 743#ifdef Bad_float_h 744 745#ifdef IEEE_Arith 746#define DBL_DIG 15 747#define DBL_MAX_10_EXP 308 748#define DBL_MAX_EXP 1024 749#define FLT_RADIX 2 750#endif /*IEEE_Arith*/ 751 752#ifdef IBM 753#define DBL_DIG 16 754#define DBL_MAX_10_EXP 75 755#define DBL_MAX_EXP 63 756#define FLT_RADIX 16 757#define DBL_MAX 7.2370055773322621e+75 758#endif 759 760#ifdef VAX 761#define DBL_DIG 16 762#define DBL_MAX_10_EXP 38 763#define DBL_MAX_EXP 127 764#define FLT_RADIX 2 765#define DBL_MAX 1.7014118346046923e+38 766#endif 767 768#ifndef LONG_MAX 769#define LONG_MAX 2147483647 770#endif 771 772#else /* ifndef Bad_float_h */ 773#include "float.h" 774#endif /* Bad_float_h */ 775 776#ifndef __MATH_H__ 777#include "math.h" 778#endif 779 780#ifdef __cplusplus 781extern "C" { 782#if 0 783} /* satisfy cc-mode */ 784#endif 785#endif 786 787#if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + defined(IBM) != 1 788Exactly one of IEEE_LITTLE_ENDIAN, IEEE_BIG_ENDIAN, VAX, or IBM should be defined. 789#endif 790 791typedef union { double d; ULong L[2]; } U; 792 793#ifdef YES_ALIAS 794typedef double double_u; 795# define dval(x) (x) 796# ifdef IEEE_LITTLE_ENDIAN 797# define word0(x) (((ULong *)&(x))[1]) 798# define word1(x) (((ULong *)&(x))[0]) 799# else 800# define word0(x) (((ULong *)&(x))[0]) 801# define word1(x) (((ULong *)&(x))[1]) 802# endif 803#else 804typedef U double_u; 805# ifdef IEEE_LITTLE_ENDIAN 806# define word0(x) ((x).L[1]) 807# define word1(x) ((x).L[0]) 808# else 809# define word0(x) ((x).L[0]) 810# define word1(x) ((x).L[1]) 811# endif 812# define dval(x) ((x).d) 813#endif 814 815/* The following definition of Storeinc is appropriate for MIPS processors. 816 * An alternative that might be better on some machines is 817 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 818 */ 819#if defined(IEEE_LITTLE_ENDIAN) + defined(VAX) + defined(__arm__) 820#define Storeinc(a,b,c) (((unsigned short *)(a))[1] = (unsigned short)(b), \ 821((unsigned short *)(a))[0] = (unsigned short)(c), (a)++) 822#else 823#define Storeinc(a,b,c) (((unsigned short *)(a))[0] = (unsigned short)(b), \ 824((unsigned short *)(a))[1] = (unsigned short)(c), (a)++) 825#endif 826 827/* #define P DBL_MANT_DIG */ 828/* Ten_pmax = floor(P*log(2)/log(5)) */ 829/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 830/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 831/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 832 833#ifdef IEEE_Arith 834#define Exp_shift 20 835#define Exp_shift1 20 836#define Exp_msk1 0x100000 837#define Exp_msk11 0x100000 838#define Exp_mask 0x7ff00000 839#define P 53 840#define Bias 1023 841#define Emin (-1022) 842#define Exp_1 0x3ff00000 843#define Exp_11 0x3ff00000 844#define Ebits 11 845#define Frac_mask 0xfffff 846#define Frac_mask1 0xfffff 847#define Ten_pmax 22 848#define Bletch 0x10 849#define Bndry_mask 0xfffff 850#define Bndry_mask1 0xfffff 851#define LSB 1 852#define Sign_bit 0x80000000 853#define Log2P 1 854#define Tiny0 0 855#define Tiny1 1 856#define Quick_max 14 857#define Int_max 14 858#ifndef NO_IEEE_Scale 859#define Avoid_Underflow 860#ifdef Flush_Denorm /* debugging option */ 861#undef Sudden_Underflow 862#endif 863#endif 864 865#ifndef Flt_Rounds 866#ifdef FLT_ROUNDS 867#define Flt_Rounds FLT_ROUNDS 868#else 869#define Flt_Rounds 1 870#endif 871#endif /*Flt_Rounds*/ 872 873#ifdef Honor_FLT_ROUNDS 874#define Rounding rounding 875#undef Check_FLT_ROUNDS 876#define Check_FLT_ROUNDS 877#else 878#define Rounding Flt_Rounds 879#endif 880 881#else /* ifndef IEEE_Arith */ 882#undef Check_FLT_ROUNDS 883#undef Honor_FLT_ROUNDS 884#undef SET_INEXACT 885#undef Sudden_Underflow 886#define Sudden_Underflow 887#ifdef IBM 888#undef Flt_Rounds 889#define Flt_Rounds 0 890#define Exp_shift 24 891#define Exp_shift1 24 892#define Exp_msk1 0x1000000 893#define Exp_msk11 0x1000000 894#define Exp_mask 0x7f000000 895#define P 14 896#define Bias 65 897#define Exp_1 0x41000000 898#define Exp_11 0x41000000 899#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 900#define Frac_mask 0xffffff 901#define Frac_mask1 0xffffff 902#define Bletch 4 903#define Ten_pmax 22 904#define Bndry_mask 0xefffff 905#define Bndry_mask1 0xffffff 906#define LSB 1 907#define Sign_bit 0x80000000 908#define Log2P 4 909#define Tiny0 0x100000 910#define Tiny1 0 911#define Quick_max 14 912#define Int_max 15 913#else /* VAX */ 914#undef Flt_Rounds 915#define Flt_Rounds 1 916#define Exp_shift 23 917#define Exp_shift1 7 918#define Exp_msk1 0x80 919#define Exp_msk11 0x800000 920#define Exp_mask 0x7f80 921#define P 56 922#define Bias 129 923#define Exp_1 0x40800000 924#define Exp_11 0x4080 925#define Ebits 8 926#define Frac_mask 0x7fffff 927#define Frac_mask1 0xffff007f 928#define Ten_pmax 24 929#define Bletch 2 930#define Bndry_mask 0xffff007f 931#define Bndry_mask1 0xffff007f 932#define LSB 0x10000 933#define Sign_bit 0x8000 934#define Log2P 1 935#define Tiny0 0x80 936#define Tiny1 0 937#define Quick_max 15 938#define Int_max 15 939#endif /* IBM, VAX */ 940#endif /* IEEE_Arith */ 941 942#ifndef IEEE_Arith 943#define ROUND_BIASED 944#endif 945 946#ifdef RND_PRODQUOT 947#define rounded_product(a,b) ((a) = rnd_prod((a), (b))) 948#define rounded_quotient(a,b) ((a) = rnd_quot((a), (b))) 949extern double rnd_prod(double, double), rnd_quot(double, double); 950#else 951#define rounded_product(a,b) ((a) *= (b)) 952#define rounded_quotient(a,b) ((a) /= (b)) 953#endif 954 955#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 956#define Big1 0xffffffff 957 958#ifndef Pack_32 959#define Pack_32 960#endif 961 962#define FFFFFFFF 0xffffffffUL 963 964#ifdef NO_LONG_LONG 965#undef ULLong 966#ifdef Just_16 967#undef Pack_32 968/* When Pack_32 is not defined, we store 16 bits per 32-bit Long. 969 * This makes some inner loops simpler and sometimes saves work 970 * during multiplications, but it often seems to make things slightly 971 * slower. Hence the default is now to store 32 bits per Long. 972 */ 973#endif 974#else /* long long available */ 975#ifndef Llong 976#define Llong long long 977#endif 978#ifndef ULLong 979#define ULLong unsigned Llong 980#endif 981#endif /* NO_LONG_LONG */ 982 983#define MULTIPLE_THREADS 1 984 985#ifndef MULTIPLE_THREADS 986#define ACQUIRE_DTOA_LOCK(n) /*nothing*/ 987#define FREE_DTOA_LOCK(n) /*nothing*/ 988#else 989#define ACQUIRE_DTOA_LOCK(n) /*unused right now*/ 990#define FREE_DTOA_LOCK(n) /*unused right now*/ 991#endif 992 993#define Kmax 15 994 995struct Bigint { 996 struct Bigint *next; 997 int k, maxwds, sign, wds; 998 ULong x[1]; 999}; 1000 1001typedef struct Bigint Bigint; 1002 1003static Bigint *freelist[Kmax+1]; 1004 1005static Bigint * 1006Balloc(int k) 1007{ 1008 int x; 1009 Bigint *rv; 1010#ifndef Omit_Private_Memory 1011 size_t len; 1012#endif 1013 1014 ACQUIRE_DTOA_LOCK(0); 1015 if (k <= Kmax && (rv = freelist[k]) != 0) { 1016 freelist[k] = rv->next; 1017 } 1018 else { 1019 x = 1 << k; 1020#ifdef Omit_Private_Memory 1021 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); 1022#else 1023 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) 1024 /sizeof(double); 1025 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) { 1026 rv = (Bigint*)pmem_next; 1027 pmem_next += len; 1028 } 1029 else 1030 rv = (Bigint*)MALLOC(len*sizeof(double)); 1031#endif 1032 rv->k = k; 1033 rv->maxwds = x; 1034 } 1035 FREE_DTOA_LOCK(0); 1036 rv->sign = rv->wds = 0; 1037 return rv; 1038} 1039 1040static void 1041Bfree(Bigint *v) 1042{ 1043 if (v) { 1044 if (v->k > Kmax) { 1045 FREE(v); 1046 return; 1047 } 1048 ACQUIRE_DTOA_LOCK(0); 1049 v->next = freelist[v->k]; 1050 freelist[v->k] = v; 1051 FREE_DTOA_LOCK(0); 1052 } 1053} 1054 1055#define Bcopy(x,y) memcpy((char *)&(x)->sign, (char *)&(y)->sign, \ 1056(y)->wds*sizeof(Long) + 2*sizeof(int)) 1057 1058static Bigint * 1059multadd(Bigint *b, int m, int a) /* multiply by m and add a */ 1060{ 1061 int i, wds; 1062 ULong *x; 1063#ifdef ULLong 1064 ULLong carry, y; 1065#else 1066 ULong carry, y; 1067#ifdef Pack_32 1068 ULong xi, z; 1069#endif 1070#endif 1071 Bigint *b1; 1072 1073 wds = b->wds; 1074 x = b->x; 1075 i = 0; 1076 carry = a; 1077 do { 1078#ifdef ULLong 1079 y = *x * (ULLong)m + carry; 1080 carry = y >> 32; 1081 *x++ = (ULong)(y & FFFFFFFF); 1082#else 1083#ifdef Pack_32 1084 xi = *x; 1085 y = (xi & 0xffff) * m + carry; 1086 z = (xi >> 16) * m + (y >> 16); 1087 carry = z >> 16; 1088 *x++ = (z << 16) + (y & 0xffff); 1089#else 1090 y = *x * m + carry; 1091 carry = y >> 16; 1092 *x++ = y & 0xffff; 1093#endif 1094#endif 1095 } while (++i < wds); 1096 if (carry) { 1097 if (wds >= b->maxwds) { 1098 b1 = Balloc(b->k+1); 1099 Bcopy(b1, b); 1100 Bfree(b); 1101 b = b1; 1102 } 1103 b->x[wds++] = (ULong)carry; 1104 b->wds = wds; 1105 } 1106 return b; 1107} 1108 1109static Bigint * 1110s2b(const char *s, int nd0, int nd, ULong y9) 1111{ 1112 Bigint *b; 1113 int i, k; 1114 Long x, y; 1115 1116 x = (nd + 8) / 9; 1117 for (k = 0, y = 1; x > y; y <<= 1, k++) ; 1118#ifdef Pack_32 1119 b = Balloc(k); 1120 b->x[0] = y9; 1121 b->wds = 1; 1122#else 1123 b = Balloc(k+1); 1124 b->x[0] = y9 & 0xffff; 1125 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; 1126#endif 1127 1128 i = 9; 1129 if (9 < nd0) { 1130 s += 9; 1131 do { 1132 b = multadd(b, 10, *s++ - '0'); 1133 } while (++i < nd0); 1134 s++; 1135 } 1136 else 1137 s += 10; 1138 for (; i < nd; i++) 1139 b = multadd(b, 10, *s++ - '0'); 1140 return b; 1141} 1142 1143static int 1144hi0bits(register ULong x) 1145{ 1146 register int k = 0; 1147 1148 if (!(x & 0xffff0000)) { 1149 k = 16; 1150 x <<= 16; 1151 } 1152 if (!(x & 0xff000000)) { 1153 k += 8; 1154 x <<= 8; 1155 } 1156 if (!(x & 0xf0000000)) { 1157 k += 4; 1158 x <<= 4; 1159 } 1160 if (!(x & 0xc0000000)) { 1161 k += 2; 1162 x <<= 2; 1163 } 1164 if (!(x & 0x80000000)) { 1165 k++; 1166 if (!(x & 0x40000000)) 1167 return 32; 1168 } 1169 return k; 1170} 1171 1172static int 1173lo0bits(ULong *y) 1174{ 1175 register int k; 1176 register ULong x = *y; 1177 1178 if (x & 7) { 1179 if (x & 1) 1180 return 0; 1181 if (x & 2) { 1182 *y = x >> 1; 1183 return 1; 1184 } 1185 *y = x >> 2; 1186 return 2; 1187 } 1188 k = 0; 1189 if (!(x & 0xffff)) { 1190 k = 16; 1191 x >>= 16; 1192 } 1193 if (!(x & 0xff)) { 1194 k += 8; 1195 x >>= 8; 1196 } 1197 if (!(x & 0xf)) { 1198 k += 4; 1199 x >>= 4; 1200 } 1201 if (!(x & 0x3)) { 1202 k += 2; 1203 x >>= 2; 1204 } 1205 if (!(x & 1)) { 1206 k++; 1207 x >>= 1; 1208 if (!x) 1209 return 32; 1210 } 1211 *y = x; 1212 return k; 1213} 1214 1215static Bigint * 1216i2b(int i) 1217{ 1218 Bigint *b; 1219 1220 b = Balloc(1); 1221 b->x[0] = i; 1222 b->wds = 1; 1223 return b; 1224} 1225 1226static Bigint * 1227mult(Bigint *a, Bigint *b) 1228{ 1229 Bigint *c; 1230 int k, wa, wb, wc; 1231 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 1232 ULong y; 1233#ifdef ULLong 1234 ULLong carry, z; 1235#else 1236 ULong carry, z; 1237#ifdef Pack_32 1238 ULong z2; 1239#endif 1240#endif 1241 1242 if (a->wds < b->wds) { 1243 c = a; 1244 a = b; 1245 b = c; 1246 } 1247 k = a->k; 1248 wa = a->wds; 1249 wb = b->wds; 1250 wc = wa + wb; 1251 if (wc > a->maxwds) 1252 k++; 1253 c = Balloc(k); 1254 for (x = c->x, xa = x + wc; x < xa; x++) 1255 *x = 0; 1256 xa = a->x; 1257 xae = xa + wa; 1258 xb = b->x; 1259 xbe = xb + wb; 1260 xc0 = c->x; 1261#ifdef ULLong 1262 for (; xb < xbe; xc0++) { 1263 if ((y = *xb++) != 0) { 1264 x = xa; 1265 xc = xc0; 1266 carry = 0; 1267 do { 1268 z = *x++ * (ULLong)y + *xc + carry; 1269 carry = z >> 32; 1270 *xc++ = (ULong)(z & FFFFFFFF); 1271 } while (x < xae); 1272 *xc = (ULong)carry; 1273 } 1274 } 1275#else 1276#ifdef Pack_32 1277 for (; xb < xbe; xb++, xc0++) { 1278 if (y = *xb & 0xffff) { 1279 x = xa; 1280 xc = xc0; 1281 carry = 0; 1282 do { 1283 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 1284 carry = z >> 16; 1285 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 1286 carry = z2 >> 16; 1287 Storeinc(xc, z2, z); 1288 } while (x < xae); 1289 *xc = (ULong)carry; 1290 } 1291 if (y = *xb >> 16) { 1292 x = xa; 1293 xc = xc0; 1294 carry = 0; 1295 z2 = *xc; 1296 do { 1297 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 1298 carry = z >> 16; 1299 Storeinc(xc, z, z2); 1300 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 1301 carry = z2 >> 16; 1302 } while (x < xae); 1303 *xc = z2; 1304 } 1305 } 1306#else 1307 for (; xb < xbe; xc0++) { 1308 if (y = *xb++) { 1309 x = xa; 1310 xc = xc0; 1311 carry = 0; 1312 do { 1313 z = *x++ * y + *xc + carry; 1314 carry = z >> 16; 1315 *xc++ = z & 0xffff; 1316 } while (x < xae); 1317 *xc = (ULong)carry; 1318 } 1319 } 1320#endif 1321#endif 1322 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 1323 c->wds = wc; 1324 return c; 1325} 1326 1327static Bigint *p5s; 1328 1329static Bigint * 1330pow5mult(Bigint *b, int k) 1331{ 1332 Bigint *b1, *p5, *p51; 1333 int i; 1334 static int p05[3] = { 5, 25, 125 }; 1335 1336 if ((i = k & 3) != 0) 1337 b = multadd(b, p05[i-1], 0); 1338 1339 if (!(k >>= 2)) 1340 return b; 1341 if (!(p5 = p5s)) { 1342 /* first time */ 1343#ifdef MULTIPLE_THREADS 1344 ACQUIRE_DTOA_LOCK(1); 1345 if (!(p5 = p5s)) { 1346 p5 = p5s = i2b(625); 1347 p5->next = 0; 1348 } 1349 FREE_DTOA_LOCK(1); 1350#else 1351 p5 = p5s = i2b(625); 1352 p5->next = 0; 1353#endif 1354 } 1355 for (;;) { 1356 if (k & 1) { 1357 b1 = mult(b, p5); 1358 Bfree(b); 1359 b = b1; 1360 } 1361 if (!(k >>= 1)) 1362 break; 1363 if (!(p51 = p5->next)) { 1364#ifdef MULTIPLE_THREADS 1365 ACQUIRE_DTOA_LOCK(1); 1366 if (!(p51 = p5->next)) { 1367 p51 = p5->next = mult(p5,p5); 1368 p51->next = 0; 1369 } 1370 FREE_DTOA_LOCK(1); 1371#else 1372 p51 = p5->next = mult(p5,p5); 1373 p51->next = 0; 1374#endif 1375 } 1376 p5 = p51; 1377 } 1378 return b; 1379} 1380 1381static Bigint * 1382lshift(Bigint *b, int k) 1383{ 1384 int i, k1, n, n1; 1385 Bigint *b1; 1386 ULong *x, *x1, *xe, z; 1387 1388#ifdef Pack_32 1389 n = k >> 5; 1390#else 1391 n = k >> 4; 1392#endif 1393 k1 = b->k; 1394 n1 = n + b->wds + 1; 1395 for (i = b->maxwds; n1 > i; i <<= 1) 1396 k1++; 1397 b1 = Balloc(k1); 1398 x1 = b1->x; 1399 for (i = 0; i < n; i++) 1400 *x1++ = 0; 1401 x = b->x; 1402 xe = x + b->wds; 1403#ifdef Pack_32 1404 if (k &= 0x1f) { 1405 k1 = 32 - k; 1406 z = 0; 1407 do { 1408 *x1++ = *x << k | z; 1409 z = *x++ >> k1; 1410 } while (x < xe); 1411 if ((*x1 = z) != 0) 1412 ++n1; 1413 } 1414#else 1415 if (k &= 0xf) { 1416 k1 = 16 - k; 1417 z = 0; 1418 do { 1419 *x1++ = *x << k & 0xffff | z; 1420 z = *x++ >> k1; 1421 } while (x < xe); 1422 if (*x1 = z) 1423 ++n1; 1424 } 1425#endif 1426 else 1427 do { 1428 *x1++ = *x++; 1429 } while (x < xe); 1430 b1->wds = n1 - 1; 1431 Bfree(b); 1432 return b1; 1433} 1434 1435static int 1436cmp(Bigint *a, Bigint *b) 1437{ 1438 ULong *xa, *xa0, *xb, *xb0; 1439 int i, j; 1440 1441 i = a->wds; 1442 j = b->wds; 1443#ifdef DEBUG 1444 if (i > 1 && !a->x[i-1]) 1445 Bug("cmp called with a->x[a->wds-1] == 0"); 1446 if (j > 1 && !b->x[j-1]) 1447 Bug("cmp called with b->x[b->wds-1] == 0"); 1448#endif 1449 if (i -= j) 1450 return i; 1451 xa0 = a->x; 1452 xa = xa0 + j; 1453 xb0 = b->x; 1454 xb = xb0 + j; 1455 for (;;) { 1456 if (*--xa != *--xb) 1457 return *xa < *xb ? -1 : 1; 1458 if (xa <= xa0) 1459 break; 1460 } 1461 return 0; 1462} 1463 1464static Bigint * 1465diff(Bigint *a, Bigint *b) 1466{ 1467 Bigint *c; 1468 int i, wa, wb; 1469 ULong *xa, *xae, *xb, *xbe, *xc; 1470#ifdef ULLong 1471 ULLong borrow, y; 1472#else 1473 ULong borrow, y; 1474#ifdef Pack_32 1475 ULong z; 1476#endif 1477#endif 1478 1479 i = cmp(a,b); 1480 if (!i) { 1481 c = Balloc(0); 1482 c->wds = 1; 1483 c->x[0] = 0; 1484 return c; 1485 } 1486 if (i < 0) { 1487 c = a; 1488 a = b; 1489 b = c; 1490 i = 1; 1491 } 1492 else 1493 i = 0; 1494 c = Balloc(a->k); 1495 c->sign = i; 1496 wa = a->wds; 1497 xa = a->x; 1498 xae = xa + wa; 1499 wb = b->wds; 1500 xb = b->x; 1501 xbe = xb + wb; 1502 xc = c->x; 1503 borrow = 0; 1504#ifdef ULLong 1505 do { 1506 y = (ULLong)*xa++ - *xb++ - borrow; 1507 borrow = y >> 32 & (ULong)1; 1508 *xc++ = (ULong)(y & FFFFFFFF); 1509 } while (xb < xbe); 1510 while (xa < xae) { 1511 y = *xa++ - borrow; 1512 borrow = y >> 32 & (ULong)1; 1513 *xc++ = (ULong)(y & FFFFFFFF); 1514 } 1515#else 1516#ifdef Pack_32 1517 do { 1518 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 1519 borrow = (y & 0x10000) >> 16; 1520 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 1521 borrow = (z & 0x10000) >> 16; 1522 Storeinc(xc, z, y); 1523 } while (xb < xbe); 1524 while (xa < xae) { 1525 y = (*xa & 0xffff) - borrow; 1526 borrow = (y & 0x10000) >> 16; 1527 z = (*xa++ >> 16) - borrow; 1528 borrow = (z & 0x10000) >> 16; 1529 Storeinc(xc, z, y); 1530 } 1531#else 1532 do { 1533 y = *xa++ - *xb++ - borrow; 1534 borrow = (y & 0x10000) >> 16; 1535 *xc++ = y & 0xffff; 1536 } while (xb < xbe); 1537 while (xa < xae) { 1538 y = *xa++ - borrow; 1539 borrow = (y & 0x10000) >> 16; 1540 *xc++ = y & 0xffff; 1541 } 1542#endif 1543#endif 1544 while (!*--xc) 1545 wa--; 1546 c->wds = wa; 1547 return c; 1548} 1549 1550static double 1551ulp(double x_) 1552{ 1553 register Long L; 1554 double_u x, a; 1555 dval(x) = x_; 1556 1557 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 1558#ifndef Avoid_Underflow 1559#ifndef Sudden_Underflow 1560 if (L > 0) { 1561#endif 1562#endif 1563#ifdef IBM 1564 L |= Exp_msk1 >> 4; 1565#endif 1566 word0(a) = L; 1567 word1(a) = 0; 1568#ifndef Avoid_Underflow 1569#ifndef Sudden_Underflow 1570 } 1571 else { 1572 L = -L >> Exp_shift; 1573 if (L < Exp_shift) { 1574 word0(a) = 0x80000 >> L; 1575 word1(a) = 0; 1576 } 1577 else { 1578 word0(a) = 0; 1579 L -= Exp_shift; 1580 word1(a) = L >= 31 ? 1 : 1 << 31 - L; 1581 } 1582 } 1583#endif 1584#endif 1585 return dval(a); 1586} 1587 1588static double 1589b2d(Bigint *a, int *e) 1590{ 1591 ULong *xa, *xa0, w, y, z; 1592 int k; 1593 double_u d; 1594#ifdef VAX 1595 ULong d0, d1; 1596#else 1597#define d0 word0(d) 1598#define d1 word1(d) 1599#endif 1600 1601 xa0 = a->x; 1602 xa = xa0 + a->wds; 1603 y = *--xa; 1604#ifdef DEBUG 1605 if (!y) Bug("zero y in b2d"); 1606#endif 1607 k = hi0bits(y); 1608 *e = 32 - k; 1609#ifdef Pack_32 1610 if (k < Ebits) { 1611 d0 = Exp_1 | y >> (Ebits - k); 1612 w = xa > xa0 ? *--xa : 0; 1613 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k); 1614 goto ret_d; 1615 } 1616 z = xa > xa0 ? *--xa : 0; 1617 if (k -= Ebits) { 1618 d0 = Exp_1 | y << k | z >> (32 - k); 1619 y = xa > xa0 ? *--xa : 0; 1620 d1 = z << k | y >> (32 - k); 1621 } 1622 else { 1623 d0 = Exp_1 | y; 1624 d1 = z; 1625 } 1626#else 1627 if (k < Ebits + 16) { 1628 z = xa > xa0 ? *--xa : 0; 1629 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 1630 w = xa > xa0 ? *--xa : 0; 1631 y = xa > xa0 ? *--xa : 0; 1632 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 1633 goto ret_d; 1634 } 1635 z = xa > xa0 ? *--xa : 0; 1636 w = xa > xa0 ? *--xa : 0; 1637 k -= Ebits + 16; 1638 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 1639 y = xa > xa0 ? *--xa : 0; 1640 d1 = w << k + 16 | y << k; 1641#endif 1642ret_d: 1643#ifdef VAX 1644 word0(d) = d0 >> 16 | d0 << 16; 1645 word1(d) = d1 >> 16 | d1 << 16; 1646#else 1647#undef d0 1648#undef d1 1649#endif 1650 return dval(d); 1651} 1652 1653static Bigint * 1654d2b(double d_, int *e, int *bits) 1655{ 1656 double_u d; 1657 Bigint *b; 1658 int de, k; 1659 ULong *x, y, z; 1660#ifndef Sudden_Underflow 1661 int i; 1662#endif 1663#ifdef VAX 1664 ULong d0, d1; 1665#endif 1666 dval(d) = d_; 1667#ifdef VAX 1668 d0 = word0(d) >> 16 | word0(d) << 16; 1669 d1 = word1(d) >> 16 | word1(d) << 16; 1670#else 1671#define d0 word0(d) 1672#define d1 word1(d) 1673#endif 1674 1675#ifdef Pack_32 1676 b = Balloc(1); 1677#else 1678 b = Balloc(2); 1679#endif 1680 x = b->x; 1681 1682 z = d0 & Frac_mask; 1683 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 1684#ifdef Sudden_Underflow 1685 de = (int)(d0 >> Exp_shift); 1686#ifndef IBM 1687 z |= Exp_msk11; 1688#endif 1689#else 1690 if ((de = (int)(d0 >> Exp_shift)) != 0) 1691 z |= Exp_msk1; 1692#endif 1693#ifdef Pack_32 1694 if ((y = d1) != 0) { 1695 if ((k = lo0bits(&y)) != 0) { 1696 x[0] = y | z << (32 - k); 1697 z >>= k; 1698 } 1699 else 1700 x[0] = y; 1701#ifndef Sudden_Underflow 1702 i = 1703#endif 1704 b->wds = (x[1] = z) ? 2 : 1; 1705 } 1706 else { 1707#ifdef DEBUG 1708 if (!z) 1709 Bug("Zero passed to d2b"); 1710#endif 1711 k = lo0bits(&z); 1712 x[0] = z; 1713#ifndef Sudden_Underflow 1714 i = 1715#endif 1716 b->wds = 1; 1717 k += 32; 1718 } 1719#else 1720 if (y = d1) { 1721 if (k = lo0bits(&y)) 1722 if (k >= 16) { 1723 x[0] = y | z << 32 - k & 0xffff; 1724 x[1] = z >> k - 16 & 0xffff; 1725 x[2] = z >> k; 1726 i = 2; 1727 } 1728 else { 1729 x[0] = y & 0xffff; 1730 x[1] = y >> 16 | z << 16 - k & 0xffff; 1731 x[2] = z >> k & 0xffff; 1732 x[3] = z >> k+16; 1733 i = 3; 1734 } 1735 else { 1736 x[0] = y & 0xffff; 1737 x[1] = y >> 16; 1738 x[2] = z & 0xffff; 1739 x[3] = z >> 16; 1740 i = 3; 1741 } 1742 } 1743 else { 1744#ifdef DEBUG 1745 if (!z) 1746 Bug("Zero passed to d2b"); 1747#endif 1748 k = lo0bits(&z); 1749 if (k >= 16) { 1750 x[0] = z; 1751 i = 0; 1752 } 1753 else { 1754 x[0] = z & 0xffff; 1755 x[1] = z >> 16; 1756 i = 1; 1757 } 1758 k += 32; 1759 } 1760 while (!x[i]) 1761 --i; 1762 b->wds = i + 1; 1763#endif 1764#ifndef Sudden_Underflow 1765 if (de) { 1766#endif 1767#ifdef IBM 1768 *e = (de - Bias - (P-1) << 2) + k; 1769 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 1770#else 1771 *e = de - Bias - (P-1) + k; 1772 *bits = P - k; 1773#endif 1774#ifndef Sudden_Underflow 1775 } 1776 else { 1777 *e = de - Bias - (P-1) + 1 + k; 1778#ifdef Pack_32 1779 *bits = 32*i - hi0bits(x[i-1]); 1780#else 1781 *bits = (i+2)*16 - hi0bits(x[i]); 1782#endif 1783 } 1784#endif 1785 return b; 1786} 1787#undef d0 1788#undef d1 1789 1790static double 1791ratio(Bigint *a, Bigint *b) 1792{ 1793 double_u da, db; 1794 int k, ka, kb; 1795 1796 dval(da) = b2d(a, &ka); 1797 dval(db) = b2d(b, &kb); 1798#ifdef Pack_32 1799 k = ka - kb + 32*(a->wds - b->wds); 1800#else 1801 k = ka - kb + 16*(a->wds - b->wds); 1802#endif 1803#ifdef IBM 1804 if (k > 0) { 1805 word0(da) += (k >> 2)*Exp_msk1; 1806 if (k &= 3) 1807 dval(da) *= 1 << k; 1808 } 1809 else { 1810 k = -k; 1811 word0(db) += (k >> 2)*Exp_msk1; 1812 if (k &= 3) 1813 dval(db) *= 1 << k; 1814 } 1815#else 1816 if (k > 0) 1817 word0(da) += k*Exp_msk1; 1818 else { 1819 k = -k; 1820 word0(db) += k*Exp_msk1; 1821 } 1822#endif 1823 return dval(da) / dval(db); 1824} 1825 1826static const double 1827tens[] = { 1828 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1829 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1830 1e20, 1e21, 1e22 1831#ifdef VAX 1832 , 1e23, 1e24 1833#endif 1834}; 1835 1836static const double 1837#ifdef IEEE_Arith 1838bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 1839static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1840#ifdef Avoid_Underflow 1841 9007199254740992.*9007199254740992.e-256 1842 /* = 2^106 * 1e-53 */ 1843#else 1844 1e-256 1845#endif 1846}; 1847/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 1848/* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 1849#define Scale_Bit 0x10 1850#define n_bigtens 5 1851#else 1852#ifdef IBM 1853bigtens[] = { 1e16, 1e32, 1e64 }; 1854static const double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 1855#define n_bigtens 3 1856#else 1857bigtens[] = { 1e16, 1e32 }; 1858static const double tinytens[] = { 1e-16, 1e-32 }; 1859#define n_bigtens 2 1860#endif 1861#endif 1862 1863#ifndef IEEE_Arith 1864#undef INFNAN_CHECK 1865#endif 1866 1867#ifdef INFNAN_CHECK 1868 1869#ifndef NAN_WORD0 1870#define NAN_WORD0 0x7ff80000 1871#endif 1872 1873#ifndef NAN_WORD1 1874#define NAN_WORD1 0 1875#endif 1876 1877static int 1878match(const char **sp, char *t) 1879{ 1880 int c, d; 1881 const char *s = *sp; 1882 1883 while (d = *t++) { 1884 if ((c = *++s) >= 'A' && c <= 'Z') 1885 c += 'a' - 'A'; 1886 if (c != d) 1887 return 0; 1888 } 1889 *sp = s + 1; 1890 return 1; 1891} 1892 1893#ifndef No_Hex_NaN 1894static void 1895hexnan(double *rvp, const char **sp) 1896{ 1897 ULong c, x[2]; 1898 const char *s; 1899 int havedig, udx0, xshift; 1900 1901 x[0] = x[1] = 0; 1902 havedig = xshift = 0; 1903 udx0 = 1; 1904 s = *sp; 1905 while (c = *(const unsigned char*)++s) { 1906 if (c >= '0' && c <= '9') 1907 c -= '0'; 1908 else if (c >= 'a' && c <= 'f') 1909 c += 10 - 'a'; 1910 else if (c >= 'A' && c <= 'F') 1911 c += 10 - 'A'; 1912 else if (c <= ' ') { 1913 if (udx0 && havedig) { 1914 udx0 = 0; 1915 xshift = 1; 1916 } 1917 continue; 1918 } 1919 else if (/*(*/ c == ')' && havedig) { 1920 *sp = s + 1; 1921 break; 1922 } 1923 else 1924 return; /* invalid form: don't change *sp */ 1925 havedig = 1; 1926 if (xshift) { 1927 xshift = 0; 1928 x[0] = x[1]; 1929 x[1] = 0; 1930 } 1931 if (udx0) 1932 x[0] = (x[0] << 4) | (x[1] >> 28); 1933 x[1] = (x[1] << 4) | c; 1934 } 1935 if ((x[0] &= 0xfffff) || x[1]) { 1936 word0(*rvp) = Exp_mask | x[0]; 1937 word1(*rvp) = x[1]; 1938 } 1939} 1940#endif /*No_Hex_NaN*/ 1941#endif /* INFNAN_CHECK */ 1942 1943double 1944ruby_strtod(const char *s00, char **se) 1945{ 1946#ifdef Avoid_Underflow 1947 int scale; 1948#endif 1949 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, 1950 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 1951 const char *s, *s0, *s1; 1952 double aadj, adj; 1953 double_u aadj1, rv, rv0; 1954 Long L; 1955 ULong y, z; 1956 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; 1957#ifdef SET_INEXACT 1958 int inexact, oldinexact; 1959#endif 1960#ifdef Honor_FLT_ROUNDS 1961 int rounding; 1962#endif 1963#ifdef USE_LOCALE 1964 const char *s2; 1965#endif 1966 1967 errno = 0; 1968 sign = nz0 = nz = 0; 1969 dval(rv) = 0.; 1970 for (s = s00;;s++) 1971 switch (*s) { 1972 case '-': 1973 sign = 1; 1974 /* no break */ 1975 case '+': 1976 if (*++s) 1977 goto break2; 1978 /* no break */ 1979 case 0: 1980 goto ret0; 1981 case '\t': 1982 case '\n': 1983 case '\v': 1984 case '\f': 1985 case '\r': 1986 case ' ': 1987 continue; 1988 default: 1989 goto break2; 1990 } 1991break2: 1992 if (*s == '0') { 1993 if (s[1] == 'x' || s[1] == 'X') { 1994 static const char hexdigit[] = "0123456789abcdef0123456789ABCDEF"; 1995 s0 = ++s; 1996 adj = 0; 1997 aadj = 1.0; 1998 nd0 = -4; 1999 2000 if (!*++s || !(s1 = strchr(hexdigit, *s))) goto ret0; 2001 if (*s == '0') { 2002 while (*++s == '0'); 2003 s1 = strchr(hexdigit, *s); 2004 } 2005 if (s1 != NULL) { 2006 do { 2007 adj += aadj * ((s1 - hexdigit) & 15); 2008 nd0 += 4; 2009 aadj /= 16; 2010 } while (*++s && (s1 = strchr(hexdigit, *s))); 2011 } 2012 2013 if (*s == '.') { 2014 dsign = 1; 2015 if (!*++s || !(s1 = strchr(hexdigit, *s))) goto ret0; 2016 if (nd0 < 0) { 2017 while (*s == '0') { 2018 s++; 2019 nd0 -= 4; 2020 } 2021 } 2022 for (; *s && (s1 = strchr(hexdigit, *s)); ++s) { 2023 adj += aadj * ((s1 - hexdigit) & 15); 2024 if ((aadj /= 16) == 0.0) { 2025 while (strchr(hexdigit, *++s)); 2026 break; 2027 } 2028 } 2029 } 2030 else { 2031 dsign = 0; 2032 } 2033 2034 if (*s == 'P' || *s == 'p') { 2035 dsign = 0x2C - *++s; /* +: 2B, -: 2D */ 2036 if (abs(dsign) == 1) s++; 2037 else dsign = 1; 2038 2039 nd = 0; 2040 c = *s; 2041 if (c < '0' || '9' < c) goto ret0; 2042 do { 2043 nd *= 10; 2044 nd += c; 2045 nd -= '0'; 2046 c = *++s; 2047 /* Float("0x0."+("0"*267)+"1fp2095") */ 2048 if (nd + dsign * nd0 > 2095) { 2049 while ('0' <= c && c <= '9') c = *++s; 2050 break; 2051 } 2052 } while ('0' <= c && c <= '9'); 2053 nd0 += nd * dsign; 2054 } 2055 else { 2056 if (dsign) goto ret0; 2057 } 2058 dval(rv) = ldexp(adj, nd0); 2059 goto ret; 2060 } 2061 nz0 = 1; 2062 while (*++s == '0') ; 2063 if (!*s) 2064 goto ret; 2065 } 2066 s0 = s; 2067 y = z = 0; 2068 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 2069 if (nd < 9) 2070 y = 10*y + c - '0'; 2071 else if (nd < 16) 2072 z = 10*z + c - '0'; 2073 nd0 = nd; 2074#ifdef USE_LOCALE 2075 s1 = localeconv()->decimal_point; 2076 if (c == *s1) { 2077 c = '.'; 2078 if (*++s1) { 2079 s2 = s; 2080 for (;;) { 2081 if (*++s2 != *s1) { 2082 c = 0; 2083 break; 2084 } 2085 if (!*++s1) { 2086 s = s2; 2087 break; 2088 } 2089 } 2090 } 2091 } 2092#endif 2093 if (c == '.') { 2094 if (!ISDIGIT(s[1])) 2095 goto dig_done; 2096 c = *++s; 2097 if (!nd) { 2098 for (; c == '0'; c = *++s) 2099 nz++; 2100 if (c > '0' && c <= '9') { 2101 s0 = s; 2102 nf += nz; 2103 nz = 0; 2104 goto have_dig; 2105 } 2106 goto dig_done; 2107 } 2108 for (; c >= '0' && c <= '9'; c = *++s) { 2109have_dig: 2110 nz++; 2111 if (nf > DBL_DIG * 4) continue; 2112 if (c -= '0') { 2113 nf += nz; 2114 for (i = 1; i < nz; i++) 2115 if (nd++ < 9) 2116 y *= 10; 2117 else if (nd <= DBL_DIG + 1) 2118 z *= 10; 2119 if (nd++ < 9) 2120 y = 10*y + c; 2121 else if (nd <= DBL_DIG + 1) 2122 z = 10*z + c; 2123 nz = 0; 2124 } 2125 } 2126 } 2127dig_done: 2128 e = 0; 2129 if (c == 'e' || c == 'E') { 2130 if (!nd && !nz && !nz0) { 2131 goto ret0; 2132 } 2133 s00 = s; 2134 esign = 0; 2135 switch (c = *++s) { 2136 case '-': 2137 esign = 1; 2138 case '+': 2139 c = *++s; 2140 } 2141 if (c >= '0' && c <= '9') { 2142 while (c == '0') 2143 c = *++s; 2144 if (c > '0' && c <= '9') { 2145 L = c - '0'; 2146 s1 = s; 2147 while ((c = *++s) >= '0' && c <= '9') 2148 L = 10*L + c - '0'; 2149 if (s - s1 > 8 || L > 19999) 2150 /* Avoid confusion from exponents 2151 * so large that e might overflow. 2152 */ 2153 e = 19999; /* safe for 16 bit ints */ 2154 else 2155 e = (int)L; 2156 if (esign) 2157 e = -e; 2158 } 2159 else 2160 e = 0; 2161 } 2162 else 2163 s = s00; 2164 } 2165 if (!nd) { 2166 if (!nz && !nz0) { 2167#ifdef INFNAN_CHECK 2168 /* Check for Nan and Infinity */ 2169 switch (c) { 2170 case 'i': 2171 case 'I': 2172 if (match(&s,"nf")) { 2173 --s; 2174 if (!match(&s,"inity")) 2175 ++s; 2176 word0(rv) = 0x7ff00000; 2177 word1(rv) = 0; 2178 goto ret; 2179 } 2180 break; 2181 case 'n': 2182 case 'N': 2183 if (match(&s, "an")) { 2184 word0(rv) = NAN_WORD0; 2185 word1(rv) = NAN_WORD1; 2186#ifndef No_Hex_NaN 2187 if (*s == '(') /*)*/ 2188 hexnan(&rv, &s); 2189#endif 2190 goto ret; 2191 } 2192 } 2193#endif /* INFNAN_CHECK */ 2194ret0: 2195 s = s00; 2196 sign = 0; 2197 } 2198 goto ret; 2199 } 2200 e1 = e -= nf; 2201 2202 /* Now we have nd0 digits, starting at s0, followed by a 2203 * decimal point, followed by nd-nd0 digits. The number we're 2204 * after is the integer represented by those digits times 2205 * 10**e */ 2206 2207 if (!nd0) 2208 nd0 = nd; 2209 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 2210 dval(rv) = y; 2211 if (k > 9) { 2212#ifdef SET_INEXACT 2213 if (k > DBL_DIG) 2214 oldinexact = get_inexact(); 2215#endif 2216 dval(rv) = tens[k - 9] * dval(rv) + z; 2217 } 2218 bd0 = bb = bd = bs = delta = 0; 2219 if (nd <= DBL_DIG 2220#ifndef RND_PRODQUOT 2221#ifndef Honor_FLT_ROUNDS 2222 && Flt_Rounds == 1 2223#endif 2224#endif 2225 ) { 2226 if (!e) 2227 goto ret; 2228 if (e > 0) { 2229 if (e <= Ten_pmax) { 2230#ifdef VAX 2231 goto vax_ovfl_check; 2232#else 2233#ifdef Honor_FLT_ROUNDS 2234 /* round correctly FLT_ROUNDS = 2 or 3 */ 2235 if (sign) { 2236 dval(rv) = -dval(rv); 2237 sign = 0; 2238 } 2239#endif 2240 /* rv = */ rounded_product(dval(rv), tens[e]); 2241 goto ret; 2242#endif 2243 } 2244 i = DBL_DIG - nd; 2245 if (e <= Ten_pmax + i) { 2246 /* A fancier test would sometimes let us do 2247 * this for larger i values. 2248 */ 2249#ifdef Honor_FLT_ROUNDS 2250 /* round correctly FLT_ROUNDS = 2 or 3 */ 2251 if (sign) { 2252 dval(rv) = -dval(rv); 2253 sign = 0; 2254 } 2255#endif 2256 e -= i; 2257 dval(rv) *= tens[i]; 2258#ifdef VAX 2259 /* VAX exponent range is so narrow we must 2260 * worry about overflow here... 2261 */ 2262vax_ovfl_check: 2263 word0(rv) -= P*Exp_msk1; 2264 /* rv = */ rounded_product(dval(rv), tens[e]); 2265 if ((word0(rv) & Exp_mask) 2266 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 2267 goto ovfl; 2268 word0(rv) += P*Exp_msk1; 2269#else 2270 /* rv = */ rounded_product(dval(rv), tens[e]); 2271#endif 2272 goto ret; 2273 } 2274 } 2275#ifndef Inaccurate_Divide 2276 else if (e >= -Ten_pmax) { 2277#ifdef Honor_FLT_ROUNDS 2278 /* round correctly FLT_ROUNDS = 2 or 3 */ 2279 if (sign) { 2280 dval(rv) = -dval(rv); 2281 sign = 0; 2282 } 2283#endif 2284 /* rv = */ rounded_quotient(dval(rv), tens[-e]); 2285 goto ret; 2286 } 2287#endif 2288 } 2289 e1 += nd - k; 2290 2291#ifdef IEEE_Arith 2292#ifdef SET_INEXACT 2293 inexact = 1; 2294 if (k <= DBL_DIG) 2295 oldinexact = get_inexact(); 2296#endif 2297#ifdef Avoid_Underflow 2298 scale = 0; 2299#endif 2300#ifdef Honor_FLT_ROUNDS 2301 if ((rounding = Flt_Rounds) >= 2) { 2302 if (sign) 2303 rounding = rounding == 2 ? 0 : 2; 2304 else 2305 if (rounding != 2) 2306 rounding = 0; 2307 } 2308#endif 2309#endif /*IEEE_Arith*/ 2310 2311 /* Get starting approximation = rv * 10**e1 */ 2312 2313 if (e1 > 0) { 2314 if ((i = e1 & 15) != 0) 2315 dval(rv) *= tens[i]; 2316 if (e1 &= ~15) { 2317 if (e1 > DBL_MAX_10_EXP) { 2318ovfl: 2319#ifndef NO_ERRNO 2320 errno = ERANGE; 2321#endif 2322 /* Can't trust HUGE_VAL */ 2323#ifdef IEEE_Arith 2324#ifdef Honor_FLT_ROUNDS 2325 switch (rounding) { 2326 case 0: /* toward 0 */ 2327 case 3: /* toward -infinity */ 2328 word0(rv) = Big0; 2329 word1(rv) = Big1; 2330 break; 2331 default: 2332 word0(rv) = Exp_mask; 2333 word1(rv) = 0; 2334 } 2335#else /*Honor_FLT_ROUNDS*/ 2336 word0(rv) = Exp_mask; 2337 word1(rv) = 0; 2338#endif /*Honor_FLT_ROUNDS*/ 2339#ifdef SET_INEXACT 2340 /* set overflow bit */ 2341 dval(rv0) = 1e300; 2342 dval(rv0) *= dval(rv0); 2343#endif 2344#else /*IEEE_Arith*/ 2345 word0(rv) = Big0; 2346 word1(rv) = Big1; 2347#endif /*IEEE_Arith*/ 2348 if (bd0) 2349 goto retfree; 2350 goto ret; 2351 } 2352 e1 >>= 4; 2353 for (j = 0; e1 > 1; j++, e1 >>= 1) 2354 if (e1 & 1) 2355 dval(rv) *= bigtens[j]; 2356 /* The last multiplication could overflow. */ 2357 word0(rv) -= P*Exp_msk1; 2358 dval(rv) *= bigtens[j]; 2359 if ((z = word0(rv) & Exp_mask) 2360 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 2361 goto ovfl; 2362 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 2363 /* set to largest number */ 2364 /* (Can't trust DBL_MAX) */ 2365 word0(rv) = Big0; 2366 word1(rv) = Big1; 2367 } 2368 else 2369 word0(rv) += P*Exp_msk1; 2370 } 2371 } 2372 else if (e1 < 0) { 2373 e1 = -e1; 2374 if ((i = e1 & 15) != 0) 2375 dval(rv) /= tens[i]; 2376 if (e1 >>= 4) { 2377 if (e1 >= 1 << n_bigtens) 2378 goto undfl; 2379#ifdef Avoid_Underflow 2380 if (e1 & Scale_Bit) 2381 scale = 2*P; 2382 for (j = 0; e1 > 0; j++, e1 >>= 1) 2383 if (e1 & 1) 2384 dval(rv) *= tinytens[j]; 2385 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask) 2386 >> Exp_shift)) > 0) { 2387 /* scaled rv is denormal; zap j low bits */ 2388 if (j >= 32) { 2389 word1(rv) = 0; 2390 if (j >= 53) 2391 word0(rv) = (P+2)*Exp_msk1; 2392 else 2393 word0(rv) &= 0xffffffff << (j-32); 2394 } 2395 else 2396 word1(rv) &= 0xffffffff << j; 2397 } 2398#else 2399 for (j = 0; e1 > 1; j++, e1 >>= 1) 2400 if (e1 & 1) 2401 dval(rv) *= tinytens[j]; 2402 /* The last multiplication could underflow. */ 2403 dval(rv0) = dval(rv); 2404 dval(rv) *= tinytens[j]; 2405 if (!dval(rv)) { 2406 dval(rv) = 2.*dval(rv0); 2407 dval(rv) *= tinytens[j]; 2408#endif 2409 if (!dval(rv)) { 2410undfl: 2411 dval(rv) = 0.; 2412#ifndef NO_ERRNO 2413 errno = ERANGE; 2414#endif 2415 if (bd0) 2416 goto retfree; 2417 goto ret; 2418 } 2419#ifndef Avoid_Underflow 2420 word0(rv) = Tiny0; 2421 word1(rv) = Tiny1; 2422 /* The refinement below will clean 2423 * this approximation up. 2424 */ 2425 } 2426#endif 2427 } 2428 } 2429 2430 /* Now the hard part -- adjusting rv to the correct value.*/ 2431 2432 /* Put digits into bd: true value = bd * 10^e */ 2433 2434 bd0 = s2b(s0, nd0, nd, y); 2435 2436 for (;;) { 2437 bd = Balloc(bd0->k); 2438 Bcopy(bd, bd0); 2439 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 2440 bs = i2b(1); 2441 2442 if (e >= 0) { 2443 bb2 = bb5 = 0; 2444 bd2 = bd5 = e; 2445 } 2446 else { 2447 bb2 = bb5 = -e; 2448 bd2 = bd5 = 0; 2449 } 2450 if (bbe >= 0) 2451 bb2 += bbe; 2452 else 2453 bd2 -= bbe; 2454 bs2 = bb2; 2455#ifdef Honor_FLT_ROUNDS 2456 if (rounding != 1) 2457 bs2++; 2458#endif 2459#ifdef Avoid_Underflow 2460 j = bbe - scale; 2461 i = j + bbbits - 1; /* logb(rv) */ 2462 if (i < Emin) /* denormal */ 2463 j += P - Emin; 2464 else 2465 j = P + 1 - bbbits; 2466#else /*Avoid_Underflow*/ 2467#ifdef Sudden_Underflow 2468#ifdef IBM 2469 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 2470#else 2471 j = P + 1 - bbbits; 2472#endif 2473#else /*Sudden_Underflow*/ 2474 j = bbe; 2475 i = j + bbbits - 1; /* logb(rv) */ 2476 if (i < Emin) /* denormal */ 2477 j += P - Emin; 2478 else 2479 j = P + 1 - bbbits; 2480#endif /*Sudden_Underflow*/ 2481#endif /*Avoid_Underflow*/ 2482 bb2 += j; 2483 bd2 += j; 2484#ifdef Avoid_Underflow 2485 bd2 += scale; 2486#endif 2487 i = bb2 < bd2 ? bb2 : bd2; 2488 if (i > bs2) 2489 i = bs2; 2490 if (i > 0) { 2491 bb2 -= i; 2492 bd2 -= i; 2493 bs2 -= i; 2494 } 2495 if (bb5 > 0) { 2496 bs = pow5mult(bs, bb5); 2497 bb1 = mult(bs, bb); 2498 Bfree(bb); 2499 bb = bb1; 2500 } 2501 if (bb2 > 0) 2502 bb = lshift(bb, bb2); 2503 if (bd5 > 0) 2504 bd = pow5mult(bd, bd5); 2505 if (bd2 > 0) 2506 bd = lshift(bd, bd2); 2507 if (bs2 > 0) 2508 bs = lshift(bs, bs2); 2509 delta = diff(bb, bd); 2510 dsign = delta->sign; 2511 delta->sign = 0; 2512 i = cmp(delta, bs); 2513#ifdef Honor_FLT_ROUNDS 2514 if (rounding != 1) { 2515 if (i < 0) { 2516 /* Error is less than an ulp */ 2517 if (!delta->x[0] && delta->wds <= 1) { 2518 /* exact */ 2519#ifdef SET_INEXACT 2520 inexact = 0; 2521#endif 2522 break; 2523 } 2524 if (rounding) { 2525 if (dsign) { 2526 adj = 1.; 2527 goto apply_adj; 2528 } 2529 } 2530 else if (!dsign) { 2531 adj = -1.; 2532 if (!word1(rv) 2533 && !(word0(rv) & Frac_mask)) { 2534 y = word0(rv) & Exp_mask; 2535#ifdef Avoid_Underflow 2536 if (!scale || y > 2*P*Exp_msk1) 2537#else 2538 if (y) 2539#endif 2540 { 2541 delta = lshift(delta,Log2P); 2542 if (cmp(delta, bs) <= 0) 2543 adj = -0.5; 2544 } 2545 } 2546apply_adj: 2547#ifdef Avoid_Underflow 2548 if (scale && (y = word0(rv) & Exp_mask) 2549 <= 2*P*Exp_msk1) 2550 word0(adj) += (2*P+1)*Exp_msk1 - y; 2551#else 2552#ifdef Sudden_Underflow 2553 if ((word0(rv) & Exp_mask) <= 2554 P*Exp_msk1) { 2555 word0(rv) += P*Exp_msk1; 2556 dval(rv) += adj*ulp(dval(rv)); 2557 word0(rv) -= P*Exp_msk1; 2558 } 2559 else 2560#endif /*Sudden_Underflow*/ 2561#endif /*Avoid_Underflow*/ 2562 dval(rv) += adj*ulp(dval(rv)); 2563 } 2564 break; 2565 } 2566 adj = ratio(delta, bs); 2567 if (adj < 1.) 2568 adj = 1.; 2569 if (adj <= 0x7ffffffe) { 2570 /* adj = rounding ? ceil(adj) : floor(adj); */ 2571 y = adj; 2572 if (y != adj) { 2573 if (!((rounding>>1) ^ dsign)) 2574 y++; 2575 adj = y; 2576 } 2577 } 2578#ifdef Avoid_Underflow 2579 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 2580 word0(adj) += (2*P+1)*Exp_msk1 - y; 2581#else 2582#ifdef Sudden_Underflow 2583 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 2584 word0(rv) += P*Exp_msk1; 2585 adj *= ulp(dval(rv)); 2586 if (dsign) 2587 dval(rv) += adj; 2588 else 2589 dval(rv) -= adj; 2590 word0(rv) -= P*Exp_msk1; 2591 goto cont; 2592 } 2593#endif /*Sudden_Underflow*/ 2594#endif /*Avoid_Underflow*/ 2595 adj *= ulp(dval(rv)); 2596 if (dsign) 2597 dval(rv) += adj; 2598 else 2599 dval(rv) -= adj; 2600 goto cont; 2601 } 2602#endif /*Honor_FLT_ROUNDS*/ 2603 2604 if (i < 0) { 2605 /* Error is less than half an ulp -- check for 2606 * special case of mantissa a power of two. 2607 */ 2608 if (dsign || word1(rv) || word0(rv) & Bndry_mask 2609#ifdef IEEE_Arith 2610#ifdef Avoid_Underflow 2611 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 2612#else 2613 || (word0(rv) & Exp_mask) <= Exp_msk1 2614#endif 2615#endif 2616 ) { 2617#ifdef SET_INEXACT 2618 if (!delta->x[0] && delta->wds <= 1) 2619 inexact = 0; 2620#endif 2621 break; 2622 } 2623 if (!delta->x[0] && delta->wds <= 1) { 2624 /* exact result */ 2625#ifdef SET_INEXACT 2626 inexact = 0; 2627#endif 2628 break; 2629 } 2630 delta = lshift(delta,Log2P); 2631 if (cmp(delta, bs) > 0) 2632 goto drop_down; 2633 break; 2634 } 2635 if (i == 0) { 2636 /* exactly half-way between */ 2637 if (dsign) { 2638 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 2639 && word1(rv) == ( 2640#ifdef Avoid_Underflow 2641 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 2642 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 2643#endif 2644 0xffffffff)) { 2645 /*boundary case -- increment exponent*/ 2646 word0(rv) = (word0(rv) & Exp_mask) 2647 + Exp_msk1 2648#ifdef IBM 2649 | Exp_msk1 >> 4 2650#endif 2651 ; 2652 word1(rv) = 0; 2653#ifdef Avoid_Underflow 2654 dsign = 0; 2655#endif 2656 break; 2657 } 2658 } 2659 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 2660drop_down: 2661 /* boundary case -- decrement exponent */ 2662#ifdef Sudden_Underflow /*{{*/ 2663 L = word0(rv) & Exp_mask; 2664#ifdef IBM 2665 if (L < Exp_msk1) 2666#else 2667#ifdef Avoid_Underflow 2668 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 2669#else 2670 if (L <= Exp_msk1) 2671#endif /*Avoid_Underflow*/ 2672#endif /*IBM*/ 2673 goto undfl; 2674 L -= Exp_msk1; 2675#else /*Sudden_Underflow}{*/ 2676#ifdef Avoid_Underflow 2677 if (scale) { 2678 L = word0(rv) & Exp_mask; 2679 if (L <= (2*P+1)*Exp_msk1) { 2680 if (L > (P+2)*Exp_msk1) 2681 /* round even ==> */ 2682 /* accept rv */ 2683 break; 2684 /* rv = smallest denormal */ 2685 goto undfl; 2686 } 2687 } 2688#endif /*Avoid_Underflow*/ 2689 L = (word0(rv) & Exp_mask) - Exp_msk1; 2690#endif /*Sudden_Underflow}}*/ 2691 word0(rv) = L | Bndry_mask1; 2692 word1(rv) = 0xffffffff; 2693#ifdef IBM 2694 goto cont; 2695#else 2696 break; 2697#endif 2698 } 2699#ifndef ROUND_BIASED 2700 if (!(word1(rv) & LSB)) 2701 break; 2702#endif 2703 if (dsign) 2704 dval(rv) += ulp(dval(rv)); 2705#ifndef ROUND_BIASED 2706 else { 2707 dval(rv) -= ulp(dval(rv)); 2708#ifndef Sudden_Underflow 2709 if (!dval(rv)) 2710 goto undfl; 2711#endif 2712 } 2713#ifdef Avoid_Underflow 2714 dsign = 1 - dsign; 2715#endif 2716#endif 2717 break; 2718 } 2719 if ((aadj = ratio(delta, bs)) <= 2.) { 2720 if (dsign) 2721 aadj = dval(aadj1) = 1.; 2722 else if (word1(rv) || word0(rv) & Bndry_mask) { 2723#ifndef Sudden_Underflow 2724 if (word1(rv) == Tiny1 && !word0(rv)) 2725 goto undfl; 2726#endif 2727 aadj = 1.; 2728 dval(aadj1) = -1.; 2729 } 2730 else { 2731 /* special case -- power of FLT_RADIX to be */ 2732 /* rounded down... */ 2733 2734 if (aadj < 2./FLT_RADIX) 2735 aadj = 1./FLT_RADIX; 2736 else 2737 aadj *= 0.5; 2738 dval(aadj1) = -aadj; 2739 } 2740 } 2741 else { 2742 aadj *= 0.5; 2743 dval(aadj1) = dsign ? aadj : -aadj; 2744#ifdef Check_FLT_ROUNDS 2745 switch (Rounding) { 2746 case 2: /* towards +infinity */ 2747 dval(aadj1) -= 0.5; 2748 break; 2749 case 0: /* towards 0 */ 2750 case 3: /* towards -infinity */ 2751 dval(aadj1) += 0.5; 2752 } 2753#else 2754 if (Flt_Rounds == 0) 2755 dval(aadj1) += 0.5; 2756#endif /*Check_FLT_ROUNDS*/ 2757 } 2758 y = word0(rv) & Exp_mask; 2759 2760 /* Check for overflow */ 2761 2762 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 2763 dval(rv0) = dval(rv); 2764 word0(rv) -= P*Exp_msk1; 2765 adj = dval(aadj1) * ulp(dval(rv)); 2766 dval(rv) += adj; 2767 if ((word0(rv) & Exp_mask) >= 2768 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 2769 if (word0(rv0) == Big0 && word1(rv0) == Big1) 2770 goto ovfl; 2771 word0(rv) = Big0; 2772 word1(rv) = Big1; 2773 goto cont; 2774 } 2775 else 2776 word0(rv) += P*Exp_msk1; 2777 } 2778 else { 2779#ifdef Avoid_Underflow 2780 if (scale && y <= 2*P*Exp_msk1) { 2781 if (aadj <= 0x7fffffff) { 2782 if ((z = (int)aadj) <= 0) 2783 z = 1; 2784 aadj = z; 2785 dval(aadj1) = dsign ? aadj : -aadj; 2786 } 2787 word0(aadj1) += (2*P+1)*Exp_msk1 - y; 2788 } 2789 adj = dval(aadj1) * ulp(dval(rv)); 2790 dval(rv) += adj; 2791#else 2792#ifdef Sudden_Underflow 2793 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 2794 dval(rv0) = dval(rv); 2795 word0(rv) += P*Exp_msk1; 2796 adj = dval(aadj1) * ulp(dval(rv)); 2797 dval(rv) += adj; 2798#ifdef IBM 2799 if ((word0(rv) & Exp_mask) < P*Exp_msk1) 2800#else 2801 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 2802#endif 2803 { 2804 if (word0(rv0) == Tiny0 && word1(rv0) == Tiny1) 2805 goto undfl; 2806 word0(rv) = Tiny0; 2807 word1(rv) = Tiny1; 2808 goto cont; 2809 } 2810 else 2811 word0(rv) -= P*Exp_msk1; 2812 } 2813 else { 2814 adj = dval(aadj1) * ulp(dval(rv)); 2815 dval(rv) += adj; 2816 } 2817#else /*Sudden_Underflow*/ 2818 /* Compute adj so that the IEEE rounding rules will 2819 * correctly round rv + adj in some half-way cases. 2820 * If rv * ulp(rv) is denormalized (i.e., 2821 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 2822 * trouble from bits lost to denormalization; 2823 * example: 1.2e-307 . 2824 */ 2825 if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 2826 dval(aadj1) = (double)(int)(aadj + 0.5); 2827 if (!dsign) 2828 dval(aadj1) = -dval(aadj1); 2829 } 2830 adj = dval(aadj1) * ulp(dval(rv)); 2831 dval(rv) += adj; 2832#endif /*Sudden_Underflow*/ 2833#endif /*Avoid_Underflow*/ 2834 } 2835 z = word0(rv) & Exp_mask; 2836#ifndef SET_INEXACT 2837#ifdef Avoid_Underflow 2838 if (!scale) 2839#endif 2840 if (y == z) { 2841 /* Can we stop now? */ 2842 L = (Long)aadj; 2843 aadj -= L; 2844 /* The tolerances below are conservative. */ 2845 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 2846 if (aadj < .4999999 || aadj > .5000001) 2847 break; 2848 } 2849 else if (aadj < .4999999/FLT_RADIX) 2850 break; 2851 } 2852#endif 2853cont: 2854 Bfree(bb); 2855 Bfree(bd); 2856 Bfree(bs); 2857 Bfree(delta); 2858 } 2859#ifdef SET_INEXACT 2860 if (inexact) { 2861 if (!oldinexact) { 2862 word0(rv0) = Exp_1 + (70 << Exp_shift); 2863 word1(rv0) = 0; 2864 dval(rv0) += 1.; 2865 } 2866 } 2867 else if (!oldinexact) 2868 clear_inexact(); 2869#endif 2870#ifdef Avoid_Underflow 2871 if (scale) { 2872 word0(rv0) = Exp_1 - 2*P*Exp_msk1; 2873 word1(rv0) = 0; 2874 dval(rv) *= dval(rv0); 2875#ifndef NO_ERRNO 2876 /* try to avoid the bug of testing an 8087 register value */ 2877 if (word0(rv) == 0 && word1(rv) == 0) 2878 errno = ERANGE; 2879#endif 2880 } 2881#endif /* Avoid_Underflow */ 2882#ifdef SET_INEXACT 2883 if (inexact && !(word0(rv) & Exp_mask)) { 2884 /* set underflow bit */ 2885 dval(rv0) = 1e-300; 2886 dval(rv0) *= dval(rv0); 2887 } 2888#endif 2889retfree: 2890 Bfree(bb); 2891 Bfree(bd); 2892 Bfree(bs); 2893 Bfree(bd0); 2894 Bfree(delta); 2895ret: 2896 if (se) 2897 *se = (char *)s; 2898 return sign ? -dval(rv) : dval(rv); 2899} 2900 2901static int 2902quorem(Bigint *b, Bigint *S) 2903{ 2904 int n; 2905 ULong *bx, *bxe, q, *sx, *sxe; 2906#ifdef ULLong 2907 ULLong borrow, carry, y, ys; 2908#else 2909 ULong borrow, carry, y, ys; 2910#ifdef Pack_32 2911 ULong si, z, zs; 2912#endif 2913#endif 2914 2915 n = S->wds; 2916#ifdef DEBUG 2917 /*debug*/ if (b->wds > n) 2918 /*debug*/ Bug("oversize b in quorem"); 2919#endif 2920 if (b->wds < n) 2921 return 0; 2922 sx = S->x; 2923 sxe = sx + --n; 2924 bx = b->x; 2925 bxe = bx + n; 2926 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 2927#ifdef DEBUG 2928 /*debug*/ if (q > 9) 2929 /*debug*/ Bug("oversized quotient in quorem"); 2930#endif 2931 if (q) { 2932 borrow = 0; 2933 carry = 0; 2934 do { 2935#ifdef ULLong 2936 ys = *sx++ * (ULLong)q + carry; 2937 carry = ys >> 32; 2938 y = *bx - (ys & FFFFFFFF) - borrow; 2939 borrow = y >> 32 & (ULong)1; 2940 *bx++ = (ULong)(y & FFFFFFFF); 2941#else 2942#ifdef Pack_32 2943 si = *sx++; 2944 ys = (si & 0xffff) * q + carry; 2945 zs = (si >> 16) * q + (ys >> 16); 2946 carry = zs >> 16; 2947 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 2948 borrow = (y & 0x10000) >> 16; 2949 z = (*bx >> 16) - (zs & 0xffff) - borrow; 2950 borrow = (z & 0x10000) >> 16; 2951 Storeinc(bx, z, y); 2952#else 2953 ys = *sx++ * q + carry; 2954 carry = ys >> 16; 2955 y = *bx - (ys & 0xffff) - borrow; 2956 borrow = (y & 0x10000) >> 16; 2957 *bx++ = y & 0xffff; 2958#endif 2959#endif 2960 } while (sx <= sxe); 2961 if (!*bxe) { 2962 bx = b->x; 2963 while (--bxe > bx && !*bxe) 2964 --n; 2965 b->wds = n; 2966 } 2967 } 2968 if (cmp(b, S) >= 0) { 2969 q++; 2970 borrow = 0; 2971 carry = 0; 2972 bx = b->x; 2973 sx = S->x; 2974 do { 2975#ifdef ULLong 2976 ys = *sx++ + carry; 2977 carry = ys >> 32; 2978 y = *bx - (ys & FFFFFFFF) - borrow; 2979 borrow = y >> 32 & (ULong)1; 2980 *bx++ = (ULong)(y & FFFFFFFF); 2981#else 2982#ifdef Pack_32 2983 si = *sx++; 2984 ys = (si & 0xffff) + carry; 2985 zs = (si >> 16) + (ys >> 16); 2986 carry = zs >> 16; 2987 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 2988 borrow = (y & 0x10000) >> 16; 2989 z = (*bx >> 16) - (zs & 0xffff) - borrow; 2990 borrow = (z & 0x10000) >> 16; 2991 Storeinc(bx, z, y); 2992#else 2993 ys = *sx++ + carry; 2994 carry = ys >> 16; 2995 y = *bx - (ys & 0xffff) - borrow; 2996 borrow = (y & 0x10000) >> 16; 2997 *bx++ = y & 0xffff; 2998#endif 2999#endif 3000 } while (sx <= sxe); 3001 bx = b->x; 3002 bxe = bx + n; 3003 if (!*bxe) { 3004 while (--bxe > bx && !*bxe) 3005 --n; 3006 b->wds = n; 3007 } 3008 } 3009 return q; 3010} 3011 3012#ifndef MULTIPLE_THREADS 3013static char *dtoa_result; 3014#endif 3015 3016#ifndef MULTIPLE_THREADS 3017static char * 3018rv_alloc(int i) 3019{ 3020 return dtoa_result = xmalloc(i); 3021} 3022#else 3023#define rv_alloc(i) xmalloc(i) 3024#endif 3025 3026static char * 3027nrv_alloc(const char *s, char **rve, size_t n) 3028{ 3029 char *rv, *t; 3030 3031 t = rv = rv_alloc(n); 3032 while ((*t = *s++) != 0) t++; 3033 if (rve) 3034 *rve = t; 3035 return rv; 3036} 3037 3038#define rv_strdup(s, rve) nrv_alloc((s), (rve), strlen(s)+1) 3039 3040#ifndef MULTIPLE_THREADS 3041/* freedtoa(s) must be used to free values s returned by dtoa 3042 * when MULTIPLE_THREADS is #defined. It should be used in all cases, 3043 * but for consistency with earlier versions of dtoa, it is optional 3044 * when MULTIPLE_THREADS is not defined. 3045 */ 3046 3047static void 3048freedtoa(char *s) 3049{ 3050 xfree(s); 3051} 3052#endif 3053 3054static const char INFSTR[] = "Infinity"; 3055static const char NANSTR[] = "NaN"; 3056static const char ZEROSTR[] = "0"; 3057 3058/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 3059 * 3060 * Inspired by "How to Print Floating-Point Numbers Accurately" by 3061 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. 3062 * 3063 * Modifications: 3064 * 1. Rather than iterating, we use a simple numeric overestimate 3065 * to determine k = floor(log10(d)). We scale relevant 3066 * quantities using O(log2(k)) rather than O(k) multiplications. 3067 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 3068 * try to generate digits strictly left to right. Instead, we 3069 * compute with fewer bits and propagate the carry if necessary 3070 * when rounding the final digit up. This is often faster. 3071 * 3. Under the assumption that input will be rounded nearest, 3072 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 3073 * That is, we allow equality in stopping tests when the 3074 * round-nearest rule will give the same floating-point value 3075 * as would satisfaction of the stopping test with strict 3076 * inequality. 3077 * 4. We remove common factors of powers of 2 from relevant 3078 * quantities. 3079 * 5. When converting floating-point integers less than 1e16, 3080 * we use floating-point arithmetic rather than resorting 3081 * to multiple-precision integers. 3082 * 6. When asked to produce fewer than 15 digits, we first try 3083 * to get by with floating-point arithmetic; we resort to 3084 * multiple-precision integer arithmetic only if we cannot 3085 * guarantee that the floating-point calculation has given 3086 * the correctly rounded result. For k requested digits and 3087 * "uniformly" distributed input, the probability is 3088 * something like 10^(k-15) that we must resort to the Long 3089 * calculation. 3090 */ 3091 3092char * 3093ruby_dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve) 3094{ 3095 /* Arguments ndigits, decpt, sign are similar to those 3096 of ecvt and fcvt; trailing zeros are suppressed from 3097 the returned string. If not null, *rve is set to point 3098 to the end of the return value. If d is +-Infinity or NaN, 3099 then *decpt is set to 9999. 3100 3101 mode: 3102 0 ==> shortest string that yields d when read in 3103 and rounded to nearest. 3104 1 ==> like 0, but with Steele & White stopping rule; 3105 e.g. with IEEE P754 arithmetic , mode 0 gives 3106 1e23 whereas mode 1 gives 9.999999999999999e22. 3107 2 ==> max(1,ndigits) significant digits. This gives a 3108 return value similar to that of ecvt, except 3109 that trailing zeros are suppressed. 3110 3 ==> through ndigits past the decimal point. This 3111 gives a return value similar to that from fcvt, 3112 except that trailing zeros are suppressed, and 3113 ndigits can be negative. 3114 4,5 ==> similar to 2 and 3, respectively, but (in 3115 round-nearest mode) with the tests of mode 0 to 3116 possibly return a shorter string that rounds to d. 3117 With IEEE arithmetic and compilation with 3118 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same 3119 as modes 2 and 3 when FLT_ROUNDS != 1. 3120 6-9 ==> Debugging modes similar to mode - 4: don't try 3121 fast floating-point estimate (if applicable). 3122 3123 Values of mode other than 0-9 are treated as mode 0. 3124 3125 Sufficient space is allocated to the return value 3126 to hold the suppressed trailing zeros. 3127 */ 3128 3129 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, 3130 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, 3131 spec_case, try_quick; 3132 Long L; 3133#ifndef Sudden_Underflow 3134 int denorm; 3135 ULong x; 3136#endif 3137 Bigint *b, *b1, *delta, *mlo = 0, *mhi = 0, *S; 3138 double ds; 3139 double_u d, d2, eps; 3140 char *s, *s0; 3141#ifdef Honor_FLT_ROUNDS 3142 int rounding; 3143#endif 3144#ifdef SET_INEXACT 3145 int inexact, oldinexact; 3146#endif 3147 3148 dval(d) = d_; 3149 3150#ifndef MULTIPLE_THREADS 3151 if (dtoa_result) { 3152 freedtoa(dtoa_result); 3153 dtoa_result = 0; 3154 } 3155#endif 3156 3157 if (word0(d) & Sign_bit) { 3158 /* set sign for everything, including 0's and NaNs */ 3159 *sign = 1; 3160 word0(d) &= ~Sign_bit; /* clear sign bit */ 3161 } 3162 else 3163 *sign = 0; 3164 3165#if defined(IEEE_Arith) + defined(VAX) 3166#ifdef IEEE_Arith 3167 if ((word0(d) & Exp_mask) == Exp_mask) 3168#else 3169 if (word0(d) == 0x8000) 3170#endif 3171 { 3172 /* Infinity or NaN */ 3173 *decpt = 9999; 3174#ifdef IEEE_Arith 3175 if (!word1(d) && !(word0(d) & 0xfffff)) 3176 return rv_strdup(INFSTR, rve); 3177#endif 3178 return rv_strdup(NANSTR, rve); 3179 } 3180#endif 3181#ifdef IBM 3182 dval(d) += 0; /* normalize */ 3183#endif 3184 if (!dval(d)) { 3185 *decpt = 1; 3186 return rv_strdup(ZEROSTR, rve); 3187 } 3188 3189#ifdef SET_INEXACT 3190 try_quick = oldinexact = get_inexact(); 3191 inexact = 1; 3192#endif 3193#ifdef Honor_FLT_ROUNDS 3194 if ((rounding = Flt_Rounds) >= 2) { 3195 if (*sign) 3196 rounding = rounding == 2 ? 0 : 2; 3197 else 3198 if (rounding != 2) 3199 rounding = 0; 3200 } 3201#endif 3202 3203 b = d2b(dval(d), &be, &bbits); 3204#ifdef Sudden_Underflow 3205 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); 3206#else 3207 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) { 3208#endif 3209 dval(d2) = dval(d); 3210 word0(d2) &= Frac_mask1; 3211 word0(d2) |= Exp_11; 3212#ifdef IBM 3213 if (j = 11 - hi0bits(word0(d2) & Frac_mask)) 3214 dval(d2) /= 1 << j; 3215#endif 3216 3217 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 3218 * log10(x) = log(x) / log(10) 3219 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 3220 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 3221 * 3222 * This suggests computing an approximation k to log10(d) by 3223 * 3224 * k = (i - Bias)*0.301029995663981 3225 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 3226 * 3227 * We want k to be too large rather than too small. 3228 * The error in the first-order Taylor series approximation 3229 * is in our favor, so we just round up the constant enough 3230 * to compensate for any error in the multiplication of 3231 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 3232 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 3233 * adding 1e-13 to the constant term more than suffices. 3234 * Hence we adjust the constant term to 0.1760912590558. 3235 * (We could get a more accurate k by invoking log10, 3236 * but this is probably not worthwhile.) 3237 */ 3238 3239 i -= Bias; 3240#ifdef IBM 3241 i <<= 2; 3242 i += j; 3243#endif 3244#ifndef Sudden_Underflow 3245 denorm = 0; 3246 } 3247 else { 3248 /* d is denormalized */ 3249 3250 i = bbits + be + (Bias + (P-1) - 1); 3251 x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32) 3252 : word1(d) << (32 - i); 3253 dval(d2) = x; 3254 word0(d2) -= 31*Exp_msk1; /* adjust exponent */ 3255 i -= (Bias + (P-1) - 1) + 1; 3256 denorm = 1; 3257 } 3258#endif 3259 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; 3260 k = (int)ds; 3261 if (ds < 0. && ds != k) 3262 k--; /* want k = floor(ds) */ 3263 k_check = 1; 3264 if (k >= 0 && k <= Ten_pmax) { 3265 if (dval(d) < tens[k]) 3266 k--; 3267 k_check = 0; 3268 } 3269 j = bbits - i - 1; 3270 if (j >= 0) { 3271 b2 = 0; 3272 s2 = j; 3273 } 3274 else { 3275 b2 = -j; 3276 s2 = 0; 3277 } 3278 if (k >= 0) { 3279 b5 = 0; 3280 s5 = k; 3281 s2 += k; 3282 } 3283 else { 3284 b2 -= k; 3285 b5 = -k; 3286 s5 = 0; 3287 } 3288 if (mode < 0 || mode > 9) 3289 mode = 0; 3290 3291#ifndef SET_INEXACT 3292#ifdef Check_FLT_ROUNDS 3293 try_quick = Rounding == 1; 3294#else 3295 try_quick = 1; 3296#endif 3297#endif /*SET_INEXACT*/ 3298 3299 if (mode > 5) { 3300 mode -= 4; 3301 try_quick = 0; 3302 } 3303 leftright = 1; 3304 ilim = ilim1 = -1; 3305 switch (mode) { 3306 case 0: 3307 case 1: 3308 i = 18; 3309 ndigits = 0; 3310 break; 3311 case 2: 3312 leftright = 0; 3313 /* no break */ 3314 case 4: 3315 if (ndigits <= 0) 3316 ndigits = 1; 3317 ilim = ilim1 = i = ndigits; 3318 break; 3319 case 3: 3320 leftright = 0; 3321 /* no break */ 3322 case 5: 3323 i = ndigits + k + 1; 3324 ilim = i; 3325 ilim1 = i - 1; 3326 if (i <= 0) 3327 i = 1; 3328 } 3329 s = s0 = rv_alloc(i+1); 3330 3331#ifdef Honor_FLT_ROUNDS 3332 if (mode > 1 && rounding != 1) 3333 leftright = 0; 3334#endif 3335 3336 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 3337 3338 /* Try to get by with floating-point arithmetic. */ 3339 3340 i = 0; 3341 dval(d2) = dval(d); 3342 k0 = k; 3343 ilim0 = ilim; 3344 ieps = 2; /* conservative */ 3345 if (k > 0) { 3346 ds = tens[k&0xf]; 3347 j = k >> 4; 3348 if (j & Bletch) { 3349 /* prevent overflows */ 3350 j &= Bletch - 1; 3351 dval(d) /= bigtens[n_bigtens-1]; 3352 ieps++; 3353 } 3354 for (; j; j >>= 1, i++) 3355 if (j & 1) { 3356 ieps++; 3357 ds *= bigtens[i]; 3358 } 3359 dval(d) /= ds; 3360 } 3361 else if ((j1 = -k) != 0) { 3362 dval(d) *= tens[j1 & 0xf]; 3363 for (j = j1 >> 4; j; j >>= 1, i++) 3364 if (j & 1) { 3365 ieps++; 3366 dval(d) *= bigtens[i]; 3367 } 3368 } 3369 if (k_check && dval(d) < 1. && ilim > 0) { 3370 if (ilim1 <= 0) 3371 goto fast_failed; 3372 ilim = ilim1; 3373 k--; 3374 dval(d) *= 10.; 3375 ieps++; 3376 } 3377 dval(eps) = ieps*dval(d) + 7.; 3378 word0(eps) -= (P-1)*Exp_msk1; 3379 if (ilim == 0) { 3380 S = mhi = 0; 3381 dval(d) -= 5.; 3382 if (dval(d) > dval(eps)) 3383 goto one_digit; 3384 if (dval(d) < -dval(eps)) 3385 goto no_digits; 3386 goto fast_failed; 3387 } 3388#ifndef No_leftright 3389 if (leftright) { 3390 /* Use Steele & White method of only 3391 * generating digits needed. 3392 */ 3393 dval(eps) = 0.5/tens[ilim-1] - dval(eps); 3394 for (i = 0;;) { 3395 L = (int)dval(d); 3396 dval(d) -= L; 3397 *s++ = '0' + (int)L; 3398 if (dval(d) < dval(eps)) 3399 goto ret1; 3400 if (1. - dval(d) < dval(eps)) 3401 goto bump_up; 3402 if (++i >= ilim) 3403 break; 3404 dval(eps) *= 10.; 3405 dval(d) *= 10.; 3406 } 3407 } 3408 else { 3409#endif 3410 /* Generate ilim digits, then fix them up. */ 3411 dval(eps) *= tens[ilim-1]; 3412 for (i = 1;; i++, dval(d) *= 10.) { 3413 L = (Long)(dval(d)); 3414 if (!(dval(d) -= L)) 3415 ilim = i; 3416 *s++ = '0' + (int)L; 3417 if (i == ilim) { 3418 if (dval(d) > 0.5 + dval(eps)) 3419 goto bump_up; 3420 else if (dval(d) < 0.5 - dval(eps)) { 3421 while (*--s == '0') ; 3422 s++; 3423 goto ret1; 3424 } 3425 break; 3426 } 3427 } 3428#ifndef No_leftright 3429 } 3430#endif 3431fast_failed: 3432 s = s0; 3433 dval(d) = dval(d2); 3434 k = k0; 3435 ilim = ilim0; 3436 } 3437 3438 /* Do we have a "small" integer? */ 3439 3440 if (be >= 0 && k <= Int_max) { 3441 /* Yes. */ 3442 ds = tens[k]; 3443 if (ndigits < 0 && ilim <= 0) { 3444 S = mhi = 0; 3445 if (ilim < 0 || dval(d) <= 5*ds) 3446 goto no_digits; 3447 goto one_digit; 3448 } 3449 for (i = 1;; i++, dval(d) *= 10.) { 3450 L = (Long)(dval(d) / ds); 3451 dval(d) -= L*ds; 3452#ifdef Check_FLT_ROUNDS 3453 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 3454 if (dval(d) < 0) { 3455 L--; 3456 dval(d) += ds; 3457 } 3458#endif 3459 *s++ = '0' + (int)L; 3460 if (!dval(d)) { 3461#ifdef SET_INEXACT 3462 inexact = 0; 3463#endif 3464 break; 3465 } 3466 if (i == ilim) { 3467#ifdef Honor_FLT_ROUNDS 3468 if (mode > 1) 3469 switch (rounding) { 3470 case 0: goto ret1; 3471 case 2: goto bump_up; 3472 } 3473#endif 3474 dval(d) += dval(d); 3475 if (dval(d) > ds || (dval(d) == ds && (L & 1))) { 3476bump_up: 3477 while (*--s == '9') 3478 if (s == s0) { 3479 k++; 3480 *s = '0'; 3481 break; 3482 } 3483 ++*s++; 3484 } 3485 break; 3486 } 3487 } 3488 goto ret1; 3489 } 3490 3491 m2 = b2; 3492 m5 = b5; 3493 if (leftright) { 3494 i = 3495#ifndef Sudden_Underflow 3496 denorm ? be + (Bias + (P-1) - 1 + 1) : 3497#endif 3498#ifdef IBM 3499 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); 3500#else 3501 1 + P - bbits; 3502#endif 3503 b2 += i; 3504 s2 += i; 3505 mhi = i2b(1); 3506 } 3507 if (m2 > 0 && s2 > 0) { 3508 i = m2 < s2 ? m2 : s2; 3509 b2 -= i; 3510 m2 -= i; 3511 s2 -= i; 3512 } 3513 if (b5 > 0) { 3514 if (leftright) { 3515 if (m5 > 0) { 3516 mhi = pow5mult(mhi, m5); 3517 b1 = mult(mhi, b); 3518 Bfree(b); 3519 b = b1; 3520 } 3521 if ((j = b5 - m5) != 0) 3522 b = pow5mult(b, j); 3523 } 3524 else 3525 b = pow5mult(b, b5); 3526 } 3527 S = i2b(1); 3528 if (s5 > 0) 3529 S = pow5mult(S, s5); 3530 3531 /* Check for special case that d is a normalized power of 2. */ 3532 3533 spec_case = 0; 3534 if ((mode < 2 || leftright) 3535#ifdef Honor_FLT_ROUNDS 3536 && rounding == 1 3537#endif 3538 ) { 3539 if (!word1(d) && !(word0(d) & Bndry_mask) 3540#ifndef Sudden_Underflow 3541 && word0(d) & (Exp_mask & ~Exp_msk1) 3542#endif 3543 ) { 3544 /* The special case */ 3545 b2 += Log2P; 3546 s2 += Log2P; 3547 spec_case = 1; 3548 } 3549 } 3550 3551 /* Arrange for convenient computation of quotients: 3552 * shift left if necessary so divisor has 4 leading 0 bits. 3553 * 3554 * Perhaps we should just compute leading 28 bits of S once 3555 * and for all and pass them and a shift to quorem, so it 3556 * can do shifts and ors to compute the numerator for q. 3557 */ 3558#ifdef Pack_32 3559 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0) 3560 i = 32 - i; 3561#else 3562 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) != 0) 3563 i = 16 - i; 3564#endif 3565 if (i > 4) { 3566 i -= 4; 3567 b2 += i; 3568 m2 += i; 3569 s2 += i; 3570 } 3571 else if (i < 4) { 3572 i += 28; 3573 b2 += i; 3574 m2 += i; 3575 s2 += i; 3576 } 3577 if (b2 > 0) 3578 b = lshift(b, b2); 3579 if (s2 > 0) 3580 S = lshift(S, s2); 3581 if (k_check) { 3582 if (cmp(b,S) < 0) { 3583 k--; 3584 b = multadd(b, 10, 0); /* we botched the k estimate */ 3585 if (leftright) 3586 mhi = multadd(mhi, 10, 0); 3587 ilim = ilim1; 3588 } 3589 } 3590 if (ilim <= 0 && (mode == 3 || mode == 5)) { 3591 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 3592 /* no digits, fcvt style */ 3593no_digits: 3594 k = -1 - ndigits; 3595 goto ret; 3596 } 3597one_digit: 3598 *s++ = '1'; 3599 k++; 3600 goto ret; 3601 } 3602 if (leftright) { 3603 if (m2 > 0) 3604 mhi = lshift(mhi, m2); 3605 3606 /* Compute mlo -- check for special case 3607 * that d is a normalized power of 2. 3608 */ 3609 3610 mlo = mhi; 3611 if (spec_case) { 3612 mhi = Balloc(mhi->k); 3613 Bcopy(mhi, mlo); 3614 mhi = lshift(mhi, Log2P); 3615 } 3616 3617 for (i = 1;;i++) { 3618 dig = quorem(b,S) + '0'; 3619 /* Do we yet have the shortest decimal string 3620 * that will round to d? 3621 */ 3622 j = cmp(b, mlo); 3623 delta = diff(S, mhi); 3624 j1 = delta->sign ? 1 : cmp(b, delta); 3625 Bfree(delta); 3626#ifndef ROUND_BIASED 3627 if (j1 == 0 && mode != 1 && !(word1(d) & 1) 3628#ifdef Honor_FLT_ROUNDS 3629 && rounding >= 1 3630#endif 3631 ) { 3632 if (dig == '9') 3633 goto round_9_up; 3634 if (j > 0) 3635 dig++; 3636#ifdef SET_INEXACT 3637 else if (!b->x[0] && b->wds <= 1) 3638 inexact = 0; 3639#endif 3640 *s++ = dig; 3641 goto ret; 3642 } 3643#endif 3644 if (j < 0 || (j == 0 && mode != 1 3645#ifndef ROUND_BIASED 3646 && !(word1(d) & 1) 3647#endif 3648 )) { 3649 if (!b->x[0] && b->wds <= 1) { 3650#ifdef SET_INEXACT 3651 inexact = 0; 3652#endif 3653 goto accept_dig; 3654 } 3655#ifdef Honor_FLT_ROUNDS 3656 if (mode > 1) 3657 switch (rounding) { 3658 case 0: goto accept_dig; 3659 case 2: goto keep_dig; 3660 } 3661#endif /*Honor_FLT_ROUNDS*/ 3662 if (j1 > 0) { 3663 b = lshift(b, 1); 3664 j1 = cmp(b, S); 3665 if ((j1 > 0 || (j1 == 0 && (dig & 1))) && dig++ == '9') 3666 goto round_9_up; 3667 } 3668accept_dig: 3669 *s++ = dig; 3670 goto ret; 3671 } 3672 if (j1 > 0) { 3673#ifdef Honor_FLT_ROUNDS 3674 if (!rounding) 3675 goto accept_dig; 3676#endif 3677 if (dig == '9') { /* possible if i == 1 */ 3678round_9_up: 3679 *s++ = '9'; 3680 goto roundoff; 3681 } 3682 *s++ = dig + 1; 3683 goto ret; 3684 } 3685#ifdef Honor_FLT_ROUNDS 3686keep_dig: 3687#endif 3688 *s++ = dig; 3689 if (i == ilim) 3690 break; 3691 b = multadd(b, 10, 0); 3692 if (mlo == mhi) 3693 mlo = mhi = multadd(mhi, 10, 0); 3694 else { 3695 mlo = multadd(mlo, 10, 0); 3696 mhi = multadd(mhi, 10, 0); 3697 } 3698 } 3699 } 3700 else 3701 for (i = 1;; i++) { 3702 *s++ = dig = quorem(b,S) + '0'; 3703 if (!b->x[0] && b->wds <= 1) { 3704#ifdef SET_INEXACT 3705 inexact = 0; 3706#endif 3707 goto ret; 3708 } 3709 if (i >= ilim) 3710 break; 3711 b = multadd(b, 10, 0); 3712 } 3713 3714 /* Round off last digit */ 3715 3716#ifdef Honor_FLT_ROUNDS 3717 switch (rounding) { 3718 case 0: goto trimzeros; 3719 case 2: goto roundoff; 3720 } 3721#endif 3722 b = lshift(b, 1); 3723 j = cmp(b, S); 3724 if (j > 0 || (j == 0 && (dig & 1))) { 3725 roundoff: 3726 while (*--s == '9') 3727 if (s == s0) { 3728 k++; 3729 *s++ = '1'; 3730 goto ret; 3731 } 3732 ++*s++; 3733 } 3734 else { 3735 while (*--s == '0') ; 3736 s++; 3737 } 3738ret: 3739 Bfree(S); 3740 if (mhi) { 3741 if (mlo && mlo != mhi) 3742 Bfree(mlo); 3743 Bfree(mhi); 3744 } 3745ret1: 3746#ifdef SET_INEXACT 3747 if (inexact) { 3748 if (!oldinexact) { 3749 word0(d) = Exp_1 + (70 << Exp_shift); 3750 word1(d) = 0; 3751 dval(d) += 1.; 3752 } 3753 } 3754 else if (!oldinexact) 3755 clear_inexact(); 3756#endif 3757 Bfree(b); 3758 *s = 0; 3759 *decpt = k + 1; 3760 if (rve) 3761 *rve = s; 3762 return s0; 3763} 3764 3765void 3766ruby_each_words(const char *str, void (*func)(const char*, int, void*), void *arg) 3767{ 3768 const char *end; 3769 int len; 3770 3771 if (!str) return; 3772 for (; *str; str = end) { 3773 while (ISSPACE(*str) || *str == ',') str++; 3774 if (!*str) break; 3775 end = str; 3776 while (*end && !ISSPACE(*end) && *end != ',') end++; 3777 len = (int)(end - str); /* assume no string exceeds INT_MAX */ 3778 (*func)(str, len, arg); 3779 } 3780} 3781 3782/*- 3783 * Copyright (c) 2004-2008 David Schultz <das@FreeBSD.ORG> 3784 * All rights reserved. 3785 * 3786 * Redistribution and use in source and binary forms, with or without 3787 * modification, are permitted provided that the following conditions 3788 * are met: 3789 * 1. Redistributions of source code must retain the above copyright 3790 * notice, this list of conditions and the following disclaimer. 3791 * 2. Redistributions in binary form must reproduce the above copyright 3792 * notice, this list of conditions and the following disclaimer in the 3793 * documentation and/or other materials provided with the distribution. 3794 * 3795 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 3796 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 3797 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 3798 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 3799 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 3800 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 3801 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 3802 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 3803 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 3804 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 3805 * SUCH DAMAGE. 3806 */ 3807 3808#define DBL_MANH_SIZE 20 3809#define DBL_MANL_SIZE 32 3810#define DBL_ADJ (DBL_MAX_EXP - 2) 3811#define SIGFIGS ((DBL_MANT_DIG + 3) / 4 + 1) 3812#define dexp_get(u) ((int)(word0(u) >> Exp_shift) & ~Exp_msk1) 3813#define dexp_set(u,v) (word0(u) = (((int)(word0(u)) & ~Exp_mask) | ((v) << Exp_shift))) 3814#define dmanh_get(u) ((uint32_t)(word0(u) & Frac_mask)) 3815#define dmanl_get(u) ((uint32_t)word1(u)) 3816 3817 3818/* 3819 * This procedure converts a double-precision number in IEEE format 3820 * into a string of hexadecimal digits and an exponent of 2. Its 3821 * behavior is bug-for-bug compatible with dtoa() in mode 2, with the 3822 * following exceptions: 3823 * 3824 * - An ndigits < 0 causes it to use as many digits as necessary to 3825 * represent the number exactly. 3826 * - The additional xdigs argument should point to either the string 3827 * "0123456789ABCDEF" or the string "0123456789abcdef", depending on 3828 * which case is desired. 3829 * - This routine does not repeat dtoa's mistake of setting decpt 3830 * to 9999 in the case of an infinity or NaN. INT_MAX is used 3831 * for this purpose instead. 3832 * 3833 * Note that the C99 standard does not specify what the leading digit 3834 * should be for non-zero numbers. For instance, 0x1.3p3 is the same 3835 * as 0x2.6p2 is the same as 0x4.cp3. This implementation always makes 3836 * the leading digit a 1. This ensures that the exponent printed is the 3837 * actual base-2 exponent, i.e., ilogb(d). 3838 * 3839 * Inputs: d, xdigs, ndigits 3840 * Outputs: decpt, sign, rve 3841 */ 3842char * 3843ruby_hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign, 3844 char **rve) 3845{ 3846 U u; 3847 char *s, *s0; 3848 int bufsize; 3849 uint32_t manh, manl; 3850 3851 u.d = d; 3852 if (word0(u) & Sign_bit) { 3853 /* set sign for everything, including 0's and NaNs */ 3854 *sign = 1; 3855 word0(u) &= ~Sign_bit; /* clear sign bit */ 3856 } 3857 else 3858 *sign = 0; 3859 3860 if (isinf(d)) { /* FP_INFINITE */ 3861 *decpt = INT_MAX; 3862 return rv_strdup(INFSTR, rve); 3863 } 3864 else if (isnan(d)) { /* FP_NAN */ 3865 *decpt = INT_MAX; 3866 return rv_strdup(NANSTR, rve); 3867 } 3868 else if (d == 0.0) { /* FP_ZERO */ 3869 *decpt = 1; 3870 return rv_strdup(ZEROSTR, rve); 3871 } 3872 else if (dexp_get(u)) { /* FP_NORMAL */ 3873 *decpt = dexp_get(u) - DBL_ADJ; 3874 } 3875 else { /* FP_SUBNORMAL */ 3876 u.d *= 5.363123171977039e+154 /* 0x1p514 */; 3877 *decpt = dexp_get(u) - (514 + DBL_ADJ); 3878 } 3879 3880 if (ndigits == 0) /* dtoa() compatibility */ 3881 ndigits = 1; 3882 3883 /* 3884 * If ndigits < 0, we are expected to auto-size, so we allocate 3885 * enough space for all the digits. 3886 */ 3887 bufsize = (ndigits > 0) ? ndigits : SIGFIGS; 3888 s0 = rv_alloc(bufsize+1); 3889 3890 /* Round to the desired number of digits. */ 3891 if (SIGFIGS > ndigits && ndigits > 0) { 3892 float redux = 1.0f; 3893 volatile double d; 3894 int offset = 4 * ndigits + DBL_MAX_EXP - 4 - DBL_MANT_DIG; 3895 dexp_set(u, offset); 3896 d = u.d; 3897 d += redux; 3898 d -= redux; 3899 u.d = d; 3900 *decpt += dexp_get(u) - offset; 3901 } 3902 3903 manh = dmanh_get(u); 3904 manl = dmanl_get(u); 3905 *s0 = '1'; 3906 for (s = s0 + 1; s < s0 + bufsize; s++) { 3907 *s = xdigs[(manh >> (DBL_MANH_SIZE - 4)) & 0xf]; 3908 manh = (manh << 4) | (manl >> (DBL_MANL_SIZE - 4)); 3909 manl <<= 4; 3910 } 3911 3912 /* If ndigits < 0, we are expected to auto-size the precision. */ 3913 if (ndigits < 0) { 3914 for (ndigits = SIGFIGS; s0[ndigits - 1] == '0'; ndigits--) 3915 ; 3916 } 3917 3918 s = s0 + ndigits; 3919 *s = '\0'; 3920 if (rve != NULL) 3921 *rve = s; 3922 return (s0); 3923} 3924 3925#ifdef __cplusplus 3926#if 0 3927{ /* satisfy cc-mode */ 3928#endif 3929} 3930#endif 3931