• 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: yugui $
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             aadj = -1;
02124 
02125             if  (!s[1]) {
02126                 rb_warn("malformed value for Float(): %s. Ruby 1.9.3 for later will raise an ArgumentError for the value.", s00);
02127             }
02128             while (*++s && (s1 = strchr(hexdigit, *s))) {
02129                 adj *= 16;
02130                 adj += (s1 - hexdigit) & 15;
02131             }
02132 
02133             if (*s == '.') {
02134                 aadj = 1.;
02135                 while (*++s && (s1 = strchr(hexdigit, *s))) {
02136                     aadj /= 16;
02137                     adj += aadj * ((s1 - hexdigit) & 15);
02138                 }
02139             }
02140 
02141             if (*s == 'P' || *s == 'p') {
02142                 dsign = 0x2C - *++s; /* +: 2B, -: 2D */
02143                 if (abs(dsign) == 1) s++;
02144                 else dsign = 1;
02145 
02146                nd = 0;
02147                c = *s;
02148                if (c < '0' || '9' < c) goto ret0;
02149                do {
02150                     nd *= 10;
02151                     nd += c;
02152                     nd -= '0';
02153                    c = *++s;
02154                } while ('0' <= c && c <= '9');
02155                 dval(rv) = ldexp(adj, nd * dsign);
02156             }
02157             else {
02158                 if (aadj != -1) {
02159                     rb_warn("malformed value for Float(): %s. Ruby 1.9.3 for later will raise an ArgumentError for the value.", s00);
02160                 }
02161                 dval(rv) = adj;
02162             }
02163             goto ret;
02164         }
02165         nz0 = 1;
02166         while (*++s == '0') ;
02167         if (!*s)
02168             goto ret;
02169     }
02170     s0 = s;
02171     y = z = 0;
02172     for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
02173         if (nd < 9)
02174             y = 10*y + c - '0';
02175         else if (nd < 16)
02176             z = 10*z + c - '0';
02177     nd0 = nd;
02178 #ifdef USE_LOCALE
02179     s1 = localeconv()->decimal_point;
02180     if (c == *s1) {
02181         c = '.';
02182         if (*++s1) {
02183             s2 = s;
02184             for (;;) {
02185                 if (*++s2 != *s1) {
02186                     c = 0;
02187                     break;
02188                 }
02189                 if (!*++s1) {
02190                     s = s2;
02191                     break;
02192                 }
02193             }
02194         }
02195     }
02196 #endif
02197     if (c == '.') {
02198         if (!ISDIGIT(s[1]))
02199             goto dig_done;
02200         c = *++s;
02201         if (!nd) {
02202             for (; c == '0'; c = *++s)
02203                 nz++;
02204             if (c > '0' && c <= '9') {
02205                 s0 = s;
02206                 nf += nz;
02207                 nz = 0;
02208                 goto have_dig;
02209             }
02210             goto dig_done;
02211         }
02212         for (; c >= '0' && c <= '9'; c = *++s) {
02213 have_dig:
02214             nz++;
02215             if (c -= '0') {
02216                 nf += nz;
02217                 for (i = 1; i < nz; i++)
02218                     if (nd++ < 9)
02219                         y *= 10;
02220                     else if (nd <= DBL_DIG + 1)
02221                         z *= 10;
02222                 if (nd++ < 9)
02223                     y = 10*y + c;
02224                 else if (nd <= DBL_DIG + 1)
02225                     z = 10*z + c;
02226                 nz = 0;
02227             }
02228         }
02229     }
02230 dig_done:
02231     e = 0;
02232     if (c == 'e' || c == 'E') {
02233         if (!nd && !nz && !nz0) {
02234             goto ret0;
02235         }
02236         s00 = s;
02237         esign = 0;
02238         switch (c = *++s) {
02239           case '-':
02240             esign = 1;
02241           case '+':
02242             c = *++s;
02243         }
02244         if (c >= '0' && c <= '9') {
02245             while (c == '0')
02246                 c = *++s;
02247             if (c > '0' && c <= '9') {
02248                 L = c - '0';
02249                 s1 = s;
02250                 while ((c = *++s) >= '0' && c <= '9')
02251                     L = 10*L + c - '0';
02252                 if (s - s1 > 8 || L > 19999)
02253                     /* Avoid confusion from exponents
02254                      * so large that e might overflow.
02255                      */
02256                     e = 19999; /* safe for 16 bit ints */
02257                 else
02258                     e = (int)L;
02259                 if (esign)
02260                     e = -e;
02261             }
02262             else
02263                 e = 0;
02264         }
02265         else
02266             s = s00;
02267     }
02268     if (!nd) {
02269         if (!nz && !nz0) {
02270 #ifdef INFNAN_CHECK
02271             /* Check for Nan and Infinity */
02272             switch (c) {
02273               case 'i':
02274               case 'I':
02275                 if (match(&s,"nf")) {
02276                     --s;
02277                     if (!match(&s,"inity"))
02278                         ++s;
02279                     word0(rv) = 0x7ff00000;
02280                     word1(rv) = 0;
02281                     goto ret;
02282                 }
02283                 break;
02284               case 'n':
02285               case 'N':
02286                 if (match(&s, "an")) {
02287                     word0(rv) = NAN_WORD0;
02288                     word1(rv) = NAN_WORD1;
02289 #ifndef No_Hex_NaN
02290                     if (*s == '(') /*)*/
02291                         hexnan(&rv, &s);
02292 #endif
02293                     goto ret;
02294                 }
02295             }
02296 #endif /* INFNAN_CHECK */
02297 ret0:
02298             s = s00;
02299             sign = 0;
02300         }
02301         goto ret;
02302     }
02303     e1 = e -= nf;
02304 
02305     /* Now we have nd0 digits, starting at s0, followed by a
02306      * decimal point, followed by nd-nd0 digits.  The number we're
02307      * after is the integer represented by those digits times
02308      * 10**e */
02309 
02310     if (!nd0)
02311         nd0 = nd;
02312     k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
02313     dval(rv) = y;
02314     if (k > 9) {
02315 #ifdef SET_INEXACT
02316         if (k > DBL_DIG)
02317             oldinexact = get_inexact();
02318 #endif
02319         dval(rv) = tens[k - 9] * dval(rv) + z;
02320     }
02321     bd0 = bb = bd = bs = delta = 0;
02322     if (nd <= DBL_DIG
02323 #ifndef RND_PRODQUOT
02324 #ifndef Honor_FLT_ROUNDS
02325         && Flt_Rounds == 1
02326 #endif
02327 #endif
02328     ) {
02329         if (!e)
02330             goto ret;
02331         if (e > 0) {
02332             if (e <= Ten_pmax) {
02333 #ifdef VAX
02334                 goto vax_ovfl_check;
02335 #else
02336 #ifdef Honor_FLT_ROUNDS
02337                 /* round correctly FLT_ROUNDS = 2 or 3 */
02338                 if (sign) {
02339                     dval(rv) = -dval(rv);
02340                     sign = 0;
02341                 }
02342 #endif
02343                 /* rv = */ rounded_product(dval(rv), tens[e]);
02344                 goto ret;
02345 #endif
02346             }
02347             i = DBL_DIG - nd;
02348             if (e <= Ten_pmax + i) {
02349                 /* A fancier test would sometimes let us do
02350                  * this for larger i values.
02351                  */
02352 #ifdef Honor_FLT_ROUNDS
02353                 /* round correctly FLT_ROUNDS = 2 or 3 */
02354                 if (sign) {
02355                     dval(rv) = -dval(rv);
02356                     sign = 0;
02357                 }
02358 #endif
02359                 e -= i;
02360                 dval(rv) *= tens[i];
02361 #ifdef VAX
02362                 /* VAX exponent range is so narrow we must
02363                  * worry about overflow here...
02364                  */
02365 vax_ovfl_check:
02366                 word0(rv) -= P*Exp_msk1;
02367                 /* rv = */ rounded_product(dval(rv), tens[e]);
02368                 if ((word0(rv) & Exp_mask)
02369                         > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
02370                     goto ovfl;
02371                 word0(rv) += P*Exp_msk1;
02372 #else
02373                 /* rv = */ rounded_product(dval(rv), tens[e]);
02374 #endif
02375                 goto ret;
02376             }
02377         }
02378 #ifndef Inaccurate_Divide
02379         else if (e >= -Ten_pmax) {
02380 #ifdef Honor_FLT_ROUNDS
02381             /* round correctly FLT_ROUNDS = 2 or 3 */
02382             if (sign) {
02383                 dval(rv) = -dval(rv);
02384                 sign = 0;
02385             }
02386 #endif
02387             /* rv = */ rounded_quotient(dval(rv), tens[-e]);
02388             goto ret;
02389         }
02390 #endif
02391     }
02392     e1 += nd - k;
02393 
02394 #ifdef IEEE_Arith
02395 #ifdef SET_INEXACT
02396     inexact = 1;
02397     if (k <= DBL_DIG)
02398         oldinexact = get_inexact();
02399 #endif
02400 #ifdef Avoid_Underflow
02401     scale = 0;
02402 #endif
02403 #ifdef Honor_FLT_ROUNDS
02404     if ((rounding = Flt_Rounds) >= 2) {
02405         if (sign)
02406             rounding = rounding == 2 ? 0 : 2;
02407         else
02408             if (rounding != 2)
02409                 rounding = 0;
02410     }
02411 #endif
02412 #endif /*IEEE_Arith*/
02413 
02414     /* Get starting approximation = rv * 10**e1 */
02415 
02416     if (e1 > 0) {
02417         if ((i = e1 & 15) != 0)
02418             dval(rv) *= tens[i];
02419         if (e1 &= ~15) {
02420             if (e1 > DBL_MAX_10_EXP) {
02421 ovfl:
02422 #ifndef NO_ERRNO
02423                 errno = ERANGE;
02424 #endif
02425                 /* Can't trust HUGE_VAL */
02426 #ifdef IEEE_Arith
02427 #ifdef Honor_FLT_ROUNDS
02428                 switch (rounding) {
02429                   case 0: /* toward 0 */
02430                   case 3: /* toward -infinity */
02431                     word0(rv) = Big0;
02432                     word1(rv) = Big1;
02433                     break;
02434                   default:
02435                     word0(rv) = Exp_mask;
02436                     word1(rv) = 0;
02437                 }
02438 #else /*Honor_FLT_ROUNDS*/
02439                 word0(rv) = Exp_mask;
02440                 word1(rv) = 0;
02441 #endif /*Honor_FLT_ROUNDS*/
02442 #ifdef SET_INEXACT
02443                 /* set overflow bit */
02444                 dval(rv0) = 1e300;
02445                 dval(rv0) *= dval(rv0);
02446 #endif
02447 #else /*IEEE_Arith*/
02448                 word0(rv) = Big0;
02449                 word1(rv) = Big1;
02450 #endif /*IEEE_Arith*/
02451                 if (bd0)
02452                     goto retfree;
02453                 goto ret;
02454             }
02455             e1 >>= 4;
02456             for (j = 0; e1 > 1; j++, e1 >>= 1)
02457                 if (e1 & 1)
02458                     dval(rv) *= bigtens[j];
02459             /* The last multiplication could overflow. */
02460             word0(rv) -= P*Exp_msk1;
02461             dval(rv) *= bigtens[j];
02462             if ((z = word0(rv) & Exp_mask)
02463                     > Exp_msk1*(DBL_MAX_EXP+Bias-P))
02464                 goto ovfl;
02465             if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
02466                 /* set to largest number */
02467                 /* (Can't trust DBL_MAX) */
02468                 word0(rv) = Big0;
02469                 word1(rv) = Big1;
02470             }
02471             else
02472                 word0(rv) += P*Exp_msk1;
02473         }
02474     }
02475     else if (e1 < 0) {
02476         e1 = -e1;
02477         if ((i = e1 & 15) != 0)
02478             dval(rv) /= tens[i];
02479         if (e1 >>= 4) {
02480             if (e1 >= 1 << n_bigtens)
02481                 goto undfl;
02482 #ifdef Avoid_Underflow
02483             if (e1 & Scale_Bit)
02484                 scale = 2*P;
02485             for (j = 0; e1 > 0; j++, e1 >>= 1)
02486                 if (e1 & 1)
02487                     dval(rv) *= tinytens[j];
02488             if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
02489                     >> Exp_shift)) > 0) {
02490                 /* scaled rv is denormal; zap j low bits */
02491                 if (j >= 32) {
02492                     word1(rv) = 0;
02493                     if (j >= 53)
02494                         word0(rv) = (P+2)*Exp_msk1;
02495                     else
02496                         word0(rv) &= 0xffffffff << (j-32);
02497                 }
02498                 else
02499                     word1(rv) &= 0xffffffff << j;
02500             }
02501 #else
02502             for (j = 0; e1 > 1; j++, e1 >>= 1)
02503                 if (e1 & 1)
02504                     dval(rv) *= tinytens[j];
02505             /* The last multiplication could underflow. */
02506             dval(rv0) = dval(rv);
02507             dval(rv) *= tinytens[j];
02508             if (!dval(rv)) {
02509                 dval(rv) = 2.*dval(rv0);
02510                 dval(rv) *= tinytens[j];
02511 #endif
02512                 if (!dval(rv)) {
02513 undfl:
02514                     dval(rv) = 0.;
02515 #ifndef NO_ERRNO
02516                     errno = ERANGE;
02517 #endif
02518                     if (bd0)
02519                         goto retfree;
02520                     goto ret;
02521                 }
02522 #ifndef Avoid_Underflow
02523                 word0(rv) = Tiny0;
02524                 word1(rv) = Tiny1;
02525                 /* The refinement below will clean
02526                  * this approximation up.
02527                  */
02528             }
02529 #endif
02530         }
02531     }
02532 
02533     /* Now the hard part -- adjusting rv to the correct value.*/
02534 
02535     /* Put digits into bd: true value = bd * 10^e */
02536 
02537     bd0 = s2b(s0, nd0, nd, y);
02538 
02539     for (;;) {
02540         bd = Balloc(bd0->k);
02541         Bcopy(bd, bd0);
02542         bb = d2b(dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
02543         bs = i2b(1);
02544 
02545         if (e >= 0) {
02546             bb2 = bb5 = 0;
02547             bd2 = bd5 = e;
02548         }
02549         else {
02550             bb2 = bb5 = -e;
02551             bd2 = bd5 = 0;
02552         }
02553         if (bbe >= 0)
02554             bb2 += bbe;
02555         else
02556             bd2 -= bbe;
02557         bs2 = bb2;
02558 #ifdef Honor_FLT_ROUNDS
02559         if (rounding != 1)
02560             bs2++;
02561 #endif
02562 #ifdef Avoid_Underflow
02563         j = bbe - scale;
02564         i = j + bbbits - 1; /* logb(rv) */
02565         if (i < Emin)   /* denormal */
02566             j += P - Emin;
02567         else
02568             j = P + 1 - bbbits;
02569 #else /*Avoid_Underflow*/
02570 #ifdef Sudden_Underflow
02571 #ifdef IBM
02572         j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
02573 #else
02574         j = P + 1 - bbbits;
02575 #endif
02576 #else /*Sudden_Underflow*/
02577         j = bbe;
02578         i = j + bbbits - 1; /* logb(rv) */
02579         if (i < Emin)   /* denormal */
02580             j += P - Emin;
02581         else
02582             j = P + 1 - bbbits;
02583 #endif /*Sudden_Underflow*/
02584 #endif /*Avoid_Underflow*/
02585         bb2 += j;
02586         bd2 += j;
02587 #ifdef Avoid_Underflow
02588         bd2 += scale;
02589 #endif
02590         i = bb2 < bd2 ? bb2 : bd2;
02591         if (i > bs2)
02592             i = bs2;
02593         if (i > 0) {
02594             bb2 -= i;
02595             bd2 -= i;
02596             bs2 -= i;
02597         }
02598         if (bb5 > 0) {
02599             bs = pow5mult(bs, bb5);
02600             bb1 = mult(bs, bb);
02601             Bfree(bb);
02602             bb = bb1;
02603         }
02604         if (bb2 > 0)
02605             bb = lshift(bb, bb2);
02606         if (bd5 > 0)
02607             bd = pow5mult(bd, bd5);
02608         if (bd2 > 0)
02609             bd = lshift(bd, bd2);
02610         if (bs2 > 0)
02611             bs = lshift(bs, bs2);
02612         delta = diff(bb, bd);
02613         dsign = delta->sign;
02614         delta->sign = 0;
02615         i = cmp(delta, bs);
02616 #ifdef Honor_FLT_ROUNDS
02617         if (rounding != 1) {
02618             if (i < 0) {
02619                 /* Error is less than an ulp */
02620                 if (!delta->x[0] && delta->wds <= 1) {
02621                     /* exact */
02622 #ifdef SET_INEXACT
02623                     inexact = 0;
02624 #endif
02625                     break;
02626                 }
02627                 if (rounding) {
02628                     if (dsign) {
02629                         adj = 1.;
02630                         goto apply_adj;
02631                     }
02632                 }
02633                 else if (!dsign) {
02634                     adj = -1.;
02635                     if (!word1(rv)
02636                      && !(word0(rv) & Frac_mask)) {
02637                         y = word0(rv) & Exp_mask;
02638 #ifdef Avoid_Underflow
02639                         if (!scale || y > 2*P*Exp_msk1)
02640 #else
02641                         if (y)
02642 #endif
02643                         {
02644                             delta = lshift(delta,Log2P);
02645                             if (cmp(delta, bs) <= 0)
02646                                 adj = -0.5;
02647                         }
02648                     }
02649 apply_adj:
02650 #ifdef Avoid_Underflow
02651                     if (scale && (y = word0(rv) & Exp_mask)
02652                             <= 2*P*Exp_msk1)
02653                         word0(adj) += (2*P+1)*Exp_msk1 - y;
02654 #else
02655 #ifdef Sudden_Underflow
02656                     if ((word0(rv) & Exp_mask) <=
02657                             P*Exp_msk1) {
02658                         word0(rv) += P*Exp_msk1;
02659                         dval(rv) += adj*ulp(dval(rv));
02660                         word0(rv) -= P*Exp_msk1;
02661                     }
02662                     else
02663 #endif /*Sudden_Underflow*/
02664 #endif /*Avoid_Underflow*/
02665                     dval(rv) += adj*ulp(dval(rv));
02666                 }
02667                 break;
02668             }
02669             adj = ratio(delta, bs);
02670             if (adj < 1.)
02671                 adj = 1.;
02672             if (adj <= 0x7ffffffe) {
02673                 /* adj = rounding ? ceil(adj) : floor(adj); */
02674                 y = adj;
02675                 if (y != adj) {
02676                     if (!((rounding>>1) ^ dsign))
02677                         y++;
02678                     adj = y;
02679                 }
02680             }
02681 #ifdef Avoid_Underflow
02682             if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02683                 word0(adj) += (2*P+1)*Exp_msk1 - y;
02684 #else
02685 #ifdef Sudden_Underflow
02686             if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02687                 word0(rv) += P*Exp_msk1;
02688                 adj *= ulp(dval(rv));
02689                 if (dsign)
02690                     dval(rv) += adj;
02691                 else
02692                     dval(rv) -= adj;
02693                 word0(rv) -= P*Exp_msk1;
02694                 goto cont;
02695             }
02696 #endif /*Sudden_Underflow*/
02697 #endif /*Avoid_Underflow*/
02698             adj *= ulp(dval(rv));
02699             if (dsign)
02700                 dval(rv) += adj;
02701             else
02702                 dval(rv) -= adj;
02703             goto cont;
02704         }
02705 #endif /*Honor_FLT_ROUNDS*/
02706 
02707         if (i < 0) {
02708             /* Error is less than half an ulp -- check for
02709              * special case of mantissa a power of two.
02710              */
02711             if (dsign || word1(rv) || word0(rv) & Bndry_mask
02712 #ifdef IEEE_Arith
02713 #ifdef Avoid_Underflow
02714                 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
02715 #else
02716                 || (word0(rv) & Exp_mask) <= Exp_msk1
02717 #endif
02718 #endif
02719             ) {
02720 #ifdef SET_INEXACT
02721                 if (!delta->x[0] && delta->wds <= 1)
02722                     inexact = 0;
02723 #endif
02724                 break;
02725             }
02726             if (!delta->x[0] && delta->wds <= 1) {
02727                 /* exact result */
02728 #ifdef SET_INEXACT
02729                 inexact = 0;
02730 #endif
02731                 break;
02732             }
02733             delta = lshift(delta,Log2P);
02734             if (cmp(delta, bs) > 0)
02735                 goto drop_down;
02736             break;
02737         }
02738         if (i == 0) {
02739             /* exactly half-way between */
02740             if (dsign) {
02741                 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
02742                         &&  word1(rv) == (
02743 #ifdef Avoid_Underflow
02744                         (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02745                         ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
02746 #endif
02747                         0xffffffff)) {
02748                     /*boundary case -- increment exponent*/
02749                     word0(rv) = (word0(rv) & Exp_mask)
02750                                 + Exp_msk1
02751 #ifdef IBM
02752                                 | Exp_msk1 >> 4
02753 #endif
02754                     ;
02755                     word1(rv) = 0;
02756 #ifdef Avoid_Underflow
02757                     dsign = 0;
02758 #endif
02759                     break;
02760                 }
02761             }
02762             else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
02763 drop_down:
02764                 /* boundary case -- decrement exponent */
02765 #ifdef Sudden_Underflow /*{{*/
02766                 L = word0(rv) & Exp_mask;
02767 #ifdef IBM
02768                 if (L <  Exp_msk1)
02769 #else
02770 #ifdef Avoid_Underflow
02771                 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
02772 #else
02773                 if (L <= Exp_msk1)
02774 #endif /*Avoid_Underflow*/
02775 #endif /*IBM*/
02776                     goto undfl;
02777                 L -= Exp_msk1;
02778 #else /*Sudden_Underflow}{*/
02779 #ifdef Avoid_Underflow
02780                 if (scale) {
02781                     L = word0(rv) & Exp_mask;
02782                     if (L <= (2*P+1)*Exp_msk1) {
02783                         if (L > (P+2)*Exp_msk1)
02784                             /* round even ==> */
02785                             /* accept rv */
02786                             break;
02787                         /* rv = smallest denormal */
02788                         goto undfl;
02789                     }
02790                 }
02791 #endif /*Avoid_Underflow*/
02792                 L = (word0(rv) & Exp_mask) - Exp_msk1;
02793 #endif /*Sudden_Underflow}}*/
02794                 word0(rv) = L | Bndry_mask1;
02795                 word1(rv) = 0xffffffff;
02796 #ifdef IBM
02797                 goto cont;
02798 #else
02799                 break;
02800 #endif
02801             }
02802 #ifndef ROUND_BIASED
02803             if (!(word1(rv) & LSB))
02804                 break;
02805 #endif
02806             if (dsign)
02807                 dval(rv) += ulp(dval(rv));
02808 #ifndef ROUND_BIASED
02809             else {
02810                 dval(rv) -= ulp(dval(rv));
02811 #ifndef Sudden_Underflow
02812                 if (!dval(rv))
02813                     goto undfl;
02814 #endif
02815             }
02816 #ifdef Avoid_Underflow
02817             dsign = 1 - dsign;
02818 #endif
02819 #endif
02820             break;
02821         }
02822         if ((aadj = ratio(delta, bs)) <= 2.) {
02823             if (dsign)
02824                 aadj = dval(aadj1) = 1.;
02825             else if (word1(rv) || word0(rv) & Bndry_mask) {
02826 #ifndef Sudden_Underflow
02827                 if (word1(rv) == Tiny1 && !word0(rv))
02828                     goto undfl;
02829 #endif
02830                 aadj = 1.;
02831                 dval(aadj1) = -1.;
02832             }
02833             else {
02834                 /* special case -- power of FLT_RADIX to be */
02835                 /* rounded down... */
02836 
02837                 if (aadj < 2./FLT_RADIX)
02838                     aadj = 1./FLT_RADIX;
02839                 else
02840                     aadj *= 0.5;
02841                 dval(aadj1) = -aadj;
02842             }
02843         }
02844         else {
02845             aadj *= 0.5;
02846             dval(aadj1) = dsign ? aadj : -aadj;
02847 #ifdef Check_FLT_ROUNDS
02848             switch (Rounding) {
02849               case 2: /* towards +infinity */
02850                 dval(aadj1) -= 0.5;
02851                 break;
02852               case 0: /* towards 0 */
02853               case 3: /* towards -infinity */
02854                 dval(aadj1) += 0.5;
02855             }
02856 #else
02857             if (Flt_Rounds == 0)
02858                 dval(aadj1) += 0.5;
02859 #endif /*Check_FLT_ROUNDS*/
02860         }
02861         y = word0(rv) & Exp_mask;
02862 
02863         /* Check for overflow */
02864 
02865         if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
02866             dval(rv0) = dval(rv);
02867             word0(rv) -= P*Exp_msk1;
02868             adj = dval(aadj1) * ulp(dval(rv));
02869             dval(rv) += adj;
02870             if ((word0(rv) & Exp_mask) >=
02871                     Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
02872                 if (word0(rv0) == Big0 && word1(rv0) == Big1)
02873                     goto ovfl;
02874                 word0(rv) = Big0;
02875                 word1(rv) = Big1;
02876                 goto cont;
02877             }
02878             else
02879                 word0(rv) += P*Exp_msk1;
02880         }
02881         else {
02882 #ifdef Avoid_Underflow
02883             if (scale && y <= 2*P*Exp_msk1) {
02884                 if (aadj <= 0x7fffffff) {
02885                     if ((z = (int)aadj) <= 0)
02886                         z = 1;
02887                     aadj = z;
02888                     dval(aadj1) = dsign ? aadj : -aadj;
02889                 }
02890                 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
02891             }
02892             adj = dval(aadj1) * ulp(dval(rv));
02893             dval(rv) += adj;
02894 #else
02895 #ifdef Sudden_Underflow
02896             if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02897                 dval(rv0) = dval(rv);
02898                 word0(rv) += P*Exp_msk1;
02899                 adj = dval(aadj1) * ulp(dval(rv));
02900                 dval(rv) += adj;
02901 #ifdef IBM
02902                 if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
02903 #else
02904                 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
02905 #endif
02906                 {
02907                     if (word0(rv0) == Tiny0 && word1(rv0) == Tiny1)
02908                         goto undfl;
02909                     word0(rv) = Tiny0;
02910                     word1(rv) = Tiny1;
02911                     goto cont;
02912                 }
02913                 else
02914                     word0(rv) -= P*Exp_msk1;
02915             }
02916             else {
02917                 adj = dval(aadj1) * ulp(dval(rv));
02918                 dval(rv) += adj;
02919             }
02920 #else /*Sudden_Underflow*/
02921             /* Compute adj so that the IEEE rounding rules will
02922              * correctly round rv + adj in some half-way cases.
02923              * If rv * ulp(rv) is denormalized (i.e.,
02924              * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
02925              * trouble from bits lost to denormalization;
02926              * example: 1.2e-307 .
02927              */
02928             if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
02929                 dval(aadj1) = (double)(int)(aadj + 0.5);
02930                 if (!dsign)
02931                     dval(aadj1) = -dval(aadj1);
02932             }
02933             adj = dval(aadj1) * ulp(dval(rv));
02934             dval(rv) += adj;
02935 #endif /*Sudden_Underflow*/
02936 #endif /*Avoid_Underflow*/
02937         }
02938         z = word0(rv) & Exp_mask;
02939 #ifndef SET_INEXACT
02940 #ifdef Avoid_Underflow
02941         if (!scale)
02942 #endif
02943         if (y == z) {
02944             /* Can we stop now? */
02945             L = (Long)aadj;
02946             aadj -= L;
02947             /* The tolerances below are conservative. */
02948             if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
02949                 if (aadj < .4999999 || aadj > .5000001)
02950                     break;
02951             }
02952             else if (aadj < .4999999/FLT_RADIX)
02953                 break;
02954         }
02955 #endif
02956 cont:
02957         Bfree(bb);
02958         Bfree(bd);
02959         Bfree(bs);
02960         Bfree(delta);
02961     }
02962 #ifdef SET_INEXACT
02963     if (inexact) {
02964         if (!oldinexact) {
02965             word0(rv0) = Exp_1 + (70 << Exp_shift);
02966             word1(rv0) = 0;
02967             dval(rv0) += 1.;
02968         }
02969     }
02970     else if (!oldinexact)
02971         clear_inexact();
02972 #endif
02973 #ifdef Avoid_Underflow
02974     if (scale) {
02975         word0(rv0) = Exp_1 - 2*P*Exp_msk1;
02976         word1(rv0) = 0;
02977         dval(rv) *= dval(rv0);
02978 #ifndef NO_ERRNO
02979         /* try to avoid the bug of testing an 8087 register value */
02980         if (word0(rv) == 0 && word1(rv) == 0)
02981             errno = ERANGE;
02982 #endif
02983     }
02984 #endif /* Avoid_Underflow */
02985 #ifdef SET_INEXACT
02986     if (inexact && !(word0(rv) & Exp_mask)) {
02987         /* set underflow bit */
02988         dval(rv0) = 1e-300;
02989         dval(rv0) *= dval(rv0);
02990     }
02991 #endif
02992 retfree:
02993     Bfree(bb);
02994     Bfree(bd);
02995     Bfree(bs);
02996     Bfree(bd0);
02997     Bfree(delta);
02998 ret:
02999     if (se)
03000         *se = (char *)s;
03001     return sign ? -dval(rv) : dval(rv);
03002 }
03003 
03004 static int
03005 quorem(Bigint *b, Bigint *S)
03006 {
03007     int n;
03008     ULong *bx, *bxe, q, *sx, *sxe;
03009 #ifdef ULLong
03010     ULLong borrow, carry, y, ys;
03011 #else
03012     ULong borrow, carry, y, ys;
03013 #ifdef Pack_32
03014     ULong si, z, zs;
03015 #endif
03016 #endif
03017 
03018     n = S->wds;
03019 #ifdef DEBUG
03020     /*debug*/ if (b->wds > n)
03021     /*debug*/   Bug("oversize b in quorem");
03022 #endif
03023     if (b->wds < n)
03024         return 0;
03025     sx = S->x;
03026     sxe = sx + --n;
03027     bx = b->x;
03028     bxe = bx + n;
03029     q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
03030 #ifdef DEBUG
03031     /*debug*/ if (q > 9)
03032     /*debug*/   Bug("oversized quotient in quorem");
03033 #endif
03034     if (q) {
03035         borrow = 0;
03036         carry = 0;
03037         do {
03038 #ifdef ULLong
03039             ys = *sx++ * (ULLong)q + carry;
03040             carry = ys >> 32;
03041             y = *bx - (ys & FFFFFFFF) - borrow;
03042             borrow = y >> 32 & (ULong)1;
03043             *bx++ = (ULong)(y & FFFFFFFF);
03044 #else
03045 #ifdef Pack_32
03046             si = *sx++;
03047             ys = (si & 0xffff) * q + carry;
03048             zs = (si >> 16) * q + (ys >> 16);
03049             carry = zs >> 16;
03050             y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
03051             borrow = (y & 0x10000) >> 16;
03052             z = (*bx >> 16) - (zs & 0xffff) - borrow;
03053             borrow = (z & 0x10000) >> 16;
03054             Storeinc(bx, z, y);
03055 #else
03056             ys = *sx++ * q + carry;
03057             carry = ys >> 16;
03058             y = *bx - (ys & 0xffff) - borrow;
03059             borrow = (y & 0x10000) >> 16;
03060             *bx++ = y & 0xffff;
03061 #endif
03062 #endif
03063         } while (sx <= sxe);
03064         if (!*bxe) {
03065             bx = b->x;
03066             while (--bxe > bx && !*bxe)
03067                 --n;
03068             b->wds = n;
03069         }
03070     }
03071     if (cmp(b, S) >= 0) {
03072         q++;
03073         borrow = 0;
03074         carry = 0;
03075         bx = b->x;
03076         sx = S->x;
03077         do {
03078 #ifdef ULLong
03079             ys = *sx++ + carry;
03080             carry = ys >> 32;
03081             y = *bx - (ys & FFFFFFFF) - borrow;
03082             borrow = y >> 32 & (ULong)1;
03083             *bx++ = (ULong)(y & FFFFFFFF);
03084 #else
03085 #ifdef Pack_32
03086             si = *sx++;
03087             ys = (si & 0xffff) + carry;
03088             zs = (si >> 16) + (ys >> 16);
03089             carry = zs >> 16;
03090             y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
03091             borrow = (y & 0x10000) >> 16;
03092             z = (*bx >> 16) - (zs & 0xffff) - borrow;
03093             borrow = (z & 0x10000) >> 16;
03094             Storeinc(bx, z, y);
03095 #else
03096             ys = *sx++ + carry;
03097             carry = ys >> 16;
03098             y = *bx - (ys & 0xffff) - borrow;
03099             borrow = (y & 0x10000) >> 16;
03100             *bx++ = y & 0xffff;
03101 #endif
03102 #endif
03103         } while (sx <= sxe);
03104         bx = b->x;
03105         bxe = bx + n;
03106         if (!*bxe) {
03107             while (--bxe > bx && !*bxe)
03108                 --n;
03109             b->wds = n;
03110         }
03111     }
03112     return q;
03113 }
03114 
03115 #ifndef MULTIPLE_THREADS
03116 static char *dtoa_result;
03117 #endif
03118 
03119 #ifndef MULTIPLE_THREADS
03120 static char *
03121 rv_alloc(int i)
03122 {
03123     return dtoa_result = xmalloc(i);
03124 }
03125 #else
03126 #define rv_alloc(i) xmalloc(i)
03127 #endif
03128 
03129 static char *
03130 nrv_alloc(const char *s, char **rve, size_t n)
03131 {
03132     char *rv, *t;
03133 
03134     t = rv = rv_alloc(n);
03135     while ((*t = *s++) != 0) t++;
03136     if (rve)
03137         *rve = t;
03138     return rv;
03139 }
03140 
03141 #define rv_strdup(s, rve) nrv_alloc(s, rve, strlen(s)+1)
03142 
03143 #ifndef MULTIPLE_THREADS
03144 /* freedtoa(s) must be used to free values s returned by dtoa
03145  * when MULTIPLE_THREADS is #defined.  It should be used in all cases,
03146  * but for consistency with earlier versions of dtoa, it is optional
03147  * when MULTIPLE_THREADS is not defined.
03148  */
03149 
03150 static void
03151 freedtoa(char *s)
03152 {
03153     xfree(s);
03154 }
03155 #endif
03156 
03157 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
03158  *
03159  * Inspired by "How to Print Floating-Point Numbers Accurately" by
03160  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
03161  *
03162  * Modifications:
03163  *  1. Rather than iterating, we use a simple numeric overestimate
03164  *     to determine k = floor(log10(d)).  We scale relevant
03165  *     quantities using O(log2(k)) rather than O(k) multiplications.
03166  *  2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
03167  *     try to generate digits strictly left to right.  Instead, we
03168  *     compute with fewer bits and propagate the carry if necessary
03169  *     when rounding the final digit up.  This is often faster.
03170  *  3. Under the assumption that input will be rounded nearest,
03171  *     mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
03172  *     That is, we allow equality in stopping tests when the
03173  *     round-nearest rule will give the same floating-point value
03174  *     as would satisfaction of the stopping test with strict
03175  *     inequality.
03176  *  4. We remove common factors of powers of 2 from relevant
03177  *     quantities.
03178  *  5. When converting floating-point integers less than 1e16,
03179  *     we use floating-point arithmetic rather than resorting
03180  *     to multiple-precision integers.
03181  *  6. When asked to produce fewer than 15 digits, we first try
03182  *     to get by with floating-point arithmetic; we resort to
03183  *     multiple-precision integer arithmetic only if we cannot
03184  *     guarantee that the floating-point calculation has given
03185  *     the correctly rounded result.  For k requested digits and
03186  *     "uniformly" distributed input, the probability is
03187  *     something like 10^(k-15) that we must resort to the Long
03188  *     calculation.
03189  */
03190 
03191 char *
03192 ruby_dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
03193 {
03194  /* Arguments ndigits, decpt, sign are similar to those
03195     of ecvt and fcvt; trailing zeros are suppressed from
03196     the returned string.  If not null, *rve is set to point
03197     to the end of the return value.  If d is +-Infinity or NaN,
03198     then *decpt is set to 9999.
03199 
03200     mode:
03201         0 ==> shortest string that yields d when read in
03202             and rounded to nearest.
03203         1 ==> like 0, but with Steele & White stopping rule;
03204             e.g. with IEEE P754 arithmetic , mode 0 gives
03205             1e23 whereas mode 1 gives 9.999999999999999e22.
03206         2 ==> max(1,ndigits) significant digits.  This gives a
03207             return value similar to that of ecvt, except
03208             that trailing zeros are suppressed.
03209         3 ==> through ndigits past the decimal point.  This
03210             gives a return value similar to that from fcvt,
03211             except that trailing zeros are suppressed, and
03212             ndigits can be negative.
03213         4,5 ==> similar to 2 and 3, respectively, but (in
03214             round-nearest mode) with the tests of mode 0 to
03215             possibly return a shorter string that rounds to d.
03216             With IEEE arithmetic and compilation with
03217             -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
03218             as modes 2 and 3 when FLT_ROUNDS != 1.
03219         6-9 ==> Debugging modes similar to mode - 4:  don't try
03220             fast floating-point estimate (if applicable).
03221 
03222         Values of mode other than 0-9 are treated as mode 0.
03223 
03224         Sufficient space is allocated to the return value
03225         to hold the suppressed trailing zeros.
03226     */
03227 
03228     int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
03229         j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
03230         spec_case, try_quick;
03231     Long L;
03232 #ifndef Sudden_Underflow
03233     int denorm;
03234     ULong x;
03235 #endif
03236     Bigint *b, *b1, *delta, *mlo = 0, *mhi = 0, *S;
03237     double ds;
03238     double_u d, d2, eps;
03239     char *s, *s0;
03240 #ifdef Honor_FLT_ROUNDS
03241     int rounding;
03242 #endif
03243 #ifdef SET_INEXACT
03244     int inexact, oldinexact;
03245 #endif
03246 
03247     dval(d) = d_;
03248 
03249 #ifndef MULTIPLE_THREADS
03250     if (dtoa_result) {
03251         freedtoa(dtoa_result);
03252         dtoa_result = 0;
03253     }
03254 #endif
03255 
03256     if (word0(d) & Sign_bit) {
03257         /* set sign for everything, including 0's and NaNs */
03258         *sign = 1;
03259         word0(d) &= ~Sign_bit;  /* clear sign bit */
03260     }
03261     else
03262         *sign = 0;
03263 
03264 #if defined(IEEE_Arith) + defined(VAX)
03265 #ifdef IEEE_Arith
03266     if ((word0(d) & Exp_mask) == Exp_mask)
03267 #else
03268     if (word0(d)  == 0x8000)
03269 #endif
03270     {
03271         /* Infinity or NaN */
03272         *decpt = 9999;
03273 #ifdef IEEE_Arith
03274         if (!word1(d) && !(word0(d) & 0xfffff))
03275             return rv_strdup("Infinity", rve);
03276 #endif
03277         return rv_strdup("NaN", rve);
03278     }
03279 #endif
03280 #ifdef IBM
03281     dval(d) += 0; /* normalize */
03282 #endif
03283     if (!dval(d)) {
03284         *decpt = 1;
03285         return rv_strdup("0", rve);
03286     }
03287 
03288 #ifdef SET_INEXACT
03289     try_quick = oldinexact = get_inexact();
03290     inexact = 1;
03291 #endif
03292 #ifdef Honor_FLT_ROUNDS
03293     if ((rounding = Flt_Rounds) >= 2) {
03294         if (*sign)
03295             rounding = rounding == 2 ? 0 : 2;
03296         else
03297             if (rounding != 2)
03298                 rounding = 0;
03299     }
03300 #endif
03301 
03302     b = d2b(dval(d), &be, &bbits);
03303 #ifdef Sudden_Underflow
03304     i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
03305 #else
03306     if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) {
03307 #endif
03308         dval(d2) = dval(d);
03309         word0(d2) &= Frac_mask1;
03310         word0(d2) |= Exp_11;
03311 #ifdef IBM
03312         if (j = 11 - hi0bits(word0(d2) & Frac_mask))
03313             dval(d2) /= 1 << j;
03314 #endif
03315 
03316         /* log(x)   ~=~ log(1.5) + (x-1.5)/1.5
03317          * log10(x)  =  log(x) / log(10)
03318          *      ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
03319          * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
03320          *
03321          * This suggests computing an approximation k to log10(d) by
03322          *
03323          * k = (i - Bias)*0.301029995663981
03324          *  + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
03325          *
03326          * We want k to be too large rather than too small.
03327          * The error in the first-order Taylor series approximation
03328          * is in our favor, so we just round up the constant enough
03329          * to compensate for any error in the multiplication of
03330          * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
03331          * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
03332          * adding 1e-13 to the constant term more than suffices.
03333          * Hence we adjust the constant term to 0.1760912590558.
03334          * (We could get a more accurate k by invoking log10,
03335          *  but this is probably not worthwhile.)
03336          */
03337 
03338         i -= Bias;
03339 #ifdef IBM
03340         i <<= 2;
03341         i += j;
03342 #endif
03343 #ifndef Sudden_Underflow
03344         denorm = 0;
03345     }
03346     else {
03347         /* d is denormalized */
03348 
03349         i = bbits + be + (Bias + (P-1) - 1);
03350         x = i > 32  ? word0(d) << (64 - i) | word1(d) >> (i - 32)
03351             : word1(d) << (32 - i);
03352         dval(d2) = x;
03353         word0(d2) -= 31*Exp_msk1; /* adjust exponent */
03354         i -= (Bias + (P-1) - 1) + 1;
03355         denorm = 1;
03356     }
03357 #endif
03358     ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
03359     k = (int)ds;
03360     if (ds < 0. && ds != k)
03361         k--;    /* want k = floor(ds) */
03362     k_check = 1;
03363     if (k >= 0 && k <= Ten_pmax) {
03364         if (dval(d) < tens[k])
03365             k--;
03366         k_check = 0;
03367     }
03368     j = bbits - i - 1;
03369     if (j >= 0) {
03370         b2 = 0;
03371         s2 = j;
03372     }
03373     else {
03374         b2 = -j;
03375         s2 = 0;
03376     }
03377     if (k >= 0) {
03378         b5 = 0;
03379         s5 = k;
03380         s2 += k;
03381     }
03382     else {
03383         b2 -= k;
03384         b5 = -k;
03385         s5 = 0;
03386     }
03387     if (mode < 0 || mode > 9)
03388         mode = 0;
03389 
03390 #ifndef SET_INEXACT
03391 #ifdef Check_FLT_ROUNDS
03392     try_quick = Rounding == 1;
03393 #else
03394     try_quick = 1;
03395 #endif
03396 #endif /*SET_INEXACT*/
03397 
03398     if (mode > 5) {
03399         mode -= 4;
03400         try_quick = 0;
03401     }
03402     leftright = 1;
03403     ilim = ilim1 = -1;
03404     switch (mode) {
03405       case 0:
03406       case 1:
03407         i = 18;
03408         ndigits = 0;
03409         break;
03410       case 2:
03411         leftright = 0;
03412         /* no break */
03413       case 4:
03414         if (ndigits <= 0)
03415             ndigits = 1;
03416         ilim = ilim1 = i = ndigits;
03417         break;
03418       case 3:
03419         leftright = 0;
03420         /* no break */
03421       case 5:
03422         i = ndigits + k + 1;
03423         ilim = i;
03424         ilim1 = i - 1;
03425         if (i <= 0)
03426             i = 1;
03427     }
03428     s = s0 = rv_alloc(i+1);
03429 
03430 #ifdef Honor_FLT_ROUNDS
03431     if (mode > 1 && rounding != 1)
03432         leftright = 0;
03433 #endif
03434 
03435     if (ilim >= 0 && ilim <= Quick_max && try_quick) {
03436 
03437         /* Try to get by with floating-point arithmetic. */
03438 
03439         i = 0;
03440         dval(d2) = dval(d);
03441         k0 = k;
03442         ilim0 = ilim;
03443         ieps = 2; /* conservative */
03444         if (k > 0) {
03445             ds = tens[k&0xf];
03446             j = k >> 4;
03447             if (j & Bletch) {
03448                 /* prevent overflows */
03449                 j &= Bletch - 1;
03450                 dval(d) /= bigtens[n_bigtens-1];
03451                 ieps++;
03452             }
03453             for (; j; j >>= 1, i++)
03454                 if (j & 1) {
03455                     ieps++;
03456                     ds *= bigtens[i];
03457                 }
03458             dval(d) /= ds;
03459         }
03460         else if ((j1 = -k) != 0) {
03461             dval(d) *= tens[j1 & 0xf];
03462             for (j = j1 >> 4; j; j >>= 1, i++)
03463                 if (j & 1) {
03464                     ieps++;
03465                     dval(d) *= bigtens[i];
03466                 }
03467         }
03468         if (k_check && dval(d) < 1. && ilim > 0) {
03469             if (ilim1 <= 0)
03470                 goto fast_failed;
03471             ilim = ilim1;
03472             k--;
03473             dval(d) *= 10.;
03474             ieps++;
03475         }
03476         dval(eps) = ieps*dval(d) + 7.;
03477         word0(eps) -= (P-1)*Exp_msk1;
03478         if (ilim == 0) {
03479             S = mhi = 0;
03480             dval(d) -= 5.;
03481             if (dval(d) > dval(eps))
03482                 goto one_digit;
03483             if (dval(d) < -dval(eps))
03484                 goto no_digits;
03485             goto fast_failed;
03486         }
03487 #ifndef No_leftright
03488         if (leftright) {
03489             /* Use Steele & White method of only
03490              * generating digits needed.
03491              */
03492             dval(eps) = 0.5/tens[ilim-1] - dval(eps);
03493             for (i = 0;;) {
03494                 L = (int)dval(d);
03495                 dval(d) -= L;
03496                 *s++ = '0' + (int)L;
03497                 if (dval(d) < dval(eps))
03498                     goto ret1;
03499                 if (1. - dval(d) < dval(eps))
03500                     goto bump_up;
03501                 if (++i >= ilim)
03502                     break;
03503                 dval(eps) *= 10.;
03504                 dval(d) *= 10.;
03505             }
03506         }
03507         else {
03508 #endif
03509             /* Generate ilim digits, then fix them up. */
03510             dval(eps) *= tens[ilim-1];
03511             for (i = 1;; i++, dval(d) *= 10.) {
03512                 L = (Long)(dval(d));
03513                 if (!(dval(d) -= L))
03514                     ilim = i;
03515                 *s++ = '0' + (int)L;
03516                 if (i == ilim) {
03517                     if (dval(d) > 0.5 + dval(eps))
03518                         goto bump_up;
03519                     else if (dval(d) < 0.5 - dval(eps)) {
03520                         while (*--s == '0') ;
03521                         s++;
03522                         goto ret1;
03523                     }
03524                     break;
03525                 }
03526             }
03527 #ifndef No_leftright
03528         }
03529 #endif
03530 fast_failed:
03531         s = s0;
03532         dval(d) = dval(d2);
03533         k = k0;
03534         ilim = ilim0;
03535     }
03536 
03537     /* Do we have a "small" integer? */
03538 
03539     if (be >= 0 && k <= Int_max) {
03540         /* Yes. */
03541         ds = tens[k];
03542         if (ndigits < 0 && ilim <= 0) {
03543             S = mhi = 0;
03544             if (ilim < 0 || dval(d) <= 5*ds)
03545                 goto no_digits;
03546             goto one_digit;
03547         }
03548         for (i = 1;; i++, dval(d) *= 10.) {
03549             L = (Long)(dval(d) / ds);
03550             dval(d) -= L*ds;
03551 #ifdef Check_FLT_ROUNDS
03552             /* If FLT_ROUNDS == 2, L will usually be high by 1 */
03553             if (dval(d) < 0) {
03554                 L--;
03555                 dval(d) += ds;
03556             }
03557 #endif
03558             *s++ = '0' + (int)L;
03559             if (!dval(d)) {
03560 #ifdef SET_INEXACT
03561                 inexact = 0;
03562 #endif
03563                 break;
03564             }
03565             if (i == ilim) {
03566 #ifdef Honor_FLT_ROUNDS
03567                 if (mode > 1)
03568                 switch (rounding) {
03569                   case 0: goto ret1;
03570                   case 2: goto bump_up;
03571                 }
03572 #endif
03573                 dval(d) += dval(d);
03574                 if (dval(d) > ds || (dval(d) == ds && (L & 1))) {
03575 bump_up:
03576                     while (*--s == '9')
03577                         if (s == s0) {
03578                             k++;
03579                             *s = '0';
03580                             break;
03581                         }
03582                     ++*s++;
03583                 }
03584                 break;
03585             }
03586         }
03587         goto ret1;
03588     }
03589 
03590     m2 = b2;
03591     m5 = b5;
03592     if (leftright) {
03593         i =
03594 #ifndef Sudden_Underflow
03595             denorm ? be + (Bias + (P-1) - 1 + 1) :
03596 #endif
03597 #ifdef IBM
03598             1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
03599 #else
03600             1 + P - bbits;
03601 #endif
03602         b2 += i;
03603         s2 += i;
03604         mhi = i2b(1);
03605     }
03606     if (m2 > 0 && s2 > 0) {
03607         i = m2 < s2 ? m2 : s2;
03608         b2 -= i;
03609         m2 -= i;
03610         s2 -= i;
03611     }
03612     if (b5 > 0) {
03613         if (leftright) {
03614             if (m5 > 0) {
03615                 mhi = pow5mult(mhi, m5);
03616                 b1 = mult(mhi, b);
03617                 Bfree(b);
03618                 b = b1;
03619             }
03620             if ((j = b5 - m5) != 0)
03621                 b = pow5mult(b, j);
03622         }
03623         else
03624             b = pow5mult(b, b5);
03625     }
03626     S = i2b(1);
03627     if (s5 > 0)
03628         S = pow5mult(S, s5);
03629 
03630     /* Check for special case that d is a normalized power of 2. */
03631 
03632     spec_case = 0;
03633     if ((mode < 2 || leftright)
03634 #ifdef Honor_FLT_ROUNDS
03635             && rounding == 1
03636 #endif
03637     ) {
03638         if (!word1(d) && !(word0(d) & Bndry_mask)
03639 #ifndef Sudden_Underflow
03640             && word0(d) & (Exp_mask & ~Exp_msk1)
03641 #endif
03642         ) {
03643             /* The special case */
03644             b2 += Log2P;
03645             s2 += Log2P;
03646             spec_case = 1;
03647         }
03648     }
03649 
03650     /* Arrange for convenient computation of quotients:
03651      * shift left if necessary so divisor has 4 leading 0 bits.
03652      *
03653      * Perhaps we should just compute leading 28 bits of S once
03654      * and for all and pass them and a shift to quorem, so it
03655      * can do shifts and ors to compute the numerator for q.
03656      */
03657 #ifdef Pack_32
03658     if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0)
03659         i = 32 - i;
03660 #else
03661     if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) != 0)
03662         i = 16 - i;
03663 #endif
03664     if (i > 4) {
03665         i -= 4;
03666         b2 += i;
03667         m2 += i;
03668         s2 += i;
03669     }
03670     else if (i < 4) {
03671         i += 28;
03672         b2 += i;
03673         m2 += i;
03674         s2 += i;
03675     }
03676     if (b2 > 0)
03677         b = lshift(b, b2);
03678     if (s2 > 0)
03679         S = lshift(S, s2);
03680     if (k_check) {
03681         if (cmp(b,S) < 0) {
03682             k--;
03683             b = multadd(b, 10, 0);  /* we botched the k estimate */
03684             if (leftright)
03685                 mhi = multadd(mhi, 10, 0);
03686             ilim = ilim1;
03687         }
03688     }
03689     if (ilim <= 0 && (mode == 3 || mode == 5)) {
03690         if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
03691             /* no digits, fcvt style */
03692 no_digits:
03693             k = -1 - ndigits;
03694             goto ret;
03695         }
03696 one_digit:
03697         *s++ = '1';
03698         k++;
03699         goto ret;
03700     }
03701     if (leftright) {
03702         if (m2 > 0)
03703             mhi = lshift(mhi, m2);
03704 
03705         /* Compute mlo -- check for special case
03706          * that d is a normalized power of 2.
03707          */
03708 
03709         mlo = mhi;
03710         if (spec_case) {
03711             mhi = Balloc(mhi->k);
03712             Bcopy(mhi, mlo);
03713             mhi = lshift(mhi, Log2P);
03714         }
03715 
03716         for (i = 1;;i++) {
03717             dig = quorem(b,S) + '0';
03718             /* Do we yet have the shortest decimal string
03719              * that will round to d?
03720              */
03721             j = cmp(b, mlo);
03722             delta = diff(S, mhi);
03723             j1 = delta->sign ? 1 : cmp(b, delta);
03724             Bfree(delta);
03725 #ifndef ROUND_BIASED
03726             if (j1 == 0 && mode != 1 && !(word1(d) & 1)
03727 #ifdef Honor_FLT_ROUNDS
03728                 && rounding >= 1
03729 #endif
03730             ) {
03731                 if (dig == '9')
03732                     goto round_9_up;
03733                 if (j > 0)
03734                     dig++;
03735 #ifdef SET_INEXACT
03736                 else if (!b->x[0] && b->wds <= 1)
03737                     inexact = 0;
03738 #endif
03739                 *s++ = dig;
03740                 goto ret;
03741             }
03742 #endif
03743             if (j < 0 || (j == 0 && mode != 1
03744 #ifndef ROUND_BIASED
03745                 && !(word1(d) & 1)
03746 #endif
03747             )) {
03748                 if (!b->x[0] && b->wds <= 1) {
03749 #ifdef SET_INEXACT
03750                     inexact = 0;
03751 #endif
03752                     goto accept_dig;
03753                 }
03754 #ifdef Honor_FLT_ROUNDS
03755                 if (mode > 1)
03756                     switch (rounding) {
03757                       case 0: goto accept_dig;
03758                       case 2: goto keep_dig;
03759                     }
03760 #endif /*Honor_FLT_ROUNDS*/
03761                 if (j1 > 0) {
03762                     b = lshift(b, 1);
03763                     j1 = cmp(b, S);
03764                     if ((j1 > 0 || (j1 == 0 && (dig & 1))) && dig++ == '9')
03765                         goto round_9_up;
03766                 }
03767 accept_dig:
03768                 *s++ = dig;
03769                 goto ret;
03770             }
03771             if (j1 > 0) {
03772 #ifdef Honor_FLT_ROUNDS
03773                 if (!rounding)
03774                     goto accept_dig;
03775 #endif
03776                 if (dig == '9') { /* possible if i == 1 */
03777 round_9_up:
03778                     *s++ = '9';
03779                     goto roundoff;
03780                 }
03781                 *s++ = dig + 1;
03782                 goto ret;
03783             }
03784 #ifdef Honor_FLT_ROUNDS
03785 keep_dig:
03786 #endif
03787             *s++ = dig;
03788             if (i == ilim)
03789                 break;
03790             b = multadd(b, 10, 0);
03791             if (mlo == mhi)
03792                 mlo = mhi = multadd(mhi, 10, 0);
03793             else {
03794                 mlo = multadd(mlo, 10, 0);
03795                 mhi = multadd(mhi, 10, 0);
03796             }
03797         }
03798     }
03799     else
03800         for (i = 1;; i++) {
03801             *s++ = dig = quorem(b,S) + '0';
03802             if (!b->x[0] && b->wds <= 1) {
03803 #ifdef SET_INEXACT
03804                 inexact = 0;
03805 #endif
03806                 goto ret;
03807             }
03808             if (i >= ilim)
03809                 break;
03810             b = multadd(b, 10, 0);
03811         }
03812 
03813     /* Round off last digit */
03814 
03815 #ifdef Honor_FLT_ROUNDS
03816     switch (rounding) {
03817       case 0: goto trimzeros;
03818       case 2: goto roundoff;
03819     }
03820 #endif
03821     b = lshift(b, 1);
03822     j = cmp(b, S);
03823     if (j > 0 || (j == 0 && (dig & 1))) {
03824  roundoff:
03825         while (*--s == '9')
03826             if (s == s0) {
03827                 k++;
03828                 *s++ = '1';
03829                 goto ret;
03830             }
03831         ++*s++;
03832     }
03833     else {
03834         while (*--s == '0') ;
03835         s++;
03836     }
03837 ret:
03838     Bfree(S);
03839     if (mhi) {
03840         if (mlo && mlo != mhi)
03841             Bfree(mlo);
03842         Bfree(mhi);
03843     }
03844 ret1:
03845 #ifdef SET_INEXACT
03846     if (inexact) {
03847         if (!oldinexact) {
03848             word0(d) = Exp_1 + (70 << Exp_shift);
03849             word1(d) = 0;
03850             dval(d) += 1.;
03851         }
03852     }
03853     else if (!oldinexact)
03854         clear_inexact();
03855 #endif
03856     Bfree(b);
03857     *s = 0;
03858     *decpt = k + 1;
03859     if (rve)
03860         *rve = s;
03861     return s0;
03862 }
03863 
03864 void
03865 ruby_each_words(const char *str, void (*func)(const char*, int, void*), void *arg)
03866 {
03867     const char *end;
03868     int len;
03869 
03870     if (!str) return;
03871     for (; *str; str = end) {
03872         while (ISSPACE(*str) || *str == ',') str++;
03873         if (!*str) break;
03874         end = str;
03875         while (*end && !ISSPACE(*end) && *end != ',') end++;
03876         len = (int)(end - str); /* assume no string exceeds INT_MAX */
03877         (*func)(str, len, arg);
03878     }
03879 }
03880 
03881 /*-
03882  * Copyright (c) 2004-2008 David Schultz <das@FreeBSD.ORG>
03883  * All rights reserved.
03884  *
03885  * Redistribution and use in source and binary forms, with or without
03886  * modification, are permitted provided that the following conditions
03887  * are met:
03888  * 1. Redistributions of source code must retain the above copyright
03889  *    notice, this list of conditions and the following disclaimer.
03890  * 2. Redistributions in binary form must reproduce the above copyright
03891  *    notice, this list of conditions and the following disclaimer in the
03892  *    documentation and/or other materials provided with the distribution.
03893  *
03894  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
03895  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
03896  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
03897  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
03898  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
03899  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
03900  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
03901  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
03902  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
03903  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
03904  * SUCH DAMAGE.
03905  */
03906 
03907 #define DBL_MANH_SIZE   20
03908 #define DBL_MANL_SIZE   32
03909 #define INFSTR  "Infinity"
03910 #define NANSTR  "NaN"
03911 #define DBL_ADJ (DBL_MAX_EXP - 2)
03912 #define SIGFIGS ((DBL_MANT_DIG + 3) / 4 + 1)
03913 #define dexp_get(u) ((int)(word0(u) >> Exp_shift) & ~Exp_msk1)
03914 #define dexp_set(u,v) (word0(u) = (((int)(word0(u)) & ~Exp_mask) | (v << Exp_shift)))
03915 #define dmanh_get(u) ((int)(word0(u) & Frac_mask))
03916 #define dmanl_get(u) ((int)word1(u))
03917 
03918 
03919 /*
03920  * This procedure converts a double-precision number in IEEE format
03921  * into a string of hexadecimal digits and an exponent of 2.  Its
03922  * behavior is bug-for-bug compatible with dtoa() in mode 2, with the
03923  * following exceptions:
03924  *
03925  * - An ndigits < 0 causes it to use as many digits as necessary to
03926  *   represent the number exactly.
03927  * - The additional xdigs argument should point to either the string
03928  *   "0123456789ABCDEF" or the string "0123456789abcdef", depending on
03929  *   which case is desired.
03930  * - This routine does not repeat dtoa's mistake of setting decpt
03931  *   to 9999 in the case of an infinity or NaN.  INT_MAX is used
03932  *   for this purpose instead.
03933  *
03934  * Note that the C99 standard does not specify what the leading digit
03935  * should be for non-zero numbers.  For instance, 0x1.3p3 is the same
03936  * as 0x2.6p2 is the same as 0x4.cp3.  This implementation always makes
03937  * the leading digit a 1. This ensures that the exponent printed is the
03938  * actual base-2 exponent, i.e., ilogb(d).
03939  *
03940  * Inputs:      d, xdigs, ndigits
03941  * Outputs:     decpt, sign, rve
03942  */
03943 char *
03944 BSD__hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
03945     char **rve)
03946 {
03947         U u;
03948         char *s, *s0;
03949         int bufsize;
03950         uint32_t manh, manl;
03951 
03952         u.d = d;
03953         if (word0(u) & Sign_bit) {
03954             /* set sign for everything, including 0's and NaNs */
03955             *sign = 1;
03956             word0(u) &= ~Sign_bit;  /* clear sign bit */
03957         }
03958         else
03959             *sign = 0;
03960 
03961         if (isinf(d)) { /* FP_INFINITE */
03962             *decpt = INT_MAX;
03963             return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
03964         }
03965         else if (isnan(d)) { /* FP_NAN */
03966             *decpt = INT_MAX;
03967             return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
03968         }
03969         else if (d == 0.0) { /* FP_ZERO */
03970             *decpt = 1;
03971             return (nrv_alloc("0", rve, 1));
03972         }
03973         else if (dexp_get(u)) { /* FP_NORMAL */
03974             *decpt = dexp_get(u) - DBL_ADJ;
03975         }
03976         else { /* FP_SUBNORMAL */
03977             u.d *= 5.363123171977039e+154 /* 0x1p514 */;
03978             *decpt = dexp_get(u) - (514 + DBL_ADJ);
03979         }
03980 
03981         if (ndigits == 0)               /* dtoa() compatibility */
03982                 ndigits = 1;
03983 
03984         /*
03985          * If ndigits < 0, we are expected to auto-size, so we allocate
03986          * enough space for all the digits.
03987          */
03988         bufsize = (ndigits > 0) ? ndigits : SIGFIGS;
03989         s0 = rv_alloc(bufsize);
03990 
03991         /* Round to the desired number of digits. */
03992         if (SIGFIGS > ndigits && ndigits > 0) {
03993                 float redux = 1.0f;
03994                 int offset = 4 * ndigits + DBL_MAX_EXP - 4 - DBL_MANT_DIG;
03995                 dexp_set(u, offset);
03996                 u.d += redux;
03997                 u.d -= redux;
03998                 *decpt += dexp_get(u) - offset;
03999         }
04000 
04001         manh = dmanh_get(u);
04002         manl = dmanl_get(u);
04003         *s0 = '1';
04004         for (s = s0 + 1; s < s0 + bufsize; s++) {
04005                 *s = xdigs[(manh >> (DBL_MANH_SIZE - 4)) & 0xf];
04006                 manh = (manh << 4) | (manl >> (DBL_MANL_SIZE - 4));
04007                 manl <<= 4;
04008         }
04009 
04010         /* If ndigits < 0, we are expected to auto-size the precision. */
04011         if (ndigits < 0) {
04012                 for (ndigits = SIGFIGS; s0[ndigits - 1] == '0'; ndigits--)
04013                         ;
04014         }
04015 
04016         s = s0 + ndigits;
04017         *s = '\0';
04018         if (rve != NULL)
04019                 *rve = s;
04020         return (s0);
04021 }
04022 
04023 #ifdef __cplusplus
04024 #if 0
04025 {
04026 #endif
04027 }
04028 #endif
04029 

Generated on Thu Sep 8 2011 03:50:48 for Ruby by  doxygen 1.7.1