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