• Main Page
  • Modules
  • Data Structures
  • Files
  • File List
  • Globals

util.c

Go to the documentation of this file.
00001 /**********************************************************************
00002 
00003   util.c -
00004 
00005   $Author: mame $
00006   created at: Fri Mar 10 17:22:34 JST 1995
00007 
00008   Copyright (C) 1993-2008 Yukihiro Matsumoto
00009 
00010 **********************************************************************/
00011 
00012 #include "ruby/ruby.h"
00013 
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); /* less than len */
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); /* less than len */
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         /*     0  1  2  3  4  5  6  7  8  9  a  b  c  d  e  f */
00062         /*0*/ -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,-1,
00064         /*2*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
00065         /*3*/  0, 1, 2, 3, 4, 5, 6, 7, 8, 9,-1,-1,-1,-1,-1,-1,
00066         /*4*/ -1,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,
00067         /*5*/ 25,26,27,28,29,30,31,32,33,34,35,-1,-1,-1,-1,-1,
00068         /*6*/ -1,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,
00069         /*7*/ 25,26,27,28,29,30,31,32,33,34,35,-1,-1,-1,-1,-1,
00070         /*8*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
00071         /*9*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
00072         /*a*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
00073         /*b*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
00074         /*c*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
00075         /*d*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
00076         /*e*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
00077         /*f*/ -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  *  Copyright (c) 1993, Intergraph Corporation
00189  *
00190  *  You may distribute under the terms of either the GNU General Public
00191  *  License or the Artistic License, as specified in the perl README file.
00192  *
00193  *  Various Unix compatibility functions and NT specific functions.
00194  *
00195  *  Some of this code was derived from the MSDOS port(s) and the OS/2 port.
00196  *
00197  */
00198 
00199 
00200 /*
00201  * Suffix appending for in-place editing under MS-DOS and OS/2 (and now NT!).
00202  *
00203  * Here are the rules:
00204  *
00205  * Style 0:  Append the suffix exactly as standard perl would do it.
00206  *           If the filesystem groks it, use it.  (HPFS will always
00207  *           grok it.  So will NTFS. FAT will rarely accept it.)
00208  *
00209  * Style 1:  The suffix begins with a '.'.  The extension is replaced.
00210  *           If the name matches the original name, use the fallback method.
00211  *
00212  * Style 2:  The suffix is a single character, not a '.'.  Try to add the
00213  *           suffix to the following places, using the first one that works.
00214  *               [1] Append to extension.
00215  *               [2] Append to filename,
00216  *               [3] Replace end of extension,
00217  *               [4] Replace end of filename.
00218  *           If the name matches the original name, use the fallback method.
00219  *
00220  * Style 3:  Any other case:  Ignore the suffix completely and use the
00221  *           fallback method.
00222  *
00223  * Fallback method:  Change the extension to ".$$$".  If that matches the
00224  *           original name, then change the extension to ".~~~".
00225  *
00226  * If filename is more than 1000 characters long, we die a horrible
00227  * death.  Sorry.
00228  *
00229  * The filename restriction is a cheat so that we can use buf[] to store
00230  * assorted temporary goo.
00231  *
00232  * Examples, assuming style 0 failed.
00233  *
00234  * suffix = ".bak" (style 1)
00235  *                foo.bar => foo.bak
00236  *                foo.bak => foo.$$$    (fallback)
00237  *                foo.$$$ => foo.~~~    (fallback)
00238  *                makefile => makefile.bak
00239  *
00240  * suffix = "~" (style 2)
00241  *                foo.c => foo.c~
00242  *                foo.c~ => foo.c~~
00243  *                foo.c~~ => foo~.c~~
00244  *                foo~.c~~ => foo~~.c~~
00245  *                foo~~~~~.c~~ => foo~~~~~.$$$ (fallback)
00246  *
00247  *                foo.pas => foo~.pas
00248  *                makefile => makefile.~
00249  *                longname.fil => longname.fi~
00250  *                longname.fi~ => longnam~.fi~
00251  *                longnam~.fi~ => longnam~.$$$
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     /* Style 0 */
00285     rb_str_cat(str, suffix, extlen);
00286     if (valid_filename(RSTRING_PTR(str))) return;
00287 
00288     /* Fooey, style 0 failed.  Fix str before continuing. */
00289     rb_str_resize(str, slen);
00290     name = StringValueCStr(str);
00291     ext = ruby_find_extname(name, &len);
00292 
00293     if (*suffix == '.') {        /* Style 1 */
00294         if (ext) {
00295             if (strEQ(ext, suffix)) {
00296                 extlen = sizeof(suffix1) - 1; /* suffix2 must be same length */
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') {  /* Style 2 */
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 { /* Style 3:  Panic */
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     // It doesn't exist, so see if we can open it.
00350     */
00351 
00352     if ((fd = open(s, O_CREAT|O_EXCL, 0666)) >= 0) {
00353         close(fd);
00354         unlink(s);      /* don't leave it laying around */
00355         return 1;
00356     }
00357     else if (errno == EEXIST) {
00358         /* if the file exists, then it's a valid filename! */
00359         return 1;
00360     }
00361     return 0;
00362 }
00363 #endif
00364 
00365 
00366 /* mm.c */
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 /* qs6.c */
00434 /*****************************************************/
00435 /*                                                   */
00436 /*          qs6   (Quick sort function)              */
00437 /*                                                   */
00438 /* by  Tomoyuki Kawamura              1995.4.21      */
00439 /* kawamura@tokuyama.ac.jp                           */
00440 /*****************************************************/
00441 
00442 typedef struct { char *LL, *RR; } stack_node; /* Stack structure for L,l,R,r */
00443 #define PUSH(ll,rr) do { top->LL = (ll); top->RR = (rr); ++top; } while (0)  /* Push L,l,R,r */
00444 #define POP(ll,rr)  do { --top; ll = top->LL; rr = top->RR; } while (0)      /* Pop L,l,R,r */
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;             /* l,r:left,right group   m:median point */
00455   register int t, eq_l, eq_r;           /* eq_l: all items in left group are equal to S */
00456   char *L = base;                       /* left end of current region */
00457   char *R = (char*)base + size*(nel-1); /* right end of current region */
00458   size_t chklim = 63;                   /* threshold of ordering element check */
00459   stack_node stack[32], *top = stack;   /* 32 is enough for 32bit CPU */
00460   int mmkind;
00461   size_t high, low, n;
00462 
00463   if (nel <= 1) return;        /* need not to sort */
00464   mmprepare(base, size);
00465   goto start;
00466 
00467   nxt:
00468   if (stack == top) return;    /* return if stack is empty */
00469   POP(L,R);
00470 
00471   for (;;) {
00472     start:
00473     if (L + size == R) {       /* 2 elements */
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;  /* number of elements */
00479     m = l + size * (n >> 1);    /* calculate median value */
00480 
00481     if (n >= 60) {
00482       register char *m1;
00483       register char *m3;
00484       if (n >= 200) {
00485         n = size*(n>>3); /* number of bytes in splitting 8 */
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); /* number of bytes in splitting 4 */
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) {                           /*3-5-?*/
00506       if ((t = (*cmp)(m,r,d)) < 0) {                         /*3-5-7*/
00507         if (chklim && nel >= chklim) {   /* check if already ascending order */
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;                                    /*3-5-7*/
00514       }
00515       if (t > 0) {
00516         if ((*cmp)(l,r,d) <= 0) {mmswap(m,r); goto loopA;}     /*3-5-4*/
00517         mmrot3(r,m,l); goto loopA;                           /*3-5-2*/
00518       }
00519       goto loopB;                                            /*3-5-5*/
00520     }
00521 
00522     if (t > 0) {                                             /*7-5-?*/
00523       if ((t = (*cmp)(m,r,d)) > 0) {                         /*7-5-3*/
00524         if (chklim && nel >= chklim) {   /* check if already ascending order */
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;}  /* reverse region */
00529           goto nxt;
00530         }
00531         fail2: mmswap(l,r); goto loopA;                      /*7-5-3*/
00532       }
00533       if (t < 0) {
00534         if ((*cmp)(l,r,d) <= 0) {mmswap(l,m); goto loopB;}   /*7-5-8*/
00535         mmrot3(l,m,r); goto loopA;                           /*7-5-6*/
00536       }
00537       mmswap(l,r); goto loopA;                               /*7-5-5*/
00538     }
00539 
00540     if ((t = (*cmp)(m,r,d)) < 0)  {goto loopA;}              /*5-5-7*/
00541     if (t > 0) {mmswap(l,r); goto loopB;}                    /*5-5-3*/
00542 
00543     /* determining splitting type in case 5-5-5 */           /*5-5-5*/
00544     for (;;) {
00545       if ((l += size) == r)      goto nxt;                   /*5-5-5*/
00546       if (l == m) continue;
00547       if ((t = (*cmp)(l,m,d)) > 0) {mmswap(l,r); l = L; goto loopA;}/*575-5*/
00548       if (t < 0)                 {mmswap(L,l); l = L; goto loopB;}  /*535-5*/
00549     }
00550 
00551     loopA: eq_l = 1; eq_r = 1;  /* splitting type A */ /* left <= median < right */
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);    /* swap left and right */
00568     }
00569 
00570     loopB: eq_l = 1; eq_r = 1;  /* splitting type B */ /* left < median <= right */
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);    /* swap left and right */
00587     }
00588 
00589     fin:
00590     if (eq_l == 0)                         /* need to sort left side */
00591       if (eq_r == 0)                       /* need to sort right side */
00592         if (l-L < R-r) {PUSH(r,R); R = l;} /* sort left side first */
00593         else           {PUSH(L,l); L = r;} /* sort right side first */
00594       else R = l;                          /* need to sort left side only */
00595     else if (eq_r == 0) L = r;             /* need to sort right side only */
00596     else goto nxt;                         /* need not to sort both sides */
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  * The author of this software is David M. Gay.
00644  *
00645  * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
00646  *
00647  * Permission to use, copy, modify, and distribute this software for any
00648  * purpose without fee is hereby granted, provided that this entire notice
00649  * is included in all copies of any software which is or includes a copy
00650  * or modification of this software and in all copies of the supporting
00651  * documentation for such software.
00652  *
00653  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
00654  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
00655  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
00656  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
00657  *
00658  ***************************************************************/
00659 
00660 /* Please send bug reports to David M. Gay (dmg at acm dot org,
00661  * with " at " changed at "@" and " dot " changed to ".").      */
00662 
00663 /* On a machine with IEEE extended-precision registers, it is
00664  * necessary to specify double-precision (53-bit) rounding precision
00665  * before invoking strtod or dtoa.  If the machine uses (the equivalent
00666  * of) Intel 80x87 arithmetic, the call
00667  *      _control87(PC_53, MCW_PC);
00668  * does this with many compilers.  Whether this or another call is
00669  * appropriate depends on the compiler; for this to work, it may be
00670  * necessary to #include "float.h" or another system-dependent header
00671  * file.
00672  */
00673 
00674 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
00675  *
00676  * This strtod returns a nearest machine number to the input decimal
00677  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
00678  * broken by the IEEE round-even rule.  Otherwise ties are broken by
00679  * biased rounding (add half and chop).
00680  *
00681  * Inspired loosely by William D. Clinger's paper "How to Read Floating
00682  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
00683  *
00684  * Modifications:
00685  *
00686  *      1. We only require IEEE, IBM, or VAX double-precision
00687  *              arithmetic (not IEEE double-extended).
00688  *      2. We get by with floating-point arithmetic in a case that
00689  *              Clinger missed -- when we're computing d * 10^n
00690  *              for a small integer d and the integer n is not too
00691  *              much larger than 22 (the maximum integer k for which
00692  *              we can represent 10^k exactly), we may be able to
00693  *              compute (d*10^k) * 10^(e-k) with just one roundoff.
00694  *      3. Rather than a bit-at-a-time adjustment of the binary
00695  *              result in the hard case, we use floating-point
00696  *              arithmetic to determine the adjustment to within
00697  *              one bit; only in really hard cases do we need to
00698  *              compute a second residual.
00699  *      4. Because of 3., we don't need a large table of powers of 10
00700  *              for ten-to-e (just some small tables, e.g. of 10^k
00701  *              for 0 <= k <= 22).
00702  */
00703 
00704 /*
00705  * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least
00706  *      significant byte has the lowest address.
00707  * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most
00708  *      significant byte has the lowest address.
00709  * #define Long int on machines with 32-bit ints and 64-bit longs.
00710  * #define IBM for IBM mainframe-style floating-point arithmetic.
00711  * #define VAX for VAX-style floating-point arithmetic (D_floating).
00712  * #define No_leftright to omit left-right logic in fast floating-point
00713  *      computation of dtoa.
00714  * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
00715  *      and strtod and dtoa should round accordingly.
00716  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
00717  *      and Honor_FLT_ROUNDS is not #defined.
00718  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
00719  *      that use extended-precision instructions to compute rounded
00720  *      products and quotients) with IBM.
00721  * #define ROUND_BIASED for IEEE-format with biased rounding.
00722  * #define Inaccurate_Divide for IEEE-format with correctly rounded
00723  *      products but inaccurate quotients, e.g., for Intel i860.
00724  * #define NO_LONG_LONG on machines that do not have a "long long"
00725  *      integer type (of >= 64 bits).  On such machines, you can
00726  *      #define Just_16 to store 16 bits per 32-bit Long when doing
00727  *      high-precision integer arithmetic.  Whether this speeds things
00728  *      up or slows things down depends on the machine and the number
00729  *      being converted.  If long long is available and the name is
00730  *      something other than "long long", #define Llong to be the name,
00731  *      and if "unsigned Llong" does not work as an unsigned version of
00732  *      Llong, #define #ULLong to be the corresponding unsigned type.
00733  * #define KR_headers for old-style C function headers.
00734  * #define Bad_float_h if your system lacks a float.h or if it does not
00735  *      define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
00736  *      FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
00737  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
00738  *      if memory is available and otherwise does something you deem
00739  *      appropriate.  If MALLOC is undefined, malloc will be invoked
00740  *      directly -- and assumed always to succeed.
00741  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
00742  *      memory allocations from a private pool of memory when possible.
00743  *      When used, the private pool is PRIVATE_MEM bytes long:  2304 bytes,
00744  *      unless #defined to be a different length.  This default length
00745  *      suffices to get rid of MALLOC calls except for unusual cases,
00746  *      such as decimal-to-binary conversion of a very long string of
00747  *      digits.  The longest string dtoa can return is about 751 bytes
00748  *      long.  For conversions by strtod of strings of 800 digits and
00749  *      all dtoa conversions in single-threaded executions with 8-byte
00750  *      pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
00751  *      pointers, PRIVATE_MEM >= 7112 appears adequate.
00752  * #define INFNAN_CHECK on IEEE systems to cause strtod to check for
00753  *      Infinity and NaN (case insensitively).  On some systems (e.g.,
00754  *      some HP systems), it may be necessary to #define NAN_WORD0
00755  *      appropriately -- to the most significant word of a quiet NaN.
00756  *      (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
00757  *      When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
00758  *      strtod also accepts (case insensitively) strings of the form
00759  *      NaN(x), where x is a string of hexadecimal digits and spaces;
00760  *      if there is only one string of hexadecimal digits, it is taken
00761  *      for the 52 fraction bits of the resulting NaN; if there are two
00762  *      or more strings of hex digits, the first is for the high 20 bits,
00763  *      the second and subsequent for the low 32 bits, with intervening
00764  *      white space ignored; but if this results in none of the 52
00765  *      fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
00766  *      and NAN_WORD1 are used instead.
00767  * #define MULTIPLE_THREADS if the system offers preemptively scheduled
00768  *      multiple threads.  In this case, you must provide (or suitably
00769  *      #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
00770  *      by FREE_DTOA_LOCK(n) for n = 0 or 1.  (The second lock, accessed
00771  *      in pow5mult, ensures lazy evaluation of only one copy of high
00772  *      powers of 5; omitting this lock would introduce a small
00773  *      probability of wasting memory, but would otherwise be harmless.)
00774  *      You must also invoke freedtoa(s) to free the value s returned by
00775  *      dtoa.  You may do so whether or not MULTIPLE_THREADS is #defined.
00776  * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
00777  *      avoids underflows on inputs whose result does not underflow.
00778  *      If you #define NO_IEEE_Scale on a machine that uses IEEE-format
00779  *      floating-point numbers and flushes underflows to zero rather
00780  *      than implementing gradual underflow, then you must also #define
00781  *      Sudden_Underflow.
00782  * #define YES_ALIAS to permit aliasing certain double values with
00783  *      arrays of ULongs.  This leads to slightly better code with
00784  *      some compilers and was always used prior to 19990916, but it
00785  *      is not strictly legal and can cause trouble with aggressively
00786  *      optimizing compilers (e.g., gcc 2.95.1 under -O2).
00787  * #define USE_LOCALE to use the current locale's decimal_point value.
00788  * #define SET_INEXACT if IEEE arithmetic is being used and extra
00789  *      computation should be done to set the inexact flag when the
00790  *      result is inexact and avoid setting inexact when the result
00791  *      is exact.  In this case, dtoa.c must be compiled in
00792  *      an environment, perhaps provided by #include "dtoa.c" in a
00793  *      suitable wrapper, that defines two functions,
00794  *              int get_inexact(void);
00795  *              void clear_inexact(void);
00796  *      such that get_inexact() returns a nonzero value if the
00797  *      inexact bit is already set, and clear_inexact() sets the
00798  *      inexact bit to 0.  When SET_INEXACT is #defined, strtod
00799  *      also does extra computations to set the underflow and overflow
00800  *      flags when appropriate (i.e., when the result is tiny and
00801  *      inexact or when it is a numeric value rounded to +-infinity).
00802  * #define NO_ERRNO if strtod should not assign errno = ERANGE when
00803  *      the result overflows to +-Infinity or underflows to 0.
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 /*IEEE_Arith*/
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 /* ifndef Bad_float_h */
00903 #include "float.h"
00904 #endif /* Bad_float_h */
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 /* The following definition of Storeinc is appropriate for MIPS processors.
00946  * An alternative that might be better on some machines is
00947  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
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 /* #define P DBL_MANT_DIG */
00958 /* Ten_pmax = floor(P*log(2)/log(5)) */
00959 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
00960 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
00961 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
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     /* debugging option */
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 /*Flt_Rounds*/
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 /* ifndef IEEE_Arith */
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 /* exponent has 7 bits, but 8 is the right value in b2d */
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 /* VAX */
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 /* IBM, VAX */
01070 #endif /* IEEE_Arith */
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 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
01099  * This makes some inner loops simpler and sometimes saves work
01100  * during multiplications, but it often seems to make things slightly
01101  * slower.  Hence the default is now to store 32 bits per Long.
01102  */
01103 #endif
01104 #else   /* long long available */
01105 #ifndef Llong
01106 #define Llong long long
01107 #endif
01108 #ifndef ULLong
01109 #define ULLong unsigned Llong
01110 #endif
01111 #endif /* NO_LONG_LONG */
01112 
01113 #define MULTIPLE_THREADS 1
01114 
01115 #ifndef MULTIPLE_THREADS
01116 #define ACQUIRE_DTOA_LOCK(n)    /*nothing*/
01117 #define FREE_DTOA_LOCK(n)       /*nothing*/
01118 #else
01119 #define ACQUIRE_DTOA_LOCK(n)    /*unused right now*/
01120 #define FREE_DTOA_LOCK(n)       /*unused right now*/
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)   /* multiply by m and add 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         /* first time */
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;   /* clear sign bit, which we ignore */
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     /* = 2^106 * 1e-53 */
01969 #else
01970     1e-256
01971 #endif
01972 };
01973 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
01974 /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
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; /* invalid form: don't change *sp */
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 /*No_Hex_NaN*/
02067 #endif /* INFNAN_CHECK */
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             /* no break */
02101           case '+':
02102             if (*++s)
02103                 goto break2;
02104             /* no break */
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; /* +: 2B, -: 2D */
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                     /* Avoid confusion from exponents
02243                      * so large that e might overflow.
02244                      */
02245                     e = 19999; /* safe for 16 bit ints */
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             /* Check for Nan and Infinity */
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 /* INFNAN_CHECK */
02286 ret0:
02287             s = s00;
02288             sign = 0;
02289         }
02290         goto ret;
02291     }
02292     e1 = e -= nf;
02293 
02294     /* Now we have nd0 digits, starting at s0, followed by a
02295      * decimal point, followed by nd-nd0 digits.  The number we're
02296      * after is the integer represented by those digits times
02297      * 10**e */
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                 /* round correctly FLT_ROUNDS = 2 or 3 */
02327                 if (sign) {
02328                     dval(rv) = -dval(rv);
02329                     sign = 0;
02330                 }
02331 #endif
02332                 /* rv = */ rounded_product(dval(rv), tens[e]);
02333                 goto ret;
02334 #endif
02335             }
02336             i = DBL_DIG - nd;
02337             if (e <= Ten_pmax + i) {
02338                 /* A fancier test would sometimes let us do
02339                  * this for larger i values.
02340                  */
02341 #ifdef Honor_FLT_ROUNDS
02342                 /* round correctly FLT_ROUNDS = 2 or 3 */
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                 /* VAX exponent range is so narrow we must
02352                  * worry about overflow here...
02353                  */
02354 vax_ovfl_check:
02355                 word0(rv) -= P*Exp_msk1;
02356                 /* rv = */ 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                 /* rv = */ 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             /* round correctly FLT_ROUNDS = 2 or 3 */
02371             if (sign) {
02372                 dval(rv) = -dval(rv);
02373                 sign = 0;
02374             }
02375 #endif
02376             /* rv = */ 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 /*IEEE_Arith*/
02402 
02403     /* Get starting approximation = rv * 10**e1 */
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                 /* Can't trust HUGE_VAL */
02415 #ifdef IEEE_Arith
02416 #ifdef Honor_FLT_ROUNDS
02417                 switch (rounding) {
02418                   case 0: /* toward 0 */
02419                   case 3: /* toward -infinity */
02420                     word0(rv) = Big0;
02421                     word1(rv) = Big1;
02422                     break;
02423                   default:
02424                     word0(rv) = Exp_mask;
02425                     word1(rv) = 0;
02426                 }
02427 #else /*Honor_FLT_ROUNDS*/
02428                 word0(rv) = Exp_mask;
02429                 word1(rv) = 0;
02430 #endif /*Honor_FLT_ROUNDS*/
02431 #ifdef SET_INEXACT
02432                 /* set overflow bit */
02433                 dval(rv0) = 1e300;
02434                 dval(rv0) *= dval(rv0);
02435 #endif
02436 #else /*IEEE_Arith*/
02437                 word0(rv) = Big0;
02438                 word1(rv) = Big1;
02439 #endif /*IEEE_Arith*/
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             /* The last multiplication could overflow. */
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                 /* set to largest number */
02456                 /* (Can't trust DBL_MAX) */
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                 /* scaled rv is denormal; zap j low bits */
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             /* The last multiplication could underflow. */
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                 /* The refinement below will clean
02515                  * this approximation up.
02516                  */
02517             }
02518 #endif
02519         }
02520     }
02521 
02522     /* Now the hard part -- adjusting rv to the correct value.*/
02523 
02524     /* Put digits into bd: true value = bd * 10^e */
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);  /* rv = bb * 2^bbe */
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; /* logb(rv) */
02554         if (i < Emin)   /* denormal */
02555             j += P - Emin;
02556         else
02557             j = P + 1 - bbbits;
02558 #else /*Avoid_Underflow*/
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 /*Sudden_Underflow*/
02566         j = bbe;
02567         i = j + bbbits - 1; /* logb(rv) */
02568         if (i < Emin)   /* denormal */
02569             j += P - Emin;
02570         else
02571             j = P + 1 - bbbits;
02572 #endif /*Sudden_Underflow*/
02573 #endif /*Avoid_Underflow*/
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                 /* Error is less than an ulp */
02609                 if (!delta->x[0] && delta->wds <= 1) {
02610                     /* exact */
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 /*Sudden_Underflow*/
02653 #endif /*Avoid_Underflow*/
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                 /* adj = rounding ? ceil(adj) : floor(adj); */
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 /*Sudden_Underflow*/
02686 #endif /*Avoid_Underflow*/
02687             adj *= ulp(dval(rv));
02688             if (dsign)
02689                 dval(rv) += adj;
02690             else
02691                 dval(rv) -= adj;
02692             goto cont;
02693         }
02694 #endif /*Honor_FLT_ROUNDS*/
02695 
02696         if (i < 0) {
02697             /* Error is less than half an ulp -- check for
02698              * special case of mantissa a power of two.
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                 /* exact result */
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             /* exactly half-way between */
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                     /*boundary case -- increment exponent*/
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                 /* boundary case -- decrement exponent */
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 /*Avoid_Underflow*/
02764 #endif /*IBM*/
02765                     goto undfl;
02766                 L -= Exp_msk1;
02767 #else /*Sudden_Underflow}{*/
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                             /* round even ==> */
02774                             /* accept rv */
02775                             break;
02776                         /* rv = smallest denormal */
02777                         goto undfl;
02778                     }
02779                 }
02780 #endif /*Avoid_Underflow*/
02781                 L = (word0(rv) & Exp_mask) - Exp_msk1;
02782 #endif /*Sudden_Underflow}}*/
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                 /* special case -- power of FLT_RADIX to be */
02824                 /* rounded down... */
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: /* towards +infinity */
02839                 dval(aadj1) -= 0.5;
02840                 break;
02841               case 0: /* towards 0 */
02842               case 3: /* towards -infinity */
02843                 dval(aadj1) += 0.5;
02844             }
02845 #else
02846             if (Flt_Rounds == 0)
02847                 dval(aadj1) += 0.5;
02848 #endif /*Check_FLT_ROUNDS*/
02849         }
02850         y = word0(rv) & Exp_mask;
02851 
02852         /* Check for overflow */
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 /*Sudden_Underflow*/
02910             /* Compute adj so that the IEEE rounding rules will
02911              * correctly round rv + adj in some half-way cases.
02912              * If rv * ulp(rv) is denormalized (i.e.,
02913              * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
02914              * trouble from bits lost to denormalization;
02915              * example: 1.2e-307 .
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 /*Sudden_Underflow*/
02925 #endif /*Avoid_Underflow*/
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             /* Can we stop now? */
02934             L = (Long)aadj;
02935             aadj -= L;
02936             /* The tolerances below are conservative. */
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         /* try to avoid the bug of testing an 8087 register value */
02969         if (word0(rv) == 0 && word1(rv) == 0)
02970             errno = ERANGE;
02971 #endif
02972     }
02973 #endif /* Avoid_Underflow */
02974 #ifdef SET_INEXACT
02975     if (inexact && !(word0(rv) & Exp_mask)) {
02976         /* set underflow bit */
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     /*debug*/ if (b->wds > n)
03010     /*debug*/   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);  /* ensure q <= true quotient */
03019 #ifdef DEBUG
03020     /*debug*/ if (q > 9)
03021     /*debug*/   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 /* freedtoa(s) must be used to free values s returned by dtoa
03134  * when MULTIPLE_THREADS is #defined.  It should be used in all cases,
03135  * but for consistency with earlier versions of dtoa, it is optional
03136  * when MULTIPLE_THREADS is not defined.
03137  */
03138 
03139 static void
03140 freedtoa(char *s)
03141 {
03142     xfree(s);
03143 }
03144 #endif
03145 
03146 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
03147  *
03148  * Inspired by "How to Print Floating-Point Numbers Accurately" by
03149  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
03150  *
03151  * Modifications:
03152  *  1. Rather than iterating, we use a simple numeric overestimate
03153  *     to determine k = floor(log10(d)).  We scale relevant
03154  *     quantities using O(log2(k)) rather than O(k) multiplications.
03155  *  2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
03156  *     try to generate digits strictly left to right.  Instead, we
03157  *     compute with fewer bits and propagate the carry if necessary
03158  *     when rounding the final digit up.  This is often faster.
03159  *  3. Under the assumption that input will be rounded nearest,
03160  *     mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
03161  *     That is, we allow equality in stopping tests when the
03162  *     round-nearest rule will give the same floating-point value
03163  *     as would satisfaction of the stopping test with strict
03164  *     inequality.
03165  *  4. We remove common factors of powers of 2 from relevant
03166  *     quantities.
03167  *  5. When converting floating-point integers less than 1e16,
03168  *     we use floating-point arithmetic rather than resorting
03169  *     to multiple-precision integers.
03170  *  6. When asked to produce fewer than 15 digits, we first try
03171  *     to get by with floating-point arithmetic; we resort to
03172  *     multiple-precision integer arithmetic only if we cannot
03173  *     guarantee that the floating-point calculation has given
03174  *     the correctly rounded result.  For k requested digits and
03175  *     "uniformly" distributed input, the probability is
03176  *     something like 10^(k-15) that we must resort to the Long
03177  *     calculation.
03178  */
03179 
03180 char *
03181 ruby_dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
03182 {
03183  /* Arguments ndigits, decpt, sign are similar to those
03184     of ecvt and fcvt; trailing zeros are suppressed from
03185     the returned string.  If not null, *rve is set to point
03186     to the end of the return value.  If d is +-Infinity or NaN,
03187     then *decpt is set to 9999.
03188 
03189     mode:
03190         0 ==> shortest string that yields d when read in
03191             and rounded to nearest.
03192         1 ==> like 0, but with Steele & White stopping rule;
03193             e.g. with IEEE P754 arithmetic , mode 0 gives
03194             1e23 whereas mode 1 gives 9.999999999999999e22.
03195         2 ==> max(1,ndigits) significant digits.  This gives a
03196             return value similar to that of ecvt, except
03197             that trailing zeros are suppressed.
03198         3 ==> through ndigits past the decimal point.  This
03199             gives a return value similar to that from fcvt,
03200             except that trailing zeros are suppressed, and
03201             ndigits can be negative.
03202         4,5 ==> similar to 2 and 3, respectively, but (in
03203             round-nearest mode) with the tests of mode 0 to
03204             possibly return a shorter string that rounds to d.
03205             With IEEE arithmetic and compilation with
03206             -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
03207             as modes 2 and 3 when FLT_ROUNDS != 1.
03208         6-9 ==> Debugging modes similar to mode - 4:  don't try
03209             fast floating-point estimate (if applicable).
03210 
03211         Values of mode other than 0-9 are treated as mode 0.
03212 
03213         Sufficient space is allocated to the return value
03214         to hold the suppressed trailing zeros.
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         /* set sign for everything, including 0's and NaNs */
03247         *sign = 1;
03248         word0(d) &= ~Sign_bit;  /* clear 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         /* Infinity or NaN */
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; /* normalize */
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         /* log(x)   ~=~ log(1.5) + (x-1.5)/1.5
03306          * log10(x)  =  log(x) / log(10)
03307          *      ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
03308          * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
03309          *
03310          * This suggests computing an approximation k to log10(d) by
03311          *
03312          * k = (i - Bias)*0.301029995663981
03313          *  + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
03314          *
03315          * We want k to be too large rather than too small.
03316          * The error in the first-order Taylor series approximation
03317          * is in our favor, so we just round up the constant enough
03318          * to compensate for any error in the multiplication of
03319          * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
03320          * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
03321          * adding 1e-13 to the constant term more than suffices.
03322          * Hence we adjust the constant term to 0.1760912590558.
03323          * (We could get a more accurate k by invoking log10,
03324          *  but this is probably not worthwhile.)
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         /* d is denormalized */
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; /* adjust exponent */
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--;    /* want k = floor(ds) */
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 /*SET_INEXACT*/
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         /* no break */
03402       case 4:
03403         if (ndigits <= 0)
03404             ndigits = 1;
03405         ilim = ilim1 = i = ndigits;
03406         break;
03407       case 3:
03408         leftright = 0;
03409         /* no break */
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         /* Try to get by with floating-point arithmetic. */
03427 
03428         i = 0;
03429         dval(d2) = dval(d);
03430         k0 = k;
03431         ilim0 = ilim;
03432         ieps = 2; /* conservative */
03433         if (k > 0) {
03434             ds = tens[k&0xf];
03435             j = k >> 4;
03436             if (j & Bletch) {
03437                 /* prevent overflows */
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             /* Use Steele & White method of only
03479              * generating digits needed.
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             /* Generate ilim digits, then fix them up. */
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     /* Do we have a "small" integer? */
03527 
03528     if (be >= 0 && k <= Int_max) {
03529         /* Yes. */
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             /* If FLT_ROUNDS == 2, L will usually be high by 1 */
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     /* Check for special case that d is a normalized power of 2. */
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             /* The special case */
03633             b2 += Log2P;
03634             s2 += Log2P;
03635             spec_case = 1;
03636         }
03637     }
03638 
03639     /* Arrange for convenient computation of quotients:
03640      * shift left if necessary so divisor has 4 leading 0 bits.
03641      *
03642      * Perhaps we should just compute leading 28 bits of S once
03643      * and for all and pass them and a shift to quorem, so it
03644      * can do shifts and ors to compute the numerator for q.
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);  /* we botched the k estimate */
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             /* no digits, fcvt style */
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         /* Compute mlo -- check for special case
03695          * that d is a normalized power of 2.
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             /* Do we yet have the shortest decimal string
03708              * that will round to d?
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 /*Honor_FLT_ROUNDS*/
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') { /* possible if i == 1 */
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     /* Round off last digit */
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); /* assume no string exceeds INT_MAX */
03866         (*func)(str, len, arg);
03867     }
03868 }
03869 
03870 /*-
03871  * Copyright (c) 2004-2008 David Schultz <das@FreeBSD.ORG>
03872  * All rights reserved.
03873  *
03874  * Redistribution and use in source and binary forms, with or without
03875  * modification, are permitted provided that the following conditions
03876  * are met:
03877  * 1. Redistributions of source code must retain the above copyright
03878  *    notice, this list of conditions and the following disclaimer.
03879  * 2. Redistributions in binary form must reproduce the above copyright
03880  *    notice, this list of conditions and the following disclaimer in the
03881  *    documentation and/or other materials provided with the distribution.
03882  *
03883  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
03884  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
03885  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
03886  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
03887  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
03888  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
03889  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
03890  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
03891  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
03892  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
03893  * SUCH DAMAGE.
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  * This procedure converts a double-precision number in IEEE format
03910  * into a string of hexadecimal digits and an exponent of 2.  Its
03911  * behavior is bug-for-bug compatible with dtoa() in mode 2, with the
03912  * following exceptions:
03913  *
03914  * - An ndigits < 0 causes it to use as many digits as necessary to
03915  *   represent the number exactly.
03916  * - The additional xdigs argument should point to either the string
03917  *   "0123456789ABCDEF" or the string "0123456789abcdef", depending on
03918  *   which case is desired.
03919  * - This routine does not repeat dtoa's mistake of setting decpt
03920  *   to 9999 in the case of an infinity or NaN.  INT_MAX is used
03921  *   for this purpose instead.
03922  *
03923  * Note that the C99 standard does not specify what the leading digit
03924  * should be for non-zero numbers.  For instance, 0x1.3p3 is the same
03925  * as 0x2.6p2 is the same as 0x4.cp3.  This implementation always makes
03926  * the leading digit a 1. This ensures that the exponent printed is the
03927  * actual base-2 exponent, i.e., ilogb(d).
03928  *
03929  * Inputs:      d, xdigs, ndigits
03930  * Outputs:     decpt, sign, rve
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             /* set sign for everything, including 0's and NaNs */
03944             *sign = 1;
03945             word0(u) &= ~Sign_bit;  /* clear sign bit */
03946         }
03947         else
03948             *sign = 0;
03949 
03950         if (isinf(d)) { /* FP_INFINITE */
03951             *decpt = INT_MAX;
03952             return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
03953         }
03954         else if (isnan(d)) { /* FP_NAN */
03955             *decpt = INT_MAX;
03956             return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
03957         }
03958         else if (d == 0.0) { /* FP_ZERO */
03959             *decpt = 1;
03960             return (nrv_alloc("0", rve, 1));
03961         }
03962         else if (dexp_get(u)) { /* FP_NORMAL */
03963             *decpt = dexp_get(u) - DBL_ADJ;
03964         }
03965         else { /* FP_SUBNORMAL */
03966             u.d *= 5.363123171977039e+154 /* 0x1p514 */;
03967             *decpt = dexp_get(u) - (514 + DBL_ADJ);
03968         }
03969 
03970         if (ndigits == 0)               /* dtoa() compatibility */
03971                 ndigits = 1;
03972 
03973         /*
03974          * If ndigits < 0, we are expected to auto-size, so we allocate
03975          * enough space for all the digits.
03976          */
03977         bufsize = (ndigits > 0) ? ndigits : SIGFIGS;
03978         s0 = rv_alloc(bufsize);
03979 
03980         /* Round to the desired number of digits. */
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         /* If ndigits < 0, we are expected to auto-size the precision. */
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 

Generated on Wed Sep 8 2010 21:55:22 for Ruby by  doxygen 1.7.1