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