Ruby 1.9.2p136(2010-12-25revision30365)

- 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