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 aadj = -1;
02124
02125 if (!s[1]) {
02126 rb_warn("malformed value for Float(): %s. Ruby 1.9.3 for later will raise an ArgumentError for the value.", s00);
02127 }
02128 while (*++s && (s1 = strchr(hexdigit, *s))) {
02129 adj *= 16;
02130 adj += (s1 - hexdigit) & 15;
02131 }
02132
02133 if (*s == '.') {
02134 aadj = 1.;
02135 while (*++s && (s1 = strchr(hexdigit, *s))) {
02136 aadj /= 16;
02137 adj += aadj * ((s1 - hexdigit) & 15);
02138 }
02139 }
02140
02141 if (*s == 'P' || *s == 'p') {
02142 dsign = 0x2C - *++s;
02143 if (abs(dsign) == 1) s++;
02144 else dsign = 1;
02145
02146 nd = 0;
02147 c = *s;
02148 if (c < '0' || '9' < c) goto ret0;
02149 do {
02150 nd *= 10;
02151 nd += c;
02152 nd -= '0';
02153 c = *++s;
02154 } while ('0' <= c && c <= '9');
02155 dval(rv) = ldexp(adj, nd * dsign);
02156 }
02157 else {
02158 if (aadj != -1) {
02159 rb_warn("malformed value for Float(): %s. Ruby 1.9.3 for later will raise an ArgumentError for the value.", s00);
02160 }
02161 dval(rv) = adj;
02162 }
02163 goto ret;
02164 }
02165 nz0 = 1;
02166 while (*++s == '0') ;
02167 if (!*s)
02168 goto ret;
02169 }
02170 s0 = s;
02171 y = z = 0;
02172 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
02173 if (nd < 9)
02174 y = 10*y + c - '0';
02175 else if (nd < 16)
02176 z = 10*z + c - '0';
02177 nd0 = nd;
02178 #ifdef USE_LOCALE
02179 s1 = localeconv()->decimal_point;
02180 if (c == *s1) {
02181 c = '.';
02182 if (*++s1) {
02183 s2 = s;
02184 for (;;) {
02185 if (*++s2 != *s1) {
02186 c = 0;
02187 break;
02188 }
02189 if (!*++s1) {
02190 s = s2;
02191 break;
02192 }
02193 }
02194 }
02195 }
02196 #endif
02197 if (c == '.') {
02198 if (!ISDIGIT(s[1]))
02199 goto dig_done;
02200 c = *++s;
02201 if (!nd) {
02202 for (; c == '0'; c = *++s)
02203 nz++;
02204 if (c > '0' && c <= '9') {
02205 s0 = s;
02206 nf += nz;
02207 nz = 0;
02208 goto have_dig;
02209 }
02210 goto dig_done;
02211 }
02212 for (; c >= '0' && c <= '9'; c = *++s) {
02213 have_dig:
02214 nz++;
02215 if (c -= '0') {
02216 nf += nz;
02217 for (i = 1; i < nz; i++)
02218 if (nd++ < 9)
02219 y *= 10;
02220 else if (nd <= DBL_DIG + 1)
02221 z *= 10;
02222 if (nd++ < 9)
02223 y = 10*y + c;
02224 else if (nd <= DBL_DIG + 1)
02225 z = 10*z + c;
02226 nz = 0;
02227 }
02228 }
02229 }
02230 dig_done:
02231 e = 0;
02232 if (c == 'e' || c == 'E') {
02233 if (!nd && !nz && !nz0) {
02234 goto ret0;
02235 }
02236 s00 = s;
02237 esign = 0;
02238 switch (c = *++s) {
02239 case '-':
02240 esign = 1;
02241 case '+':
02242 c = *++s;
02243 }
02244 if (c >= '0' && c <= '9') {
02245 while (c == '0')
02246 c = *++s;
02247 if (c > '0' && c <= '9') {
02248 L = c - '0';
02249 s1 = s;
02250 while ((c = *++s) >= '0' && c <= '9')
02251 L = 10*L + c - '0';
02252 if (s - s1 > 8 || L > 19999)
02253
02254
02255
02256 e = 19999;
02257 else
02258 e = (int)L;
02259 if (esign)
02260 e = -e;
02261 }
02262 else
02263 e = 0;
02264 }
02265 else
02266 s = s00;
02267 }
02268 if (!nd) {
02269 if (!nz && !nz0) {
02270 #ifdef INFNAN_CHECK
02271
02272 switch (c) {
02273 case 'i':
02274 case 'I':
02275 if (match(&s,"nf")) {
02276 --s;
02277 if (!match(&s,"inity"))
02278 ++s;
02279 word0(rv) = 0x7ff00000;
02280 word1(rv) = 0;
02281 goto ret;
02282 }
02283 break;
02284 case 'n':
02285 case 'N':
02286 if (match(&s, "an")) {
02287 word0(rv) = NAN_WORD0;
02288 word1(rv) = NAN_WORD1;
02289 #ifndef No_Hex_NaN
02290 if (*s == '(')
02291 hexnan(&rv, &s);
02292 #endif
02293 goto ret;
02294 }
02295 }
02296 #endif
02297 ret0:
02298 s = s00;
02299 sign = 0;
02300 }
02301 goto ret;
02302 }
02303 e1 = e -= nf;
02304
02305
02306
02307
02308
02309
02310 if (!nd0)
02311 nd0 = nd;
02312 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
02313 dval(rv) = y;
02314 if (k > 9) {
02315 #ifdef SET_INEXACT
02316 if (k > DBL_DIG)
02317 oldinexact = get_inexact();
02318 #endif
02319 dval(rv) = tens[k - 9] * dval(rv) + z;
02320 }
02321 bd0 = bb = bd = bs = delta = 0;
02322 if (nd <= DBL_DIG
02323 #ifndef RND_PRODQUOT
02324 #ifndef Honor_FLT_ROUNDS
02325 && Flt_Rounds == 1
02326 #endif
02327 #endif
02328 ) {
02329 if (!e)
02330 goto ret;
02331 if (e > 0) {
02332 if (e <= Ten_pmax) {
02333 #ifdef VAX
02334 goto vax_ovfl_check;
02335 #else
02336 #ifdef Honor_FLT_ROUNDS
02337
02338 if (sign) {
02339 dval(rv) = -dval(rv);
02340 sign = 0;
02341 }
02342 #endif
02343 rounded_product(dval(rv), tens[e]);
02344 goto ret;
02345 #endif
02346 }
02347 i = DBL_DIG - nd;
02348 if (e <= Ten_pmax + i) {
02349
02350
02351
02352 #ifdef Honor_FLT_ROUNDS
02353
02354 if (sign) {
02355 dval(rv) = -dval(rv);
02356 sign = 0;
02357 }
02358 #endif
02359 e -= i;
02360 dval(rv) *= tens[i];
02361 #ifdef VAX
02362
02363
02364
02365 vax_ovfl_check:
02366 word0(rv) -= P*Exp_msk1;
02367 rounded_product(dval(rv), tens[e]);
02368 if ((word0(rv) & Exp_mask)
02369 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
02370 goto ovfl;
02371 word0(rv) += P*Exp_msk1;
02372 #else
02373 rounded_product(dval(rv), tens[e]);
02374 #endif
02375 goto ret;
02376 }
02377 }
02378 #ifndef Inaccurate_Divide
02379 else if (e >= -Ten_pmax) {
02380 #ifdef Honor_FLT_ROUNDS
02381
02382 if (sign) {
02383 dval(rv) = -dval(rv);
02384 sign = 0;
02385 }
02386 #endif
02387 rounded_quotient(dval(rv), tens[-e]);
02388 goto ret;
02389 }
02390 #endif
02391 }
02392 e1 += nd - k;
02393
02394 #ifdef IEEE_Arith
02395 #ifdef SET_INEXACT
02396 inexact = 1;
02397 if (k <= DBL_DIG)
02398 oldinexact = get_inexact();
02399 #endif
02400 #ifdef Avoid_Underflow
02401 scale = 0;
02402 #endif
02403 #ifdef Honor_FLT_ROUNDS
02404 if ((rounding = Flt_Rounds) >= 2) {
02405 if (sign)
02406 rounding = rounding == 2 ? 0 : 2;
02407 else
02408 if (rounding != 2)
02409 rounding = 0;
02410 }
02411 #endif
02412 #endif
02413
02414
02415
02416 if (e1 > 0) {
02417 if ((i = e1 & 15) != 0)
02418 dval(rv) *= tens[i];
02419 if (e1 &= ~15) {
02420 if (e1 > DBL_MAX_10_EXP) {
02421 ovfl:
02422 #ifndef NO_ERRNO
02423 errno = ERANGE;
02424 #endif
02425
02426 #ifdef IEEE_Arith
02427 #ifdef Honor_FLT_ROUNDS
02428 switch (rounding) {
02429 case 0:
02430 case 3:
02431 word0(rv) = Big0;
02432 word1(rv) = Big1;
02433 break;
02434 default:
02435 word0(rv) = Exp_mask;
02436 word1(rv) = 0;
02437 }
02438 #else
02439 word0(rv) = Exp_mask;
02440 word1(rv) = 0;
02441 #endif
02442 #ifdef SET_INEXACT
02443
02444 dval(rv0) = 1e300;
02445 dval(rv0) *= dval(rv0);
02446 #endif
02447 #else
02448 word0(rv) = Big0;
02449 word1(rv) = Big1;
02450 #endif
02451 if (bd0)
02452 goto retfree;
02453 goto ret;
02454 }
02455 e1 >>= 4;
02456 for (j = 0; e1 > 1; j++, e1 >>= 1)
02457 if (e1 & 1)
02458 dval(rv) *= bigtens[j];
02459
02460 word0(rv) -= P*Exp_msk1;
02461 dval(rv) *= bigtens[j];
02462 if ((z = word0(rv) & Exp_mask)
02463 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
02464 goto ovfl;
02465 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
02466
02467
02468 word0(rv) = Big0;
02469 word1(rv) = Big1;
02470 }
02471 else
02472 word0(rv) += P*Exp_msk1;
02473 }
02474 }
02475 else if (e1 < 0) {
02476 e1 = -e1;
02477 if ((i = e1 & 15) != 0)
02478 dval(rv) /= tens[i];
02479 if (e1 >>= 4) {
02480 if (e1 >= 1 << n_bigtens)
02481 goto undfl;
02482 #ifdef Avoid_Underflow
02483 if (e1 & Scale_Bit)
02484 scale = 2*P;
02485 for (j = 0; e1 > 0; j++, e1 >>= 1)
02486 if (e1 & 1)
02487 dval(rv) *= tinytens[j];
02488 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
02489 >> Exp_shift)) > 0) {
02490
02491 if (j >= 32) {
02492 word1(rv) = 0;
02493 if (j >= 53)
02494 word0(rv) = (P+2)*Exp_msk1;
02495 else
02496 word0(rv) &= 0xffffffff << (j-32);
02497 }
02498 else
02499 word1(rv) &= 0xffffffff << j;
02500 }
02501 #else
02502 for (j = 0; e1 > 1; j++, e1 >>= 1)
02503 if (e1 & 1)
02504 dval(rv) *= tinytens[j];
02505
02506 dval(rv0) = dval(rv);
02507 dval(rv) *= tinytens[j];
02508 if (!dval(rv)) {
02509 dval(rv) = 2.*dval(rv0);
02510 dval(rv) *= tinytens[j];
02511 #endif
02512 if (!dval(rv)) {
02513 undfl:
02514 dval(rv) = 0.;
02515 #ifndef NO_ERRNO
02516 errno = ERANGE;
02517 #endif
02518 if (bd0)
02519 goto retfree;
02520 goto ret;
02521 }
02522 #ifndef Avoid_Underflow
02523 word0(rv) = Tiny0;
02524 word1(rv) = Tiny1;
02525
02526
02527
02528 }
02529 #endif
02530 }
02531 }
02532
02533
02534
02535
02536
02537 bd0 = s2b(s0, nd0, nd, y);
02538
02539 for (;;) {
02540 bd = Balloc(bd0->k);
02541 Bcopy(bd, bd0);
02542 bb = d2b(dval(rv), &bbe, &bbbits);
02543 bs = i2b(1);
02544
02545 if (e >= 0) {
02546 bb2 = bb5 = 0;
02547 bd2 = bd5 = e;
02548 }
02549 else {
02550 bb2 = bb5 = -e;
02551 bd2 = bd5 = 0;
02552 }
02553 if (bbe >= 0)
02554 bb2 += bbe;
02555 else
02556 bd2 -= bbe;
02557 bs2 = bb2;
02558 #ifdef Honor_FLT_ROUNDS
02559 if (rounding != 1)
02560 bs2++;
02561 #endif
02562 #ifdef Avoid_Underflow
02563 j = bbe - scale;
02564 i = j + bbbits - 1;
02565 if (i < Emin)
02566 j += P - Emin;
02567 else
02568 j = P + 1 - bbbits;
02569 #else
02570 #ifdef Sudden_Underflow
02571 #ifdef IBM
02572 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
02573 #else
02574 j = P + 1 - bbbits;
02575 #endif
02576 #else
02577 j = bbe;
02578 i = j + bbbits - 1;
02579 if (i < Emin)
02580 j += P - Emin;
02581 else
02582 j = P + 1 - bbbits;
02583 #endif
02584 #endif
02585 bb2 += j;
02586 bd2 += j;
02587 #ifdef Avoid_Underflow
02588 bd2 += scale;
02589 #endif
02590 i = bb2 < bd2 ? bb2 : bd2;
02591 if (i > bs2)
02592 i = bs2;
02593 if (i > 0) {
02594 bb2 -= i;
02595 bd2 -= i;
02596 bs2 -= i;
02597 }
02598 if (bb5 > 0) {
02599 bs = pow5mult(bs, bb5);
02600 bb1 = mult(bs, bb);
02601 Bfree(bb);
02602 bb = bb1;
02603 }
02604 if (bb2 > 0)
02605 bb = lshift(bb, bb2);
02606 if (bd5 > 0)
02607 bd = pow5mult(bd, bd5);
02608 if (bd2 > 0)
02609 bd = lshift(bd, bd2);
02610 if (bs2 > 0)
02611 bs = lshift(bs, bs2);
02612 delta = diff(bb, bd);
02613 dsign = delta->sign;
02614 delta->sign = 0;
02615 i = cmp(delta, bs);
02616 #ifdef Honor_FLT_ROUNDS
02617 if (rounding != 1) {
02618 if (i < 0) {
02619
02620 if (!delta->x[0] && delta->wds <= 1) {
02621
02622 #ifdef SET_INEXACT
02623 inexact = 0;
02624 #endif
02625 break;
02626 }
02627 if (rounding) {
02628 if (dsign) {
02629 adj = 1.;
02630 goto apply_adj;
02631 }
02632 }
02633 else if (!dsign) {
02634 adj = -1.;
02635 if (!word1(rv)
02636 && !(word0(rv) & Frac_mask)) {
02637 y = word0(rv) & Exp_mask;
02638 #ifdef Avoid_Underflow
02639 if (!scale || y > 2*P*Exp_msk1)
02640 #else
02641 if (y)
02642 #endif
02643 {
02644 delta = lshift(delta,Log2P);
02645 if (cmp(delta, bs) <= 0)
02646 adj = -0.5;
02647 }
02648 }
02649 apply_adj:
02650 #ifdef Avoid_Underflow
02651 if (scale && (y = word0(rv) & Exp_mask)
02652 <= 2*P*Exp_msk1)
02653 word0(adj) += (2*P+1)*Exp_msk1 - y;
02654 #else
02655 #ifdef Sudden_Underflow
02656 if ((word0(rv) & Exp_mask) <=
02657 P*Exp_msk1) {
02658 word0(rv) += P*Exp_msk1;
02659 dval(rv) += adj*ulp(dval(rv));
02660 word0(rv) -= P*Exp_msk1;
02661 }
02662 else
02663 #endif
02664 #endif
02665 dval(rv) += adj*ulp(dval(rv));
02666 }
02667 break;
02668 }
02669 adj = ratio(delta, bs);
02670 if (adj < 1.)
02671 adj = 1.;
02672 if (adj <= 0x7ffffffe) {
02673
02674 y = adj;
02675 if (y != adj) {
02676 if (!((rounding>>1) ^ dsign))
02677 y++;
02678 adj = y;
02679 }
02680 }
02681 #ifdef Avoid_Underflow
02682 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02683 word0(adj) += (2*P+1)*Exp_msk1 - y;
02684 #else
02685 #ifdef Sudden_Underflow
02686 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02687 word0(rv) += P*Exp_msk1;
02688 adj *= ulp(dval(rv));
02689 if (dsign)
02690 dval(rv) += adj;
02691 else
02692 dval(rv) -= adj;
02693 word0(rv) -= P*Exp_msk1;
02694 goto cont;
02695 }
02696 #endif
02697 #endif
02698 adj *= ulp(dval(rv));
02699 if (dsign)
02700 dval(rv) += adj;
02701 else
02702 dval(rv) -= adj;
02703 goto cont;
02704 }
02705 #endif
02706
02707 if (i < 0) {
02708
02709
02710
02711 if (dsign || word1(rv) || word0(rv) & Bndry_mask
02712 #ifdef IEEE_Arith
02713 #ifdef Avoid_Underflow
02714 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
02715 #else
02716 || (word0(rv) & Exp_mask) <= Exp_msk1
02717 #endif
02718 #endif
02719 ) {
02720 #ifdef SET_INEXACT
02721 if (!delta->x[0] && delta->wds <= 1)
02722 inexact = 0;
02723 #endif
02724 break;
02725 }
02726 if (!delta->x[0] && delta->wds <= 1) {
02727
02728 #ifdef SET_INEXACT
02729 inexact = 0;
02730 #endif
02731 break;
02732 }
02733 delta = lshift(delta,Log2P);
02734 if (cmp(delta, bs) > 0)
02735 goto drop_down;
02736 break;
02737 }
02738 if (i == 0) {
02739
02740 if (dsign) {
02741 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
02742 && word1(rv) == (
02743 #ifdef Avoid_Underflow
02744 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02745 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
02746 #endif
02747 0xffffffff)) {
02748
02749 word0(rv) = (word0(rv) & Exp_mask)
02750 + Exp_msk1
02751 #ifdef IBM
02752 | Exp_msk1 >> 4
02753 #endif
02754 ;
02755 word1(rv) = 0;
02756 #ifdef Avoid_Underflow
02757 dsign = 0;
02758 #endif
02759 break;
02760 }
02761 }
02762 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
02763 drop_down:
02764
02765 #ifdef Sudden_Underflow
02766 L = word0(rv) & Exp_mask;
02767 #ifdef IBM
02768 if (L < Exp_msk1)
02769 #else
02770 #ifdef Avoid_Underflow
02771 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
02772 #else
02773 if (L <= Exp_msk1)
02774 #endif
02775 #endif
02776 goto undfl;
02777 L -= Exp_msk1;
02778 #else
02779 #ifdef Avoid_Underflow
02780 if (scale) {
02781 L = word0(rv) & Exp_mask;
02782 if (L <= (2*P+1)*Exp_msk1) {
02783 if (L > (P+2)*Exp_msk1)
02784
02785
02786 break;
02787
02788 goto undfl;
02789 }
02790 }
02791 #endif
02792 L = (word0(rv) & Exp_mask) - Exp_msk1;
02793 #endif
02794 word0(rv) = L | Bndry_mask1;
02795 word1(rv) = 0xffffffff;
02796 #ifdef IBM
02797 goto cont;
02798 #else
02799 break;
02800 #endif
02801 }
02802 #ifndef ROUND_BIASED
02803 if (!(word1(rv) & LSB))
02804 break;
02805 #endif
02806 if (dsign)
02807 dval(rv) += ulp(dval(rv));
02808 #ifndef ROUND_BIASED
02809 else {
02810 dval(rv) -= ulp(dval(rv));
02811 #ifndef Sudden_Underflow
02812 if (!dval(rv))
02813 goto undfl;
02814 #endif
02815 }
02816 #ifdef Avoid_Underflow
02817 dsign = 1 - dsign;
02818 #endif
02819 #endif
02820 break;
02821 }
02822 if ((aadj = ratio(delta, bs)) <= 2.) {
02823 if (dsign)
02824 aadj = dval(aadj1) = 1.;
02825 else if (word1(rv) || word0(rv) & Bndry_mask) {
02826 #ifndef Sudden_Underflow
02827 if (word1(rv) == Tiny1 && !word0(rv))
02828 goto undfl;
02829 #endif
02830 aadj = 1.;
02831 dval(aadj1) = -1.;
02832 }
02833 else {
02834
02835
02836
02837 if (aadj < 2./FLT_RADIX)
02838 aadj = 1./FLT_RADIX;
02839 else
02840 aadj *= 0.5;
02841 dval(aadj1) = -aadj;
02842 }
02843 }
02844 else {
02845 aadj *= 0.5;
02846 dval(aadj1) = dsign ? aadj : -aadj;
02847 #ifdef Check_FLT_ROUNDS
02848 switch (Rounding) {
02849 case 2:
02850 dval(aadj1) -= 0.5;
02851 break;
02852 case 0:
02853 case 3:
02854 dval(aadj1) += 0.5;
02855 }
02856 #else
02857 if (Flt_Rounds == 0)
02858 dval(aadj1) += 0.5;
02859 #endif
02860 }
02861 y = word0(rv) & Exp_mask;
02862
02863
02864
02865 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
02866 dval(rv0) = dval(rv);
02867 word0(rv) -= P*Exp_msk1;
02868 adj = dval(aadj1) * ulp(dval(rv));
02869 dval(rv) += adj;
02870 if ((word0(rv) & Exp_mask) >=
02871 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
02872 if (word0(rv0) == Big0 && word1(rv0) == Big1)
02873 goto ovfl;
02874 word0(rv) = Big0;
02875 word1(rv) = Big1;
02876 goto cont;
02877 }
02878 else
02879 word0(rv) += P*Exp_msk1;
02880 }
02881 else {
02882 #ifdef Avoid_Underflow
02883 if (scale && y <= 2*P*Exp_msk1) {
02884 if (aadj <= 0x7fffffff) {
02885 if ((z = (int)aadj) <= 0)
02886 z = 1;
02887 aadj = z;
02888 dval(aadj1) = dsign ? aadj : -aadj;
02889 }
02890 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
02891 }
02892 adj = dval(aadj1) * ulp(dval(rv));
02893 dval(rv) += adj;
02894 #else
02895 #ifdef Sudden_Underflow
02896 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02897 dval(rv0) = dval(rv);
02898 word0(rv) += P*Exp_msk1;
02899 adj = dval(aadj1) * ulp(dval(rv));
02900 dval(rv) += adj;
02901 #ifdef IBM
02902 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
02903 #else
02904 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
02905 #endif
02906 {
02907 if (word0(rv0) == Tiny0 && word1(rv0) == Tiny1)
02908 goto undfl;
02909 word0(rv) = Tiny0;
02910 word1(rv) = Tiny1;
02911 goto cont;
02912 }
02913 else
02914 word0(rv) -= P*Exp_msk1;
02915 }
02916 else {
02917 adj = dval(aadj1) * ulp(dval(rv));
02918 dval(rv) += adj;
02919 }
02920 #else
02921
02922
02923
02924
02925
02926
02927
02928 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
02929 dval(aadj1) = (double)(int)(aadj + 0.5);
02930 if (!dsign)
02931 dval(aadj1) = -dval(aadj1);
02932 }
02933 adj = dval(aadj1) * ulp(dval(rv));
02934 dval(rv) += adj;
02935 #endif
02936 #endif
02937 }
02938 z = word0(rv) & Exp_mask;
02939 #ifndef SET_INEXACT
02940 #ifdef Avoid_Underflow
02941 if (!scale)
02942 #endif
02943 if (y == z) {
02944
02945 L = (Long)aadj;
02946 aadj -= L;
02947
02948 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
02949 if (aadj < .4999999 || aadj > .5000001)
02950 break;
02951 }
02952 else if (aadj < .4999999/FLT_RADIX)
02953 break;
02954 }
02955 #endif
02956 cont:
02957 Bfree(bb);
02958 Bfree(bd);
02959 Bfree(bs);
02960 Bfree(delta);
02961 }
02962 #ifdef SET_INEXACT
02963 if (inexact) {
02964 if (!oldinexact) {
02965 word0(rv0) = Exp_1 + (70 << Exp_shift);
02966 word1(rv0) = 0;
02967 dval(rv0) += 1.;
02968 }
02969 }
02970 else if (!oldinexact)
02971 clear_inexact();
02972 #endif
02973 #ifdef Avoid_Underflow
02974 if (scale) {
02975 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
02976 word1(rv0) = 0;
02977 dval(rv) *= dval(rv0);
02978 #ifndef NO_ERRNO
02979
02980 if (word0(rv) == 0 && word1(rv) == 0)
02981 errno = ERANGE;
02982 #endif
02983 }
02984 #endif
02985 #ifdef SET_INEXACT
02986 if (inexact && !(word0(rv) & Exp_mask)) {
02987
02988 dval(rv0) = 1e-300;
02989 dval(rv0) *= dval(rv0);
02990 }
02991 #endif
02992 retfree:
02993 Bfree(bb);
02994 Bfree(bd);
02995 Bfree(bs);
02996 Bfree(bd0);
02997 Bfree(delta);
02998 ret:
02999 if (se)
03000 *se = (char *)s;
03001 return sign ? -dval(rv) : dval(rv);
03002 }
03003
03004 static int
03005 quorem(Bigint *b, Bigint *S)
03006 {
03007 int n;
03008 ULong *bx, *bxe, q, *sx, *sxe;
03009 #ifdef ULLong
03010 ULLong borrow, carry, y, ys;
03011 #else
03012 ULong borrow, carry, y, ys;
03013 #ifdef Pack_32
03014 ULong si, z, zs;
03015 #endif
03016 #endif
03017
03018 n = S->wds;
03019 #ifdef DEBUG
03020 if (b->wds > n)
03021 Bug("oversize b in quorem");
03022 #endif
03023 if (b->wds < n)
03024 return 0;
03025 sx = S->x;
03026 sxe = sx + --n;
03027 bx = b->x;
03028 bxe = bx + n;
03029 q = *bxe / (*sxe + 1);
03030 #ifdef DEBUG
03031 if (q > 9)
03032 Bug("oversized quotient in quorem");
03033 #endif
03034 if (q) {
03035 borrow = 0;
03036 carry = 0;
03037 do {
03038 #ifdef ULLong
03039 ys = *sx++ * (ULLong)q + carry;
03040 carry = ys >> 32;
03041 y = *bx - (ys & FFFFFFFF) - borrow;
03042 borrow = y >> 32 & (ULong)1;
03043 *bx++ = (ULong)(y & FFFFFFFF);
03044 #else
03045 #ifdef Pack_32
03046 si = *sx++;
03047 ys = (si & 0xffff) * q + carry;
03048 zs = (si >> 16) * q + (ys >> 16);
03049 carry = zs >> 16;
03050 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
03051 borrow = (y & 0x10000) >> 16;
03052 z = (*bx >> 16) - (zs & 0xffff) - borrow;
03053 borrow = (z & 0x10000) >> 16;
03054 Storeinc(bx, z, y);
03055 #else
03056 ys = *sx++ * q + carry;
03057 carry = ys >> 16;
03058 y = *bx - (ys & 0xffff) - borrow;
03059 borrow = (y & 0x10000) >> 16;
03060 *bx++ = y & 0xffff;
03061 #endif
03062 #endif
03063 } while (sx <= sxe);
03064 if (!*bxe) {
03065 bx = b->x;
03066 while (--bxe > bx && !*bxe)
03067 --n;
03068 b->wds = n;
03069 }
03070 }
03071 if (cmp(b, S) >= 0) {
03072 q++;
03073 borrow = 0;
03074 carry = 0;
03075 bx = b->x;
03076 sx = S->x;
03077 do {
03078 #ifdef ULLong
03079 ys = *sx++ + carry;
03080 carry = ys >> 32;
03081 y = *bx - (ys & FFFFFFFF) - borrow;
03082 borrow = y >> 32 & (ULong)1;
03083 *bx++ = (ULong)(y & FFFFFFFF);
03084 #else
03085 #ifdef Pack_32
03086 si = *sx++;
03087 ys = (si & 0xffff) + carry;
03088 zs = (si >> 16) + (ys >> 16);
03089 carry = zs >> 16;
03090 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
03091 borrow = (y & 0x10000) >> 16;
03092 z = (*bx >> 16) - (zs & 0xffff) - borrow;
03093 borrow = (z & 0x10000) >> 16;
03094 Storeinc(bx, z, y);
03095 #else
03096 ys = *sx++ + carry;
03097 carry = ys >> 16;
03098 y = *bx - (ys & 0xffff) - borrow;
03099 borrow = (y & 0x10000) >> 16;
03100 *bx++ = y & 0xffff;
03101 #endif
03102 #endif
03103 } while (sx <= sxe);
03104 bx = b->x;
03105 bxe = bx + n;
03106 if (!*bxe) {
03107 while (--bxe > bx && !*bxe)
03108 --n;
03109 b->wds = n;
03110 }
03111 }
03112 return q;
03113 }
03114
03115 #ifndef MULTIPLE_THREADS
03116 static char *dtoa_result;
03117 #endif
03118
03119 #ifndef MULTIPLE_THREADS
03120 static char *
03121 rv_alloc(int i)
03122 {
03123 return dtoa_result = xmalloc(i);
03124 }
03125 #else
03126 #define rv_alloc(i) xmalloc(i)
03127 #endif
03128
03129 static char *
03130 nrv_alloc(const char *s, char **rve, size_t n)
03131 {
03132 char *rv, *t;
03133
03134 t = rv = rv_alloc(n);
03135 while ((*t = *s++) != 0) t++;
03136 if (rve)
03137 *rve = t;
03138 return rv;
03139 }
03140
03141 #define rv_strdup(s, rve) nrv_alloc(s, rve, strlen(s)+1)
03142
03143 #ifndef MULTIPLE_THREADS
03144
03145
03146
03147
03148
03149
03150 static void
03151 freedtoa(char *s)
03152 {
03153 xfree(s);
03154 }
03155 #endif
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
03181
03182
03183
03184
03185
03186
03187
03188
03189
03190
03191 char *
03192 ruby_dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
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
03218
03219
03220
03221
03222
03223
03224
03225
03226
03227
03228 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
03229 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
03230 spec_case, try_quick;
03231 Long L;
03232 #ifndef Sudden_Underflow
03233 int denorm;
03234 ULong x;
03235 #endif
03236 Bigint *b, *b1, *delta, *mlo = 0, *mhi = 0, *S;
03237 double ds;
03238 double_u d, d2, eps;
03239 char *s, *s0;
03240 #ifdef Honor_FLT_ROUNDS
03241 int rounding;
03242 #endif
03243 #ifdef SET_INEXACT
03244 int inexact, oldinexact;
03245 #endif
03246
03247 dval(d) = d_;
03248
03249 #ifndef MULTIPLE_THREADS
03250 if (dtoa_result) {
03251 freedtoa(dtoa_result);
03252 dtoa_result = 0;
03253 }
03254 #endif
03255
03256 if (word0(d) & Sign_bit) {
03257
03258 *sign = 1;
03259 word0(d) &= ~Sign_bit;
03260 }
03261 else
03262 *sign = 0;
03263
03264 #if defined(IEEE_Arith) + defined(VAX)
03265 #ifdef IEEE_Arith
03266 if ((word0(d) & Exp_mask) == Exp_mask)
03267 #else
03268 if (word0(d) == 0x8000)
03269 #endif
03270 {
03271
03272 *decpt = 9999;
03273 #ifdef IEEE_Arith
03274 if (!word1(d) && !(word0(d) & 0xfffff))
03275 return rv_strdup("Infinity", rve);
03276 #endif
03277 return rv_strdup("NaN", rve);
03278 }
03279 #endif
03280 #ifdef IBM
03281 dval(d) += 0;
03282 #endif
03283 if (!dval(d)) {
03284 *decpt = 1;
03285 return rv_strdup("0", rve);
03286 }
03287
03288 #ifdef SET_INEXACT
03289 try_quick = oldinexact = get_inexact();
03290 inexact = 1;
03291 #endif
03292 #ifdef Honor_FLT_ROUNDS
03293 if ((rounding = Flt_Rounds) >= 2) {
03294 if (*sign)
03295 rounding = rounding == 2 ? 0 : 2;
03296 else
03297 if (rounding != 2)
03298 rounding = 0;
03299 }
03300 #endif
03301
03302 b = d2b(dval(d), &be, &bbits);
03303 #ifdef Sudden_Underflow
03304 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
03305 #else
03306 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) {
03307 #endif
03308 dval(d2) = dval(d);
03309 word0(d2) &= Frac_mask1;
03310 word0(d2) |= Exp_11;
03311 #ifdef IBM
03312 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
03313 dval(d2) /= 1 << j;
03314 #endif
03315
03316
03317
03318
03319
03320
03321
03322
03323
03324
03325
03326
03327
03328
03329
03330
03331
03332
03333
03334
03335
03336
03337
03338 i -= Bias;
03339 #ifdef IBM
03340 i <<= 2;
03341 i += j;
03342 #endif
03343 #ifndef Sudden_Underflow
03344 denorm = 0;
03345 }
03346 else {
03347
03348
03349 i = bbits + be + (Bias + (P-1) - 1);
03350 x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
03351 : word1(d) << (32 - i);
03352 dval(d2) = x;
03353 word0(d2) -= 31*Exp_msk1;
03354 i -= (Bias + (P-1) - 1) + 1;
03355 denorm = 1;
03356 }
03357 #endif
03358 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
03359 k = (int)ds;
03360 if (ds < 0. && ds != k)
03361 k--;
03362 k_check = 1;
03363 if (k >= 0 && k <= Ten_pmax) {
03364 if (dval(d) < tens[k])
03365 k--;
03366 k_check = 0;
03367 }
03368 j = bbits - i - 1;
03369 if (j >= 0) {
03370 b2 = 0;
03371 s2 = j;
03372 }
03373 else {
03374 b2 = -j;
03375 s2 = 0;
03376 }
03377 if (k >= 0) {
03378 b5 = 0;
03379 s5 = k;
03380 s2 += k;
03381 }
03382 else {
03383 b2 -= k;
03384 b5 = -k;
03385 s5 = 0;
03386 }
03387 if (mode < 0 || mode > 9)
03388 mode = 0;
03389
03390 #ifndef SET_INEXACT
03391 #ifdef Check_FLT_ROUNDS
03392 try_quick = Rounding == 1;
03393 #else
03394 try_quick = 1;
03395 #endif
03396 #endif
03397
03398 if (mode > 5) {
03399 mode -= 4;
03400 try_quick = 0;
03401 }
03402 leftright = 1;
03403 ilim = ilim1 = -1;
03404 switch (mode) {
03405 case 0:
03406 case 1:
03407 i = 18;
03408 ndigits = 0;
03409 break;
03410 case 2:
03411 leftright = 0;
03412
03413 case 4:
03414 if (ndigits <= 0)
03415 ndigits = 1;
03416 ilim = ilim1 = i = ndigits;
03417 break;
03418 case 3:
03419 leftright = 0;
03420
03421 case 5:
03422 i = ndigits + k + 1;
03423 ilim = i;
03424 ilim1 = i - 1;
03425 if (i <= 0)
03426 i = 1;
03427 }
03428 s = s0 = rv_alloc(i+1);
03429
03430 #ifdef Honor_FLT_ROUNDS
03431 if (mode > 1 && rounding != 1)
03432 leftright = 0;
03433 #endif
03434
03435 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
03436
03437
03438
03439 i = 0;
03440 dval(d2) = dval(d);
03441 k0 = k;
03442 ilim0 = ilim;
03443 ieps = 2;
03444 if (k > 0) {
03445 ds = tens[k&0xf];
03446 j = k >> 4;
03447 if (j & Bletch) {
03448
03449 j &= Bletch - 1;
03450 dval(d) /= bigtens[n_bigtens-1];
03451 ieps++;
03452 }
03453 for (; j; j >>= 1, i++)
03454 if (j & 1) {
03455 ieps++;
03456 ds *= bigtens[i];
03457 }
03458 dval(d) /= ds;
03459 }
03460 else if ((j1 = -k) != 0) {
03461 dval(d) *= tens[j1 & 0xf];
03462 for (j = j1 >> 4; j; j >>= 1, i++)
03463 if (j & 1) {
03464 ieps++;
03465 dval(d) *= bigtens[i];
03466 }
03467 }
03468 if (k_check && dval(d) < 1. && ilim > 0) {
03469 if (ilim1 <= 0)
03470 goto fast_failed;
03471 ilim = ilim1;
03472 k--;
03473 dval(d) *= 10.;
03474 ieps++;
03475 }
03476 dval(eps) = ieps*dval(d) + 7.;
03477 word0(eps) -= (P-1)*Exp_msk1;
03478 if (ilim == 0) {
03479 S = mhi = 0;
03480 dval(d) -= 5.;
03481 if (dval(d) > dval(eps))
03482 goto one_digit;
03483 if (dval(d) < -dval(eps))
03484 goto no_digits;
03485 goto fast_failed;
03486 }
03487 #ifndef No_leftright
03488 if (leftright) {
03489
03490
03491
03492 dval(eps) = 0.5/tens[ilim-1] - dval(eps);
03493 for (i = 0;;) {
03494 L = (int)dval(d);
03495 dval(d) -= L;
03496 *s++ = '0' + (int)L;
03497 if (dval(d) < dval(eps))
03498 goto ret1;
03499 if (1. - dval(d) < dval(eps))
03500 goto bump_up;
03501 if (++i >= ilim)
03502 break;
03503 dval(eps) *= 10.;
03504 dval(d) *= 10.;
03505 }
03506 }
03507 else {
03508 #endif
03509
03510 dval(eps) *= tens[ilim-1];
03511 for (i = 1;; i++, dval(d) *= 10.) {
03512 L = (Long)(dval(d));
03513 if (!(dval(d) -= L))
03514 ilim = i;
03515 *s++ = '0' + (int)L;
03516 if (i == ilim) {
03517 if (dval(d) > 0.5 + dval(eps))
03518 goto bump_up;
03519 else if (dval(d) < 0.5 - dval(eps)) {
03520 while (*--s == '0') ;
03521 s++;
03522 goto ret1;
03523 }
03524 break;
03525 }
03526 }
03527 #ifndef No_leftright
03528 }
03529 #endif
03530 fast_failed:
03531 s = s0;
03532 dval(d) = dval(d2);
03533 k = k0;
03534 ilim = ilim0;
03535 }
03536
03537
03538
03539 if (be >= 0 && k <= Int_max) {
03540
03541 ds = tens[k];
03542 if (ndigits < 0 && ilim <= 0) {
03543 S = mhi = 0;
03544 if (ilim < 0 || dval(d) <= 5*ds)
03545 goto no_digits;
03546 goto one_digit;
03547 }
03548 for (i = 1;; i++, dval(d) *= 10.) {
03549 L = (Long)(dval(d) / ds);
03550 dval(d) -= L*ds;
03551 #ifdef Check_FLT_ROUNDS
03552
03553 if (dval(d) < 0) {
03554 L--;
03555 dval(d) += ds;
03556 }
03557 #endif
03558 *s++ = '0' + (int)L;
03559 if (!dval(d)) {
03560 #ifdef SET_INEXACT
03561 inexact = 0;
03562 #endif
03563 break;
03564 }
03565 if (i == ilim) {
03566 #ifdef Honor_FLT_ROUNDS
03567 if (mode > 1)
03568 switch (rounding) {
03569 case 0: goto ret1;
03570 case 2: goto bump_up;
03571 }
03572 #endif
03573 dval(d) += dval(d);
03574 if (dval(d) > ds || (dval(d) == ds && (L & 1))) {
03575 bump_up:
03576 while (*--s == '9')
03577 if (s == s0) {
03578 k++;
03579 *s = '0';
03580 break;
03581 }
03582 ++*s++;
03583 }
03584 break;
03585 }
03586 }
03587 goto ret1;
03588 }
03589
03590 m2 = b2;
03591 m5 = b5;
03592 if (leftright) {
03593 i =
03594 #ifndef Sudden_Underflow
03595 denorm ? be + (Bias + (P-1) - 1 + 1) :
03596 #endif
03597 #ifdef IBM
03598 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
03599 #else
03600 1 + P - bbits;
03601 #endif
03602 b2 += i;
03603 s2 += i;
03604 mhi = i2b(1);
03605 }
03606 if (m2 > 0 && s2 > 0) {
03607 i = m2 < s2 ? m2 : s2;
03608 b2 -= i;
03609 m2 -= i;
03610 s2 -= i;
03611 }
03612 if (b5 > 0) {
03613 if (leftright) {
03614 if (m5 > 0) {
03615 mhi = pow5mult(mhi, m5);
03616 b1 = mult(mhi, b);
03617 Bfree(b);
03618 b = b1;
03619 }
03620 if ((j = b5 - m5) != 0)
03621 b = pow5mult(b, j);
03622 }
03623 else
03624 b = pow5mult(b, b5);
03625 }
03626 S = i2b(1);
03627 if (s5 > 0)
03628 S = pow5mult(S, s5);
03629
03630
03631
03632 spec_case = 0;
03633 if ((mode < 2 || leftright)
03634 #ifdef Honor_FLT_ROUNDS
03635 && rounding == 1
03636 #endif
03637 ) {
03638 if (!word1(d) && !(word0(d) & Bndry_mask)
03639 #ifndef Sudden_Underflow
03640 && word0(d) & (Exp_mask & ~Exp_msk1)
03641 #endif
03642 ) {
03643
03644 b2 += Log2P;
03645 s2 += Log2P;
03646 spec_case = 1;
03647 }
03648 }
03649
03650
03651
03652
03653
03654
03655
03656
03657 #ifdef Pack_32
03658 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0)
03659 i = 32 - i;
03660 #else
03661 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) != 0)
03662 i = 16 - i;
03663 #endif
03664 if (i > 4) {
03665 i -= 4;
03666 b2 += i;
03667 m2 += i;
03668 s2 += i;
03669 }
03670 else if (i < 4) {
03671 i += 28;
03672 b2 += i;
03673 m2 += i;
03674 s2 += i;
03675 }
03676 if (b2 > 0)
03677 b = lshift(b, b2);
03678 if (s2 > 0)
03679 S = lshift(S, s2);
03680 if (k_check) {
03681 if (cmp(b,S) < 0) {
03682 k--;
03683 b = multadd(b, 10, 0);
03684 if (leftright)
03685 mhi = multadd(mhi, 10, 0);
03686 ilim = ilim1;
03687 }
03688 }
03689 if (ilim <= 0 && (mode == 3 || mode == 5)) {
03690 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
03691
03692 no_digits:
03693 k = -1 - ndigits;
03694 goto ret;
03695 }
03696 one_digit:
03697 *s++ = '1';
03698 k++;
03699 goto ret;
03700 }
03701 if (leftright) {
03702 if (m2 > 0)
03703 mhi = lshift(mhi, m2);
03704
03705
03706
03707
03708
03709 mlo = mhi;
03710 if (spec_case) {
03711 mhi = Balloc(mhi->k);
03712 Bcopy(mhi, mlo);
03713 mhi = lshift(mhi, Log2P);
03714 }
03715
03716 for (i = 1;;i++) {
03717 dig = quorem(b,S) + '0';
03718
03719
03720
03721 j = cmp(b, mlo);
03722 delta = diff(S, mhi);
03723 j1 = delta->sign ? 1 : cmp(b, delta);
03724 Bfree(delta);
03725 #ifndef ROUND_BIASED
03726 if (j1 == 0 && mode != 1 && !(word1(d) & 1)
03727 #ifdef Honor_FLT_ROUNDS
03728 && rounding >= 1
03729 #endif
03730 ) {
03731 if (dig == '9')
03732 goto round_9_up;
03733 if (j > 0)
03734 dig++;
03735 #ifdef SET_INEXACT
03736 else if (!b->x[0] && b->wds <= 1)
03737 inexact = 0;
03738 #endif
03739 *s++ = dig;
03740 goto ret;
03741 }
03742 #endif
03743 if (j < 0 || (j == 0 && mode != 1
03744 #ifndef ROUND_BIASED
03745 && !(word1(d) & 1)
03746 #endif
03747 )) {
03748 if (!b->x[0] && b->wds <= 1) {
03749 #ifdef SET_INEXACT
03750 inexact = 0;
03751 #endif
03752 goto accept_dig;
03753 }
03754 #ifdef Honor_FLT_ROUNDS
03755 if (mode > 1)
03756 switch (rounding) {
03757 case 0: goto accept_dig;
03758 case 2: goto keep_dig;
03759 }
03760 #endif
03761 if (j1 > 0) {
03762 b = lshift(b, 1);
03763 j1 = cmp(b, S);
03764 if ((j1 > 0 || (j1 == 0 && (dig & 1))) && dig++ == '9')
03765 goto round_9_up;
03766 }
03767 accept_dig:
03768 *s++ = dig;
03769 goto ret;
03770 }
03771 if (j1 > 0) {
03772 #ifdef Honor_FLT_ROUNDS
03773 if (!rounding)
03774 goto accept_dig;
03775 #endif
03776 if (dig == '9') {
03777 round_9_up:
03778 *s++ = '9';
03779 goto roundoff;
03780 }
03781 *s++ = dig + 1;
03782 goto ret;
03783 }
03784 #ifdef Honor_FLT_ROUNDS
03785 keep_dig:
03786 #endif
03787 *s++ = dig;
03788 if (i == ilim)
03789 break;
03790 b = multadd(b, 10, 0);
03791 if (mlo == mhi)
03792 mlo = mhi = multadd(mhi, 10, 0);
03793 else {
03794 mlo = multadd(mlo, 10, 0);
03795 mhi = multadd(mhi, 10, 0);
03796 }
03797 }
03798 }
03799 else
03800 for (i = 1;; i++) {
03801 *s++ = dig = quorem(b,S) + '0';
03802 if (!b->x[0] && b->wds <= 1) {
03803 #ifdef SET_INEXACT
03804 inexact = 0;
03805 #endif
03806 goto ret;
03807 }
03808 if (i >= ilim)
03809 break;
03810 b = multadd(b, 10, 0);
03811 }
03812
03813
03814
03815 #ifdef Honor_FLT_ROUNDS
03816 switch (rounding) {
03817 case 0: goto trimzeros;
03818 case 2: goto roundoff;
03819 }
03820 #endif
03821 b = lshift(b, 1);
03822 j = cmp(b, S);
03823 if (j > 0 || (j == 0 && (dig & 1))) {
03824 roundoff:
03825 while (*--s == '9')
03826 if (s == s0) {
03827 k++;
03828 *s++ = '1';
03829 goto ret;
03830 }
03831 ++*s++;
03832 }
03833 else {
03834 while (*--s == '0') ;
03835 s++;
03836 }
03837 ret:
03838 Bfree(S);
03839 if (mhi) {
03840 if (mlo && mlo != mhi)
03841 Bfree(mlo);
03842 Bfree(mhi);
03843 }
03844 ret1:
03845 #ifdef SET_INEXACT
03846 if (inexact) {
03847 if (!oldinexact) {
03848 word0(d) = Exp_1 + (70 << Exp_shift);
03849 word1(d) = 0;
03850 dval(d) += 1.;
03851 }
03852 }
03853 else if (!oldinexact)
03854 clear_inexact();
03855 #endif
03856 Bfree(b);
03857 *s = 0;
03858 *decpt = k + 1;
03859 if (rve)
03860 *rve = s;
03861 return s0;
03862 }
03863
03864 void
03865 ruby_each_words(const char *str, void (*func)(const char*, int, void*), void *arg)
03866 {
03867 const char *end;
03868 int len;
03869
03870 if (!str) return;
03871 for (; *str; str = end) {
03872 while (ISSPACE(*str) || *str == ',') str++;
03873 if (!*str) break;
03874 end = str;
03875 while (*end && !ISSPACE(*end) && *end != ',') end++;
03876 len = (int)(end - str);
03877 (*func)(str, len, arg);
03878 }
03879 }
03880
03881
03882
03883
03884
03885
03886
03887
03888
03889
03890
03891
03892
03893
03894
03895
03896
03897
03898
03899
03900
03901
03902
03903
03904
03905
03906
03907 #define DBL_MANH_SIZE 20
03908 #define DBL_MANL_SIZE 32
03909 #define INFSTR "Infinity"
03910 #define NANSTR "NaN"
03911 #define DBL_ADJ (DBL_MAX_EXP - 2)
03912 #define SIGFIGS ((DBL_MANT_DIG + 3) / 4 + 1)
03913 #define dexp_get(u) ((int)(word0(u) >> Exp_shift) & ~Exp_msk1)
03914 #define dexp_set(u,v) (word0(u) = (((int)(word0(u)) & ~Exp_mask) | (v << Exp_shift)))
03915 #define dmanh_get(u) ((int)(word0(u) & Frac_mask))
03916 #define dmanl_get(u) ((int)word1(u))
03917
03918
03919
03920
03921
03922
03923
03924
03925
03926
03927
03928
03929
03930
03931
03932
03933
03934
03935
03936
03937
03938
03939
03940
03941
03942
03943 char *
03944 BSD__hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
03945 char **rve)
03946 {
03947 U u;
03948 char *s, *s0;
03949 int bufsize;
03950 uint32_t manh, manl;
03951
03952 u.d = d;
03953 if (word0(u) & Sign_bit) {
03954
03955 *sign = 1;
03956 word0(u) &= ~Sign_bit;
03957 }
03958 else
03959 *sign = 0;
03960
03961 if (isinf(d)) {
03962 *decpt = INT_MAX;
03963 return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
03964 }
03965 else if (isnan(d)) {
03966 *decpt = INT_MAX;
03967 return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
03968 }
03969 else if (d == 0.0) {
03970 *decpt = 1;
03971 return (nrv_alloc("0", rve, 1));
03972 }
03973 else if (dexp_get(u)) {
03974 *decpt = dexp_get(u) - DBL_ADJ;
03975 }
03976 else {
03977 u.d *= 5.363123171977039e+154 ;
03978 *decpt = dexp_get(u) - (514 + DBL_ADJ);
03979 }
03980
03981 if (ndigits == 0)
03982 ndigits = 1;
03983
03984
03985
03986
03987
03988 bufsize = (ndigits > 0) ? ndigits : SIGFIGS;
03989 s0 = rv_alloc(bufsize);
03990
03991
03992 if (SIGFIGS > ndigits && ndigits > 0) {
03993 float redux = 1.0f;
03994 int offset = 4 * ndigits + DBL_MAX_EXP - 4 - DBL_MANT_DIG;
03995 dexp_set(u, offset);
03996 u.d += redux;
03997 u.d -= redux;
03998 *decpt += dexp_get(u) - offset;
03999 }
04000
04001 manh = dmanh_get(u);
04002 manl = dmanl_get(u);
04003 *s0 = '1';
04004 for (s = s0 + 1; s < s0 + bufsize; s++) {
04005 *s = xdigs[(manh >> (DBL_MANH_SIZE - 4)) & 0xf];
04006 manh = (manh << 4) | (manl >> (DBL_MANL_SIZE - 4));
04007 manl <<= 4;
04008 }
04009
04010
04011 if (ndigits < 0) {
04012 for (ndigits = SIGFIGS; s0[ndigits - 1] == '0'; ndigits--)
04013 ;
04014 }
04015
04016 s = s0 + ndigits;
04017 *s = '\0';
04018 if (rve != NULL)
04019 *rve = s;
04020 return (s0);
04021 }
04022
04023 #ifdef __cplusplus
04024 #if 0
04025 {
04026 #endif
04027 }
04028 #endif
04029