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