Ruby  2.1.3p242(2014-09-19revision47630)
util.c
Go to the documentation of this file.
1 /**********************************************************************
2 
3  util.c -
4 
5  $Author: nobu $
6  created at: Fri Mar 10 17:22:34 JST 1995
7 
8  Copyright (C) 1993-2008 Yukihiro Matsumoto
9 
10 **********************************************************************/
11 
12 #include "ruby/ruby.h"
13 #include "internal.h"
14 
15 #include <ctype.h>
16 #include <stdio.h>
17 #include <errno.h>
18 #include <math.h>
19 #include <float.h>
20 
21 #ifdef _WIN32
22 #include "missing/file.h"
23 #endif
24 
25 #include "ruby/util.h"
26 
27 unsigned long
28 ruby_scan_oct(const char *start, size_t len, size_t *retlen)
29 {
30  register const char *s = start;
31  register unsigned long retval = 0;
32 
33  while (len-- && *s >= '0' && *s <= '7') {
34  retval <<= 3;
35  retval |= *s++ - '0';
36  }
37  *retlen = (int)(s - start); /* less than len */
38  return retval;
39 }
40 
41 unsigned long
42 ruby_scan_hex(const char *start, size_t len, size_t *retlen)
43 {
44  static const char hexdigit[] = "0123456789abcdef0123456789ABCDEF";
45  register const char *s = start;
46  register unsigned long retval = 0;
47  const char *tmp;
48 
49  while (len-- && *s && (tmp = strchr(hexdigit, *s))) {
50  retval <<= 4;
51  retval |= (tmp - hexdigit) & 15;
52  s++;
53  }
54  *retlen = (int)(s - start); /* less than len */
55  return retval;
56 }
57 
58 const signed char ruby_digit36_to_number_table[] = {
59  /* 0 1 2 3 4 5 6 7 8 9 a b c d e f */
60  /*0*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
61  /*1*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
62  /*2*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
63  /*3*/ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,-1,-1,-1,-1,-1,-1,
64  /*4*/ -1,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,
65  /*5*/ 25,26,27,28,29,30,31,32,33,34,35,-1,-1,-1,-1,-1,
66  /*6*/ -1,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,
67  /*7*/ 25,26,27,28,29,30,31,32,33,34,35,-1,-1,-1,-1,-1,
68  /*8*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
69  /*9*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
70  /*a*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
71  /*b*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
72  /*c*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
73  /*d*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
74  /*e*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
75  /*f*/ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
76 };
77 
78 static unsigned long
79 scan_digits(const char *str, int base, size_t *retlen, int *overflow)
80 {
81 
82  const char *start = str;
83  unsigned long ret = 0, x;
84  unsigned long mul_overflow = (~(unsigned long)0) / base;
85  int c;
86  *overflow = 0;
87 
88  while ((c = (unsigned char)*str++) != '\0') {
90  if (d == -1 || base <= d) {
91  *retlen = (str-1) - start;
92  return ret;
93  }
94  if (mul_overflow < ret)
95  *overflow = 1;
96  ret *= base;
97  x = ret;
98  ret += d;
99  if (ret < x)
100  *overflow = 1;
101  }
102  *retlen = (str-1) - start;
103  return ret;
104 }
105 
106 unsigned long
107 ruby_strtoul(const char *str, char **endptr, int base)
108 {
109  int c, b, overflow;
110  int sign = 0;
111  size_t len;
112  unsigned long ret;
113  const char *subject_found = str;
114 
115  if (base == 1 || 36 < base) {
116  errno = EINVAL;
117  return 0;
118  }
119 
120  while ((c = *str) && ISSPACE(c))
121  str++;
122 
123  if (c == '+') {
124  sign = 1;
125  str++;
126  }
127  else if (c == '-') {
128  sign = -1;
129  str++;
130  }
131 
132  if (str[0] == '0') {
133  subject_found = str+1;
134  if (base == 0 || base == 16) {
135  if (str[1] == 'x' || str[1] == 'X') {
136  b = 16;
137  str += 2;
138  }
139  else {
140  b = base == 0 ? 8 : 16;
141  str++;
142  }
143  }
144  else {
145  b = base;
146  str++;
147  }
148  }
149  else {
150  b = base == 0 ? 10 : base;
151  }
152 
153  ret = scan_digits(str, b, &len, &overflow);
154 
155  if (0 < len)
156  subject_found = str+len;
157 
158  if (endptr)
159  *endptr = (char*)subject_found;
160 
161  if (overflow) {
162  errno = ERANGE;
163  return ULONG_MAX;
164  }
165 
166  if (sign < 0) {
167  ret = (unsigned long)(-(long)ret);
168  return ret;
169  }
170  else {
171  return ret;
172  }
173 }
174 
175 #include <sys/types.h>
176 #include <sys/stat.h>
177 #ifdef HAVE_UNISTD_H
178 #include <unistd.h>
179 #endif
180 #if defined(HAVE_FCNTL_H)
181 #include <fcntl.h>
182 #endif
183 
184 #ifndef S_ISDIR
185 # define S_ISDIR(m) (((m) & S_IFMT) == S_IFDIR)
186 #endif
187 
188 
189 /* mm.c */
190 
191 #define mmtype long
192 #define mmcount (16 / SIZEOF_LONG)
193 #define A ((mmtype*)a)
194 #define B ((mmtype*)b)
195 #define C ((mmtype*)c)
196 #define D ((mmtype*)d)
197 
198 #define mmstep (sizeof(mmtype) * mmcount)
199 #define mmprepare(base, size) do {\
200  if (((VALUE)(base) % sizeof(mmtype)) == 0 && ((size) % sizeof(mmtype)) == 0) \
201  if ((size) >= mmstep) mmkind = 1;\
202  else mmkind = 0;\
203  else mmkind = -1;\
204  high = ((size) / mmstep) * mmstep;\
205  low = ((size) % mmstep);\
206 } while (0)\
207 
208 #define mmarg mmkind, size, high, low
209 #define mmargdecl int mmkind, size_t size, size_t high, size_t low
210 
211 static void mmswap_(register char *a, register char *b, mmargdecl)
212 {
213  if (a == b) return;
214  if (mmkind >= 0) {
215  register mmtype s;
216 #if mmcount > 1
217  if (mmkind > 0) {
218  register char *t = a + high;
219  do {
220  s = A[0]; A[0] = B[0]; B[0] = s;
221  s = A[1]; A[1] = B[1]; B[1] = s;
222 #if mmcount > 2
223  s = A[2]; A[2] = B[2]; B[2] = s;
224 #if mmcount > 3
225  s = A[3]; A[3] = B[3]; B[3] = s;
226 #endif
227 #endif
228  a += mmstep; b += mmstep;
229  } while (a < t);
230  }
231 #endif
232  if (low != 0) { s = A[0]; A[0] = B[0]; B[0] = s;
233 #if mmcount > 2
234  if (low >= 2 * sizeof(mmtype)) { s = A[1]; A[1] = B[1]; B[1] = s;
235 #if mmcount > 3
236  if (low >= 3 * sizeof(mmtype)) {s = A[2]; A[2] = B[2]; B[2] = s;}
237 #endif
238  }
239 #endif
240  }
241  }
242  else {
243  register char *t = a + size, s;
244  do {s = *a; *a++ = *b; *b++ = s;} while (a < t);
245  }
246 }
247 #define mmswap(a,b) mmswap_((a),(b),mmarg)
248 
249 /* a, b, c = b, c, a */
250 static void mmrot3_(register char *a, register char *b, register char *c, mmargdecl)
251 {
252  if (mmkind >= 0) {
253  register mmtype s;
254 #if mmcount > 1
255  if (mmkind > 0) {
256  register char *t = a + high;
257  do {
258  s = A[0]; A[0] = B[0]; B[0] = C[0]; C[0] = s;
259  s = A[1]; A[1] = B[1]; B[1] = C[1]; C[1] = s;
260 #if mmcount > 2
261  s = A[2]; A[2] = B[2]; B[2] = C[2]; C[2] = s;
262 #if mmcount > 3
263  s = A[3]; A[3] = B[3]; B[3] = C[3]; C[3] = s;
264 #endif
265 #endif
266  a += mmstep; b += mmstep; c += mmstep;
267  } while (a < t);
268  }
269 #endif
270  if (low != 0) { s = A[0]; A[0] = B[0]; B[0] = C[0]; C[0] = s;
271 #if mmcount > 2
272  if (low >= 2 * sizeof(mmtype)) { s = A[1]; A[1] = B[1]; B[1] = C[1]; C[1] = s;
273 #if mmcount > 3
274  if (low == 3 * sizeof(mmtype)) {s = A[2]; A[2] = B[2]; B[2] = C[2]; C[2] = s;}
275 #endif
276  }
277 #endif
278  }
279  }
280  else {
281  register char *t = a + size, s;
282  do {s = *a; *a++ = *b; *b++ = *c; *c++ = s;} while (a < t);
283  }
284 }
285 #define mmrot3(a,b,c) mmrot3_((a),(b),(c),mmarg)
286 
287 /* qs6.c */
288 /*****************************************************/
289 /* */
290 /* qs6 (Quick sort function) */
291 /* */
292 /* by Tomoyuki Kawamura 1995.4.21 */
293 /* kawamura@tokuyama.ac.jp */
294 /*****************************************************/
295 
296 typedef struct { char *LL, *RR; } stack_node; /* Stack structure for L,l,R,r */
297 #define PUSH(ll,rr) do { top->LL = (ll); top->RR = (rr); ++top; } while (0) /* Push L,l,R,r */
298 #define POP(ll,rr) do { --top; (ll) = top->LL; (rr) = top->RR; } while (0) /* Pop L,l,R,r */
299 
300 #define med3(a,b,c) ((*cmp)((a),(b),d)<0 ? \
301  ((*cmp)((b),(c),d)<0 ? (b) : ((*cmp)((a),(c),d)<0 ? (c) : (a))) : \
302  ((*cmp)((b),(c),d)>0 ? (b) : ((*cmp)((a),(c),d)<0 ? (a) : (c))))
303 
304 typedef int (cmpfunc_t)(const void*, const void*, void*);
305 void
306 ruby_qsort(void* base, const size_t nel, const size_t size, cmpfunc_t *cmp, void *d)
307 {
308  register char *l, *r, *m; /* l,r:left,right group m:median point */
309  register int t, eq_l, eq_r; /* eq_l: all items in left group are equal to S */
310  char *L = base; /* left end of current region */
311  char *R = (char*)base + size*(nel-1); /* right end of current region */
312  size_t chklim = 63; /* threshold of ordering element check */
313  enum {size_bits = sizeof(size) * CHAR_BIT};
314  stack_node stack[size_bits]; /* enough for size_t size */
315  stack_node *top = stack;
316  int mmkind;
317  size_t high, low, n;
318 
319  if (nel <= 1) return; /* need not to sort */
320  mmprepare(base, size);
321  goto start;
322 
323  nxt:
324  if (stack == top) return; /* return if stack is empty */
325  POP(L,R);
326 
327  for (;;) {
328  start:
329  if (L + size == R) { /* 2 elements */
330  if ((*cmp)(L,R,d) > 0) mmswap(L,R); goto nxt;
331  }
332 
333  l = L; r = R;
334  n = (r - l + size) / size; /* number of elements */
335  m = l + size * (n >> 1); /* calculate median value */
336 
337  if (n >= 60) {
338  register char *m1;
339  register char *m3;
340  if (n >= 200) {
341  n = size*(n>>3); /* number of bytes in splitting 8 */
342  {
343  register char *p1 = l + n;
344  register char *p2 = p1 + n;
345  register char *p3 = p2 + n;
346  m1 = med3(p1, p2, p3);
347  p1 = m + n;
348  p2 = p1 + n;
349  p3 = p2 + n;
350  m3 = med3(p1, p2, p3);
351  }
352  }
353  else {
354  n = size*(n>>2); /* number of bytes in splitting 4 */
355  m1 = l + n;
356  m3 = m + n;
357  }
358  m = med3(m1, m, m3);
359  }
360 
361  if ((t = (*cmp)(l,m,d)) < 0) { /*3-5-?*/
362  if ((t = (*cmp)(m,r,d)) < 0) { /*3-5-7*/
363  if (chklim && nel >= chklim) { /* check if already ascending order */
364  char *p;
365  chklim = 0;
366  for (p=l; p<r; p+=size) if ((*cmp)(p,p+size,d) > 0) goto fail;
367  goto nxt;
368  }
369  fail: goto loopA; /*3-5-7*/
370  }
371  if (t > 0) {
372  if ((*cmp)(l,r,d) <= 0) {mmswap(m,r); goto loopA;} /*3-5-4*/
373  mmrot3(r,m,l); goto loopA; /*3-5-2*/
374  }
375  goto loopB; /*3-5-5*/
376  }
377 
378  if (t > 0) { /*7-5-?*/
379  if ((t = (*cmp)(m,r,d)) > 0) { /*7-5-3*/
380  if (chklim && nel >= chklim) { /* check if already ascending order */
381  char *p;
382  chklim = 0;
383  for (p=l; p<r; p+=size) if ((*cmp)(p,p+size,d) < 0) goto fail2;
384  while (l<r) {mmswap(l,r); l+=size; r-=size;} /* reverse region */
385  goto nxt;
386  }
387  fail2: mmswap(l,r); goto loopA; /*7-5-3*/
388  }
389  if (t < 0) {
390  if ((*cmp)(l,r,d) <= 0) {mmswap(l,m); goto loopB;} /*7-5-8*/
391  mmrot3(l,m,r); goto loopA; /*7-5-6*/
392  }
393  mmswap(l,r); goto loopA; /*7-5-5*/
394  }
395 
396  if ((t = (*cmp)(m,r,d)) < 0) {goto loopA;} /*5-5-7*/
397  if (t > 0) {mmswap(l,r); goto loopB;} /*5-5-3*/
398 
399  /* determining splitting type in case 5-5-5 */ /*5-5-5*/
400  for (;;) {
401  if ((l += size) == r) goto nxt; /*5-5-5*/
402  if (l == m) continue;
403  if ((t = (*cmp)(l,m,d)) > 0) {mmswap(l,r); l = L; goto loopA;}/*575-5*/
404  if (t < 0) {mmswap(L,l); l = L; goto loopB;} /*535-5*/
405  }
406 
407  loopA: eq_l = 1; eq_r = 1; /* splitting type A */ /* left <= median < right */
408  for (;;) {
409  for (;;) {
410  if ((l += size) == r)
411  {l -= size; if (l != m) mmswap(m,l); l -= size; goto fin;}
412  if (l == m) continue;
413  if ((t = (*cmp)(l,m,d)) > 0) {eq_r = 0; break;}
414  if (t < 0) eq_l = 0;
415  }
416  for (;;) {
417  if (l == (r -= size))
418  {l -= size; if (l != m) mmswap(m,l); l -= size; goto fin;}
419  if (r == m) {m = l; break;}
420  if ((t = (*cmp)(r,m,d)) < 0) {eq_l = 0; break;}
421  if (t == 0) break;
422  }
423  mmswap(l,r); /* swap left and right */
424  }
425 
426  loopB: eq_l = 1; eq_r = 1; /* splitting type B */ /* left < median <= right */
427  for (;;) {
428  for (;;) {
429  if (l == (r -= size))
430  {r += size; if (r != m) mmswap(r,m); r += size; goto fin;}
431  if (r == m) continue;
432  if ((t = (*cmp)(r,m,d)) < 0) {eq_l = 0; break;}
433  if (t > 0) eq_r = 0;
434  }
435  for (;;) {
436  if ((l += size) == r)
437  {r += size; if (r != m) mmswap(r,m); r += size; goto fin;}
438  if (l == m) {m = r; break;}
439  if ((t = (*cmp)(l,m,d)) > 0) {eq_r = 0; break;}
440  if (t == 0) break;
441  }
442  mmswap(l,r); /* swap left and right */
443  }
444 
445  fin:
446  if (eq_l == 0) /* need to sort left side */
447  if (eq_r == 0) /* need to sort right side */
448  if (l-L < R-r) {PUSH(r,R); R = l;} /* sort left side first */
449  else {PUSH(L,l); L = r;} /* sort right side first */
450  else R = l; /* need to sort left side only */
451  else if (eq_r == 0) L = r; /* need to sort right side only */
452  else goto nxt; /* need not to sort both sides */
453  }
454 }
455 
456 char *
457 ruby_strdup(const char *str)
458 {
459  char *tmp;
460  size_t len = strlen(str) + 1;
461 
462  tmp = xmalloc(len);
463  memcpy(tmp, str, len);
464 
465  return tmp;
466 }
467 
468 #ifdef __native_client__
469 char *
470 ruby_getcwd(void)
471 {
472  char *buf = xmalloc(2);
473  strcpy(buf, ".");
474  return buf;
475 }
476 #else
477 char *
479 {
480 #ifdef HAVE_GETCWD
481  int size = 200;
482  char *buf = xmalloc(size);
483 
484  while (!getcwd(buf, size)) {
485  if (errno != ERANGE) {
486  xfree(buf);
487  rb_sys_fail("getcwd");
488  }
489  size *= 2;
490  buf = xrealloc(buf, size);
491  }
492 #else
493 # ifndef PATH_MAX
494 # define PATH_MAX 8192
495 # endif
496  char *buf = xmalloc(PATH_MAX+1);
497 
498  if (!getwd(buf)) {
499  xfree(buf);
500  rb_sys_fail("getwd");
501  }
502 #endif
503  return buf;
504 }
505 #endif
506 
507 /****************************************************************
508  *
509  * The author of this software is David M. Gay.
510  *
511  * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
512  *
513  * Permission to use, copy, modify, and distribute this software for any
514  * purpose without fee is hereby granted, provided that this entire notice
515  * is included in all copies of any software which is or includes a copy
516  * or modification of this software and in all copies of the supporting
517  * documentation for such software.
518  *
519  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
520  * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
521  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
522  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
523  *
524  ***************************************************************/
525 
526 /* Please send bug reports to David M. Gay (dmg at acm dot org,
527  * with " at " changed at "@" and " dot " changed to "."). */
528 
529 /* On a machine with IEEE extended-precision registers, it is
530  * necessary to specify double-precision (53-bit) rounding precision
531  * before invoking strtod or dtoa. If the machine uses (the equivalent
532  * of) Intel 80x87 arithmetic, the call
533  * _control87(PC_53, MCW_PC);
534  * does this with many compilers. Whether this or another call is
535  * appropriate depends on the compiler; for this to work, it may be
536  * necessary to #include "float.h" or another system-dependent header
537  * file.
538  */
539 
540 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
541  *
542  * This strtod returns a nearest machine number to the input decimal
543  * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
544  * broken by the IEEE round-even rule. Otherwise ties are broken by
545  * biased rounding (add half and chop).
546  *
547  * Inspired loosely by William D. Clinger's paper "How to Read Floating
548  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
549  *
550  * Modifications:
551  *
552  * 1. We only require IEEE, IBM, or VAX double-precision
553  * arithmetic (not IEEE double-extended).
554  * 2. We get by with floating-point arithmetic in a case that
555  * Clinger missed -- when we're computing d * 10^n
556  * for a small integer d and the integer n is not too
557  * much larger than 22 (the maximum integer k for which
558  * we can represent 10^k exactly), we may be able to
559  * compute (d*10^k) * 10^(e-k) with just one roundoff.
560  * 3. Rather than a bit-at-a-time adjustment of the binary
561  * result in the hard case, we use floating-point
562  * arithmetic to determine the adjustment to within
563  * one bit; only in really hard cases do we need to
564  * compute a second residual.
565  * 4. Because of 3., we don't need a large table of powers of 10
566  * for ten-to-e (just some small tables, e.g. of 10^k
567  * for 0 <= k <= 22).
568  */
569 
570 /*
571  * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least
572  * significant byte has the lowest address.
573  * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most
574  * significant byte has the lowest address.
575  * #define Long int on machines with 32-bit ints and 64-bit longs.
576  * #define IBM for IBM mainframe-style floating-point arithmetic.
577  * #define VAX for VAX-style floating-point arithmetic (D_floating).
578  * #define No_leftright to omit left-right logic in fast floating-point
579  * computation of dtoa.
580  * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
581  * and strtod and dtoa should round accordingly.
582  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
583  * and Honor_FLT_ROUNDS is not #defined.
584  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
585  * that use extended-precision instructions to compute rounded
586  * products and quotients) with IBM.
587  * #define ROUND_BIASED for IEEE-format with biased rounding.
588  * #define Inaccurate_Divide for IEEE-format with correctly rounded
589  * products but inaccurate quotients, e.g., for Intel i860.
590  * #define NO_LONG_LONG on machines that do not have a "long long"
591  * integer type (of >= 64 bits). On such machines, you can
592  * #define Just_16 to store 16 bits per 32-bit Long when doing
593  * high-precision integer arithmetic. Whether this speeds things
594  * up or slows things down depends on the machine and the number
595  * being converted. If long long is available and the name is
596  * something other than "long long", #define Llong to be the name,
597  * and if "unsigned Llong" does not work as an unsigned version of
598  * Llong, #define #ULLong to be the corresponding unsigned type.
599  * #define KR_headers for old-style C function headers.
600  * #define Bad_float_h if your system lacks a float.h or if it does not
601  * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
602  * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
603  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
604  * if memory is available and otherwise does something you deem
605  * appropriate. If MALLOC is undefined, malloc will be invoked
606  * directly -- and assumed always to succeed.
607  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
608  * memory allocations from a private pool of memory when possible.
609  * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
610  * unless #defined to be a different length. This default length
611  * suffices to get rid of MALLOC calls except for unusual cases,
612  * such as decimal-to-binary conversion of a very long string of
613  * digits. The longest string dtoa can return is about 751 bytes
614  * long. For conversions by strtod of strings of 800 digits and
615  * all dtoa conversions in single-threaded executions with 8-byte
616  * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
617  * pointers, PRIVATE_MEM >= 7112 appears adequate.
618  * #define INFNAN_CHECK on IEEE systems to cause strtod to check for
619  * Infinity and NaN (case insensitively). On some systems (e.g.,
620  * some HP systems), it may be necessary to #define NAN_WORD0
621  * appropriately -- to the most significant word of a quiet NaN.
622  * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
623  * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
624  * strtod also accepts (case insensitively) strings of the form
625  * NaN(x), where x is a string of hexadecimal digits and spaces;
626  * if there is only one string of hexadecimal digits, it is taken
627  * for the 52 fraction bits of the resulting NaN; if there are two
628  * or more strings of hex digits, the first is for the high 20 bits,
629  * the second and subsequent for the low 32 bits, with intervening
630  * white space ignored; but if this results in none of the 52
631  * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
632  * and NAN_WORD1 are used instead.
633  * #define MULTIPLE_THREADS if the system offers preemptively scheduled
634  * multiple threads. In this case, you must provide (or suitably
635  * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
636  * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
637  * in pow5mult, ensures lazy evaluation of only one copy of high
638  * powers of 5; omitting this lock would introduce a small
639  * probability of wasting memory, but would otherwise be harmless.)
640  * You must also invoke freedtoa(s) to free the value s returned by
641  * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
642  * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
643  * avoids underflows on inputs whose result does not underflow.
644  * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
645  * floating-point numbers and flushes underflows to zero rather
646  * than implementing gradual underflow, then you must also #define
647  * Sudden_Underflow.
648  * #define YES_ALIAS to permit aliasing certain double values with
649  * arrays of ULongs. This leads to slightly better code with
650  * some compilers and was always used prior to 19990916, but it
651  * is not strictly legal and can cause trouble with aggressively
652  * optimizing compilers (e.g., gcc 2.95.1 under -O2).
653  * #define USE_LOCALE to use the current locale's decimal_point value.
654  * #define SET_INEXACT if IEEE arithmetic is being used and extra
655  * computation should be done to set the inexact flag when the
656  * result is inexact and avoid setting inexact when the result
657  * is exact. In this case, dtoa.c must be compiled in
658  * an environment, perhaps provided by #include "dtoa.c" in a
659  * suitable wrapper, that defines two functions,
660  * int get_inexact(void);
661  * void clear_inexact(void);
662  * such that get_inexact() returns a nonzero value if the
663  * inexact bit is already set, and clear_inexact() sets the
664  * inexact bit to 0. When SET_INEXACT is #defined, strtod
665  * also does extra computations to set the underflow and overflow
666  * flags when appropriate (i.e., when the result is tiny and
667  * inexact or when it is a numeric value rounded to +-infinity).
668  * #define NO_ERRNO if strtod should not assign errno = ERANGE when
669  * the result overflows to +-Infinity or underflows to 0.
670  */
671 
672 #ifdef WORDS_BIGENDIAN
673 #define IEEE_BIG_ENDIAN
674 #else
675 #define IEEE_LITTLE_ENDIAN
676 #endif
677 
678 #ifdef __vax__
679 #define VAX
680 #undef IEEE_BIG_ENDIAN
681 #undef IEEE_LITTLE_ENDIAN
682 #endif
683 
684 #if defined(__arm__) && !defined(__VFP_FP__)
685 #define IEEE_BIG_ENDIAN
686 #undef IEEE_LITTLE_ENDIAN
687 #endif
688 
689 #undef Long
690 #undef ULong
691 
692 #if SIZEOF_INT == 4
693 #define Long int
694 #define ULong unsigned int
695 #elif SIZEOF_LONG == 4
696 #define Long long int
697 #define ULong unsigned long int
698 #endif
699 
700 #if HAVE_LONG_LONG
701 #define Llong LONG_LONG
702 #endif
703 
704 #ifdef DEBUG
705 #include "stdio.h"
706 #define Bug(x) {fprintf(stderr, "%s\n", (x)); exit(EXIT_FAILURE);}
707 #endif
708 
709 #include "stdlib.h"
710 #include "string.h"
711 
712 #ifdef USE_LOCALE
713 #include "locale.h"
714 #endif
715 
716 #ifdef MALLOC
717 extern void *MALLOC(size_t);
718 #else
719 #define MALLOC malloc
720 #endif
721 #ifdef FREE
722 extern void FREE(void*);
723 #else
724 #define FREE free
725 #endif
726 
727 #ifndef Omit_Private_Memory
728 #ifndef PRIVATE_MEM
729 #define PRIVATE_MEM 2304
730 #endif
731 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
733 #endif
734 
735 #undef IEEE_Arith
736 #undef Avoid_Underflow
737 #ifdef IEEE_BIG_ENDIAN
738 #define IEEE_Arith
739 #endif
740 #ifdef IEEE_LITTLE_ENDIAN
741 #define IEEE_Arith
742 #endif
743 
744 #ifdef Bad_float_h
745 
746 #ifdef IEEE_Arith
747 #define DBL_DIG 15
748 #define DBL_MAX_10_EXP 308
749 #define DBL_MAX_EXP 1024
750 #define FLT_RADIX 2
751 #endif /*IEEE_Arith*/
752 
753 #ifdef IBM
754 #define DBL_DIG 16
755 #define DBL_MAX_10_EXP 75
756 #define DBL_MAX_EXP 63
757 #define FLT_RADIX 16
758 #define DBL_MAX 7.2370055773322621e+75
759 #endif
760 
761 #ifdef VAX
762 #define DBL_DIG 16
763 #define DBL_MAX_10_EXP 38
764 #define DBL_MAX_EXP 127
765 #define FLT_RADIX 2
766 #define DBL_MAX 1.7014118346046923e+38
767 #endif
768 
769 #ifndef LONG_MAX
770 #define LONG_MAX 2147483647
771 #endif
772 
773 #else /* ifndef Bad_float_h */
774 #include "float.h"
775 #endif /* Bad_float_h */
776 
777 #ifndef __MATH_H__
778 #include "math.h"
779 #endif
780 
781 #ifdef __cplusplus
782 extern "C" {
783 #if 0
784 } /* satisfy cc-mode */
785 #endif
786 #endif
787 
788 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + defined(IBM) != 1
789 Exactly one of IEEE_LITTLE_ENDIAN, IEEE_BIG_ENDIAN, VAX, or IBM should be defined.
790 #endif
791 
792 typedef union { double d; ULong L[2]; } U;
793 
794 #ifdef YES_ALIAS
795 typedef double double_u;
796 # define dval(x) (x)
797 # ifdef IEEE_LITTLE_ENDIAN
798 # define word0(x) (((ULong *)&(x))[1])
799 # define word1(x) (((ULong *)&(x))[0])
800 # else
801 # define word0(x) (((ULong *)&(x))[0])
802 # define word1(x) (((ULong *)&(x))[1])
803 # endif
804 #else
805 typedef U double_u;
806 # ifdef IEEE_LITTLE_ENDIAN
807 # define word0(x) ((x).L[1])
808 # define word1(x) ((x).L[0])
809 # else
810 # define word0(x) ((x).L[0])
811 # define word1(x) ((x).L[1])
812 # endif
813 # define dval(x) ((x).d)
814 #endif
815 
816 /* The following definition of Storeinc is appropriate for MIPS processors.
817  * An alternative that might be better on some machines is
818  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
819  */
820 #if defined(IEEE_LITTLE_ENDIAN) + defined(VAX) + defined(__arm__)
821 #define Storeinc(a,b,c) (((unsigned short *)(a))[1] = (unsigned short)(b), \
822 ((unsigned short *)(a))[0] = (unsigned short)(c), (a)++)
823 #else
824 #define Storeinc(a,b,c) (((unsigned short *)(a))[0] = (unsigned short)(b), \
825 ((unsigned short *)(a))[1] = (unsigned short)(c), (a)++)
826 #endif
827 
828 /* #define P DBL_MANT_DIG */
829 /* Ten_pmax = floor(P*log(2)/log(5)) */
830 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
831 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
832 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
833 
834 #ifdef IEEE_Arith
835 #define Exp_shift 20
836 #define Exp_shift1 20
837 #define Exp_msk1 0x100000
838 #define Exp_msk11 0x100000
839 #define Exp_mask 0x7ff00000
840 #define P 53
841 #define Bias 1023
842 #define Emin (-1022)
843 #define Exp_1 0x3ff00000
844 #define Exp_11 0x3ff00000
845 #define Ebits 11
846 #define Frac_mask 0xfffff
847 #define Frac_mask1 0xfffff
848 #define Ten_pmax 22
849 #define Bletch 0x10
850 #define Bndry_mask 0xfffff
851 #define Bndry_mask1 0xfffff
852 #define LSB 1
853 #define Sign_bit 0x80000000
854 #define Log2P 1
855 #define Tiny0 0
856 #define Tiny1 1
857 #define Quick_max 14
858 #define Int_max 14
859 #ifndef NO_IEEE_Scale
860 #define Avoid_Underflow
861 #ifdef Flush_Denorm /* debugging option */
862 #undef Sudden_Underflow
863 #endif
864 #endif
865 
866 #ifndef Flt_Rounds
867 #ifdef FLT_ROUNDS
868 #define Flt_Rounds FLT_ROUNDS
869 #else
870 #define Flt_Rounds 1
871 #endif
872 #endif /*Flt_Rounds*/
873 
874 #ifdef Honor_FLT_ROUNDS
875 #define Rounding rounding
876 #undef Check_FLT_ROUNDS
877 #define Check_FLT_ROUNDS
878 #else
879 #define Rounding Flt_Rounds
880 #endif
881 
882 #else /* ifndef IEEE_Arith */
883 #undef Check_FLT_ROUNDS
884 #undef Honor_FLT_ROUNDS
885 #undef SET_INEXACT
886 #undef Sudden_Underflow
887 #define Sudden_Underflow
888 #ifdef IBM
889 #undef Flt_Rounds
890 #define Flt_Rounds 0
891 #define Exp_shift 24
892 #define Exp_shift1 24
893 #define Exp_msk1 0x1000000
894 #define Exp_msk11 0x1000000
895 #define Exp_mask 0x7f000000
896 #define P 14
897 #define Bias 65
898 #define Exp_1 0x41000000
899 #define Exp_11 0x41000000
900 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
901 #define Frac_mask 0xffffff
902 #define Frac_mask1 0xffffff
903 #define Bletch 4
904 #define Ten_pmax 22
905 #define Bndry_mask 0xefffff
906 #define Bndry_mask1 0xffffff
907 #define LSB 1
908 #define Sign_bit 0x80000000
909 #define Log2P 4
910 #define Tiny0 0x100000
911 #define Tiny1 0
912 #define Quick_max 14
913 #define Int_max 15
914 #else /* VAX */
915 #undef Flt_Rounds
916 #define Flt_Rounds 1
917 #define Exp_shift 23
918 #define Exp_shift1 7
919 #define Exp_msk1 0x80
920 #define Exp_msk11 0x800000
921 #define Exp_mask 0x7f80
922 #define P 56
923 #define Bias 129
924 #define Exp_1 0x40800000
925 #define Exp_11 0x4080
926 #define Ebits 8
927 #define Frac_mask 0x7fffff
928 #define Frac_mask1 0xffff007f
929 #define Ten_pmax 24
930 #define Bletch 2
931 #define Bndry_mask 0xffff007f
932 #define Bndry_mask1 0xffff007f
933 #define LSB 0x10000
934 #define Sign_bit 0x8000
935 #define Log2P 1
936 #define Tiny0 0x80
937 #define Tiny1 0
938 #define Quick_max 15
939 #define Int_max 15
940 #endif /* IBM, VAX */
941 #endif /* IEEE_Arith */
942 
943 #ifndef IEEE_Arith
944 #define ROUND_BIASED
945 #endif
946 
947 #ifdef RND_PRODQUOT
948 #define rounded_product(a,b) ((a) = rnd_prod((a), (b)))
949 #define rounded_quotient(a,b) ((a) = rnd_quot((a), (b)))
950 extern double rnd_prod(double, double), rnd_quot(double, double);
951 #else
952 #define rounded_product(a,b) ((a) *= (b))
953 #define rounded_quotient(a,b) ((a) /= (b))
954 #endif
955 
956 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
957 #define Big1 0xffffffff
958 
959 #ifndef Pack_32
960 #define Pack_32
961 #endif
962 
963 #define FFFFFFFF 0xffffffffUL
964 
965 #ifdef NO_LONG_LONG
966 #undef ULLong
967 #ifdef Just_16
968 #undef Pack_32
969 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
970  * This makes some inner loops simpler and sometimes saves work
971  * during multiplications, but it often seems to make things slightly
972  * slower. Hence the default is now to store 32 bits per Long.
973  */
974 #endif
975 #else /* long long available */
976 #ifndef Llong
977 #define Llong long long
978 #endif
979 #ifndef ULLong
980 #define ULLong unsigned Llong
981 #endif
982 #endif /* NO_LONG_LONG */
983 
984 #define MULTIPLE_THREADS 1
985 
986 #ifndef MULTIPLE_THREADS
987 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
988 #define FREE_DTOA_LOCK(n) /*nothing*/
989 #else
990 #define ACQUIRE_DTOA_LOCK(n) /*unused right now*/
991 #define FREE_DTOA_LOCK(n) /*unused right now*/
992 #endif
993 
994 #define Kmax 15
995 
996 struct Bigint {
997  struct Bigint *next;
998  int k, maxwds, sign, wds;
999  ULong x[1];
1000 };
1001 
1002 typedef struct Bigint Bigint;
1003 
1004 static Bigint *freelist[Kmax+1];
1005 
1006 static Bigint *
1007 Balloc(int k)
1008 {
1009  int x;
1010  Bigint *rv;
1011 #ifndef Omit_Private_Memory
1012  size_t len;
1013 #endif
1014 
1015  ACQUIRE_DTOA_LOCK(0);
1016  if (k <= Kmax && (rv = freelist[k]) != 0) {
1017  freelist[k] = rv->next;
1018  }
1019  else {
1020  x = 1 << k;
1021 #ifdef Omit_Private_Memory
1022  rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
1023 #else
1024  len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
1025  /sizeof(double);
1026  if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
1027  rv = (Bigint*)pmem_next;
1028  pmem_next += len;
1029  }
1030  else
1031  rv = (Bigint*)MALLOC(len*sizeof(double));
1032 #endif
1033  rv->k = k;
1034  rv->maxwds = x;
1035  }
1036  FREE_DTOA_LOCK(0);
1037  rv->sign = rv->wds = 0;
1038  return rv;
1039 }
1040 
1041 static void
1043 {
1044  if (v) {
1045  if (v->k > Kmax) {
1046  FREE(v);
1047  return;
1048  }
1049  ACQUIRE_DTOA_LOCK(0);
1050  v->next = freelist[v->k];
1051  freelist[v->k] = v;
1052  FREE_DTOA_LOCK(0);
1053  }
1054 }
1055 
1056 #define Bcopy(x,y) memcpy((char *)&(x)->sign, (char *)&(y)->sign, \
1057 (y)->wds*sizeof(Long) + 2*sizeof(int))
1058 
1059 static Bigint *
1060 multadd(Bigint *b, int m, int a) /* multiply by m and add a */
1061 {
1062  int i, wds;
1063  ULong *x;
1064 #ifdef ULLong
1065  ULLong carry, y;
1066 #else
1067  ULong carry, y;
1068 #ifdef Pack_32
1069  ULong xi, z;
1070 #endif
1071 #endif
1072  Bigint *b1;
1073 
1074  wds = b->wds;
1075  x = b->x;
1076  i = 0;
1077  carry = a;
1078  do {
1079 #ifdef ULLong
1080  y = *x * (ULLong)m + carry;
1081  carry = y >> 32;
1082  *x++ = (ULong)(y & FFFFFFFF);
1083 #else
1084 #ifdef Pack_32
1085  xi = *x;
1086  y = (xi & 0xffff) * m + carry;
1087  z = (xi >> 16) * m + (y >> 16);
1088  carry = z >> 16;
1089  *x++ = (z << 16) + (y & 0xffff);
1090 #else
1091  y = *x * m + carry;
1092  carry = y >> 16;
1093  *x++ = y & 0xffff;
1094 #endif
1095 #endif
1096  } while (++i < wds);
1097  if (carry) {
1098  if (wds >= b->maxwds) {
1099  b1 = Balloc(b->k+1);
1100  Bcopy(b1, b);
1101  Bfree(b);
1102  b = b1;
1103  }
1104  b->x[wds++] = (ULong)carry;
1105  b->wds = wds;
1106  }
1107  return b;
1108 }
1109 
1110 static Bigint *
1111 s2b(const char *s, int nd0, int nd, ULong y9)
1112 {
1113  Bigint *b;
1114  int i, k;
1115  Long x, y;
1116 
1117  x = (nd + 8) / 9;
1118  for (k = 0, y = 1; x > y; y <<= 1, k++) ;
1119 #ifdef Pack_32
1120  b = Balloc(k);
1121  b->x[0] = y9;
1122  b->wds = 1;
1123 #else
1124  b = Balloc(k+1);
1125  b->x[0] = y9 & 0xffff;
1126  b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
1127 #endif
1128 
1129  i = 9;
1130  if (9 < nd0) {
1131  s += 9;
1132  do {
1133  b = multadd(b, 10, *s++ - '0');
1134  } while (++i < nd0);
1135  s++;
1136  }
1137  else
1138  s += 10;
1139  for (; i < nd; i++)
1140  b = multadd(b, 10, *s++ - '0');
1141  return b;
1142 }
1143 
1144 static int
1145 hi0bits(register ULong x)
1146 {
1147  register int k = 0;
1148 
1149  if (!(x & 0xffff0000)) {
1150  k = 16;
1151  x <<= 16;
1152  }
1153  if (!(x & 0xff000000)) {
1154  k += 8;
1155  x <<= 8;
1156  }
1157  if (!(x & 0xf0000000)) {
1158  k += 4;
1159  x <<= 4;
1160  }
1161  if (!(x & 0xc0000000)) {
1162  k += 2;
1163  x <<= 2;
1164  }
1165  if (!(x & 0x80000000)) {
1166  k++;
1167  if (!(x & 0x40000000))
1168  return 32;
1169  }
1170  return k;
1171 }
1172 
1173 static int
1174 lo0bits(ULong *y)
1175 {
1176  register int k;
1177  register ULong x = *y;
1178 
1179  if (x & 7) {
1180  if (x & 1)
1181  return 0;
1182  if (x & 2) {
1183  *y = x >> 1;
1184  return 1;
1185  }
1186  *y = x >> 2;
1187  return 2;
1188  }
1189  k = 0;
1190  if (!(x & 0xffff)) {
1191  k = 16;
1192  x >>= 16;
1193  }
1194  if (!(x & 0xff)) {
1195  k += 8;
1196  x >>= 8;
1197  }
1198  if (!(x & 0xf)) {
1199  k += 4;
1200  x >>= 4;
1201  }
1202  if (!(x & 0x3)) {
1203  k += 2;
1204  x >>= 2;
1205  }
1206  if (!(x & 1)) {
1207  k++;
1208  x >>= 1;
1209  if (!x)
1210  return 32;
1211  }
1212  *y = x;
1213  return k;
1214 }
1215 
1216 static Bigint *
1217 i2b(int i)
1218 {
1219  Bigint *b;
1220 
1221  b = Balloc(1);
1222  b->x[0] = i;
1223  b->wds = 1;
1224  return b;
1225 }
1226 
1227 static Bigint *
1229 {
1230  Bigint *c;
1231  int k, wa, wb, wc;
1232  ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
1233  ULong y;
1234 #ifdef ULLong
1235  ULLong carry, z;
1236 #else
1237  ULong carry, z;
1238 #ifdef Pack_32
1239  ULong z2;
1240 #endif
1241 #endif
1242 
1243  if (a->wds < b->wds) {
1244  c = a;
1245  a = b;
1246  b = c;
1247  }
1248  k = a->k;
1249  wa = a->wds;
1250  wb = b->wds;
1251  wc = wa + wb;
1252  if (wc > a->maxwds)
1253  k++;
1254  c = Balloc(k);
1255  for (x = c->x, xa = x + wc; x < xa; x++)
1256  *x = 0;
1257  xa = a->x;
1258  xae = xa + wa;
1259  xb = b->x;
1260  xbe = xb + wb;
1261  xc0 = c->x;
1262 #ifdef ULLong
1263  for (; xb < xbe; xc0++) {
1264  if ((y = *xb++) != 0) {
1265  x = xa;
1266  xc = xc0;
1267  carry = 0;
1268  do {
1269  z = *x++ * (ULLong)y + *xc + carry;
1270  carry = z >> 32;
1271  *xc++ = (ULong)(z & FFFFFFFF);
1272  } while (x < xae);
1273  *xc = (ULong)carry;
1274  }
1275  }
1276 #else
1277 #ifdef Pack_32
1278  for (; xb < xbe; xb++, xc0++) {
1279  if (y = *xb & 0xffff) {
1280  x = xa;
1281  xc = xc0;
1282  carry = 0;
1283  do {
1284  z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
1285  carry = z >> 16;
1286  z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
1287  carry = z2 >> 16;
1288  Storeinc(xc, z2, z);
1289  } while (x < xae);
1290  *xc = (ULong)carry;
1291  }
1292  if (y = *xb >> 16) {
1293  x = xa;
1294  xc = xc0;
1295  carry = 0;
1296  z2 = *xc;
1297  do {
1298  z = (*x & 0xffff) * y + (*xc >> 16) + carry;
1299  carry = z >> 16;
1300  Storeinc(xc, z, z2);
1301  z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
1302  carry = z2 >> 16;
1303  } while (x < xae);
1304  *xc = z2;
1305  }
1306  }
1307 #else
1308  for (; xb < xbe; xc0++) {
1309  if (y = *xb++) {
1310  x = xa;
1311  xc = xc0;
1312  carry = 0;
1313  do {
1314  z = *x++ * y + *xc + carry;
1315  carry = z >> 16;
1316  *xc++ = z & 0xffff;
1317  } while (x < xae);
1318  *xc = (ULong)carry;
1319  }
1320  }
1321 #endif
1322 #endif
1323  for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
1324  c->wds = wc;
1325  return c;
1326 }
1327 
1328 static Bigint *p5s;
1329 
1330 static Bigint *
1332 {
1333  Bigint *b1, *p5, *p51;
1334  int i;
1335  static int p05[3] = { 5, 25, 125 };
1336 
1337  if ((i = k & 3) != 0)
1338  b = multadd(b, p05[i-1], 0);
1339 
1340  if (!(k >>= 2))
1341  return b;
1342  if (!(p5 = p5s)) {
1343  /* first time */
1344 #ifdef MULTIPLE_THREADS
1345  ACQUIRE_DTOA_LOCK(1);
1346  if (!(p5 = p5s)) {
1347  p5 = p5s = i2b(625);
1348  p5->next = 0;
1349  }
1350  FREE_DTOA_LOCK(1);
1351 #else
1352  p5 = p5s = i2b(625);
1353  p5->next = 0;
1354 #endif
1355  }
1356  for (;;) {
1357  if (k & 1) {
1358  b1 = mult(b, p5);
1359  Bfree(b);
1360  b = b1;
1361  }
1362  if (!(k >>= 1))
1363  break;
1364  if (!(p51 = p5->next)) {
1365 #ifdef MULTIPLE_THREADS
1366  ACQUIRE_DTOA_LOCK(1);
1367  if (!(p51 = p5->next)) {
1368  p51 = p5->next = mult(p5,p5);
1369  p51->next = 0;
1370  }
1371  FREE_DTOA_LOCK(1);
1372 #else
1373  p51 = p5->next = mult(p5,p5);
1374  p51->next = 0;
1375 #endif
1376  }
1377  p5 = p51;
1378  }
1379  return b;
1380 }
1381 
1382 static Bigint *
1383 lshift(Bigint *b, int k)
1384 {
1385  int i, k1, n, n1;
1386  Bigint *b1;
1387  ULong *x, *x1, *xe, z;
1388 
1389 #ifdef Pack_32
1390  n = k >> 5;
1391 #else
1392  n = k >> 4;
1393 #endif
1394  k1 = b->k;
1395  n1 = n + b->wds + 1;
1396  for (i = b->maxwds; n1 > i; i <<= 1)
1397  k1++;
1398  b1 = Balloc(k1);
1399  x1 = b1->x;
1400  for (i = 0; i < n; i++)
1401  *x1++ = 0;
1402  x = b->x;
1403  xe = x + b->wds;
1404 #ifdef Pack_32
1405  if (k &= 0x1f) {
1406  k1 = 32 - k;
1407  z = 0;
1408  do {
1409  *x1++ = *x << k | z;
1410  z = *x++ >> k1;
1411  } while (x < xe);
1412  if ((*x1 = z) != 0)
1413  ++n1;
1414  }
1415 #else
1416  if (k &= 0xf) {
1417  k1 = 16 - k;
1418  z = 0;
1419  do {
1420  *x1++ = *x << k & 0xffff | z;
1421  z = *x++ >> k1;
1422  } while (x < xe);
1423  if (*x1 = z)
1424  ++n1;
1425  }
1426 #endif
1427  else
1428  do {
1429  *x1++ = *x++;
1430  } while (x < xe);
1431  b1->wds = n1 - 1;
1432  Bfree(b);
1433  return b1;
1434 }
1435 
1436 static int
1438 {
1439  ULong *xa, *xa0, *xb, *xb0;
1440  int i, j;
1441 
1442  i = a->wds;
1443  j = b->wds;
1444 #ifdef DEBUG
1445  if (i > 1 && !a->x[i-1])
1446  Bug("cmp called with a->x[a->wds-1] == 0");
1447  if (j > 1 && !b->x[j-1])
1448  Bug("cmp called with b->x[b->wds-1] == 0");
1449 #endif
1450  if (i -= j)
1451  return i;
1452  xa0 = a->x;
1453  xa = xa0 + j;
1454  xb0 = b->x;
1455  xb = xb0 + j;
1456  for (;;) {
1457  if (*--xa != *--xb)
1458  return *xa < *xb ? -1 : 1;
1459  if (xa <= xa0)
1460  break;
1461  }
1462  return 0;
1463 }
1464 
1465 static Bigint *
1467 {
1468  Bigint *c;
1469  int i, wa, wb;
1470  ULong *xa, *xae, *xb, *xbe, *xc;
1471 #ifdef ULLong
1472  ULLong borrow, y;
1473 #else
1474  ULong borrow, y;
1475 #ifdef Pack_32
1476  ULong z;
1477 #endif
1478 #endif
1479 
1480  i = cmp(a,b);
1481  if (!i) {
1482  c = Balloc(0);
1483  c->wds = 1;
1484  c->x[0] = 0;
1485  return c;
1486  }
1487  if (i < 0) {
1488  c = a;
1489  a = b;
1490  b = c;
1491  i = 1;
1492  }
1493  else
1494  i = 0;
1495  c = Balloc(a->k);
1496  c->sign = i;
1497  wa = a->wds;
1498  xa = a->x;
1499  xae = xa + wa;
1500  wb = b->wds;
1501  xb = b->x;
1502  xbe = xb + wb;
1503  xc = c->x;
1504  borrow = 0;
1505 #ifdef ULLong
1506  do {
1507  y = (ULLong)*xa++ - *xb++ - borrow;
1508  borrow = y >> 32 & (ULong)1;
1509  *xc++ = (ULong)(y & FFFFFFFF);
1510  } while (xb < xbe);
1511  while (xa < xae) {
1512  y = *xa++ - borrow;
1513  borrow = y >> 32 & (ULong)1;
1514  *xc++ = (ULong)(y & FFFFFFFF);
1515  }
1516 #else
1517 #ifdef Pack_32
1518  do {
1519  y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1520  borrow = (y & 0x10000) >> 16;
1521  z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1522  borrow = (z & 0x10000) >> 16;
1523  Storeinc(xc, z, y);
1524  } while (xb < xbe);
1525  while (xa < xae) {
1526  y = (*xa & 0xffff) - borrow;
1527  borrow = (y & 0x10000) >> 16;
1528  z = (*xa++ >> 16) - borrow;
1529  borrow = (z & 0x10000) >> 16;
1530  Storeinc(xc, z, y);
1531  }
1532 #else
1533  do {
1534  y = *xa++ - *xb++ - borrow;
1535  borrow = (y & 0x10000) >> 16;
1536  *xc++ = y & 0xffff;
1537  } while (xb < xbe);
1538  while (xa < xae) {
1539  y = *xa++ - borrow;
1540  borrow = (y & 0x10000) >> 16;
1541  *xc++ = y & 0xffff;
1542  }
1543 #endif
1544 #endif
1545  while (!*--xc)
1546  wa--;
1547  c->wds = wa;
1548  return c;
1549 }
1550 
1551 static double
1552 ulp(double x_)
1553 {
1554  register Long L;
1555  double_u x, a;
1556  dval(x) = x_;
1557 
1558  L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1559 #ifndef Avoid_Underflow
1560 #ifndef Sudden_Underflow
1561  if (L > 0) {
1562 #endif
1563 #endif
1564 #ifdef IBM
1565  L |= Exp_msk1 >> 4;
1566 #endif
1567  word0(a) = L;
1568  word1(a) = 0;
1569 #ifndef Avoid_Underflow
1570 #ifndef Sudden_Underflow
1571  }
1572  else {
1573  L = -L >> Exp_shift;
1574  if (L < Exp_shift) {
1575  word0(a) = 0x80000 >> L;
1576  word1(a) = 0;
1577  }
1578  else {
1579  word0(a) = 0;
1580  L -= Exp_shift;
1581  word1(a) = L >= 31 ? 1 : 1 << 31 - L;
1582  }
1583  }
1584 #endif
1585 #endif
1586  return dval(a);
1587 }
1588 
1589 static double
1590 b2d(Bigint *a, int *e)
1591 {
1592  ULong *xa, *xa0, w, y, z;
1593  int k;
1594  double_u d;
1595 #ifdef VAX
1596  ULong d0, d1;
1597 #else
1598 #define d0 word0(d)
1599 #define d1 word1(d)
1600 #endif
1601 
1602  xa0 = a->x;
1603  xa = xa0 + a->wds;
1604  y = *--xa;
1605 #ifdef DEBUG
1606  if (!y) Bug("zero y in b2d");
1607 #endif
1608  k = hi0bits(y);
1609  *e = 32 - k;
1610 #ifdef Pack_32
1611  if (k < Ebits) {
1612  d0 = Exp_1 | y >> (Ebits - k);
1613  w = xa > xa0 ? *--xa : 0;
1614  d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1615  goto ret_d;
1616  }
1617  z = xa > xa0 ? *--xa : 0;
1618  if (k -= Ebits) {
1619  d0 = Exp_1 | y << k | z >> (32 - k);
1620  y = xa > xa0 ? *--xa : 0;
1621  d1 = z << k | y >> (32 - k);
1622  }
1623  else {
1624  d0 = Exp_1 | y;
1625  d1 = z;
1626  }
1627 #else
1628  if (k < Ebits + 16) {
1629  z = xa > xa0 ? *--xa : 0;
1630  d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1631  w = xa > xa0 ? *--xa : 0;
1632  y = xa > xa0 ? *--xa : 0;
1633  d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1634  goto ret_d;
1635  }
1636  z = xa > xa0 ? *--xa : 0;
1637  w = xa > xa0 ? *--xa : 0;
1638  k -= Ebits + 16;
1639  d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1640  y = xa > xa0 ? *--xa : 0;
1641  d1 = w << k + 16 | y << k;
1642 #endif
1643 ret_d:
1644 #ifdef VAX
1645  word0(d) = d0 >> 16 | d0 << 16;
1646  word1(d) = d1 >> 16 | d1 << 16;
1647 #else
1648 #undef d0
1649 #undef d1
1650 #endif
1651  return dval(d);
1652 }
1653 
1654 static Bigint *
1655 d2b(double d_, int *e, int *bits)
1656 {
1657  double_u d;
1658  Bigint *b;
1659  int de, k;
1660  ULong *x, y, z;
1661 #ifndef Sudden_Underflow
1662  int i;
1663 #endif
1664 #ifdef VAX
1665  ULong d0, d1;
1666 #endif
1667  dval(d) = d_;
1668 #ifdef VAX
1669  d0 = word0(d) >> 16 | word0(d) << 16;
1670  d1 = word1(d) >> 16 | word1(d) << 16;
1671 #else
1672 #define d0 word0(d)
1673 #define d1 word1(d)
1674 #endif
1675 
1676 #ifdef Pack_32
1677  b = Balloc(1);
1678 #else
1679  b = Balloc(2);
1680 #endif
1681  x = b->x;
1682 
1683  z = d0 & Frac_mask;
1684  d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1685 #ifdef Sudden_Underflow
1686  de = (int)(d0 >> Exp_shift);
1687 #ifndef IBM
1688  z |= Exp_msk11;
1689 #endif
1690 #else
1691  if ((de = (int)(d0 >> Exp_shift)) != 0)
1692  z |= Exp_msk1;
1693 #endif
1694 #ifdef Pack_32
1695  if ((y = d1) != 0) {
1696  if ((k = lo0bits(&y)) != 0) {
1697  x[0] = y | z << (32 - k);
1698  z >>= k;
1699  }
1700  else
1701  x[0] = y;
1702 #ifndef Sudden_Underflow
1703  i =
1704 #endif
1705  b->wds = (x[1] = z) ? 2 : 1;
1706  }
1707  else {
1708 #ifdef DEBUG
1709  if (!z)
1710  Bug("Zero passed to d2b");
1711 #endif
1712  k = lo0bits(&z);
1713  x[0] = z;
1714 #ifndef Sudden_Underflow
1715  i =
1716 #endif
1717  b->wds = 1;
1718  k += 32;
1719  }
1720 #else
1721  if (y = d1) {
1722  if (k = lo0bits(&y))
1723  if (k >= 16) {
1724  x[0] = y | z << 32 - k & 0xffff;
1725  x[1] = z >> k - 16 & 0xffff;
1726  x[2] = z >> k;
1727  i = 2;
1728  }
1729  else {
1730  x[0] = y & 0xffff;
1731  x[1] = y >> 16 | z << 16 - k & 0xffff;
1732  x[2] = z >> k & 0xffff;
1733  x[3] = z >> k+16;
1734  i = 3;
1735  }
1736  else {
1737  x[0] = y & 0xffff;
1738  x[1] = y >> 16;
1739  x[2] = z & 0xffff;
1740  x[3] = z >> 16;
1741  i = 3;
1742  }
1743  }
1744  else {
1745 #ifdef DEBUG
1746  if (!z)
1747  Bug("Zero passed to d2b");
1748 #endif
1749  k = lo0bits(&z);
1750  if (k >= 16) {
1751  x[0] = z;
1752  i = 0;
1753  }
1754  else {
1755  x[0] = z & 0xffff;
1756  x[1] = z >> 16;
1757  i = 1;
1758  }
1759  k += 32;
1760  }
1761  while (!x[i])
1762  --i;
1763  b->wds = i + 1;
1764 #endif
1765 #ifndef Sudden_Underflow
1766  if (de) {
1767 #endif
1768 #ifdef IBM
1769  *e = (de - Bias - (P-1) << 2) + k;
1770  *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1771 #else
1772  *e = de - Bias - (P-1) + k;
1773  *bits = P - k;
1774 #endif
1775 #ifndef Sudden_Underflow
1776  }
1777  else {
1778  *e = de - Bias - (P-1) + 1 + k;
1779 #ifdef Pack_32
1780  *bits = 32*i - hi0bits(x[i-1]);
1781 #else
1782  *bits = (i+2)*16 - hi0bits(x[i]);
1783 #endif
1784  }
1785 #endif
1786  return b;
1787 }
1788 #undef d0
1789 #undef d1
1790 
1791 static double
1793 {
1794  double_u da, db;
1795  int k, ka, kb;
1796 
1797  dval(da) = b2d(a, &ka);
1798  dval(db) = b2d(b, &kb);
1799 #ifdef Pack_32
1800  k = ka - kb + 32*(a->wds - b->wds);
1801 #else
1802  k = ka - kb + 16*(a->wds - b->wds);
1803 #endif
1804 #ifdef IBM
1805  if (k > 0) {
1806  word0(da) += (k >> 2)*Exp_msk1;
1807  if (k &= 3)
1808  dval(da) *= 1 << k;
1809  }
1810  else {
1811  k = -k;
1812  word0(db) += (k >> 2)*Exp_msk1;
1813  if (k &= 3)
1814  dval(db) *= 1 << k;
1815  }
1816 #else
1817  if (k > 0)
1818  word0(da) += k*Exp_msk1;
1819  else {
1820  k = -k;
1821  word0(db) += k*Exp_msk1;
1822  }
1823 #endif
1824  return dval(da) / dval(db);
1825 }
1826 
1827 static const double
1828 tens[] = {
1829  1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1830  1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1831  1e20, 1e21, 1e22
1832 #ifdef VAX
1833  , 1e23, 1e24
1834 #endif
1835 };
1836 
1837 static const double
1838 #ifdef IEEE_Arith
1839 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1840 static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1841 #ifdef Avoid_Underflow
1842  9007199254740992.*9007199254740992.e-256
1843  /* = 2^106 * 1e-53 */
1844 #else
1845  1e-256
1846 #endif
1847 };
1848 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1849 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1850 #define Scale_Bit 0x10
1851 #define n_bigtens 5
1852 #else
1853 #ifdef IBM
1854 bigtens[] = { 1e16, 1e32, 1e64 };
1855 static const double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1856 #define n_bigtens 3
1857 #else
1858 bigtens[] = { 1e16, 1e32 };
1859 static const double tinytens[] = { 1e-16, 1e-32 };
1860 #define n_bigtens 2
1861 #endif
1862 #endif
1863 
1864 #ifndef IEEE_Arith
1865 #undef INFNAN_CHECK
1866 #endif
1867 
1868 #ifdef INFNAN_CHECK
1869 
1870 #ifndef NAN_WORD0
1871 #define NAN_WORD0 0x7ff80000
1872 #endif
1873 
1874 #ifndef NAN_WORD1
1875 #define NAN_WORD1 0
1876 #endif
1877 
1878 static int
1879 match(const char **sp, char *t)
1880 {
1881  int c, d;
1882  const char *s = *sp;
1883 
1884  while (d = *t++) {
1885  if ((c = *++s) >= 'A' && c <= 'Z')
1886  c += 'a' - 'A';
1887  if (c != d)
1888  return 0;
1889  }
1890  *sp = s + 1;
1891  return 1;
1892 }
1893 
1894 #ifndef No_Hex_NaN
1895 static void
1896 hexnan(double *rvp, const char **sp)
1897 {
1898  ULong c, x[2];
1899  const char *s;
1900  int havedig, udx0, xshift;
1901 
1902  x[0] = x[1] = 0;
1903  havedig = xshift = 0;
1904  udx0 = 1;
1905  s = *sp;
1906  while (c = *(const unsigned char*)++s) {
1907  if (c >= '0' && c <= '9')
1908  c -= '0';
1909  else if (c >= 'a' && c <= 'f')
1910  c += 10 - 'a';
1911  else if (c >= 'A' && c <= 'F')
1912  c += 10 - 'A';
1913  else if (c <= ' ') {
1914  if (udx0 && havedig) {
1915  udx0 = 0;
1916  xshift = 1;
1917  }
1918  continue;
1919  }
1920  else if (/*(*/ c == ')' && havedig) {
1921  *sp = s + 1;
1922  break;
1923  }
1924  else
1925  return; /* invalid form: don't change *sp */
1926  havedig = 1;
1927  if (xshift) {
1928  xshift = 0;
1929  x[0] = x[1];
1930  x[1] = 0;
1931  }
1932  if (udx0)
1933  x[0] = (x[0] << 4) | (x[1] >> 28);
1934  x[1] = (x[1] << 4) | c;
1935  }
1936  if ((x[0] &= 0xfffff) || x[1]) {
1937  word0(*rvp) = Exp_mask | x[0];
1938  word1(*rvp) = x[1];
1939  }
1940 }
1941 #endif /*No_Hex_NaN*/
1942 #endif /* INFNAN_CHECK */
1943 
1944 double
1945 ruby_strtod(const char *s00, char **se)
1946 {
1947 #ifdef Avoid_Underflow
1948  int scale;
1949 #endif
1950  int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1951  e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1952  const char *s, *s0, *s1;
1953  double aadj, adj;
1954  double_u aadj1, rv, rv0;
1955  Long L;
1956  ULong y, z;
1957  Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1958 #ifdef SET_INEXACT
1959  int inexact, oldinexact;
1960 #endif
1961 #ifdef Honor_FLT_ROUNDS
1962  int rounding;
1963 #endif
1964 #ifdef USE_LOCALE
1965  const char *s2;
1966 #endif
1967 
1968  errno = 0;
1969  sign = nz0 = nz = 0;
1970  dval(rv) = 0.;
1971  for (s = s00;;s++)
1972  switch (*s) {
1973  case '-':
1974  sign = 1;
1975  /* no break */
1976  case '+':
1977  if (*++s)
1978  goto break2;
1979  /* no break */
1980  case 0:
1981  goto ret0;
1982  case '\t':
1983  case '\n':
1984  case '\v':
1985  case '\f':
1986  case '\r':
1987  case ' ':
1988  continue;
1989  default:
1990  goto break2;
1991  }
1992 break2:
1993  if (*s == '0') {
1994  if (s[1] == 'x' || s[1] == 'X') {
1995  static const char hexdigit[] = "0123456789abcdef0123456789ABCDEF";
1996  s0 = ++s;
1997  adj = 0;
1998  aadj = 1.0;
1999  nd0 = -4;
2000 
2001  if (!*++s || !(s1 = strchr(hexdigit, *s))) goto ret0;
2002  if (*s == '0') {
2003  while (*++s == '0');
2004  s1 = strchr(hexdigit, *s);
2005  }
2006  if (s1 != NULL) {
2007  do {
2008  adj += aadj * ((s1 - hexdigit) & 15);
2009  nd0 += 4;
2010  aadj /= 16;
2011  } while (*++s && (s1 = strchr(hexdigit, *s)));
2012  }
2013 
2014  if (*s == '.') {
2015  dsign = 1;
2016  if (!*++s || !(s1 = strchr(hexdigit, *s))) goto ret0;
2017  if (nd0 < 0) {
2018  while (*s == '0') {
2019  s++;
2020  nd0 -= 4;
2021  }
2022  }
2023  for (; *s && (s1 = strchr(hexdigit, *s)); ++s) {
2024  adj += aadj * ((s1 - hexdigit) & 15);
2025  if ((aadj /= 16) == 0.0) {
2026  while (strchr(hexdigit, *++s));
2027  break;
2028  }
2029  }
2030  }
2031  else {
2032  dsign = 0;
2033  }
2034 
2035  if (*s == 'P' || *s == 'p') {
2036  dsign = 0x2C - *++s; /* +: 2B, -: 2D */
2037  if (abs(dsign) == 1) s++;
2038  else dsign = 1;
2039 
2040  nd = 0;
2041  c = *s;
2042  if (c < '0' || '9' < c) goto ret0;
2043  do {
2044  nd *= 10;
2045  nd += c;
2046  nd -= '0';
2047  c = *++s;
2048  /* Float("0x0."+("0"*267)+"1fp2095") */
2049  if (nd + dsign * nd0 > 2095) {
2050  while ('0' <= c && c <= '9') c = *++s;
2051  break;
2052  }
2053  } while ('0' <= c && c <= '9');
2054  nd0 += nd * dsign;
2055  }
2056  else {
2057  if (dsign) goto ret0;
2058  }
2059  dval(rv) = ldexp(adj, nd0);
2060  goto ret;
2061  }
2062  nz0 = 1;
2063  while (*++s == '0') ;
2064  if (!*s)
2065  goto ret;
2066  }
2067  s0 = s;
2068  y = z = 0;
2069  for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
2070  if (nd < 9)
2071  y = 10*y + c - '0';
2072  else if (nd < 16)
2073  z = 10*z + c - '0';
2074  nd0 = nd;
2075 #ifdef USE_LOCALE
2076  s1 = localeconv()->decimal_point;
2077  if (c == *s1) {
2078  c = '.';
2079  if (*++s1) {
2080  s2 = s;
2081  for (;;) {
2082  if (*++s2 != *s1) {
2083  c = 0;
2084  break;
2085  }
2086  if (!*++s1) {
2087  s = s2;
2088  break;
2089  }
2090  }
2091  }
2092  }
2093 #endif
2094  if (c == '.') {
2095  if (!ISDIGIT(s[1]))
2096  goto dig_done;
2097  c = *++s;
2098  if (!nd) {
2099  for (; c == '0'; c = *++s)
2100  nz++;
2101  if (c > '0' && c <= '9') {
2102  s0 = s;
2103  nf += nz;
2104  nz = 0;
2105  goto have_dig;
2106  }
2107  goto dig_done;
2108  }
2109  for (; c >= '0' && c <= '9'; c = *++s) {
2110 have_dig:
2111  nz++;
2112  if (nf > DBL_DIG * 4) continue;
2113  if (c -= '0') {
2114  nf += nz;
2115  for (i = 1; i < nz; i++)
2116  if (nd++ < 9)
2117  y *= 10;
2118  else if (nd <= DBL_DIG + 1)
2119  z *= 10;
2120  if (nd++ < 9)
2121  y = 10*y + c;
2122  else if (nd <= DBL_DIG + 1)
2123  z = 10*z + c;
2124  nz = 0;
2125  }
2126  }
2127  }
2128 dig_done:
2129  e = 0;
2130  if (c == 'e' || c == 'E') {
2131  if (!nd && !nz && !nz0) {
2132  goto ret0;
2133  }
2134  s00 = s;
2135  esign = 0;
2136  switch (c = *++s) {
2137  case '-':
2138  esign = 1;
2139  case '+':
2140  c = *++s;
2141  }
2142  if (c >= '0' && c <= '9') {
2143  while (c == '0')
2144  c = *++s;
2145  if (c > '0' && c <= '9') {
2146  L = c - '0';
2147  s1 = s;
2148  while ((c = *++s) >= '0' && c <= '9')
2149  L = 10*L + c - '0';
2150  if (s - s1 > 8 || L > 19999)
2151  /* Avoid confusion from exponents
2152  * so large that e might overflow.
2153  */
2154  e = 19999; /* safe for 16 bit ints */
2155  else
2156  e = (int)L;
2157  if (esign)
2158  e = -e;
2159  }
2160  else
2161  e = 0;
2162  }
2163  else
2164  s = s00;
2165  }
2166  if (!nd) {
2167  if (!nz && !nz0) {
2168 #ifdef INFNAN_CHECK
2169  /* Check for Nan and Infinity */
2170  switch (c) {
2171  case 'i':
2172  case 'I':
2173  if (match(&s,"nf")) {
2174  --s;
2175  if (!match(&s,"inity"))
2176  ++s;
2177  word0(rv) = 0x7ff00000;
2178  word1(rv) = 0;
2179  goto ret;
2180  }
2181  break;
2182  case 'n':
2183  case 'N':
2184  if (match(&s, "an")) {
2185  word0(rv) = NAN_WORD0;
2186  word1(rv) = NAN_WORD1;
2187 #ifndef No_Hex_NaN
2188  if (*s == '(') /*)*/
2189  hexnan(&rv, &s);
2190 #endif
2191  goto ret;
2192  }
2193  }
2194 #endif /* INFNAN_CHECK */
2195 ret0:
2196  s = s00;
2197  sign = 0;
2198  }
2199  goto ret;
2200  }
2201  e1 = e -= nf;
2202 
2203  /* Now we have nd0 digits, starting at s0, followed by a
2204  * decimal point, followed by nd-nd0 digits. The number we're
2205  * after is the integer represented by those digits times
2206  * 10**e */
2207 
2208  if (!nd0)
2209  nd0 = nd;
2210  k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
2211  dval(rv) = y;
2212  if (k > 9) {
2213 #ifdef SET_INEXACT
2214  if (k > DBL_DIG)
2215  oldinexact = get_inexact();
2216 #endif
2217  dval(rv) = tens[k - 9] * dval(rv) + z;
2218  }
2219  bd0 = bb = bd = bs = delta = 0;
2220  if (nd <= DBL_DIG
2221 #ifndef RND_PRODQUOT
2222 #ifndef Honor_FLT_ROUNDS
2223  && Flt_Rounds == 1
2224 #endif
2225 #endif
2226  ) {
2227  if (!e)
2228  goto ret;
2229  if (e > 0) {
2230  if (e <= Ten_pmax) {
2231 #ifdef VAX
2232  goto vax_ovfl_check;
2233 #else
2234 #ifdef Honor_FLT_ROUNDS
2235  /* round correctly FLT_ROUNDS = 2 or 3 */
2236  if (sign) {
2237  dval(rv) = -dval(rv);
2238  sign = 0;
2239  }
2240 #endif
2241  /* rv = */ rounded_product(dval(rv), tens[e]);
2242  goto ret;
2243 #endif
2244  }
2245  i = DBL_DIG - nd;
2246  if (e <= Ten_pmax + i) {
2247  /* A fancier test would sometimes let us do
2248  * this for larger i values.
2249  */
2250 #ifdef Honor_FLT_ROUNDS
2251  /* round correctly FLT_ROUNDS = 2 or 3 */
2252  if (sign) {
2253  dval(rv) = -dval(rv);
2254  sign = 0;
2255  }
2256 #endif
2257  e -= i;
2258  dval(rv) *= tens[i];
2259 #ifdef VAX
2260  /* VAX exponent range is so narrow we must
2261  * worry about overflow here...
2262  */
2263 vax_ovfl_check:
2264  word0(rv) -= P*Exp_msk1;
2265  /* rv = */ rounded_product(dval(rv), tens[e]);
2266  if ((word0(rv) & Exp_mask)
2267  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
2268  goto ovfl;
2269  word0(rv) += P*Exp_msk1;
2270 #else
2271  /* rv = */ rounded_product(dval(rv), tens[e]);
2272 #endif
2273  goto ret;
2274  }
2275  }
2276 #ifndef Inaccurate_Divide
2277  else if (e >= -Ten_pmax) {
2278 #ifdef Honor_FLT_ROUNDS
2279  /* round correctly FLT_ROUNDS = 2 or 3 */
2280  if (sign) {
2281  dval(rv) = -dval(rv);
2282  sign = 0;
2283  }
2284 #endif
2285  /* rv = */ rounded_quotient(dval(rv), tens[-e]);
2286  goto ret;
2287  }
2288 #endif
2289  }
2290  e1 += nd - k;
2291 
2292 #ifdef IEEE_Arith
2293 #ifdef SET_INEXACT
2294  inexact = 1;
2295  if (k <= DBL_DIG)
2296  oldinexact = get_inexact();
2297 #endif
2298 #ifdef Avoid_Underflow
2299  scale = 0;
2300 #endif
2301 #ifdef Honor_FLT_ROUNDS
2302  if ((rounding = Flt_Rounds) >= 2) {
2303  if (sign)
2304  rounding = rounding == 2 ? 0 : 2;
2305  else
2306  if (rounding != 2)
2307  rounding = 0;
2308  }
2309 #endif
2310 #endif /*IEEE_Arith*/
2311 
2312  /* Get starting approximation = rv * 10**e1 */
2313 
2314  if (e1 > 0) {
2315  if ((i = e1 & 15) != 0)
2316  dval(rv) *= tens[i];
2317  if (e1 &= ~15) {
2318  if (e1 > DBL_MAX_10_EXP) {
2319 ovfl:
2320 #ifndef NO_ERRNO
2321  errno = ERANGE;
2322 #endif
2323  /* Can't trust HUGE_VAL */
2324 #ifdef IEEE_Arith
2325 #ifdef Honor_FLT_ROUNDS
2326  switch (rounding) {
2327  case 0: /* toward 0 */
2328  case 3: /* toward -infinity */
2329  word0(rv) = Big0;
2330  word1(rv) = Big1;
2331  break;
2332  default:
2333  word0(rv) = Exp_mask;
2334  word1(rv) = 0;
2335  }
2336 #else /*Honor_FLT_ROUNDS*/
2337  word0(rv) = Exp_mask;
2338  word1(rv) = 0;
2339 #endif /*Honor_FLT_ROUNDS*/
2340 #ifdef SET_INEXACT
2341  /* set overflow bit */
2342  dval(rv0) = 1e300;
2343  dval(rv0) *= dval(rv0);
2344 #endif
2345 #else /*IEEE_Arith*/
2346  word0(rv) = Big0;
2347  word1(rv) = Big1;
2348 #endif /*IEEE_Arith*/
2349  if (bd0)
2350  goto retfree;
2351  goto ret;
2352  }
2353  e1 >>= 4;
2354  for (j = 0; e1 > 1; j++, e1 >>= 1)
2355  if (e1 & 1)
2356  dval(rv) *= bigtens[j];
2357  /* The last multiplication could overflow. */
2358  word0(rv) -= P*Exp_msk1;
2359  dval(rv) *= bigtens[j];
2360  if ((z = word0(rv) & Exp_mask)
2361  > Exp_msk1*(DBL_MAX_EXP+Bias-P))
2362  goto ovfl;
2363  if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
2364  /* set to largest number */
2365  /* (Can't trust DBL_MAX) */
2366  word0(rv) = Big0;
2367  word1(rv) = Big1;
2368  }
2369  else
2370  word0(rv) += P*Exp_msk1;
2371  }
2372  }
2373  else if (e1 < 0) {
2374  e1 = -e1;
2375  if ((i = e1 & 15) != 0)
2376  dval(rv) /= tens[i];
2377  if (e1 >>= 4) {
2378  if (e1 >= 1 << n_bigtens)
2379  goto undfl;
2380 #ifdef Avoid_Underflow
2381  if (e1 & Scale_Bit)
2382  scale = 2*P;
2383  for (j = 0; e1 > 0; j++, e1 >>= 1)
2384  if (e1 & 1)
2385  dval(rv) *= tinytens[j];
2386  if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
2387  >> Exp_shift)) > 0) {
2388  /* scaled rv is denormal; zap j low bits */
2389  if (j >= 32) {
2390  word1(rv) = 0;
2391  if (j >= 53)
2392  word0(rv) = (P+2)*Exp_msk1;
2393  else
2394  word0(rv) &= 0xffffffff << (j-32);
2395  }
2396  else
2397  word1(rv) &= 0xffffffff << j;
2398  }
2399 #else
2400  for (j = 0; e1 > 1; j++, e1 >>= 1)
2401  if (e1 & 1)
2402  dval(rv) *= tinytens[j];
2403  /* The last multiplication could underflow. */
2404  dval(rv0) = dval(rv);
2405  dval(rv) *= tinytens[j];
2406  if (!dval(rv)) {
2407  dval(rv) = 2.*dval(rv0);
2408  dval(rv) *= tinytens[j];
2409 #endif
2410  if (!dval(rv)) {
2411 undfl:
2412  dval(rv) = 0.;
2413 #ifndef NO_ERRNO
2414  errno = ERANGE;
2415 #endif
2416  if (bd0)
2417  goto retfree;
2418  goto ret;
2419  }
2420 #ifndef Avoid_Underflow
2421  word0(rv) = Tiny0;
2422  word1(rv) = Tiny1;
2423  /* The refinement below will clean
2424  * this approximation up.
2425  */
2426  }
2427 #endif
2428  }
2429  }
2430 
2431  /* Now the hard part -- adjusting rv to the correct value.*/
2432 
2433  /* Put digits into bd: true value = bd * 10^e */
2434 
2435  bd0 = s2b(s0, nd0, nd, y);
2436 
2437  for (;;) {
2438  bd = Balloc(bd0->k);
2439  Bcopy(bd, bd0);
2440  bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
2441  bs = i2b(1);
2442 
2443  if (e >= 0) {
2444  bb2 = bb5 = 0;
2445  bd2 = bd5 = e;
2446  }
2447  else {
2448  bb2 = bb5 = -e;
2449  bd2 = bd5 = 0;
2450  }
2451  if (bbe >= 0)
2452  bb2 += bbe;
2453  else
2454  bd2 -= bbe;
2455  bs2 = bb2;
2456 #ifdef Honor_FLT_ROUNDS
2457  if (rounding != 1)
2458  bs2++;
2459 #endif
2460 #ifdef Avoid_Underflow
2461  j = bbe - scale;
2462  i = j + bbbits - 1; /* logb(rv) */
2463  if (i < Emin) /* denormal */
2464  j += P - Emin;
2465  else
2466  j = P + 1 - bbbits;
2467 #else /*Avoid_Underflow*/
2468 #ifdef Sudden_Underflow
2469 #ifdef IBM
2470  j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
2471 #else
2472  j = P + 1 - bbbits;
2473 #endif
2474 #else /*Sudden_Underflow*/
2475  j = bbe;
2476  i = j + bbbits - 1; /* logb(rv) */
2477  if (i < Emin) /* denormal */
2478  j += P - Emin;
2479  else
2480  j = P + 1 - bbbits;
2481 #endif /*Sudden_Underflow*/
2482 #endif /*Avoid_Underflow*/
2483  bb2 += j;
2484  bd2 += j;
2485 #ifdef Avoid_Underflow
2486  bd2 += scale;
2487 #endif
2488  i = bb2 < bd2 ? bb2 : bd2;
2489  if (i > bs2)
2490  i = bs2;
2491  if (i > 0) {
2492  bb2 -= i;
2493  bd2 -= i;
2494  bs2 -= i;
2495  }
2496  if (bb5 > 0) {
2497  bs = pow5mult(bs, bb5);
2498  bb1 = mult(bs, bb);
2499  Bfree(bb);
2500  bb = bb1;
2501  }
2502  if (bb2 > 0)
2503  bb = lshift(bb, bb2);
2504  if (bd5 > 0)
2505  bd = pow5mult(bd, bd5);
2506  if (bd2 > 0)
2507  bd = lshift(bd, bd2);
2508  if (bs2 > 0)
2509  bs = lshift(bs, bs2);
2510  delta = diff(bb, bd);
2511  dsign = delta->sign;
2512  delta->sign = 0;
2513  i = cmp(delta, bs);
2514 #ifdef Honor_FLT_ROUNDS
2515  if (rounding != 1) {
2516  if (i < 0) {
2517  /* Error is less than an ulp */
2518  if (!delta->x[0] && delta->wds <= 1) {
2519  /* exact */
2520 #ifdef SET_INEXACT
2521  inexact = 0;
2522 #endif
2523  break;
2524  }
2525  if (rounding) {
2526  if (dsign) {
2527  adj = 1.;
2528  goto apply_adj;
2529  }
2530  }
2531  else if (!dsign) {
2532  adj = -1.;
2533  if (!word1(rv)
2534  && !(word0(rv) & Frac_mask)) {
2535  y = word0(rv) & Exp_mask;
2536 #ifdef Avoid_Underflow
2537  if (!scale || y > 2*P*Exp_msk1)
2538 #else
2539  if (y)
2540 #endif
2541  {
2542  delta = lshift(delta,Log2P);
2543  if (cmp(delta, bs) <= 0)
2544  adj = -0.5;
2545  }
2546  }
2547 apply_adj:
2548 #ifdef Avoid_Underflow
2549  if (scale && (y = word0(rv) & Exp_mask)
2550  <= 2*P*Exp_msk1)
2551  word0(adj) += (2*P+1)*Exp_msk1 - y;
2552 #else
2553 #ifdef Sudden_Underflow
2554  if ((word0(rv) & Exp_mask) <=
2555  P*Exp_msk1) {
2556  word0(rv) += P*Exp_msk1;
2557  dval(rv) += adj*ulp(dval(rv));
2558  word0(rv) -= P*Exp_msk1;
2559  }
2560  else
2561 #endif /*Sudden_Underflow*/
2562 #endif /*Avoid_Underflow*/
2563  dval(rv) += adj*ulp(dval(rv));
2564  }
2565  break;
2566  }
2567  adj = ratio(delta, bs);
2568  if (adj < 1.)
2569  adj = 1.;
2570  if (adj <= 0x7ffffffe) {
2571  /* adj = rounding ? ceil(adj) : floor(adj); */
2572  y = adj;
2573  if (y != adj) {
2574  if (!((rounding>>1) ^ dsign))
2575  y++;
2576  adj = y;
2577  }
2578  }
2579 #ifdef Avoid_Underflow
2580  if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2581  word0(adj) += (2*P+1)*Exp_msk1 - y;
2582 #else
2583 #ifdef Sudden_Underflow
2584  if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2585  word0(rv) += P*Exp_msk1;
2586  adj *= ulp(dval(rv));
2587  if (dsign)
2588  dval(rv) += adj;
2589  else
2590  dval(rv) -= adj;
2591  word0(rv) -= P*Exp_msk1;
2592  goto cont;
2593  }
2594 #endif /*Sudden_Underflow*/
2595 #endif /*Avoid_Underflow*/
2596  adj *= ulp(dval(rv));
2597  if (dsign)
2598  dval(rv) += adj;
2599  else
2600  dval(rv) -= adj;
2601  goto cont;
2602  }
2603 #endif /*Honor_FLT_ROUNDS*/
2604 
2605  if (i < 0) {
2606  /* Error is less than half an ulp -- check for
2607  * special case of mantissa a power of two.
2608  */
2609  if (dsign || word1(rv) || word0(rv) & Bndry_mask
2610 #ifdef IEEE_Arith
2611 #ifdef Avoid_Underflow
2612  || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2613 #else
2614  || (word0(rv) & Exp_mask) <= Exp_msk1
2615 #endif
2616 #endif
2617  ) {
2618 #ifdef SET_INEXACT
2619  if (!delta->x[0] && delta->wds <= 1)
2620  inexact = 0;
2621 #endif
2622  break;
2623  }
2624  if (!delta->x[0] && delta->wds <= 1) {
2625  /* exact result */
2626 #ifdef SET_INEXACT
2627  inexact = 0;
2628 #endif
2629  break;
2630  }
2631  delta = lshift(delta,Log2P);
2632  if (cmp(delta, bs) > 0)
2633  goto drop_down;
2634  break;
2635  }
2636  if (i == 0) {
2637  /* exactly half-way between */
2638  if (dsign) {
2639  if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2640  && word1(rv) == (
2641 #ifdef Avoid_Underflow
2642  (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2643  ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2644 #endif
2645  0xffffffff)) {
2646  /*boundary case -- increment exponent*/
2647  word0(rv) = (word0(rv) & Exp_mask)
2648  + Exp_msk1
2649 #ifdef IBM
2650  | Exp_msk1 >> 4
2651 #endif
2652  ;
2653  word1(rv) = 0;
2654 #ifdef Avoid_Underflow
2655  dsign = 0;
2656 #endif
2657  break;
2658  }
2659  }
2660  else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2661 drop_down:
2662  /* boundary case -- decrement exponent */
2663 #ifdef Sudden_Underflow /*{{*/
2664  L = word0(rv) & Exp_mask;
2665 #ifdef IBM
2666  if (L < Exp_msk1)
2667 #else
2668 #ifdef Avoid_Underflow
2669  if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
2670 #else
2671  if (L <= Exp_msk1)
2672 #endif /*Avoid_Underflow*/
2673 #endif /*IBM*/
2674  goto undfl;
2675  L -= Exp_msk1;
2676 #else /*Sudden_Underflow}{*/
2677 #ifdef Avoid_Underflow
2678  if (scale) {
2679  L = word0(rv) & Exp_mask;
2680  if (L <= (2*P+1)*Exp_msk1) {
2681  if (L > (P+2)*Exp_msk1)
2682  /* round even ==> */
2683  /* accept rv */
2684  break;
2685  /* rv = smallest denormal */
2686  goto undfl;
2687  }
2688  }
2689 #endif /*Avoid_Underflow*/
2690  L = (word0(rv) & Exp_mask) - Exp_msk1;
2691 #endif /*Sudden_Underflow}}*/
2692  word0(rv) = L | Bndry_mask1;
2693  word1(rv) = 0xffffffff;
2694 #ifdef IBM
2695  goto cont;
2696 #else
2697  break;
2698 #endif
2699  }
2700 #ifndef ROUND_BIASED
2701  if (!(word1(rv) & LSB))
2702  break;
2703 #endif
2704  if (dsign)
2705  dval(rv) += ulp(dval(rv));
2706 #ifndef ROUND_BIASED
2707  else {
2708  dval(rv) -= ulp(dval(rv));
2709 #ifndef Sudden_Underflow
2710  if (!dval(rv))
2711  goto undfl;
2712 #endif
2713  }
2714 #ifdef Avoid_Underflow
2715  dsign = 1 - dsign;
2716 #endif
2717 #endif
2718  break;
2719  }
2720  if ((aadj = ratio(delta, bs)) <= 2.) {
2721  if (dsign)
2722  aadj = dval(aadj1) = 1.;
2723  else if (word1(rv) || word0(rv) & Bndry_mask) {
2724 #ifndef Sudden_Underflow
2725  if (word1(rv) == Tiny1 && !word0(rv))
2726  goto undfl;
2727 #endif
2728  aadj = 1.;
2729  dval(aadj1) = -1.;
2730  }
2731  else {
2732  /* special case -- power of FLT_RADIX to be */
2733  /* rounded down... */
2734 
2735  if (aadj < 2./FLT_RADIX)
2736  aadj = 1./FLT_RADIX;
2737  else
2738  aadj *= 0.5;
2739  dval(aadj1) = -aadj;
2740  }
2741  }
2742  else {
2743  aadj *= 0.5;
2744  dval(aadj1) = dsign ? aadj : -aadj;
2745 #ifdef Check_FLT_ROUNDS
2746  switch (Rounding) {
2747  case 2: /* towards +infinity */
2748  dval(aadj1) -= 0.5;
2749  break;
2750  case 0: /* towards 0 */
2751  case 3: /* towards -infinity */
2752  dval(aadj1) += 0.5;
2753  }
2754 #else
2755  if (Flt_Rounds == 0)
2756  dval(aadj1) += 0.5;
2757 #endif /*Check_FLT_ROUNDS*/
2758  }
2759  y = word0(rv) & Exp_mask;
2760 
2761  /* Check for overflow */
2762 
2763  if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2764  dval(rv0) = dval(rv);
2765  word0(rv) -= P*Exp_msk1;
2766  adj = dval(aadj1) * ulp(dval(rv));
2767  dval(rv) += adj;
2768  if ((word0(rv) & Exp_mask) >=
2769  Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2770  if (word0(rv0) == Big0 && word1(rv0) == Big1)
2771  goto ovfl;
2772  word0(rv) = Big0;
2773  word1(rv) = Big1;
2774  goto cont;
2775  }
2776  else
2777  word0(rv) += P*Exp_msk1;
2778  }
2779  else {
2780 #ifdef Avoid_Underflow
2781  if (scale && y <= 2*P*Exp_msk1) {
2782  if (aadj <= 0x7fffffff) {
2783  if ((z = (int)aadj) <= 0)
2784  z = 1;
2785  aadj = z;
2786  dval(aadj1) = dsign ? aadj : -aadj;
2787  }
2788  word0(aadj1) += (2*P+1)*Exp_msk1 - y;
2789  }
2790  adj = dval(aadj1) * ulp(dval(rv));
2791  dval(rv) += adj;
2792 #else
2793 #ifdef Sudden_Underflow
2794  if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2795  dval(rv0) = dval(rv);
2796  word0(rv) += P*Exp_msk1;
2797  adj = dval(aadj1) * ulp(dval(rv));
2798  dval(rv) += adj;
2799 #ifdef IBM
2800  if ((word0(rv) & Exp_mask) < P*Exp_msk1)
2801 #else
2802  if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2803 #endif
2804  {
2805  if (word0(rv0) == Tiny0 && word1(rv0) == Tiny1)
2806  goto undfl;
2807  word0(rv) = Tiny0;
2808  word1(rv) = Tiny1;
2809  goto cont;
2810  }
2811  else
2812  word0(rv) -= P*Exp_msk1;
2813  }
2814  else {
2815  adj = dval(aadj1) * ulp(dval(rv));
2816  dval(rv) += adj;
2817  }
2818 #else /*Sudden_Underflow*/
2819  /* Compute adj so that the IEEE rounding rules will
2820  * correctly round rv + adj in some half-way cases.
2821  * If rv * ulp(rv) is denormalized (i.e.,
2822  * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2823  * trouble from bits lost to denormalization;
2824  * example: 1.2e-307 .
2825  */
2826  if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
2827  dval(aadj1) = (double)(int)(aadj + 0.5);
2828  if (!dsign)
2829  dval(aadj1) = -dval(aadj1);
2830  }
2831  adj = dval(aadj1) * ulp(dval(rv));
2832  dval(rv) += adj;
2833 #endif /*Sudden_Underflow*/
2834 #endif /*Avoid_Underflow*/
2835  }
2836  z = word0(rv) & Exp_mask;
2837 #ifndef SET_INEXACT
2838 #ifdef Avoid_Underflow
2839  if (!scale)
2840 #endif
2841  if (y == z) {
2842  /* Can we stop now? */
2843  L = (Long)aadj;
2844  aadj -= L;
2845  /* The tolerances below are conservative. */
2846  if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2847  if (aadj < .4999999 || aadj > .5000001)
2848  break;
2849  }
2850  else if (aadj < .4999999/FLT_RADIX)
2851  break;
2852  }
2853 #endif
2854 cont:
2855  Bfree(bb);
2856  Bfree(bd);
2857  Bfree(bs);
2858  Bfree(delta);
2859  }
2860 #ifdef SET_INEXACT
2861  if (inexact) {
2862  if (!oldinexact) {
2863  word0(rv0) = Exp_1 + (70 << Exp_shift);
2864  word1(rv0) = 0;
2865  dval(rv0) += 1.;
2866  }
2867  }
2868  else if (!oldinexact)
2869  clear_inexact();
2870 #endif
2871 #ifdef Avoid_Underflow
2872  if (scale) {
2873  word0(rv0) = Exp_1 - 2*P*Exp_msk1;
2874  word1(rv0) = 0;
2875  dval(rv) *= dval(rv0);
2876 #ifndef NO_ERRNO
2877  /* try to avoid the bug of testing an 8087 register value */
2878  if (word0(rv) == 0 && word1(rv) == 0)
2879  errno = ERANGE;
2880 #endif
2881  }
2882 #endif /* Avoid_Underflow */
2883 #ifdef SET_INEXACT
2884  if (inexact && !(word0(rv) & Exp_mask)) {
2885  /* set underflow bit */
2886  dval(rv0) = 1e-300;
2887  dval(rv0) *= dval(rv0);
2888  }
2889 #endif
2890 retfree:
2891  Bfree(bb);
2892  Bfree(bd);
2893  Bfree(bs);
2894  Bfree(bd0);
2895  Bfree(delta);
2896 ret:
2897  if (se)
2898  *se = (char *)s;
2899  return sign ? -dval(rv) : dval(rv);
2900 }
2901 
2902 static int
2904 {
2905  int n;
2906  ULong *bx, *bxe, q, *sx, *sxe;
2907 #ifdef ULLong
2908  ULLong borrow, carry, y, ys;
2909 #else
2910  ULong borrow, carry, y, ys;
2911 #ifdef Pack_32
2912  ULong si, z, zs;
2913 #endif
2914 #endif
2915 
2916  n = S->wds;
2917 #ifdef DEBUG
2918  /*debug*/ if (b->wds > n)
2919  /*debug*/ Bug("oversize b in quorem");
2920 #endif
2921  if (b->wds < n)
2922  return 0;
2923  sx = S->x;
2924  sxe = sx + --n;
2925  bx = b->x;
2926  bxe = bx + n;
2927  q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
2928 #ifdef DEBUG
2929  /*debug*/ if (q > 9)
2930  /*debug*/ Bug("oversized quotient in quorem");
2931 #endif
2932  if (q) {
2933  borrow = 0;
2934  carry = 0;
2935  do {
2936 #ifdef ULLong
2937  ys = *sx++ * (ULLong)q + carry;
2938  carry = ys >> 32;
2939  y = *bx - (ys & FFFFFFFF) - borrow;
2940  borrow = y >> 32 & (ULong)1;
2941  *bx++ = (ULong)(y & FFFFFFFF);
2942 #else
2943 #ifdef Pack_32
2944  si = *sx++;
2945  ys = (si & 0xffff) * q + carry;
2946  zs = (si >> 16) * q + (ys >> 16);
2947  carry = zs >> 16;
2948  y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2949  borrow = (y & 0x10000) >> 16;
2950  z = (*bx >> 16) - (zs & 0xffff) - borrow;
2951  borrow = (z & 0x10000) >> 16;
2952  Storeinc(bx, z, y);
2953 #else
2954  ys = *sx++ * q + carry;
2955  carry = ys >> 16;
2956  y = *bx - (ys & 0xffff) - borrow;
2957  borrow = (y & 0x10000) >> 16;
2958  *bx++ = y & 0xffff;
2959 #endif
2960 #endif
2961  } while (sx <= sxe);
2962  if (!*bxe) {
2963  bx = b->x;
2964  while (--bxe > bx && !*bxe)
2965  --n;
2966  b->wds = n;
2967  }
2968  }
2969  if (cmp(b, S) >= 0) {
2970  q++;
2971  borrow = 0;
2972  carry = 0;
2973  bx = b->x;
2974  sx = S->x;
2975  do {
2976 #ifdef ULLong
2977  ys = *sx++ + carry;
2978  carry = ys >> 32;
2979  y = *bx - (ys & FFFFFFFF) - borrow;
2980  borrow = y >> 32 & (ULong)1;
2981  *bx++ = (ULong)(y & FFFFFFFF);
2982 #else
2983 #ifdef Pack_32
2984  si = *sx++;
2985  ys = (si & 0xffff) + carry;
2986  zs = (si >> 16) + (ys >> 16);
2987  carry = zs >> 16;
2988  y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2989  borrow = (y & 0x10000) >> 16;
2990  z = (*bx >> 16) - (zs & 0xffff) - borrow;
2991  borrow = (z & 0x10000) >> 16;
2992  Storeinc(bx, z, y);
2993 #else
2994  ys = *sx++ + carry;
2995  carry = ys >> 16;
2996  y = *bx - (ys & 0xffff) - borrow;
2997  borrow = (y & 0x10000) >> 16;
2998  *bx++ = y & 0xffff;
2999 #endif
3000 #endif
3001  } while (sx <= sxe);
3002  bx = b->x;
3003  bxe = bx + n;
3004  if (!*bxe) {
3005  while (--bxe > bx && !*bxe)
3006  --n;
3007  b->wds = n;
3008  }
3009  }
3010  return q;
3011 }
3012 
3013 #ifndef MULTIPLE_THREADS
3014 static char *dtoa_result;
3015 #endif
3016 
3017 #ifndef MULTIPLE_THREADS
3018 static char *
3019 rv_alloc(int i)
3020 {
3021  return dtoa_result = xmalloc(i);
3022 }
3023 #else
3024 #define rv_alloc(i) xmalloc(i)
3025 #endif
3026 
3027 static char *
3028 nrv_alloc(const char *s, char **rve, size_t n)
3029 {
3030  char *rv, *t;
3031 
3032  t = rv = rv_alloc(n);
3033  while ((*t = *s++) != 0) t++;
3034  if (rve)
3035  *rve = t;
3036  return rv;
3037 }
3038 
3039 #define rv_strdup(s, rve) nrv_alloc((s), (rve), strlen(s)+1)
3040 
3041 #ifndef MULTIPLE_THREADS
3042 /* freedtoa(s) must be used to free values s returned by dtoa
3043  * when MULTIPLE_THREADS is #defined. It should be used in all cases,
3044  * but for consistency with earlier versions of dtoa, it is optional
3045  * when MULTIPLE_THREADS is not defined.
3046  */
3047 
3048 static void
3049 freedtoa(char *s)
3050 {
3051  xfree(s);
3052 }
3053 #endif
3054 
3055 static const char INFSTR[] = "Infinity";
3056 static const char NANSTR[] = "NaN";
3057 static const char ZEROSTR[] = "0";
3058 
3059 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
3060  *
3061  * Inspired by "How to Print Floating-Point Numbers Accurately" by
3062  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
3063  *
3064  * Modifications:
3065  * 1. Rather than iterating, we use a simple numeric overestimate
3066  * to determine k = floor(log10(d)). We scale relevant
3067  * quantities using O(log2(k)) rather than O(k) multiplications.
3068  * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
3069  * try to generate digits strictly left to right. Instead, we
3070  * compute with fewer bits and propagate the carry if necessary
3071  * when rounding the final digit up. This is often faster.
3072  * 3. Under the assumption that input will be rounded nearest,
3073  * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
3074  * That is, we allow equality in stopping tests when the
3075  * round-nearest rule will give the same floating-point value
3076  * as would satisfaction of the stopping test with strict
3077  * inequality.
3078  * 4. We remove common factors of powers of 2 from relevant
3079  * quantities.
3080  * 5. When converting floating-point integers less than 1e16,
3081  * we use floating-point arithmetic rather than resorting
3082  * to multiple-precision integers.
3083  * 6. When asked to produce fewer than 15 digits, we first try
3084  * to get by with floating-point arithmetic; we resort to
3085  * multiple-precision integer arithmetic only if we cannot
3086  * guarantee that the floating-point calculation has given
3087  * the correctly rounded result. For k requested digits and
3088  * "uniformly" distributed input, the probability is
3089  * something like 10^(k-15) that we must resort to the Long
3090  * calculation.
3091  */
3092 
3093 char *
3094 ruby_dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
3095 {
3096  /* Arguments ndigits, decpt, sign are similar to those
3097  of ecvt and fcvt; trailing zeros are suppressed from
3098  the returned string. If not null, *rve is set to point
3099  to the end of the return value. If d is +-Infinity or NaN,
3100  then *decpt is set to 9999.
3101 
3102  mode:
3103  0 ==> shortest string that yields d when read in
3104  and rounded to nearest.
3105  1 ==> like 0, but with Steele & White stopping rule;
3106  e.g. with IEEE P754 arithmetic , mode 0 gives
3107  1e23 whereas mode 1 gives 9.999999999999999e22.
3108  2 ==> max(1,ndigits) significant digits. This gives a
3109  return value similar to that of ecvt, except
3110  that trailing zeros are suppressed.
3111  3 ==> through ndigits past the decimal point. This
3112  gives a return value similar to that from fcvt,
3113  except that trailing zeros are suppressed, and
3114  ndigits can be negative.
3115  4,5 ==> similar to 2 and 3, respectively, but (in
3116  round-nearest mode) with the tests of mode 0 to
3117  possibly return a shorter string that rounds to d.
3118  With IEEE arithmetic and compilation with
3119  -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
3120  as modes 2 and 3 when FLT_ROUNDS != 1.
3121  6-9 ==> Debugging modes similar to mode - 4: don't try
3122  fast floating-point estimate (if applicable).
3123 
3124  Values of mode other than 0-9 are treated as mode 0.
3125 
3126  Sufficient space is allocated to the return value
3127  to hold the suppressed trailing zeros.
3128  */
3129 
3130  int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
3131  j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
3132  spec_case, try_quick;
3133  Long L;
3134 #ifndef Sudden_Underflow
3135  int denorm;
3136  ULong x;
3137 #endif
3138  Bigint *b, *b1, *delta, *mlo = 0, *mhi = 0, *S;
3139  double ds;
3140  double_u d, d2, eps;
3141  char *s, *s0;
3142 #ifdef Honor_FLT_ROUNDS
3143  int rounding;
3144 #endif
3145 #ifdef SET_INEXACT
3146  int inexact, oldinexact;
3147 #endif
3148 
3149  dval(d) = d_;
3150 
3151 #ifndef MULTIPLE_THREADS
3152  if (dtoa_result) {
3153  freedtoa(dtoa_result);
3154  dtoa_result = 0;
3155  }
3156 #endif
3157 
3158  if (word0(d) & Sign_bit) {
3159  /* set sign for everything, including 0's and NaNs */
3160  *sign = 1;
3161  word0(d) &= ~Sign_bit; /* clear sign bit */
3162  }
3163  else
3164  *sign = 0;
3165 
3166 #if defined(IEEE_Arith) + defined(VAX)
3167 #ifdef IEEE_Arith
3168  if ((word0(d) & Exp_mask) == Exp_mask)
3169 #else
3170  if (word0(d) == 0x8000)
3171 #endif
3172  {
3173  /* Infinity or NaN */
3174  *decpt = 9999;
3175 #ifdef IEEE_Arith
3176  if (!word1(d) && !(word0(d) & 0xfffff))
3177  return rv_strdup(INFSTR, rve);
3178 #endif
3179  return rv_strdup(NANSTR, rve);
3180  }
3181 #endif
3182 #ifdef IBM
3183  dval(d) += 0; /* normalize */
3184 #endif
3185  if (!dval(d)) {
3186  *decpt = 1;
3187  return rv_strdup(ZEROSTR, rve);
3188  }
3189 
3190 #ifdef SET_INEXACT
3191  try_quick = oldinexact = get_inexact();
3192  inexact = 1;
3193 #endif
3194 #ifdef Honor_FLT_ROUNDS
3195  if ((rounding = Flt_Rounds) >= 2) {
3196  if (*sign)
3197  rounding = rounding == 2 ? 0 : 2;
3198  else
3199  if (rounding != 2)
3200  rounding = 0;
3201  }
3202 #endif
3203 
3204  b = d2b(dval(d), &be, &bbits);
3205 #ifdef Sudden_Underflow
3206  i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
3207 #else
3208  if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) {
3209 #endif
3210  dval(d2) = dval(d);
3211  word0(d2) &= Frac_mask1;
3212  word0(d2) |= Exp_11;
3213 #ifdef IBM
3214  if (j = 11 - hi0bits(word0(d2) & Frac_mask))
3215  dval(d2) /= 1 << j;
3216 #endif
3217 
3218  /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
3219  * log10(x) = log(x) / log(10)
3220  * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
3221  * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
3222  *
3223  * This suggests computing an approximation k to log10(d) by
3224  *
3225  * k = (i - Bias)*0.301029995663981
3226  * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
3227  *
3228  * We want k to be too large rather than too small.
3229  * The error in the first-order Taylor series approximation
3230  * is in our favor, so we just round up the constant enough
3231  * to compensate for any error in the multiplication of
3232  * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
3233  * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
3234  * adding 1e-13 to the constant term more than suffices.
3235  * Hence we adjust the constant term to 0.1760912590558.
3236  * (We could get a more accurate k by invoking log10,
3237  * but this is probably not worthwhile.)
3238  */
3239 
3240  i -= Bias;
3241 #ifdef IBM
3242  i <<= 2;
3243  i += j;
3244 #endif
3245 #ifndef Sudden_Underflow
3246  denorm = 0;
3247  }
3248  else {
3249  /* d is denormalized */
3250 
3251  i = bbits + be + (Bias + (P-1) - 1);
3252  x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
3253  : word1(d) << (32 - i);
3254  dval(d2) = x;
3255  word0(d2) -= 31*Exp_msk1; /* adjust exponent */
3256  i -= (Bias + (P-1) - 1) + 1;
3257  denorm = 1;
3258  }
3259 #endif
3260  ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
3261  k = (int)ds;
3262  if (ds < 0. && ds != k)
3263  k--; /* want k = floor(ds) */
3264  k_check = 1;
3265  if (k >= 0 && k <= Ten_pmax) {
3266  if (dval(d) < tens[k])
3267  k--;
3268  k_check = 0;
3269  }
3270  j = bbits - i - 1;
3271  if (j >= 0) {
3272  b2 = 0;
3273  s2 = j;
3274  }
3275  else {
3276  b2 = -j;
3277  s2 = 0;
3278  }
3279  if (k >= 0) {
3280  b5 = 0;
3281  s5 = k;
3282  s2 += k;
3283  }
3284  else {
3285  b2 -= k;
3286  b5 = -k;
3287  s5 = 0;
3288  }
3289  if (mode < 0 || mode > 9)
3290  mode = 0;
3291 
3292 #ifndef SET_INEXACT
3293 #ifdef Check_FLT_ROUNDS
3294  try_quick = Rounding == 1;
3295 #else
3296  try_quick = 1;
3297 #endif
3298 #endif /*SET_INEXACT*/
3299 
3300  if (mode > 5) {
3301  mode -= 4;
3302  try_quick = 0;
3303  }
3304  leftright = 1;
3305  ilim = ilim1 = -1;
3306  switch (mode) {
3307  case 0:
3308  case 1:
3309  i = 18;
3310  ndigits = 0;
3311  break;
3312  case 2:
3313  leftright = 0;
3314  /* no break */
3315  case 4:
3316  if (ndigits <= 0)
3317  ndigits = 1;
3318  ilim = ilim1 = i = ndigits;
3319  break;
3320  case 3:
3321  leftright = 0;
3322  /* no break */
3323  case 5:
3324  i = ndigits + k + 1;
3325  ilim = i;
3326  ilim1 = i - 1;
3327  if (i <= 0)
3328  i = 1;
3329  }
3330  s = s0 = rv_alloc(i+1);
3331 
3332 #ifdef Honor_FLT_ROUNDS
3333  if (mode > 1 && rounding != 1)
3334  leftright = 0;
3335 #endif
3336 
3337  if (ilim >= 0 && ilim <= Quick_max && try_quick) {
3338 
3339  /* Try to get by with floating-point arithmetic. */
3340 
3341  i = 0;
3342  dval(d2) = dval(d);
3343  k0 = k;
3344  ilim0 = ilim;
3345  ieps = 2; /* conservative */
3346  if (k > 0) {
3347  ds = tens[k&0xf];
3348  j = k >> 4;
3349  if (j & Bletch) {
3350  /* prevent overflows */
3351  j &= Bletch - 1;
3352  dval(d) /= bigtens[n_bigtens-1];
3353  ieps++;
3354  }
3355  for (; j; j >>= 1, i++)
3356  if (j & 1) {
3357  ieps++;
3358  ds *= bigtens[i];
3359  }
3360  dval(d) /= ds;
3361  }
3362  else if ((j1 = -k) != 0) {
3363  dval(d) *= tens[j1 & 0xf];
3364  for (j = j1 >> 4; j; j >>= 1, i++)
3365  if (j & 1) {
3366  ieps++;
3367  dval(d) *= bigtens[i];
3368  }
3369  }
3370  if (k_check && dval(d) < 1. && ilim > 0) {
3371  if (ilim1 <= 0)
3372  goto fast_failed;
3373  ilim = ilim1;
3374  k--;
3375  dval(d) *= 10.;
3376  ieps++;
3377  }
3378  dval(eps) = ieps*dval(d) + 7.;
3379  word0(eps) -= (P-1)*Exp_msk1;
3380  if (ilim == 0) {
3381  S = mhi = 0;
3382  dval(d) -= 5.;
3383  if (dval(d) > dval(eps))
3384  goto one_digit;
3385  if (dval(d) < -dval(eps))
3386  goto no_digits;
3387  goto fast_failed;
3388  }
3389 #ifndef No_leftright
3390  if (leftright) {
3391  /* Use Steele & White method of only
3392  * generating digits needed.
3393  */
3394  dval(eps) = 0.5/tens[ilim-1] - dval(eps);
3395  for (i = 0;;) {
3396  L = (int)dval(d);
3397  dval(d) -= L;
3398  *s++ = '0' + (int)L;
3399  if (dval(d) < dval(eps))
3400  goto ret1;
3401  if (1. - dval(d) < dval(eps))
3402  goto bump_up;
3403  if (++i >= ilim)
3404  break;
3405  dval(eps) *= 10.;
3406  dval(d) *= 10.;
3407  }
3408  }
3409  else {
3410 #endif
3411  /* Generate ilim digits, then fix them up. */
3412  dval(eps) *= tens[ilim-1];
3413  for (i = 1;; i++, dval(d) *= 10.) {
3414  L = (Long)(dval(d));
3415  if (!(dval(d) -= L))
3416  ilim = i;
3417  *s++ = '0' + (int)L;
3418  if (i == ilim) {
3419  if (dval(d) > 0.5 + dval(eps))
3420  goto bump_up;
3421  else if (dval(d) < 0.5 - dval(eps)) {
3422  while (*--s == '0') ;
3423  s++;
3424  goto ret1;
3425  }
3426  break;
3427  }
3428  }
3429 #ifndef No_leftright
3430  }
3431 #endif
3432 fast_failed:
3433  s = s0;
3434  dval(d) = dval(d2);
3435  k = k0;
3436  ilim = ilim0;
3437  }
3438 
3439  /* Do we have a "small" integer? */
3440 
3441  if (be >= 0 && k <= Int_max) {
3442  /* Yes. */
3443  ds = tens[k];
3444  if (ndigits < 0 && ilim <= 0) {
3445  S = mhi = 0;
3446  if (ilim < 0 || dval(d) <= 5*ds)
3447  goto no_digits;
3448  goto one_digit;
3449  }
3450  for (i = 1;; i++, dval(d) *= 10.) {
3451  L = (Long)(dval(d) / ds);
3452  dval(d) -= L*ds;
3453 #ifdef Check_FLT_ROUNDS
3454  /* If FLT_ROUNDS == 2, L will usually be high by 1 */
3455  if (dval(d) < 0) {
3456  L--;
3457  dval(d) += ds;
3458  }
3459 #endif
3460  *s++ = '0' + (int)L;
3461  if (!dval(d)) {
3462 #ifdef SET_INEXACT
3463  inexact = 0;
3464 #endif
3465  break;
3466  }
3467  if (i == ilim) {
3468 #ifdef Honor_FLT_ROUNDS
3469  if (mode > 1)
3470  switch (rounding) {
3471  case 0: goto ret1;
3472  case 2: goto bump_up;
3473  }
3474 #endif
3475  dval(d) += dval(d);
3476  if (dval(d) > ds || (dval(d) == ds && (L & 1))) {
3477 bump_up:
3478  while (*--s == '9')
3479  if (s == s0) {
3480  k++;
3481  *s = '0';
3482  break;
3483  }
3484  ++*s++;
3485  }
3486  break;
3487  }
3488  }
3489  goto ret1;
3490  }
3491 
3492  m2 = b2;
3493  m5 = b5;
3494  if (leftright) {
3495  i =
3496 #ifndef Sudden_Underflow
3497  denorm ? be + (Bias + (P-1) - 1 + 1) :
3498 #endif
3499 #ifdef IBM
3500  1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
3501 #else
3502  1 + P - bbits;
3503 #endif
3504  b2 += i;
3505  s2 += i;
3506  mhi = i2b(1);
3507  }
3508  if (m2 > 0 && s2 > 0) {
3509  i = m2 < s2 ? m2 : s2;
3510  b2 -= i;
3511  m2 -= i;
3512  s2 -= i;
3513  }
3514  if (b5 > 0) {
3515  if (leftright) {
3516  if (m5 > 0) {
3517  mhi = pow5mult(mhi, m5);
3518  b1 = mult(mhi, b);
3519  Bfree(b);
3520  b = b1;
3521  }
3522  if ((j = b5 - m5) != 0)
3523  b = pow5mult(b, j);
3524  }
3525  else
3526  b = pow5mult(b, b5);
3527  }
3528  S = i2b(1);
3529  if (s5 > 0)
3530  S = pow5mult(S, s5);
3531 
3532  /* Check for special case that d is a normalized power of 2. */
3533 
3534  spec_case = 0;
3535  if ((mode < 2 || leftright)
3536 #ifdef Honor_FLT_ROUNDS
3537  && rounding == 1
3538 #endif
3539  ) {
3540  if (!word1(d) && !(word0(d) & Bndry_mask)
3541 #ifndef Sudden_Underflow
3542  && word0(d) & (Exp_mask & ~Exp_msk1)
3543 #endif
3544  ) {
3545  /* The special case */
3546  b2 += Log2P;
3547  s2 += Log2P;
3548  spec_case = 1;
3549  }
3550  }
3551 
3552  /* Arrange for convenient computation of quotients:
3553  * shift left if necessary so divisor has 4 leading 0 bits.
3554  *
3555  * Perhaps we should just compute leading 28 bits of S once
3556  * and for all and pass them and a shift to quorem, so it
3557  * can do shifts and ors to compute the numerator for q.
3558  */
3559 #ifdef Pack_32
3560  if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0)
3561  i = 32 - i;
3562 #else
3563  if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) != 0)
3564  i = 16 - i;
3565 #endif
3566  if (i > 4) {
3567  i -= 4;
3568  b2 += i;
3569  m2 += i;
3570  s2 += i;
3571  }
3572  else if (i < 4) {
3573  i += 28;
3574  b2 += i;
3575  m2 += i;
3576  s2 += i;
3577  }
3578  if (b2 > 0)
3579  b = lshift(b, b2);
3580  if (s2 > 0)
3581  S = lshift(S, s2);
3582  if (k_check) {
3583  if (cmp(b,S) < 0) {
3584  k--;
3585  b = multadd(b, 10, 0); /* we botched the k estimate */
3586  if (leftright)
3587  mhi = multadd(mhi, 10, 0);
3588  ilim = ilim1;
3589  }
3590  }
3591  if (ilim <= 0 && (mode == 3 || mode == 5)) {
3592  if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
3593  /* no digits, fcvt style */
3594 no_digits:
3595  k = -1 - ndigits;
3596  goto ret;
3597  }
3598 one_digit:
3599  *s++ = '1';
3600  k++;
3601  goto ret;
3602  }
3603  if (leftright) {
3604  if (m2 > 0)
3605  mhi = lshift(mhi, m2);
3606 
3607  /* Compute mlo -- check for special case
3608  * that d is a normalized power of 2.
3609  */
3610 
3611  mlo = mhi;
3612  if (spec_case) {
3613  mhi = Balloc(mhi->k);
3614  Bcopy(mhi, mlo);
3615  mhi = lshift(mhi, Log2P);
3616  }
3617 
3618  for (i = 1;;i++) {
3619  dig = quorem(b,S) + '0';
3620  /* Do we yet have the shortest decimal string
3621  * that will round to d?
3622  */
3623  j = cmp(b, mlo);
3624  delta = diff(S, mhi);
3625  j1 = delta->sign ? 1 : cmp(b, delta);
3626  Bfree(delta);
3627 #ifndef ROUND_BIASED
3628  if (j1 == 0 && mode != 1 && !(word1(d) & 1)
3629 #ifdef Honor_FLT_ROUNDS
3630  && rounding >= 1
3631 #endif
3632  ) {
3633  if (dig == '9')
3634  goto round_9_up;
3635  if (j > 0)
3636  dig++;
3637 #ifdef SET_INEXACT
3638  else if (!b->x[0] && b->wds <= 1)
3639  inexact = 0;
3640 #endif
3641  *s++ = dig;
3642  goto ret;
3643  }
3644 #endif
3645  if (j < 0 || (j == 0 && mode != 1
3646 #ifndef ROUND_BIASED
3647  && !(word1(d) & 1)
3648 #endif
3649  )) {
3650  if (!b->x[0] && b->wds <= 1) {
3651 #ifdef SET_INEXACT
3652  inexact = 0;
3653 #endif
3654  goto accept_dig;
3655  }
3656 #ifdef Honor_FLT_ROUNDS
3657  if (mode > 1)
3658  switch (rounding) {
3659  case 0: goto accept_dig;
3660  case 2: goto keep_dig;
3661  }
3662 #endif /*Honor_FLT_ROUNDS*/
3663  if (j1 > 0) {
3664  b = lshift(b, 1);
3665  j1 = cmp(b, S);
3666  if ((j1 > 0 || (j1 == 0 && (dig & 1))) && dig++ == '9')
3667  goto round_9_up;
3668  }
3669 accept_dig:
3670  *s++ = dig;
3671  goto ret;
3672  }
3673  if (j1 > 0) {
3674 #ifdef Honor_FLT_ROUNDS
3675  if (!rounding)
3676  goto accept_dig;
3677 #endif
3678  if (dig == '9') { /* possible if i == 1 */
3679 round_9_up:
3680  *s++ = '9';
3681  goto roundoff;
3682  }
3683  *s++ = dig + 1;
3684  goto ret;
3685  }
3686 #ifdef Honor_FLT_ROUNDS
3687 keep_dig:
3688 #endif
3689  *s++ = dig;
3690  if (i == ilim)
3691  break;
3692  b = multadd(b, 10, 0);
3693  if (mlo == mhi)
3694  mlo = mhi = multadd(mhi, 10, 0);
3695  else {
3696  mlo = multadd(mlo, 10, 0);
3697  mhi = multadd(mhi, 10, 0);
3698  }
3699  }
3700  }
3701  else
3702  for (i = 1;; i++) {
3703  *s++ = dig = quorem(b,S) + '0';
3704  if (!b->x[0] && b->wds <= 1) {
3705 #ifdef SET_INEXACT
3706  inexact = 0;
3707 #endif
3708  goto ret;
3709  }
3710  if (i >= ilim)
3711  break;
3712  b = multadd(b, 10, 0);
3713  }
3714 
3715  /* Round off last digit */
3716 
3717 #ifdef Honor_FLT_ROUNDS
3718  switch (rounding) {
3719  case 0: goto trimzeros;
3720  case 2: goto roundoff;
3721  }
3722 #endif
3723  b = lshift(b, 1);
3724  j = cmp(b, S);
3725  if (j > 0 || (j == 0 && (dig & 1))) {
3726  roundoff:
3727  while (*--s == '9')
3728  if (s == s0) {
3729  k++;
3730  *s++ = '1';
3731  goto ret;
3732  }
3733  ++*s++;
3734  }
3735  else {
3736  while (*--s == '0') ;
3737  s++;
3738  }
3739 ret:
3740  Bfree(S);
3741  if (mhi) {
3742  if (mlo && mlo != mhi)
3743  Bfree(mlo);
3744  Bfree(mhi);
3745  }
3746 ret1:
3747 #ifdef SET_INEXACT
3748  if (inexact) {
3749  if (!oldinexact) {
3750  word0(d) = Exp_1 + (70 << Exp_shift);
3751  word1(d) = 0;
3752  dval(d) += 1.;
3753  }
3754  }
3755  else if (!oldinexact)
3756  clear_inexact();
3757 #endif
3758  Bfree(b);
3759  *s = 0;
3760  *decpt = k + 1;
3761  if (rve)
3762  *rve = s;
3763  return s0;
3764 }
3765 
3766 void
3767 ruby_each_words(const char *str, void (*func)(const char*, int, void*), void *arg)
3768 {
3769  const char *end;
3770  int len;
3771 
3772  if (!str) return;
3773  for (; *str; str = end) {
3774  while (ISSPACE(*str) || *str == ',') str++;
3775  if (!*str) break;
3776  end = str;
3777  while (*end && !ISSPACE(*end) && *end != ',') end++;
3778  len = (int)(end - str); /* assume no string exceeds INT_MAX */
3779  (*func)(str, len, arg);
3780  }
3781 }
3782 
3783 /*-
3784  * Copyright (c) 2004-2008 David Schultz <das@FreeBSD.ORG>
3785  * All rights reserved.
3786  *
3787  * Redistribution and use in source and binary forms, with or without
3788  * modification, are permitted provided that the following conditions
3789  * are met:
3790  * 1. Redistributions of source code must retain the above copyright
3791  * notice, this list of conditions and the following disclaimer.
3792  * 2. Redistributions in binary form must reproduce the above copyright
3793  * notice, this list of conditions and the following disclaimer in the
3794  * documentation and/or other materials provided with the distribution.
3795  *
3796  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
3797  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
3798  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
3799  * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
3800  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
3801  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
3802  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
3803  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
3804  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
3805  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
3806  * SUCH DAMAGE.
3807  */
3808 
3809 #define DBL_MANH_SIZE 20
3810 #define DBL_MANL_SIZE 32
3811 #define DBL_ADJ (DBL_MAX_EXP - 2)
3812 #define SIGFIGS ((DBL_MANT_DIG + 3) / 4 + 1)
3813 #define dexp_get(u) ((int)(word0(u) >> Exp_shift) & ~Exp_msk1)
3814 #define dexp_set(u,v) (word0(u) = (((int)(word0(u)) & ~Exp_mask) | ((v) << Exp_shift)))
3815 #define dmanh_get(u) ((uint32_t)(word0(u) & Frac_mask))
3816 #define dmanl_get(u) ((uint32_t)word1(u))
3817 
3818 
3819 /*
3820  * This procedure converts a double-precision number in IEEE format
3821  * into a string of hexadecimal digits and an exponent of 2. Its
3822  * behavior is bug-for-bug compatible with dtoa() in mode 2, with the
3823  * following exceptions:
3824  *
3825  * - An ndigits < 0 causes it to use as many digits as necessary to
3826  * represent the number exactly.
3827  * - The additional xdigs argument should point to either the string
3828  * "0123456789ABCDEF" or the string "0123456789abcdef", depending on
3829  * which case is desired.
3830  * - This routine does not repeat dtoa's mistake of setting decpt
3831  * to 9999 in the case of an infinity or NaN. INT_MAX is used
3832  * for this purpose instead.
3833  *
3834  * Note that the C99 standard does not specify what the leading digit
3835  * should be for non-zero numbers. For instance, 0x1.3p3 is the same
3836  * as 0x2.6p2 is the same as 0x4.cp3. This implementation always makes
3837  * the leading digit a 1. This ensures that the exponent printed is the
3838  * actual base-2 exponent, i.e., ilogb(d).
3839  *
3840  * Inputs: d, xdigs, ndigits
3841  * Outputs: decpt, sign, rve
3842  */
3843 char *
3844 ruby_hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
3845  char **rve)
3846 {
3847  U u;
3848  char *s, *s0;
3849  int bufsize;
3850  uint32_t manh, manl;
3851 
3852  u.d = d;
3853  if (word0(u) & Sign_bit) {
3854  /* set sign for everything, including 0's and NaNs */
3855  *sign = 1;
3856  word0(u) &= ~Sign_bit; /* clear sign bit */
3857  }
3858  else
3859  *sign = 0;
3860 
3861  if (isinf(d)) { /* FP_INFINITE */
3862  *decpt = INT_MAX;
3863  return rv_strdup(INFSTR, rve);
3864  }
3865  else if (isnan(d)) { /* FP_NAN */
3866  *decpt = INT_MAX;
3867  return rv_strdup(NANSTR, rve);
3868  }
3869  else if (d == 0.0) { /* FP_ZERO */
3870  *decpt = 1;
3871  return rv_strdup(ZEROSTR, rve);
3872  }
3873  else if (dexp_get(u)) { /* FP_NORMAL */
3874  *decpt = dexp_get(u) - DBL_ADJ;
3875  }
3876  else { /* FP_SUBNORMAL */
3877  u.d *= 5.363123171977039e+154 /* 0x1p514 */;
3878  *decpt = dexp_get(u) - (514 + DBL_ADJ);
3879  }
3880 
3881  if (ndigits == 0) /* dtoa() compatibility */
3882  ndigits = 1;
3883 
3884  /*
3885  * If ndigits < 0, we are expected to auto-size, so we allocate
3886  * enough space for all the digits.
3887  */
3888  bufsize = (ndigits > 0) ? ndigits : SIGFIGS;
3889  s0 = rv_alloc(bufsize+1);
3890 
3891  /* Round to the desired number of digits. */
3892  if (SIGFIGS > ndigits && ndigits > 0) {
3893  float redux = 1.0f;
3894  int offset = 4 * ndigits + DBL_MAX_EXP - 4 - DBL_MANT_DIG;
3895  dexp_set(u, offset);
3896  u.d += redux;
3897  u.d -= redux;
3898  *decpt += dexp_get(u) - offset;
3899  }
3900 
3901  manh = dmanh_get(u);
3902  manl = dmanl_get(u);
3903  *s0 = '1';
3904  for (s = s0 + 1; s < s0 + bufsize; s++) {
3905  *s = xdigs[(manh >> (DBL_MANH_SIZE - 4)) & 0xf];
3906  manh = (manh << 4) | (manl >> (DBL_MANL_SIZE - 4));
3907  manl <<= 4;
3908  }
3909 
3910  /* If ndigits < 0, we are expected to auto-size the precision. */
3911  if (ndigits < 0) {
3912  for (ndigits = SIGFIGS; s0[ndigits - 1] == '0'; ndigits--)
3913  ;
3914  }
3915 
3916  s = s0 + ndigits;
3917  *s = '\0';
3918  if (rve != NULL)
3919  *rve = s;
3920  return (s0);
3921 }
3922 
3923 #ifdef __cplusplus
3924 #if 0
3925 { /* satisfy cc-mode */
3926 #endif
3927 }
3928 #endif
#define d0
#define Sign_bit
Definition: util.c:853
#define mmstep
Definition: util.c:198
#define dexp_set(u, v)
Definition: util.c:3814
#define ISDIGIT(c)
Definition: ruby.h:1775
static int lo0bits(ULong *y)
Definition: util.c:1174
#define FLT_RADIX
Definition: numeric.c:43
#define Big1
Definition: util.c:957
#define R(b, x)
Definition: sha2.c:203
size_t strlen(const char *)
#define ACQUIRE_DTOA_LOCK(n)
Definition: util.c:990
#define rv_alloc(i)
Definition: util.c:3024
#define d1
#define P
Definition: util.c:840
static Bigint * Balloc(int k)
Definition: util.c:1007
#define rv_strdup(s, rve)
Definition: util.c:3039
#define DBL_DIG
Definition: numeric.c:67
#define dmanh_get(u)
Definition: util.c:3815
#define PATH_MAX
static Bigint * pow5mult(Bigint *b, int k)
Definition: util.c:1331
#define med3(a, b, c)
Definition: util.c:300
#define Exp_1
Definition: util.c:843
unsigned long ruby_scan_hex(const char *start, size_t len, size_t *retlen)
Definition: util.c:42
static const char ZEROSTR[]
Definition: util.c:3057
#define dmanl_get(u)
Definition: util.c:3816
#define Kmax
Definition: util.c:994
char * RR
Definition: util.c:296
#define Quick_max
Definition: util.c:857
#define DBL_MANL_SIZE
Definition: util.c:3810
#define Tiny0
Definition: util.c:855
int sign
Definition: util.c:998
#define Bletch
Definition: util.c:849
static void mmrot3_(register char *a, register char *b, register char *c, mmargdecl)
Definition: util.c:250
#define C
Definition: util.c:195
static double ulp(double x_)
Definition: util.c:1552
#define Storeinc(a, b, c)
Definition: util.c:821
struct Bigint Bigint
Definition: util.c:1002
#define Int_max
Definition: util.c:858
#define Ebits
Definition: util.c:845
#define PRIVATE_mem
Definition: util.c:731
#define S(s)
#define A
Definition: util.c:193
static int hi0bits(register ULong x)
Definition: util.c:1145
static Bigint * lshift(Bigint *b, int k)
Definition: util.c:1383
Definition: util.c:996
static const char NANSTR[]
Definition: util.c:3056
static double private_mem[PRIVATE_mem]
Definition: util.c:732
#define DBL_MAX_10_EXP
Definition: numeric.c:64
int( cmpfunc_t)(const void *, const void *, void *)
Definition: util.c:304
#define LSB
Definition: util.c:852
ULong x[1]
Definition: util.c:999
#define Emin
Definition: util.c:842
static char * nrv_alloc(const char *s, char **rve, size_t n)
Definition: util.c:3028
#define fail()
#define dexp_get(u)
Definition: util.c:3813
static unsigned long scan_digits(const char *str, int base, size_t *retlen, int *overflow)
Definition: util.c:79
#define mmargdecl
Definition: util.c:209
#define B
Definition: util.c:194
#define Exp_mask
Definition: util.c:839
#define Rounding
Definition: util.c:879
static int quorem(Bigint *b, Bigint *S)
Definition: util.c:2903
#define DBL_MANH_SIZE
Definition: util.c:3809
static double one(void)
Definition: isinf.c:52
#define n_bigtens
Definition: util.c:1851
#define SIGFIGS
Definition: util.c:3812
#define POP(ll, rr)
Definition: util.c:298
#define Bndry_mask1
Definition: util.c:851
RUBY_EXTERN int isinf(double)
Definition: isinf.c:56
int maxwds
Definition: util.c:998
double d
Definition: util.c:792
static Bigint * s2b(const char *s, int nd0, int nd, ULong y9)
Definition: util.c:1111
static double * pmem_next
Definition: util.c:732
Definition: util.c:792
double ruby_strtod(const char *s00, char **se)
Definition: util.c:1945
#define Scale_Bit
Definition: util.c:1850
#define IEEE_Arith
Definition: util.c:741
SSL_METHOD *(* func)(void)
Definition: ossl_ssl.c:113
int errno
static const double tens[]
Definition: util.c:1828
#define rounded_product(a, b)
Definition: util.c:952
int wds
Definition: util.c:998
static const char INFSTR[]
Definition: util.c:3055
void ruby_qsort(void *base, const size_t nel, const size_t size, cmpfunc_t *cmp, void *d)
Definition: util.c:306
#define Avoid_Underflow
Definition: util.c:860
unsigned char buf[MIME_BUF_SIZE]
Definition: nkf.c:4308
#define no_digits()
static const double bigtens[]
Definition: util.c:1839
static int cmp(Bigint *a, Bigint *b)
Definition: util.c:1437
const signed char ruby_digit36_to_number_table[]
Definition: util.c:58
#define Log2P
Definition: util.c:854
#define rounded_quotient(a, b)
Definition: util.c:953
static Bigint * multadd(Bigint *b, int m, int a)
Definition: util.c:1060
char * strchr(char *, char)
static double b2d(Bigint *a, int *e)
Definition: util.c:1590
#define Ten_pmax
Definition: util.c:848
#define mmswap(a, b)
Definition: util.c:247
#define isnan(x)
Definition: win32.h:376
void rb_sys_fail(const char *mesg)
Definition: error.c:1973
#define MALLOC
Definition: util.c:719
#define DBL_MANT_DIG
Definition: acosh.c:19
#define Tiny1
Definition: util.c:856
#define ULLong
Definition: util.c:980
#define word1(x)
Definition: util.c:808
#define CHAR_BIT
Definition: ruby.h:198
#define DBL_MAX_EXP
Definition: numeric.c:58
unsigned long ruby_scan_oct(const char *start, size_t len, size_t *retlen)
Definition: util.c:28
static Bigint * p5s
Definition: util.c:1328
unsigned int uint32_t
Definition: sha2.h:101
#define mmtype
Definition: util.c:191
#define Bndry_mask
Definition: util.c:850
unsigned int top
Definition: nkf.c:4309
#define IEEE_LITTLE_ENDIAN
Definition: util.c:675
#define Exp_shift1
Definition: util.c:836
int size
Definition: encoding.c:49
#define DBL_ADJ
Definition: util.c:3811
static Bigint * d2b(double d_, int *e, int *bits)
Definition: util.c:1655
#define xmalloc
Definition: defines.h:108
#define PUSH(ll, rr)
Definition: util.c:297
char * ruby_hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign, char **rve)
Definition: util.c:3844
#define Exp_msk11
Definition: util.c:838
U double_u
Definition: util.c:805
#define Frac_mask1
Definition: util.c:847
#define FREE_DTOA_LOCK(n)
Definition: util.c:991
#define mmprepare(base, size)
Definition: util.c:199
static void mmswap_(register char *a, register char *b, mmargdecl)
Definition: util.c:211
#define mmrot3(a, b, c)
Definition: util.c:285
static Bigint * diff(Bigint *a, Bigint *b)
Definition: util.c:1466
#define Big0
Definition: util.c:956
struct Bigint * next
Definition: util.c:997
#define Exp_msk1
Definition: util.c:837
static double ratio(Bigint *a, Bigint *b)
Definition: util.c:1792
int k
Definition: util.c:998
#define xrealloc
Definition: defines.h:111
static void Bfree(Bigint *v)
Definition: util.c:1042
#define FFFFFFFF
Definition: util.c:963
unsigned long ruby_strtoul(const char *str, char **endptr, int base)
Definition: util.c:107
#define Flt_Rounds
Definition: util.c:870
char * ruby_getcwd(void)
Definition: util.c:478
void ruby_each_words(const char *str, void(*func)(const char *, int, void *), void *arg)
Definition: util.c:3767
void void xfree(void *)
static const double tinytens[]
Definition: util.c:1840
#define Bcopy(x, y)
Definition: util.c:1056
#define word0(x)
Definition: util.c:807
#define NULL
Definition: _sdbm.c:103
static Bigint * freelist[Kmax+1]
Definition: util.c:1004
static int match(VALUE str, VALUE pat, VALUE hash, int(*cb)(VALUE, VALUE))
Definition: date_parse.c:273
#define Exp_11
Definition: util.c:844
#define dval(x)
Definition: util.c:813
char * ruby_strdup(const char *str)
Definition: util.c:457
static Bigint * i2b(int i)
Definition: util.c:1217
static Bigint * mult(Bigint *a, Bigint *b)
Definition: util.c:1228
#define Bias
Definition: util.c:841
#define ISSPACE(c)
Definition: ruby.h:1770
#define Frac_mask
Definition: util.c:846
#define Exp_shift
Definition: util.c:835
char * ruby_dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
Definition: util.c:3094
#define FREE
Definition: util.c:724