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