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