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